Как найти сопряженную задачу

Once the solution of the adjoint problem is found, the functional may be represented solely by the second integral in Eq. (13.58), as follows

(13.60)F(x0)=∫0Lδ(x−x0)QC⁎(x)dx=QC⁎(x0)

The minimization of F(x0) is now trivial, as C⁎ decays exponentially downstream of the target point. The optimum location for x0 would be at the downstream boundary, but for all practical purposes, equally good solutions can be found much closer to the target point by adopting a small, but finite limit for the concentration. More importantly, only one solution of the adjoint problem is necessary, compared to N solutions of the main problem in order to arrive at the same conclusion.

It is clear that the main and adjoint problems have identical characteristic directions. At the upstream boundary of the main problem, a C+ characteristic travels away from the boundary in the forward time direction. Along this path, as dt becomes infinitesimally small, the compatibility relation reads

(14.29)dq−(u+c)dh=0

Similarly, at the downstream boundary a C− characteristic travels away from the boundary in the forward time direction. Along this direction, the corresponding compatibility relation is as follows

(14.30)dq−(u−c)dh=0

Therefore, substitution of (14.29) for δh(0,t) and (14.30) for δh(L,t) into the expression for the first variation of the functional, yields

(14.31)δF=−∫0T{[(ϕ+(u−c)ψ)δq]x=L−[(ϕ+(u+c)ψ)δq]x=0}dt

This ensures that perturbations generated within the channel travel unobstructed through the boundaries, which is precisely the behavior that is desired at open computational boundaries. Now, recalling that the adjoint problem solution evaluates the system response to individual, instantaneous system perturbations, either δq(0,tp) or δq(L,tp) is non-zero at the perturbation time tp, but not both. Therefore, the sensitivity relations for upstream and downstream perturbations of discharge are given by

(14.32)δFδq(0,tp)=[ϕ+(u+c)ψ]x=0,t=tp

and

(14.33)δFδq(L,tp)=−[ϕ+(u−c)ψ]x=L,t=tp

The expressions, A1=ϕ+(u+c)ψ and A2=ϕ+(u−c)ψ are representative of the characteristic variables of the adjoint problem, representing sensitivity signals traveling towards the boundaries in the reverse time direction.

Similar expressions corresponding to the sensitivity of the functional F to perturbations in depth are derived by substituting (14.29) for δq(0,t) and (14.30) for δq(L,t) into the functional variation. This yields

(14.34)δFδh(0,tp)=[(ϕ+(u+c)ψ)(u+c)]x=0,t=tp

and

(14.35)δFδh(L,tp)=−[(ϕ+(u−c)ψ)(u−c)]x=L,t=tp

In the limiting process used to obtain the adjoint problem, the following expressions are left, for the no-ablation and ablation periods, respectively:

(13a)ΔSQτ=−∫τ=0τ0λ1τΔQτdτ

(13b)ΔSQτ=−∫τ=τ0τfλBτ,τΔQτdτ

By assuming that the unknown function Q(τ) belongs to the space of square integrable functions in the time domain 0 < τ < τ_f, the directional derivative of S[Q(τ)] in the direction of the perturbation ΔQ(τ) can be written as [4]

(14)ΔSQτ=−∫τ=0τf∇SQτΔQτdτ

Therefore, by comparing Eqs. (13a,b) with Eq. (14), the following equations result for the gradient of the functional:

(15a,b)∇SQτ=−λ1τ,fortheno-ablationperiod−λBτ,τ,fortheablationperiod

The concept of least squares fitting has additional justification in self-adjoint problems in which an energy functional is minimized. In such cases, typical of a displacement formulation of elasticity, it can be readily shown that the minimization is equivalent to a least squares fit of the approximate stresses to the exact ones. Thus quite generally we can start from a theory given by the differential equation

(15.16)Lu=ST(ASu)=p

In the above, L is a self-adjoint operator defined by S and A (symmetric) and p are prescribed matrices of position. The minimization of an energy functionalΠ defined as

(15.17)Π=12∫ΩSuTASudΩ-∫ΩuTpdΩ

gives at an absolute minimum the exact solution u=u¯ and is equivalent to minimization of another functional Π★ defined as

(15.18)Π★=12∫ΩSu-u¯TASu-u¯dΩ

The above quadratic functional [Eq. (15.17)] arises in all linear self-adjoint problems.

The proof of the Herrmann theorem is as follows. The variation of Π defined in Eq. (15.17) gives, at u=u¯ (the exact solution),

δΠ=12∫ΩSδuTASu¯dΩ+12∫ΩSu¯TASδudΩ-∫ΩδuTpdΩ=0

or as A is symmetric

δΠ=∫ΩSδuTASu¯dΩ-∫ΩδuTpdΩ=0

in which δu is any arbitrary variation. Thus we can select

δu=u

and

∫ΩSuTASu¯dΩ-∫ΩuTpdΩ=0

Subtracting the above from Eq. (15.17) and noting the symmetry of the A matrix, we can write

(15.19)Π=12∫ΩSu-u¯TASu-u¯dΩ-12∫ΩSu¯TASu¯dΩ

where the last term is not subject to variation. Thus

(15.20)Π★=Π+constant

and its stationarity is equivalent to the stationarity of Π.

It follows directly from the Herrmann theorem that, for one dimension and by a well-known property of the Gauss-Legendre quadrature points, if the approximate gradients are defined by a polynomial of degree p-1, where p is the degree of the polynomial used for the unknown function u, then stresses taken at these quadrature points must be superconvergent. The single point at the center of an element integrates precisely all linear functions passing through that point and, hence, if the stresses are exact to the linear form they will be exact at that point of integration. For any higher order polynomial of order p, the Gauss-Legendre points numbering p will also provide points of superconvergent sampling. We see this from Fig. 15.5 directly. Here we indicate one-, two-, and three-point Gauss-Legendre quadrature showing why exact results are recovered there for gradients and stresses.

Figure 15.5. The integration property of Gauss points: p=1,p=2, and p=3 which guarantees superconvergence.

For points based on rectangles and products of polynomial functions it is clear that the exact integration points will exist at the product points as shown in Fig. 15.6 for various rectangular elements assuming that the weighting matrix A is diagonal. In the same figure we show some triangles that appear to be “good” but are not necessarily superconvergent sampling points. Though we find that superconvergent points do not exist in triangles, the points shown in Fig. 15.6 are optimal. In Fig. 15.6 we contrast these points with the minimum number of quadrature points necessary for obtaining an accurate (though not always stable) stiffness representation and find these to be almost coincident at all times.

Figure 15.6. Optimal superconvergent sampling and minimum integration points for some C0 elements.

In Fig. 15.7 representing an analysis of a cantilever beam by four rectangular quadratic serendipity elements we see how well the stresses sampled at superconvergent points behave compared to the overall stress pattern computed in each element.

Figure 15.7. Cantilever beam with four quadratic (Q8) elements. Stress sampling at cubic order (2×2) Gauss points with extrapolation to nodes.

As a first step in the investigation of the properties of equation (68), the adjoint problem, equation (66), is multiplied with an arbitrary fuel composition, subject only to criticality and reloading constraints, as given in Section 3.3.1. Solving for γc∞ yields

(69)γC∞=(N’,C_*w*)(N’,w*).

Revolving the vectors in equation (69) turns the matrix in the numerator into the respective adjoint and thus back to the real matrix (C^** = C). This gives

(70)γC∞=(w*,C_N’)(w*,N’)=(w*,[PL−]N’)(w*,N’).

In the same way as γc∞ in equation (65) was independent of the worth factors, γc∞ in equation (70) is independent of the fuel composition (provided it satisfies the criticality and reloading constraints).

The accuracy that can be expected from the γ-estimation formula based on breeding worth weighting (equation (68) or equation (70)) is investigated in the following. For this, estimates and exact values need to be distinguished, notationally. Let w^* in equation (66) be the ‘exact’ value, corresponding to the matrix C, and wˆ* the approximate value, used for the γ-estimate:

(71)γCest=(w˚*,C_N’)(w˚*,N’).

Note, the w^* in equation (68) is meant to be a ‘typical value’; it therefore corresponds to the wˆ* in equation (71); whereas w^* in equation (70) solves the adjoint problem of the C matrix in the same equation:

(72)C_*w*=γC∞w*.

The Ns for that problem are obtained by solving

(73)C_NC=γC∞NC.

The deviations from the correct values are given by

(74)w*=w˙*+Δw*.

(75)NC=N′+ΔN.

Inserting equations (74) and (75) into equation (71) and applying some algebra leads to the following equation

(76)γCest=γC∞+(Δw˚*,[C_−γC∞]ΔN)(w˚*,N’)

Thus, by applying the breeding worth weighting, the growth-rate estimate becomes stationary, in the sense that the deviation between γcest and γc∞ is of second order. Again, the estimate becomes exact when either Δw^* or ΔN are zero. The stationarity of equation (68) is demonstrated numerically in Section 6 below.

The equation of motion associated with the problem defined by (2.45) is

(2.47)EI∂4w∂x4+P∂2w∂x2+m∂2w∂t2=0,

with the boundary conditions as given in (2.46). This clearly represents free motions of the system; hence, of interest are the eigenfrequencies and the corresponding eigenfunctions and how they vary with P (or its nondimensional counterpart, PL2/EI). This can be done by direct application of the Galerkin method with wN=∑jϕj(x)qj(t), in which the cantilever-beam eigenfunctions (2.27) are used as comparison functions, since they satisfy boundary conditions (2.46), which are identical to (2.23). In this way, one obtains an equation similar to (2.35), i.e.

(2.48)[M]{q¨}+[K]{q}={0},

but with only [M] being diagonal, while [K] is nondiagonal. In fact, the elements of [K] are

krj=EIλr4Lδrj+P∫0Lϕrϕj″dx,

the prime denoting differentiation with respect to x.

Suppose now that this system is subjected also to a distributed force, F(x,t), so that the equation of motion is

(2.49)EI∂4w∂x4+P∂2w∂x2+m∂2w∂t2=F(x,t);

(2.50)[M]{q¨}+[K]{q}={Q}.

If this had been a self-adjoint conservative system, matrices [M] and [K] in equation (2.50) would both be symmetric. For the problem at hand, however, the system is non-self-adjoint, as remarked earlier, and hence [K] is asymmetric, by virtue of the fact that ∫0Lϕrϕj″dx≠∫0Lϕjϕr″dx. Hence, the decoupling procedure leading to equation (2.15) should be adopted.

Before proceeding further, however, it is useful to transform equation (2.49) into dimensionless form, which also serves to introduce the kind of dimensionless terms appearing frequently in the following chapters. Hence, defining

(2.51)ξ=x/L,η=w/L,τ=(EI/mL4)12t,P=PL2/EI,f=FL3/EI,ω=(EI/mL4)-12Ω

and taking, as a concrete example, f=f0ξsinωfτ — a triangularly distributed load along the beam, as in Figure 2.2, — substitution into (2.49) yields

(2.52)η⁗+Pη″+η¨=f0ξsinωfτ,

in which primes and overdots denote, respectively, partial differentiation with respect to ξ and τ. The discretized form of (2.52) is

(2.53)[I]{q¨}+[K]{q}={Q}sinωfτ;

[I] is the unit matrix, and the elements of [K] and {Q} are

(2.54)kij=λi4δij+P∫01ϕiϕj″dξ,Qi=∫01f0ξϕidξ,

in which the ϕi≡ϕi(ξ), the dimensionless version of (2.27). The decoupled equation, corresponding to equation (2.15), is

(2.55){y¨}+[Λ]{y}=[A]-1{Q}sinωfτ={Ψ}sinωfτ,

where [Λ] is the diagonal matrix of the eigenvalues; the solution therefore is

(2.56)yk=αkcosΛk1/2τ+βksinΛk1/2τ+Ψk/Λk-ωf2sinωfτ,k=1,2,…,N.

Numerical results for the case of P=1,f0=7,ωf=0.6 are shown in Figure 2.3: (a) for αk=βk=0, i.e. showing only the particular solution, and (b,c) for η(1,0)=0.15,η̇(1,0)=1.5. The dimensionless natural frequencies, obtained with N=4, are found to be ω1=3.64,ω2=21.73,w3=61.32 and ω4=120.5;ωf is chosen to be far below all of them.

Figure 2.3. Solutions to equation (2.52) showing η(1,τ) versus τ, for P=1,f0=7 and ωf=0.6: (a) the particular solution alone [i.e. αk=βk=0 in equation (2.56)], which would correspond to the steady-state solution if damping were included; (b) full solution for N=1; (c) full solution for N=2 (—) and N=4 (- – -) on an expanded scale of τ.

The fact that the response in Figure 2.3(b,c) is quasiperiodic is most apparent in the phase plane, as shown for example in Figure 2.4. It is seen that the response evolves by winding itself around a torus, the projection of which is shown in the figure, instead of tracing a planar curve, as would be the case for periodic motion.

Figure 2.4. The response of Figure 2.3(b) plotted in the phase plane: the dimensionless tip velocity, η̇(1,τ), versus displacement, η(1,τ).

There is another, more general method for obtaining the response of such a system, specific to non-self-adjoint problems (Washizu 1966, 1968Washizu, 1966Washizu, 1968; Anderson 1972). This begins with the determination of the adjoint problem.^† If the eigenfunctions of the homogeneous form of equation (2.47), i.e. of the compressively loaded beam, are χi(ξ) and those of the adjoint problem ψi(x), the adjoint problem is defined through the new operators L∗ and M∗, such that

(2.57)∫Dψ(ξ)L[χ(ξ)]dD-∫Dχ(ξ)L∗[ψ(ξ)]dD=C[χ(ξ),ψ(ξ)]D,

in which it is required that the so-called concomitant, C, vanish. A similar expression for M should be satisfied, but since M here is a scalar, we immediately have M∗=M. In the nondimensional notation used here, D=[0,1] and ξ=x/L. This problem has in fact been solved by Chen (1987), but it is not difficult to reproduce the results. One finds L∗=L, but a new set of boundary conditions for the adjoint problem, namely

(2.58)ψ(0)=ψ′(0)=0,ψ″(1)+Pψ(1)=0,ψ‴(1)+Pψ′(1)=0.

Solving the two eigenvalue problems, one obtains^‡

(2.59)χ(ξ)=A1sinpξ+A2cospξ+A3sinhqξ+A4coshqξ,ψ(ξ)=B1sinpξ+B2cospξ+B3sinhqξ+B4coshqξ,p=14P2+λ12+12P12,q=14P2+λ12-12P12,A1=1,A2=-(p2sinp+pqsinhq)/(p2cosp+q2coshq),A3=-p/q,A4=-A2,B1=1,B2=-(P-p2)sinp-(p/q)(P+q2)sinhq÷(P-p2)cosp-(P+q2)coshq,B3=-p/q,B4=-B2.

The characteristic equation is

P2+2λ(1+cospcoshq)+Pλsinpsinhq=0,

and it is the same for both problems; hence, so are the eigenvalues.

The essence of this method is that it achieves direct decoupling of the equations of motion via the so-called biorthogonality of the initial and adjoint eigenfunctions, viz.

(2.60)∫01ψrL[χj]dξ=krjδrj,∫01ψrM[χj]dξ=mrjδrj.

By introducing ηN=Σχj(ξ)qj(τ) into equation (2.52), then multiplying by ψr(ξ) and integrating over D, the system is decoupled in a single operation, by virtue of relations (2.60), yielding

(2.61)mjq¨j+kjqj=fjsinωfτ,j=1,2,…,N.

The equation of motion associated with the problem defined by (2.45) is

(2.47)EI∂4w∂x4+P∂2w∂x2+ m∂2w∂t2=0,

with the boundary conditions as given in (2.46). This clearly represents free motions of the system; hence, of interest are the eigenfrequencies and the corresponding eigenfunctions and how they vary with P (or its nondimensional counterpart, PL²/EI). This can be done by direct application of the Galerkin method with
wN=∑jφj(x)qj(t), in which the cantilever-beam eigenfunctions (2.27) are used as comparison functions, since they satisfy boundary conditions (2.46), which are identical to (2.23). In this way, one obtains an equation similar to (2.35), i.e.

(2.48)[M] {q¨}+[K] {q}={0},

but with only [M] being diagonal, while [K] is nondiagonal. In fact, the elements of [K] are

krj= EIλr4Lδrj+P∫L0φrφ″j    dx,

the prime denoting differentiation with respect to x.

Suppose now that this system is subjected also to a distributed force, F(x, t), so that the equation of motion is

(2.49)EI∂4w∂x4+P∂2w∂x2+ m∂2w∂t2=F(x,t);

(2.50)[M] {q¨}+[K] {q}={Q}⋅

If this had been a self-adjoint conservative system, matrices [M] and [K] in equation (2.50) would both be symmetric. For the problem at hand, however, the system is non-self-adjoint, as remarked earlier, and hence [K] is asymmetric, by virtue of the fact that
∫L0φrφ″jdx≠ ∫L0φjφ″rdx⋅Hence, the decoupling procedure leading to equation (2.15) should be adopted.

Before proceeding further, however, it is useful to transform equation (2.49) into dimensionless form, which serves to introduce the kind of dimensionless terms appearing frequently in the following chapters. Hence, defining

(2.51)ξ = x/L,      η=w/L,       τ=(EI/mL4)1/2t,P=PL2/EI ,      f=FL3/EI,   ω=(EI/mL4)−1/2Ω

and taking, as a concrete example,
f=f0ξ  sin  (ωfτ) —representing a triangularly distributed load along the beam, as shown in Figure 2.2 —substitution into (2.49) yields

(2.52)η‴′+Pη″+η¨= f0ξ sin(ωfτ),

in which primes and overdots denote, respectively, partial differentiation with respect to ξ and τ. The discretized form of (2.52) is

(2.53)[I] {q¨}+[K] {q}={Q }sin(ωfτ),

and the elements of [K] and [Q] are

(2.54)kij=λi4δij+P∫01φiφ″j dξ,      Qi=∫01f0ξφi dξ,

in which the ϕ_i ≡ ϕ_i(ξ), the dimensionless version of (2.27). The decoupled equation, corresponding to equation (2.15), is

(2.55){y¨}+[Λ {y} = [A] −1{Q}sin(ωfτ)={ψ}sin(ωfτ),

in which [Λ] is the diagonal matrix of the eigenvalues; the solution therefore is

(2.56)yk=αkcos Λk1/2τ+βksin Λk1/2τ+ [ψk/(Λk−wf2) sin(ωfτ),   k=  1,2,⋯,N⋅

Numerical results for the case of P = 1, f₀ = 7, ω_f = 0.6 are shown in Figure 2.3: (a) for α_k = β_k = 0, i.e. showing only the particular solution, and (b,c) for η(1, 0) = 0.15, η˙(1, 0) = 1.5. The dimensionless natural frequencies, obtained with N = 4, are found to be ω₁ = 3.64, ω₂ = 21.73, ω₃ = 61.32 and ω₄ = 120.5; ω_f is chosen to be far below all of them.

Figure 2.3. Solutions to equation (2.52) showing η(1, τ) versus τ, for P = 1, f₀ = 7 and ω_f = 0.6: (a) the particular solution alone [i.e. α_k = β_k = 0 in equation (2.56)], which would correspond to the steady-state solution if damping were included; (b) full solution for N = 1; (c) full solution for N = 2 (—) and N = 4 (· · ·) on an expanded scale of τ.

The fact that the response in Figure 2.3(b,c) is quasiperiodic is most apparent in the phase plane, as shown for example in Figure 2.4. It is seen that the response evolves by winding itself around a torus, the projection of which is shown in the figure, instead of tracing a planar curve, as would be the case for periodic motion.

Figure 2.4. The response of Figure 2.3(b) plotted in the phase plane: the dimensionless tip velocity, η(1, τ), versus displacement, η(1, τ).

There is another, more general method for obtaining the response of such a system, specific to non-self-adjoint problems (Washizu 1966, 1968; Anderson 1972). This begins with the determination of the adjoint problem.^†
If the eigenfunctions of the homogeneous form of equation (2.47), i.e. of the compressively loaded beam, are χ_i(ξ) and those of the adjoint problem ψ_i(ξ), the adjoint problem is defined through the new operators L^* and M^*, such that

(2.57)∫Dψ(ξ)L[ℵ(ξ) dD  −   ∫Dℵ(ξ)L*[ψ(ξ) dD  =  C[ℵ(ξ),ψ(ξ) D 

in which it is required that the so-called concomitant, C, vanish. A similar expression for M should be satisfied, but since for the problem at hand M is a scalar, we immediatelyhave
M^* =
M. In the nondimensional notation used here,
D = [0, 1] and ξ = x/L. This problem has in fact been solved by Chen (1987), but it is not difficult to reproduce the results. One finds
L^* =
L, but a new set of boundary conditions for the adjoint problem, namely

(2.58)ψ(0)=ψ′(0)=0,    ψ″(1)+ψP(1)=0,        ψ‴(1)+ψ′(1)=0⋅

Solving the two eigenvalue problems, one obtains

(2.59)X(ξ)=A1 sin pξ+A2cos pξ+A3sin qξ+A4cosh qξ,ψ(ξ)=B1 sin pξ+B2cos pξ+B3sin qξ+B4cosh qξ,      p=[(14P2+λ)1/2+12P 1/2,   q=[(14P2+λ)1/2−12P 1/2,†    A1=1,   A2=−(p2sin p+pqsinq)/(p2cos p+q2coshq),    A3=−p/q,   A4=−A2,    B1=1,        B2=[(P−p2)sinp−(p/q)(P+q2)sinq] [(P−p2)cosp−(P+q2)coshq]     B3=−p/q,      B4=−B2.

^‡
The characteristic equation is

P2+2λ(1+cos p cosh q)  + P λsin p sinh q=0,

and it is the same for both problems; hence, so are the eigenvalues.

The essence of this method is that it achieves direct decoupling of the equations of motion via the so-called biorthogonality of the initial and adjoint eigenfunctions, viz.

(2.60)∫01ψrL[ℵj dξ=krjδrj,        ∫01ψrM[ℵj dξ=mrjδrj⋅

By introducing
ηN= ∑  ℵj(ξ)qj(τ) into equation(2.52), then multiplying by ψ_r(ξ) and integrating over
D, the system is decoupled in a single operation, by virtue of relations (2.60), yielding

(2.61)mjq¨j+kjqj=fj sin ωfτ,    j=1,2,⋯,N⋅

The methods of estimating errors and adaptive refinement which are described in this and the previous chapter constitute a very important tool for practical application of finite element methods. The range of applications is large and we have only touched here upon the relatively simple range of linear elasticity and similar self-adjoint problems. We would like to reiterate that in other classes of problems many different norms or measures of error can be used and that for some problems the energy norm is not in fact “natural.” A good example of this is given by problems of high-speed gas flow, where very steep gradients (shocks) can develop. The formulation of such problems is complex [31], but this is not necessary for the present argument. Applications of error estimates and adaptive refinement in this and many more areas are reported in Refs. [32, 33] and the reader is referred to these publications for additional details. pFor problems in fluid mechanics discussed in Ref. [31] and similarly for problems of strain localization in plastic softening discussed in Ref. [34] no global norms can be used effectively. In such situations it is convenient to base the refinement on the value of the maximum curvatures developed by the solution of u. On occasion an elongation of the elements will be used to refine the mesh appropriately. Such anisotropic mesh adaptation is especially effective when it is based on the recovery type of error estimators [35–38]. Figure 16.23 shows a typical problem of shock capturing solved adaptively.

Figure 16.23. Directional mesh refinement. Gas flow past a circular cylinder: Mach number 3. Third refinement mesh 709 nodes (1348 elements): (a) Local mesh and (b) Pressure coefficients.

The equation adjoint to the integrodifferential form of the Boltzmann equation can be expressed as

(55)−Ω⋅▿φ*(r,E,Ω)+Σ(r,E,Ω)φ*(r,E,Ω)=∫E’dE’∫Ω’dΩ’K(r,E;E’,Ω;Ω’)φ*(r,E’,Ω’)+ΣR(r,E),

where ϕ^* is the adjoint function and Σ_R is the response function which serves as the source for the adjoint problem. The Boltzmann equation is often referred to as the “forward” equation and the adjoint as the “backward” equation. It can be shown that by combining the two equations (forward and adjoint) and integrating over all space. a response of interest, R, can be calculated using either

(56)R=∫rdr∫EdE∫ΩdΩ∑R(r,E)φr,E,Ω

or

(57)R=∫rdr∫EdE∫ΩdΩφ*(r,E,Ω)Q(r,E,Ω).