0X2&7683/93 $6.00+ .lO Pergamon Press Ltd
hr. J. Solids StrucruresVol. 30, No. 9, pp. IISI1212, 1993 Printed in Great Britain
A SYSTEMATIC APPROACH TO SHAPE SENSITIVITY ANALYSIS!DANIEL A. TORTORELLI Department of Mechanical and Industrial Engineering and Department of Theoretical and Applied Mechanics, University of Illinois, Urbana, IL 61801, U.S.A.
and ZIXIAN WANG Department of Theoretical and Applied Mechanics, University of Illinois, Urbana, IL 61801, U.S.A. (Received 16 March 1992; in revisedform 21 September 1992)
AbstractA fully detailed development of the domain parameterization method is presented for shape sensitivity analysis. It is shown that equivalent results are obtained from the material derivative method. In fact, the material derivative method may be viewed as a special case of the domain parameterization method which occurs when the reference configuration coincides with the body configuration. The method is illustrated for the Laplace problem in which explicit shape sensitivities are derived by the adjoint and direct differentiation methods. Both finite element and boundary element applications are discussed. The similarities between this approach and the isoparametric finite/boundary element method are transparent. In the finite element approach, it is shown that the sensitivity integrals may be transformed to the boundary (as is commonly done in the material derivative method) for the adjoint method, however, this does not seem possible for the direct differentiation method. Finally, in the boundary element approach, the sensitivities do not require the differentiation of the fundamental solutions.
1. INTRODUCTION
In “conventional” design sensitivity analyses, derivatives of response measures are evaluated with respect to design parameters which describe the physical and material properties of the domain and the loads and boundary conditions. For example, the derivative of the temperature at a point is determined with respect to the thermal conductivity of an isotropic homogeneous body. To derive these sensitivities, the variational equations which govern the system are differentiated with respect to the design parameters. Similarly, in a “shape” sensitivity analysis, the response derivatives are evaluated with respect to shape parameters which describe the geometry of the domain. Difficulties are encountered here when the governing variational equations are differentiated with respect to the shape parameters. In effect, we are differentiating integrals with respect to their region of integration. The problem is analogous to that encountered in Eulerian (spatial) motion descriptions where a control volume is introduced. In the following, a systematic approach for shape sensitivity analysis is presented. The foundation of the approach is domain parameterization in which a reference configuration is introduced. All quantities which are defined over the body configuration are transformed to field quantities defined on the fixed reference configuration. This approach is commonly employed in finite elasticity to transform the governing equations to the undeformed configuration and in the isoparametric finite/boundary element methods to transform the governing equations to the parent element. To obtain the sensitivities, the mapping which relates the reference and body configurations is varied. The adjoint and direct differentiation approaches for sensitivity analysis may then be applied to derive the explicit sensitivities. Both of these approaches are utilized to derive shape sensitivities of the Laplace problem. It is shown that the material derivative method for shape sensitivity analysis may be viewed as a special case of the domain parameterization method which arises when the reference and body configurations coincide. Both finite element and boundary element analysis applications are discussed ; and in the finite element technique, it is shown that for the 7 Presented at the AIAA 31st Aerospace Sciences Meeting and Exhibit, Reno, NV, 1993, 1181
1182
D. A. TORTORELLIand Z. WANG
adjoint derivation, the domain integrals may be transformed to the boundary as is done in the so called “boundary” approach of the material derivative method. However, the ability to convert domain integrals to boundary integrals in the direct approach does not seem plausible. The boundary element sensitivity expressions are naturally derived and do not require the differentiation of the fundamental solutions. Domain paramete~zations~ there referred to as the image of a fixed domain, were briefly discussed by Cea (1981) as early as 1981. However, the method did not receive much attention, as the material derivative (or speed method), was primarily used to derive shape sensitivities throughout the eighties [for example see Cea (1981a,b), Zolesio (1981), Choi and Huag (1983), Dems and Mroz (1984), Haug et al. (1986), Arora et al. (1991) and Petryk and Mroz (1986)]. In 1987 Haber (1987) used the domain parameterization method and a mixed form of the mutual energy principle to evaluate explicit adjoint shape sensitivities for a linear elastostatic continuum. Since then, numerous other applications of the domain parameterization method for shape sensitivity analysis have appeared [cf. Cardoso and Arora (1988), Phelan et al. (1989), Tortorelli and Haber (1989), Tortorelli et al. (1990), Tortorelli et al. (1990, 1991), Tortorelli et al. (1991), Arora et al. (1991) and Phelan and Haber (1989)j. Comparisons between these methods are denoted in Tortorelli et ai. (1991) and Arora et al. (1991). All of these studies utilized the finite element method to evaluate the system performance and the sensitivities. Shape sensitivities have also been derived using the boundary element method, via both the adjoint variable approach (Merit, 1988 ; Park and Yoo, 1988; Choi and Kwak, 1988a,b; Aithal and Saigal, 1990; Mota Soares et al., 1984) and the direct differentiation method (Barone and Yang, 1988, 1989 ; Zhang and Mukherjee, 1991; Zhao, 1991). These boundary element analyses relied on either the material derivative technique or differentiation of the discretized equations to evaluate the sensitivities. In the following, a detailed formulation of the domain parameterization method is presented. First, a review of some necessary continuum mechanics concepts and results are noted, e.g. the body and reference configurations, norm, transpose, tensor product, trace, derivative, gradient, divergence, Laplacian, divergence theorem, composite map, change of variable theorem, Jacobian and its derivatives, etc. Most of the material is self contained, although we do cite the texts by Gurtin (1981) and Ciarlet (1988) for derivations of some well known relations. At this point, we proceed to transform quantities from the body to the reference configuration and take their variations with respect to design changes. Next, we derive explicit shape sensitivities for the Laplace problem and discuss finite element and boundary element approaches for evaluating the sensitivities. The direct differentiation and adjoint methods are used to derive the sensitivities. Then we compare the results of the domain parameterization and the material derivative methods. Finally an example problem is provided.
2. CONTINUUM
MECHANICS
REVIEW
An attempt has been made to use consistent notation throughout this article. Direct notation is employed with uppercase bold faced Latin symbols denoting secondorder tensor fields, lower case bold faced Latin symbols denoting vector fields, and lower case Greek symbols denoting scalar fields, with the exceptions of the design variation vector cp and reference map function x. Lower case and upper case openface symbols denote points in the body configuration b and reference configuration 3, respectively. Calligraphic letters represent sets. The inner product, trace, tensor product, transpose, inverse, inverse transpose, determinant, norm, gradient and divergence operators are denoted by *, tr, @, T,  ‘, T, det, I/ 11,V and div. The x delineates the set product and the caretdenotes referential held quantities, defined through the composition 0. Lin (Js, @) denotes the space of linear operators from & to
[email protected] %?(&,@) denotes the space of n times continuously differentiable functions from d to
[email protected] and lower case sans serif symbols are used to denote vectors in the reference and body configurations, respectively. In our examples we use indicial notation where subscripts denote components with respect to a Cartesian coordinate
A systematic approach to shape sensitivity analysis
1183
system with basis vectors e and E in b and B, respectively. Finally, D denotes differentiation with respect to the subscripted argument. Where possible, the following results have been specialized to hasten the derivations in the following sections. 2.1. The body and reference configurations and reference map In the ensuing analysis, we consider a body comprised of particles. We identify the particIes by the places which they occupy in two configurations, the body conjiguration and the reference conjiguration. The body configuration is denoted by b c 6” where B” is an ndimensional Euclidean point space with the associated vector space V’. The body configuration is assumed to be a regular region in b” with outward unit normal vector n E T’ on the bounding surface ab. Points in the body configuration are denoted by x E b and are defined by their coordinates relative to the rectangular Cartesian coordinate system (o; e,}, i = 1, n. In general, lower case symbols are reserved for quantities associated with the body configuration. The reference configuration is denoted by B c 8” where gN is an Ndimensional Euclidean point space with the associated vector space V”‘. The reference configuration is also assumed to be a regular region with outward unit normal vector N E VN on the bounding surface aB. Points in the reference configuration are denoted by X E B and are defined by their components relative to the rectangular Cartesian coordinate system (0 ; Ei}, i = 1, N. Uppercase symbols are reserved for quantities associated with the reference configuration. We require that 1 < N < n Q 3. We define a referential function x : B + 6” which maps the reference configuration onto the body configuration. This mapping may be thought of as the deformation which takes the undeformed reference configuration into the deformed configuration in elasticity. The present approach is also consistent with isoparametric finite and boundary element formulations where the reference configuration is the parent element. We insist that the mapping x have the properties of a proper deformation, i.e. it is invertible, differentiable, and orientation preserving. The Laplace problem which we study in the sequel is defined over the body configuration. The body configuration, in turn, is defined on the reference configuration through the mapping x. Thus, the results of our analysis clearly depend on the mapping x. The objective of the shape sensitivity analysis is to evaluate this dependence. It is imperative, then, to quantify the referential map’s dependence on the design, to wit we redefine the map as x0 : B x gM + 6”. This definition is analogous to that of a motion in elasticity where ~~~~ replaces time. SYMis a subset of the Mdimensional real number space B?’ (i.e. gM c W”), and its elements are the design vectors @ which contain the parameters that define the geometry. Note that @ could be replaced by an infinite dimensional function space, however, we have selected a finite dimensional set to be consistent with numerical discretization schemes such as the finite and boundary element methods. Now the body configuration is expressed b. =
[email protected](B) to note its dependence on the design @. Other quantities defined over bg e.g. x and n, are not subscripted with @, however we emphasize that they too are functions of the design @. For notational convenience, we define x(X, @) =
[email protected](X) and write x = x(X, 0). Then ba = x(B, @) and we define the design trajectory as the set of pairs Y = {(x, a) : x E
[email protected], @~9~). Again, analogies may be drawn between the definition of the design trajectory and the trajectory used in continuum mechanics. Here the design @ replaces time. In essence, the design trajectory is the set of body configurations which is spanned by the design space. Two examples of the mapping are illustrated in Figs 1 and 2. In the first example (Cea, 1981), N = n = 2 and the design @ is an element of an infinite dimensional function space 9. The reference domain is the set of ordered pairs B = {(X1, X,) ~9’ : 0 6 X, 6 1 and OdX2< 1);andthemapx: Bx9 + g2 is subsequently defined by the components Ml =
x,,
x2
[email protected]@;(l),
=
(1)
1184
D. A. TORTORELLIand Z. WANG
1 Fig. 1. Mapping from 2dimensions
L to 2dimensions.
where the position vectors of x and X are given by X:i’=,wiei and Xi’=, XiEi, respectively. In the second example the isoparametric map of a threedimensional continuum boundary element is used. The reference domain, or parent element, is the set of ordered pairs 1
(2) where the position vectors of x and X are given respectively by x = C/=, w,ei and X = Cf= 1XiEi. In the above, X: denotes the ith component of the Jth node and Y’: B + W is the Zth element shape (interpolation) function. Here N = 2 while n = 3 and the design parameter vector is the element node coordinate vector @ = (w ], xi, xi, xi, xi, w:, . . . , WY,xi, xi) E 924. This latter finite dimensional example will be referenced occasionally throughout the article. 2.2. The norm, tensor product, transpose and trace The norm, tensor product, transpose and trace yield unique quantities defined through their respective operations. The trace, transpose and inner product also possess several properties which we draw on for the ensuing sensitivity analyses. The definitions are given in direct notation. Where deemed appropriate complementary indicial notation definitions are supplied ; we recall our restriction to orthonormal basis vectors. Given the vector a E V’, then the norm ]]a]]EW is defined as [cf. Gurtin (1981)]
Fig. 2. Mapping from 2dimensions
to 3dimensions.
A systematic approach to shape sensitivity analysis
Ml = (
llall = (a*a)‘12,
!,aiui>“2>
1185
(3)
where * denotes the vector inner product operator [cf. Gurtin (I 98 l)], The tensor product of the vectors a, b E P defines the tensor a @ b E Lin (“lm, ly”) such that ([email protected] b)c = (bc)a,
(a @ b)ij = aibj,
(4)
for all c E Y” [cf. Gurtin (198 I)]. Given the vectors a, be V’ and the secondorder tensor A E Lin (V”, Y”), then the ?~ff~sp~~eATc Lin (V”, P) is defined via [cf. Gurtin (1981), p. 3] aAb = ATa*b.
(5)
Given the tensor A E Lin (V”, V”“) and the vectors a, be V’, then the trace tr (A) E 9l is defined through the relation [cf. Gurtin (1981)] tr([email protected])
= ab,
tr(A) = i Ai<
(6)
i= 1
The inner product of the tensors A, BE Lin (V”, V’) is defined through
AB = tr(ATB),
A*B = i i [email protected],, in 1 f= t
(7)
[cf. Gurtin (198 l)]. Given the tensors A, B, C E Lin (V, V’), it may be shown that the transpose, trace and inner product yield the following relations which appear in Gurtin (198 1) : tr(A) = A1, tr (A) = tr (AT), tr (AB) = tr (BA), A * (BC) = (BTA)  C = (ACT)  B, (AB)T = BTAT,
(8)
where I is the identity tensor [cf. Gurtin (198 I)]. 2.3, The derivative, gradient, variation, divergence, Laplacian and partial derivative In our derivations we require the derivative, second derivative, partial derivative, variation, gradient, divergence and Laplacian of scalar, vector and tensor fields. For completeness we choose to review these definitions. Throughout this section A, a and c1denote tensor, vector and scalar valued functions on [email protected]; b, P and 6 denote tensor, vector and scalar valued functions on B x ~23~; and the notation of the previous section is enforced. The following represent elements of their noted sets u, v E Vn, U E TN, and o, E 9”. In all cases we assume that the functions are su~ciently smooth to perform the indicated operations. The derivatives are defined on the open interiors of the sets, denoted by (“) ; we assume that extensions to the set boundaries a( ) are possible. Given the function 01: bg + 9, the derivatiue at x is the linear operator DC+) E Lin (V, W) such that
1186
D.
A. TORT~RELU and
Z.
WANG
Da(ru)u = a(x + u) a(w) + o(u),
(9)
where o(u) tends to zero faster than u. Thus, Da : b& + Lin (Y”, B?). In essence, the derivative is a linear map which approximates the difference a(x +u) CL(X) for small j/uII. If the derivative exists and is continuous for all XE&, then we say that a belongs to the set of continuously differentiable functions, i.e. a E U’(&, 9). [See Gurtin (1981) for additional discussion.] Similar definitions apply to vector and tensor fields. If a maps into the vector space Y”, i.e. if a E%?‘(b&V”), then the derivative defines a secondorder tensor, i.e. Da(x) E Lin (Y’, Y”“). Where we recall that Y’ is the vector field associated with be. Considering the derivative Da : &, + Lin (Ike”,a) as a function, we define the second ~r~v~~ive (Ciarlet, 1988) D 2a : 6% + Lin (Y”, Lin (Y”, a)) through D “a(x) = DDa(x).
(10)
If the second derivative exists and is continuous for all x E b;Bthen we say that CIbelongs to the set of twice continuously differentiable functions, i.e. aE%2(b& W). It may be shown that D*a(w) : Y” x V” P ~2% is a bilinear operator such that (Ciarlet, 1988) D *a(x)(u, v) = (DDa(w)u)v.
(11)
Similar definitions apply to vector and tensor fields. So the second derivative approximates the difference of the derivative, i.e. Dafw +v)u Da(x)u for small ]lul] and ]lvll . We also note the symmetry of the second derivative (Ciarlet, 1988) D*a(x)(u,v)
= D*a(x)(v,u).
(12)
If a E %P(&,, B) then a has the v~ri~~io~ or dire~~~~n~~derivative 6a which is defined through (Ciarlet, 1988) Ga(w; u) = Da(w)u,
(13)
for all u E OP. Thus, for the differentiable function a, the variation Ga(x ; u) E W equals the derivative acting on the increment u.t For every linear scalar valued function, i.e. linear functional #Olin (Y”,W), the representation theorem for linear forms ensures the existence of a unique vector a E T” such that t&u) = a * u for any vector u E Y” (Gurtin, 1981). Armed with this theorem and noting that Da(x) E Lin (P, a), we define the gradient Va : 6; + V” through Va(x) *u = Da(x)u
(14)
[see Gurtin (1981)] The divergence div a : b& ) 92 of the vector field a~%‘(&& +““) is defined as (Gurtin, 1981) diva(x)
= tr Da(x).
(15)
The divergence div A : bg , Y” of a differentiable tensor field A E V’(b&, Lin (YR, Y)) is defined through (Gurtin, 1981)
7 We note that the variation is actually defined through the limit process, is not necessarily a linear operator, and may exist even if the derivative does not exist (Ciarlet, 1988). In the forthcoming analysis we limit ourselves to differentiable functions, so the above, over restrictive, definition suffices.
I187
A systematic approach to shape sensitivity analysis
A(x) *u = div (AT(x)u)
(16)
for all constant vectors u E Y”. The Laplacian Au : bi + W of a twice differentiable scalar field a E V2(b& W) is defined as (Gurtin, 1981) Aa
= div Vu(x).
(17)
Next we define the partial derivative and its properties for the referential fields. As our referential functions are defined over the set product of the domain B and the design field JP”, i.e. B x @‘, we require the use of the partial derivative. Formally, given the differentiable function aiE V’(B’ x PM, a), we define the partial derivatives with respect to X and CD(the first and second arguments) respectively, by (Ciarlet, 1988) D,di(X,@)U = Doi(X,cD)(U,O), [email protected],
wcp
=
wx,
W(O,
cp>,
(18)
hence, D,8 : B” x @’ + Lin (P’, a) and D20i: B” x .GSM + Lin (g”, W) (Ciarlet, 1988). The linearity of the derivative then yields (Ciarlet, 1988) D&X, @)(U, cp) = Dloi(X, @)U+D&(X,
@)cp.
(19)
Similar definitions apply to vector and tensor valued functions. We define the function D, D2$ : B” x WM + Lin (V”, Lin (g”, W)) in a similar fashion to the second derivative, and it has the property that
DAW,
@)cp)U,
W(U cp)= (DAW,
(20)
so that D, D2&(X, @) approximates the difference D2&(X +U, O)cp  D&X, @)cp for small [lull and (Icp(I. Finally, like the second derivative, DID23 is bilinear and exhibits the symmetry
D,DzV& @NJ, rp) = DzD,WC W(cp,v).
(21)
Again, similar definitions apply to vector and tensor valued functions. The following partial variation and partial gradient are used in the sequel 82W%
@
; cp)
=
D245W.,
Wcp,
V,oi(X,@)*U = D,&(X,@)U. For the variation,
(22)
similar definitions apply to vector and tensor valued functions. The
partial divergences and partial Laplacian are defined through
div, ii(X, Q) = tr D 1ii(X, a), div, .h(X, CD)* u = div 1(A’(.&‘, @)u), A,di(X,@) = div,V,B(X,@), where oi~~‘(B”~9’~,~),
PE%‘~(B”xP~,P),
(23)
and AE%?‘(B~xS?‘~,L~~(Y~,Y”)).
2.4. Elementary rules of dflerentiation In this section, we present some elementary rules of differentiation. agreements of the previous sections remain in effect.
All notational
D. A. TORTORELLI and Z. WANG
1188
Fig. 3. Composite
function.
Consider the differentiable maps CI,BE%‘@,, W), a, beGf’(b& ^Y”), and A, BE+?’ (B& Lin (V’, Ilr”)) then it may be shown by the product rule of differentiation (Gurtin, 1981) that D(cf(x)a(x))u
= Da(w)ua(x) +a(x)Da(w)u,
D (a(w)A(x))u = Dcr(w)uA(x) + a(x)DA(x)u, D(A(w)a(x))u
= DA(w)ua(w) +A(x)Da(x)u,
D(A(x)B(x))u
= DA(x)uB(x)
+A(x)DB(x)u.
(24)
Using the above definitions for CI,a, b and A, the following relations may be verified (Gurtin, 1981) : div (cc(x)a(x)) = a(x) div a(x) + a(x) *Vor(x), div (AT(m
div (a(w
= A(w) Da(m) + a(x) * div A(x), = a(x) div A(x) +A(x)Va(x),
V(a(x) * b(x)) = (Db(m))Ta(x) + (Da(x))Tb(x), V(~(Mw)
= B(ww)
+mVlw.
(25)
2.5. The composition, chain rule and Jacobian
To transform the governing equations from the body to the reference configuration we use composite functions. Temporarily excluding the design from the domain, we consider the composite function (see Fig. 3) di: B + 9 defined by oi = a 0 x so that WQ = 4x(W),
(26)
where u~W’(b~, 9%‘)and x E%‘(B, 8”) are both differentiable. Then by the chain rule of dljerentiation (Gurtin, 198 1) for any U E VN, Doi(X)U = &(x(X))
0 Dx(X)U
= Mx(W)Px(WU),
(27)
where X E B”. Similar results are obtained for vector and tensor valued functions on b,.
1189
A systematic approach to shape sensitivity analysis
Fig. 4. Composite
function.
Now
customize the for our 8: Bx~~+%?, on the a: Y x = x(X, Cp) the composite like function (see Fig. 4)
where we the design’s and x: BxC~~+~‘. Then for
oi.(X,@) = a(x(X Q), @)
(28)
is differentiable. To evaluate the derivatives, we define the functions g, : B+b”,g, b”, h,: B+b”, and hZ: P”j+&“: Sl(W
= x(X
:CP+
@),
92(@) = x(X, 9, h,(X)
= Q,
hd’W = @,
and their derivatives
where we use eqn (19) and note the independence between X and Q to derive the above. Using the above relations and eqns (19) (22) and (27) we obtain the following relations :
[email protected],WJ
WC%
= D,a(s,(X),h,(X))(Dg,(X)U)+D*a(g,(X),h,(X))(Dhl(%)U) = [email protected],
W, ‘WD,[email protected],
WJ)
= VdxW,
9,
WJ,
@)cp = Dl~(g2W)~ = DdxW, = VdxW,@L
‘J9 *D,xW, [email protected])
9,
[email protected],
+D&d%
WWWWW)
Wcp> [email protected],@),
@)rp
~).D2x(~,~)cp+V2a(x(X,~),~).cp.
(29)
Again, drawing analogies from continuum mechanics, we note that D&X, cD)cpresembles the material derivative for spatial fields, where time and @D are interchanged (Gurtin, 1981). In the sequel, the caretalways indicates composite functions defined as above [eqn (28)J. Since all of our composite functions use x, for convenience, we define the Jacobian J: Bxg”+Lin(YN,Vn) as
1190
D. A. TORTORELLI and Z. WANG
J(%@) = D,x([email protected])
(30)
and the determinant of the Jacobian as J(X, a) = det J(X, Q),
(31)
where det denotes the determinant (Guttin, 1981). The determinant only makes sense if A’ = n. However, this does not present a problem as we may always redefine the map so that N = n (DoCarmo, 1976). For our example problem [eqn (2)] one possible modification yields x : (B x {92>) x 924 + 8” :
Y’(X,,~,wl,
Xl = i I=
x2 =
X3=
I
; Y'(X,,X,)& I=,
(32)
~y'(x,,%*)&.+X3,
I= I
where X3 is equated to zero and (X,, X2) E B c 6%‘. Likewise, any other functions defined on B, are redefined on (Bx {W}), for the example illustrated in Fig. 2, oi(X , , % 2, X 3) = oi(X , , X *). This modified example map x will be used in the sequel. Note that problems for which n # N may also be solved using manifold theory as described in Rousselet (199 1). For our example [eqn (32)],
,~,v,w~,,x,)d
[J(XWl =
; v,‘P(x,,x2)x:
I= I
2 I=
2
v*~“w,,x*)x:
0
V*~'(~,,X2)~12
0
I
I= I
3
(33)
where cp = (~x~,~x:,SX:,~X:,BX:,BX:, perturbations of the node coordinates.
. . . ,6x:, 6x;, Sx$ for each element and represents
2.6. The divergence theorem and change of variable theorem Here we review some integral theorems, namely the divergence theorem which is used in variational formulations and the change of variable theorem which is used to transform integrals from the body to the reference configuration. Consider the differentiable functions c(: .F + 93, a : F + V’, and A: F + Lin (Vn, Vn), then by the divergence theorem (Gurtin, 1981) :
1191
A systematic approach to shape sensitivity analysis
J
a(x, cD)n(x, Cp)da, =
db,
[email protected], @) dr,, s 6.v
J J
div, a(w, @) dv,, a(w, 0) * n(x, @) da, = sa&l bo
J
and
div, A(x) dv,, A(x)n(x) da, = abo b,
(34)
where n(x, @) is the normal vector at the point x to the body surface &, ; and da, and du, represent differential elements in ab,,, and bo. Given x: B x W’ ) P and the above definition for u, then by the change of variable theorem (Gurtin, 1981) :
J
[email protected], a) dv, = x(&W
@I, @)J(X, @) J J&(X, dvw =
B
oi(X, @)J(X, @) dvw,
(35)
B
where dvX is an elemental volume of B. Recall from Section 2.1 that b. = x(B, @) is the image of B under x for the design @. Thus, this equation equates integrals in b. to B. Similarly,
where dax is an elemental volume of 8B and K(X, @) = J(X, B) 11 JmT(X, @)N(X) 11(cf. Section 3.4). Finally note that the previous two relations are applicable to vector and tensor valued functions.
2.7. Special derivatives Eventually, we require the design variations of the Jacobian J, its inverse J ‘, their transposes, and its determinant J. Some of these results require the condition that N = n, which we assume holds henceforth. However, this is not a restriction, as we may always modify the map so that N = n as exemplified above [eqn (32)]. Of course, for the inverse J ’ to exist we assume that J is nonsingular, i.e. its determinant is nonzero and/or the columns of the matrix J are linearly independent, this condition is ensured if the mapping x is invertible as previously discussed. For J : B x &’ + Lin (V”‘, V) we may readily verify that its design variation is given by the spatial gradient of the design variation of x, i.e. h2J = D,~*x. Indeed, by the symmetry of the partial derivative, the arbitrariness U, and the definitions of the variation and Jacobian [cf. eqns (20), (22,), (30)] we arrive at
(37) For our example, a direct calculation on eqn (33) yields us
30:9D
1192
D. A. TORT~RELLIand Z. WANG
bJW,@; cp)=0,[email protected],@;
cp)
(38)
By invoking the chain rule [eqn (27)], the linearity of the trace, and the above result [eqn (37)], it is easily verified that
6,J’(X, ‘J’; cp)= W&W,
@; cp).
(39)
This result of differentiation of the transpose also appears in Gurtin (198 1). To determine the design variation of J ’ we use the relation J ‘J = I where  ’ denotes the inverse operator (Gurtin, 1981). Differentiation of this relation by the product rule, the chain rule, the definition of the variation, the previous result [eqns (24,), (27,), (22,) and (37)] the invertibility of J, and the fact that I is constant gives J ‘(X, @)J(X, a,) = I, D,J‘(X,cD)cpJ(X,~)+J‘(X,~)D,J(X,Q)cp
= 0,
D&‘(X,cD)cp
= J‘(X,cD)D,J(X,~)cpJ‘(X,~),
which leads to iS2J‘(X,@; cp) = J‘(X,@)D,&x(X,@;
cp)J‘(X,0).
A similar result appears in Gurtin (1981) for the differentiation Likewise, using the identity JmTJT = I we‘obtain &J=(X,@;
cp) = JT(X,Q)D,62~T(X,Q);
(40)
of tensor inverses.
(p)JT(X,Q),
(41)
where the notation JT = (J‘)T is enforced (Gurtin, 1981). To determine the design variation of J we first consider the map 71: Lin (V”, VN)+ such that rr3 = det A.
w
(42)
Then, it may be shown [see Gurtin (1981) or see the alternative approach provided in the Appendix] that Dn3(A)W = det A tr (WA ‘),
(43)
where A, W Olin (VN, VN). Using the chain rule [eqn (27,)] where we equate J(X, QD)= n3 0J(X, @) so that A = J(X, 0) and W = D2J(X, @)cpin the above, we obtain
= det J(X, a) tr (D*J(%, @)qJ
’ (X, a)),
which upon combining eqns (22 ,), (8 ,), (37) and (44), yields the result
(44)
1193
A systematic approach to shape sensitivity analysis WW,
0;
cp) = 4X
9) tr (D1bx(X
@; cp)J
= J(X,@)D,&x(X,O;
3. TRANSFORMATION
‘(X
W)
&JT(X,@).
TO THE REFERENCE
(45)
CONFIGURATION
In order to transform our problem to the reference domain we must convert gradients, the divergence, the Laplacian, normal vectors, and elemental areas from the body configuration to the reference domain. To accomplish this, we repeatedly apply the chain rule and Piola’s identity div, (JJ=) = 0 (Ciarlet, 1988). We also evaluate the design variation of some of these quantities. 3.1. Transformation of gradients Consider the composite value scalar map ai: B x QM + 9 defined by eqn (28) where a : .T + W and x : B x @” +6” are defined in the usual manner. Then from the chain rule and the definitions of the gradient, Jacobian and transpose [eqns (29,), (22,), (30) and (5)J Vloi(X,@)*U = V,a(x(X,@),@).J(X,@)U = J’K
WV,~(x(X,
@),a)
.U,
so that JT(X~)V,oi(X,cP)
= V,a(x(X,@),6),
(46)
where we used the invertibility of J and the arbitrariness of U. Next, we consider the vector value composite map B : B x @’ + Y” defined in an analogous manner to oi. Then by the chain rule, the definition of J, invertibility of J, and the arbitrariness of U [eqns (29 i), (30)] respectively, we obtain DiP(X,cP)U = D,a(x(X,@),@)J(X,@)U so that D,a(X,cD)J‘(X,(D)
= D,a(@[email protected]),@).
(47)
3.2. Transformation of the divergence Consider the equation div, a(x, @) = 0
for
(x, a) E Y,
(48)
where a : F + V’. For example, this represents the energy equation for a steadystate conduction system with no source where a is the heat flux vector and Q, represents the current design. We may write for any b& c [email protected]
s
div, a(x, @) dv, = 0,
bi
which after an application of the change of variable theorem [eqn (35)] becomes
(49)
1194
D. A. TORTORELLIand Z. WANG
div t a(x(X, a), (D)J(X, #) dv, =
s B’
s ~0’)
div, a(x, @) du, = 0,
(50)
where x(B’) = b&. Then by the localization theorem (Gurtin, 198l)t div, a(x(X, CD),4D)J(X, @) = 0
for
(X, @) E (B x 9”).
(50
Unfortunately, this equation is not expressed in terms of referential field quantities. To this end, we perform the following transformations :
where we use eqns (23 II, (47), (7), (25,), and Piola’s identity div, (J(X, #)JpT(X, @)) = 0 (see Appendix for details). Thus, we now have the equivalent equation defined solely by referential quantities divl (J(X, Qp)J ‘(X, @)ii(X, @)) = 0
for
(X, a) E (B x &Y”),
(531
where the fields are related through the composition a(x(X, @), a,) = [email protected], Q). 3.3. Trff~~~or~u~~o~of the Laplace equutio~ Consider the Laplace equation A,a(w,dD) = 0
for
(x,@)E$“,
(541
where tl : F + $8. Proceeding as before and using the definition of the Laplacian [eqn (23,)] we express for any b&, c ba
s
div, V,cl(w,Q)dv, VbP
= 0,
(55)
which after an application of the change of variable theorem [eqn (3511 and localization theorem becomes divr V,a(x(X,cD),Q))J(X,@)
= 0
for
(X,Q)E([email protected]‘).
(561
We now commence to transform this equation by sequentially applying eqns (23,), (8 ,), JTJmT= I, (84), (47), (25,), Piola’s identity div, (J(X, a))JT(X, @)) = 0, and (46) :
div, V,dxtX @,),WV,
W =
tr t~~V~atxt~,~)))~t~, @I
= JT(X, @)JT(X, Q) *D,V,&$K, =
@), (Ib)J(X, ai)
J(X,~)JTIX,QI)*D,[‘C7,a(xIX,QD),cD)l
= div, (J(X, @D)J‘(X, #D)V,(r(x(X, (D), @)) div,(J(X,dP)JT(X,Q))*V,a(~(%,@),@) = divx (~(~,~)J‘(~,#)JT(~,~)V~~(~,~~).
(57)
For clarity we used the subscripts 1xand X to denote partial differentiation with respect to fTo use the localization theorem we require div, al to be continuous and B’ to be open.
1195
A systematic approach to shape sensitivity anatysis
the first argument for functions on .F and B x @, respectively. Thus, we now express the equivalent expression for eqn (54) defined solely in terms of referential quantities : div, (J(X,~)J‘(X,QD)JT(X,~)V,oi(X,~))
= 0
where the fields are again related through the imposition
for
(X,@)E([email protected]‘),
(58)
a(@I, @), a) = c?i(X,clp).
3.4. Nanson’s relation Combining the results of Section 3.2 [eqns (50) and (52)] yields
J
J
div, a(w, ip) dv, = XfB’f
div, (~(~,~)J‘(~,~)~(~,~))d~~
(59)
B
which after an application of the divergence theorem [eqn (34,)] becomes
J
a(x, #)  n(x, @) da, = J(X, @)J ‘(X, @)[email protected], @) *N(X) dux. JC3B ax(w)
(60)
Now by the arbitrariness of a(X(X, @)) = ri(X, CD)and the definition of the transpose [eqn (5)] we are led to Nansoni relation :
n([email protected],@Ida,(x,ej= J(X, (D)JvT(X, at)N(X) da,.
(61)
This relation appears in Gurtin (1981). Taking the norm of both sides and using the normalized normal vector n, gives [cf. Ciarlet (1988)]
K(X, Ur) then is a metric that relates surface differentials between the reference and body configurations [cf. eqn (36)]. Upon combining the previous two results we obtain
ntx(X Qp),@l= JTtX, ~)N(~)/tl JTWWfW>
It
= J(X, @)J‘(X, @)N(X)/K(%, #),
(63)
which is an expression relating normal vectors between the two configurations. To derive the sensitivities we require the design variation of K. Multiple applications of eqns (3), (24), (27), (45), (41), symmet~ of the inner product, and eqns (62) and (63) give
wm,~;
p) =D,J(X,~P)~~~~J=(X,~)N(X)~I+~J(X,~)J~~X,~)N(X) [email protected],Wd’+JfW/Il J=tWWW> + ~~(~,~)D~J=(~,~)~N(~) . J‘PW)N(~),‘II
J=(X WWQ It
= JW,@)tr(D,[email protected]; J(X,UQJ‘(X, ItJ=(X,
II
cp)Jl(X,~))IIJT(X,cP)N(X)lI
@)N(X) *J=(X, QI)D,&x~(X, a; ‘p)JT(X, @)N(X)/
W&W II (64)
1196
D. A. TORT~RELLIand Z. WANG 4. SHAPE SENSITIVITIES
FOR LAPLACE’S
PROBLEM
In the following we transform Laplace’s problem to the reference configuration and draw analogies to the isoparametric map used in the finite and boundary element methods. We also introduce the concept of a response functional and derive shape sensitivities via the direct and adjoint methods. 4.1. Laplace’s equation and the isoparametric formulations Consider the boundaryvalue problem governed by Laplace’s equation for
(x, Cp)E 9,
a(~,@) = olp(x,@)
for
(w,*) Ea$;,
V,cc(x, Qp)* n(x, @) = qP(x, a,)
for
(w, Qp)E&Y&
Ahlcr(x,@) = 0
(65)
whereaFZ = ((~,a); x~db*,, @~9~), aFZ = ((~,a); x~db~,,@&“‘). This could represent a steadystate conduction system with no source. Then, we have the prescribed tem~rature ap on 8% and the prescribed flux qp on &FZ. As previously mentioned, to derive shape sensitivities we transform the above to the reference configuration, which from eqns (58), (46) and (63) becomes div, (J(X,@)J‘(X,@)JT(X,@)V,B(X,@))
= 0
d(X, @) = 3(X, to)
for
(X,@)EBXW’
for
(W, cll) o 899, x gM
J‘(X,~P)V,~(X,@)XF(X,~)J~(X,Q)N(X) = @(X, @)K(X, CD) for
(X, @) E dB, x P,
(66)
where dP(X, CD)= a”(x(X, @), CD)and &‘(W,4) = #Q(W, co>,iD) are composite functions. We now solve the above equations for the referential response a and VI& Once these quantities are determined a and Va may be evaluated via eqns (28) and (46), i.e. a(X(X,@),@) = oi(X,@)and V,a(x(X,@),QD) = JT(X,Q))V,di(%,@). A comment here, regarding the functions oip and 4” seems appropriate. These are functions of a reduced domain, as they are defined on the surface 8B and not within the body B, and therefore have no normal derivatives. This has r~fi~tions when we take the gradient of BP in later sections, e.g. eqn (76,). Note that oi: B + W whereas dp : 8B + S?, thus their spatial gradients which are defined on different domains are not equal, V,ai # V ,oi{. Similar remarks regarding Qpalso apply. The above equation may be solved by using the isoparametric finite and boundary element methods. Indeed, if we multiply eqn (66,) by an arbitrary function fi: B x gM + 9 and integrate over B we obtain
s B
idiv,
(JJ‘JTVIB)dvx
= 0,
(67)
where we suppressed the ar~ents for conciseness; a practice which we continue henceforth. Now application of eqns (25 ,) and (34,) to the above gives
s
8BfiJ‘JTVIOi*N
da,IB
V,fi*JJ‘JTV,Oidv,
= 0.
(68)
To derive the finite element equation, we restrict oi to the space of square integrable functions which satisfies eqn (66,); and we restrict 1 to the space of square integrable functions which equal zero on 8B,. Then we determine the oi from the admissable space of functions which satisfies
A systematic approach to shape sensitivity analysis
1197
I
(69)
ae &*Kdax 2
=
BJTV,6*JTV,&Id~x s
for all admissible f We derive the above from the definition of the transpose [eqn (5)], eqn (66& and the restriction that fi = 0 on aB,. The finite e’lement method systematically develops the basis functions for the Galerkin method (Becker et al., 1981) which is used to obtain the approximate solution. We use the left and righthand sides of the above to form, in a piecewise manner, the load vector and symmetric stiffness matrix, respectively. In this piecewise approach, the map x and Jacobian J are defined locally over elements (subregions) whose images collectively approximate b*, In this way, the (typically numerical) integrations are performed individually over the elements and then summed. To be precise, in the isoparametric finite element method, we approximate oi and 1 with the same interpolation functions which define x. For our example [eqn (32)], oi = Xc,“=. , Y’oi’and V ,d = C,“= , V ,Y’oi’ where the unknown (to be determined) parameter oi’ is the value of & at the node with coordinates x’. These results are consistent with those appearing in Becker et al. (1981). To derive the boundary element equation we start with eqn (68) and again use eqns (25,), (34,), (5), (66*) and (66,) to obtain
1
oiPJTV,iS JJTN
dax +
aBl
l
s 84
lgpKdax
= 
s
&JTV,~*JJ*N
daX
a&
A7~2~ JJTN da%
s dB1
sB
oidivi (JJ‘JPTVifi)dvX.
(70)
Next we choose x(X, a) = n&(X, @)) where 1: b,,, + 8 is the design independent fundamental solution [cf. Carey and Oden (1983)] so that the volume integral becomes f$(Y, @) for V E f3B. Hence, the above reduces to
1
oiPJTVIx= JJeTN dax +
8Bl

@*Kdax
s as2
s
IJJ~v~~*JJTN aB1
=
daX +
&J‘V,&JJTN s aBz fai(V,@)
dax for
VEaB.
(71)
The left and righthand sides of the above are used to form the load vector and nonsymmetric stiffness matrix, respectively. These results are consistent with those appearing in Carey and Oden (1983). Again, as in the finite elment method, the map x and Jacobian J are defined piecewise over the elements (subsurfaces). 4.2. The response functional To gage the performance of a system we may evaluate a general response functional of the form
G(a)
=
s
g(oz,V,a*n,x)da,. 8b.a
(72)
Clearly, we expect G {#) to vary as
1198
D. A. TORTORELLI and 2. WANG
The above functional is defined in terms of the body’s response quantities. However, we are not directly evaluating a and V,a during our analysis (so that we may derive the shape sensitivities and be consistent with the isoparametric formulations). Instead, we invoke eqns (28), (36), (46) and (63) and the definition of x to write the equivalent functional over the reference configuration
G(Q) =
s aB
@(oi,JmTV,oi* JJTN/K,X)KdaX,
(73)
where 8 is defined on the set product of the V’ continuously differentiable functions. S~cifically,~: ~‘(Bx~~,~)x~i(Bx~~,W)x~‘(Bx~M,~~)~~. The objective of a sensitivity analysis is to determine the design variation of G, i.e. we wish to evaluate GG(cD; p) which upon noting the independence of X from Q, and applying eqns (13), (19), (20), (22,), (22,) (24), (27,) and (5) becomes
where the identities J* = J*T = JJ‘JmT
and r&J* = c?,J*~ = 62JJ1JT+JS2J1JT+ Note that since x, its derivatives, and 3 are explicitly defined, then all of the terms in eqn (74) are readily evaluated with the exception of [email protected] and V,&&. These terms are implicitly defined by the design # through the boundary value problem [eqn (66)]. In the following sensitivity analyses, we will show that it is possible to explicitly evaluate 6G. JJJ ‘C&J‘i’are defined for convenience.
In the direct differentiation method we ma&ate && and V&Z. Once these quantities are determined, we may determine 6G directly from eqn (74) or the response variation via eqn (76). We compute the referential response variations by differentiating eqn (66) with respect to the design, whereupon [using eqns (20), (24) and (27.2)] div, (~~J*V,~+J*V~~~~)
= 0
Bz&= JTV,oiP*6~~+S#’ V,&0i.J*N+V,Oi.&J*N
= JTV,~P.62~K+S2qPK+(iP62K
for
(X, #) E B x ?F,
for
(X, aP) E aB, x ZP,
for
(X,@)EC%X~~, (75)
where it is henceforth understood that all occurrences of CC and qp are evaluated at (w,@) = (x(X,@),@). Now we solve the above for &,oi and V,@ and substitute these quantities into eqn (74) to obtain the explicit sensitivity. To solve this equation we may use either the finite or boundary element methods. Once we know &A!,we may evaluate &a. Indeed by using the chain rule [eqn (29,)] we obtain the relation
so that
where we used eqns (22,), (22,) and (46). The above represents the design variation of a
1199
A systematic approach to shape sensitivity analysis
at fixed position x in the body configuration. A similar relation for V&a may also be developed. To obtain the finite and boundary element equations, we proceed as before. First we multiply eqn (75 ,) by j, integrate over B, and apply eqns (2.5,) and (34,)
s aB
IJJ*V,~~~.N dax +
V,;i+ J*v,S~~~ dvx,
(77)
s LI
where the explicit and implicit quantities appear on the left and righthand sides, respectively. Now, for the finite element equation, we restrict &oi to the space of square integrable functions which satisfy eqn (75,) ; and we restrict 1 to the space of square integrable functions which equal zero on dB,. Then, we are to determine the && from the space of admissible functions which satisfies

J
V&&J*V,&
dvX =
J
J=Vi& JTVIG&Jdv,
(78)
B
B
for all admissible x. We derived the above from eqn (77) by using eqns (5) and (75& the symmetry of J*, and the restriction that ;i = 0 on c?B,. We could have also derived this by direct differentiation of eqn (69). Note that this equation resembles eqn (69) in that the known terms are on the lefthand side and the unknown quantities are on the riothand side. In the finite element analysis, the lefthand side forms the load vector while the righthand side forms the stiffness matrix. Further, note that the righthand sides of eqns (69) and (78) are similar. In fact, they both define the same finite element stiffness matrix. This implies that the decomposed stiffness matrix which is used to evaluate the primal response oi can also be used to evaluate the pseudo response 6& Thus, the pseudo response is e~ciently computed by merely perfo~ing a pseudo load vector assembly correspon~ng to the lefthand side of eqn (78) and a back substitution into the existing decomposed stiffness matrix. Of course, we compute Vi&a from &(7iin the usual finite element manner. Once we have determined these quantities we may evaluate SG via eqn (74). For the boundary element approach, we begin with eqn (77), apply eqns (251) and (34,) to the righthand side and use eqns (75,) and (75,) to obtain (after isolating the implicit and explicit vacations)
J
=:

J
~JTJ*VISzoi* N dax + S#J*V,x* 8% ae,
s B
C&Sdiv, (J*V,fi) dvX.
N daX
(79)
This equation may be solved via the boundary element, however we now have a load term which is defined in the domain, to which the following transformations are applied :
D. A. TORTORELLI and
=
2. WANG
V’;i*(JD,82~.JTJ‘JTJJ‘D,6z~J’JT sB JJ‘JTD,62xTJT)V,Bdvx
= B[JJTV,~*JT(D,(JTV,~~)T6Z~~JJTV,~*JT(D,(JTV’~))Tb~~ s +JTV,fJTV,OiD,d2~.JJT(D,(JTV,~)))T~Z~.JJ‘JTV,a (D,(JTV,oi))T62~.JJ‘JTV,fiD,~2~TJTV’~.JJ‘JTV,~  D,62xTJTV,f* JJ‘JTV,b]
dvx
= B[D’(JTV,~*JTV,B)*JJ‘82~+JTV’~*JTV’Oidiv, s
(JJ,S?,x)
I= dB[JTV,fi.JTV,diS2~JTV,ti.~2~JTV,~JTV,/Z.62xJTV,8] f
lJJeTN dax+
[JTV,oi*62Xdiv, (JJ‘JTV,i)+JTV,fi sB
*Szx div’ (JJ‘JTV,ti)]
dvx,
(80)
where we systematically applied the definition of C&J*and eqns (4Q, (40), (41), (5), &), (25,), (25,), Piola’s identity div, (JJsT) = 0, (25 ,) and (34,). The above proof also requires the symmetry of JeT(Di (J‘V , x)i>‘, indeed, eqn (8,) and repeated use of eqn (47) yields (JT(D,(JrV,x))T)T = Dx(JTVx~)J’ = D,(V,~). So symmetry is proved as D,(V,n) is symmetric by eqn (12). Again, for clarity, the subscripts X and x denote differentiation of fields on F and B x QMwith respect to the first argument. Substituting eqn (80) into eqn (79) gives
&&J*V,&*N daxf 88,
[JTV,fi* JTV,oi&~
JTV,h*62~JTV,~
s as
JTV,~*~~xJTV,~]*JJTN
dux
[JTV’&*6,xdiv,
(JJ‘JTV’$
sB +JTV,fi*82~div’
F&~~K+~~~~K)

5
fi*V,82&N
(JJ‘JTV,Q]dvx+
dax 
da,f
JB,

Szoi sB
div, (J*V,advx.
li(JTV,~p~62~K s as2
aB (JTV,oiP.82~+SZ~P)J*V,fi.N s I &aiJ*V,&N dax f Jh
dax =
1201
A systematic approach to shape sensitivity analysis
Now, we use eqn (66,) and the property remaining domain integrals, whereby
s
&2J*V,di* N dax 
as [J‘V,fi.
JTVloid2~
solution to eliminate the
JTVIS*62~JTVi~
s
JBI
JTV,~.62~JTV,oi].JJTN +
of the fundamental
daxtJ~TV,oi(Y,~).82~(V,Q,;
~(JTV,~P~~2~K+6zqPK+~P62K)dax
s aB2
+62~P)J*Vlfi*N  @,oi(Y,
a)
da, = for
s JBI
cp)
JB (JTV12S62~
s,
;i~*v,8~oi. N daX +
~5~0iJ*V,f*N daX s JB2
Y E dB.
(82)
When incorporating the boundary element method, this equation may be solved in the same manner as eqn (70) so that no additional stiffness matrix assemblies or decompositions are required. As in the finite element case, we only require the pseudo load vector formation and back substitution to evaluate the pseudo response. Note that we could have derived this equation directly from eqn (70) with the use of eqn (80). Had we instead differentiated eqn (71) we would require the design derivatives of the fundamental solution l, as fi is design dependent even though 1 is not. In essence the difference between differentiating eqn (70) and (71) has to do with the fact that in eqn (70), lis any admissible function, therefore we can use S,i as it too will be admissible. Thus, when we differentiate eqn (70), the coefficients of the 6 2i terms cancel by virtue of eqn (70). On the other hand, if we differentiate eqn (71), the S,i terms do not cancel, and therefore must be included. The above analyses are valid for a specific design variation cp. However, we may compute each of the response gradient components (V,oi), by appropriately defining the components of cp as ‘PJ= 6, where 6ij is the Kronecker delta. Upon obtaining all of the (V,oi)i we may compute VG and subsequently 6G = VG * cp for any design perturbation. Note that the evaluation of Voirequires M pseudo analyses, M being the dimension of the design space P. 4.4. Aa’joint method
In the adjoint method we eliminate the response variations d2di and V,8$i. This is accomplished via the standard Lagrange multiplier method where we impose a subsidiary condition on G. In this case, the subsidiary condition is the governing integral equation (67) and we define the augmented functional G+ as
G+(a)
=
g(oi,JTV,oi*JJTN/K,X)Kdaxs JB
ldiv,
(JJ‘JTV,&)dvX,
(83)
s B
where fi plays the role of the Lagrange multiplier. Note that G = G+ since the augmented term is identically zero. Next, we proceed to take the design variation of G+ and note that 6G = 6G+ since the design variation of the augmented term is zero. Using eqns (74) and (79) [the latter is derived from the design variation of eqn (67)], and (67) (where &f replaces 1) we arrive at the following expression for 6G+ : 6G+ =
s
[V,@20iK+V&‘,62&*
JB
J*N +V,oi.&J*N
V,oi* J*N&K/K)
D. A. TORT~RELLIand Z. WANG
1202

a(JTV,riP*6211K+52qPK+~P62~da,+ s as2
+B&‘)J*Vr& 
N dax 
(JTV,~P*&2~ f JBI
~J*V,&&* N da,+ 6&J*Vrl*N s a& s aB78,
dax
SZoidiv, (J*Vll)doX.
(84)
sB Rearranging the above to isolate the explicit (rGE+, and implicit SC:, variations gives 6G+ = 6G$ +6G:,
(85)
where


&J*V,d* s as1 f
ae
N da, +
V,fi.&J*Vri
dv,
sB
~(JTV,BP~82~K+62qPK+~P62K)dax
z
+ as CJTVroiP*(T2X+82aP)(J*V,~*N+V,~K)dax, s, and using the sy~etry
SC: =
s
(V&;i)V,S#* JBI
(861
of J*
J*N dax +
s a&
(V,~K+J*V,&N)&oi
1
dax
6,8 div, (J*V,fi)dvx.
(87)
B
The objective now is to determine the ad&&r response fi which drives 6G: to zero. Then we can substitute this fi into eqn (84) so that 6G = 6G+ = SG,+ can be explicitly computed. Upon examination of the above, we see that to annihilate 6G: the following conditions must be attained : div, (JJ‘JTVI&
JTV,&JJTN
= 0
for
(X,icl6)EBxP,
XV,&
for
(X,4P)EaB,
for
(X,@)E~~B~X~~,
= V,&K
x2P, (88)
where we use the definition of J*. Note that with the exception of the data, this adjoint problem is identical to the primal problem [cf. eqn (66)] and hence, we may evaluate the adjoint response via either the finite or boundary element methods. So when these methods are used to evaluate the primal response, then the decomposed stiffness matrix may be used to compute the adjoint response. We merely form the adjoint load vector where we replace $’ and @”with V,& and V,& respectively ; and perform the back substitution. Substituting eqn (88,) into (86), gives
1203
A systematic approach to shape sensitivity analysis
+

s s
V2&JTV,~P*S2~K+62qPK)dax+
sB
ab
V,f&J*V18dvx
~(~~V,p~6~~K+ii~q~K+~~&K)da~
w
+
(J~v,~~P.~~x+~~c~~)(J*V,~~.N+V,~K)~~,.
(89)
s 84
This equation is best suited for finite element applications where all of the information is accurately extracted within the domain, If the boundary element method is used to evaluate the sensitivities it is advantageous to use eqn (80) to express
+
s
V2~(JTV,~P*62XK+&4P1Y)dax+
aB[JTV,& JTV,:&X s
a4
JTV,~.62~JTV,0iJTV,bi.BZ~JTV,fi].JJTN
dax
[JTV,&*62x div, (JJ‘JTV,&+JTV,;i.82X

div, (JJ‘JTV,B)]
du,
sB 
s
&TV,~P~~2~K+62qPK+~P82K)dax aB2
(90)
[(JTV,oiP*62x+62aP)(V,&K+J*V,fi.N)]dax.
+
We use eqns (66 ,) and (88 ,) to eliminate the domain integral, thus
SC,’ =
s as +
[V34*[email protected] dax +
s s s,
[V,g(V,ri*
J*NG,K/K)] da,
s aB,
V2&JTV,~P.~Z~K+62qPK)dax+
[JTV,fiSJTV,%~
a4
JTV,f~G2xJTV,oiJTV,0i*~2~JTV,~].JJTN 
dax
fi(JTV,~P.62~K+62qPK+BPBZK)da,
a4
+ aB [(JTV,0iP~62~+6zaP)(V,~K+J*V,fi*N)]dax.
(91)
This equation is best suited for boundary element applications in which accurate gradients are readily extracted over the surface [however, it may also be evaluated via a finite element analysis, although it has been shown to be less accurate than eqn (86) (Haug et al., 1986)]. Finally, note that for the boundary element adjoint sensitivity analysis, we again require no differentiation of the fundamental equation. In the adjoint approach, we are able to compute VG directly after the solution of one adjoint problem, whereas in the direct method we require the solution of M pseudo
D. A. TORTORELLI and Z. WANG
1204
problems. However, the adjoint approach requires the solution of one problem for each functional G, so if the functions out number the dimension of the design space, then the direct approach is preferred, and vice versa. 5. MATERIAL
DERIVATIVE METHOD
We may derive the material derivative method solutions directly from the domain parameterization method. In material derivative formulations 4 E 9 c 8 and the scalar # is analogous to time [cf. Haug et al. (1986), eqn (3.2.1)]. The reference configuration coincides with the body configuration so that x = X = x(X, cli). The design velocity v : B x 5SM+ Lin (sS~, 7y”) is defined as
This definition merely equates the design velocity to the design variation of x in the direction 40. Forthischoiceofx,D,x(~,~)~ = Usothat~(~,~) = l,K(X,cD) = 1, J(X,@) = I and n(x(X, a), @) = N(X, a). Substitution of these values into eqns (37), (39), (40), (41) and (45) allows us to obtain eqns (3.2.15)(3.2.18) in Haug et al. (1986), respectively, i.e. DzJ=Dlv;
DzJT =L(DlvT; D&’
= D,v;
DzJyT
=
D,yT;
Dzj=;
divv.
(93)
The simpIified equations (fil), (63) and (62) are identical to eqns (3.2.29)(3.2.31) in Haug et al. (1986). For this choice of the reference configuration, eqn (64) becomes &K=divvnsD,vTn
=vnH,
(94)
which isconsistent with eqns (3.2.34) and (3.2.95) in Haug et al. (1986), where the restriction v(X, a’; cp) = Y(X, a; (p)N(X), i.e. v is a normal velocity with magnitude V(X, bl, cp), and H = div n is the surface curvature which is applied to the rightmost equality. Equation (89) for ?I& is consistent with the equation derived via the “domain” approach of the material derivative method which appears in Haug et al. (1986). To reduce this equation, we let V$ = V,@ = 0, recall that J = 1, K = 1, J = I, v = f’n, i&K = Hv* n, and use the identity (resulting from an application of the divergence theroem on the (6&J‘JwT term)
s
V,&G,(JJ‘J‘)V,B B
dvw =

s a3
JTV,~*JTV,&.
JJmTN dax
B[JTV,~.JTV,(JTV,~.v)+JTV,~*JTV,(JTV,~~v)]Jd~~ s
to obtain
SC,’ =
s
V,AV,cxvn dax
@Hv*n da,+
Cm
s 8B 
sB
[V,l.V,(V,a*v)fV,a*V,(V,A*v)]dvx.
(95)
This is identical to eqn (24) in Park and Yoo (1988) after the appropriate simpli~cations have been made to their analysis. Equation (91) is consistent with the “boundary” version of the sensitivity equation which appears in Haug et al. (1986). To arrive at the material derivative version of this
1205
A systematic approach to shape sensitivity analysis
‘~fjjIbX, L
1 Fig. 5. Referential map.
equation, we set V,@ = 0, use J = 1, K = 1 and J = I, n = N ; and require that v = I/n, so that d2K = Hv * n. In light of these equalities, eqn (91) reduces to
6G,+=
s
8Hv.n dax
a.3
s JBI
V2~(V1a*nHvsn)daxt s JB,
s s 
V,u*vV,A.n
[V,I.V,avV,1.vV,al.n
dax
s JB
JB

dax+
[email protected](V,qP*v+&qP)da,
A(V,qP~v+82qP+qPHv~n)dax+
J4
s JBI
(V1aP*v+62aP)(VI~+VI~~n)dax.
(96) This equation is identical to eqn (23) in Merit (1988) after the appropriate are made to their analysis.
simplifications
6. EXAMPLE
In the following example, we present a simple problem with an analytical solution to validate the previous methodologies. It is henceforth understood that vectors and tensors will be identified by their matrices of components. 6.1. Mapping Consider the reference configuration via the map (see Fig. 5)
where (X,, X2), (x , , x2) and (w, L) are respectively the components of X E b2, 8:E 82 and @E&?‘. This mapping transforms the reference configuration ((X ,, X2) ; 0 < X1 < 1,0 ,
JKW The above results yield
= V,x(K’W =
(98)
1206
D. A. TORTORELLIand Z. WANG .I =
det J(X, @) = Lw,
(99)
(100)
and
J* = JJ’
JpT
zz
El[a;j=[:;j.
L
(101)
We evaluate Kfrom eqn (62). On the boundaries X2 = 0 and X2 = 1, N(X) = (0,  l)T and N(X) = (0, l)‘, hence
K(X,@)
= Lw
L.
(102)
ll[
Similarly on the boundaries X , = 0 and X , = 1, N(X) = ( 1, QT
K(X,@)
= Lw ll[
1 L
0
0
1
+1 = w. j( 0 )II
and
N(X)
=
(1, QT
(103)
W
For a design variation of cp = (6w, 6L), we obtain C&Xfrom eqns (22,)
(104) We may verify eqn (37) from (104)
V1~2x([email protected];
cp) =
6L
0
o
6w
(
(105)
)
and verify eqn (40) from (100)
r 6L L2 0
O SW ~ W2
1
(106)
.
Similarly we may verify eqn (45) from (99) 62J(X,@,cp) = L6w+wdL and from eqn (lOl), we may verify the chain rule
(107)
A systematic
approach
to shape sensitivity
1201
analysis
Vlcr.n=T,,xl W
a=0
sz 91 Q = ToLxz
s,
%
0 Fig. 6. Rectangular
domain
) Xl
L
n = T,,xl
Via.
with boundary
conditions.
&J* = 62JJ‘JT+J&J‘JT+JJ‘82JT,
(108)
On the boundaries X2 = 0, IX* = 1, we verify eqn (64) from (102) d,K(X,@;
cp) = 6L
(109)
and on the boundaries X, = 0, X, = 1, we verify eqn (64) from (103) 6,K(X,
@ ; cp>= 6w.
(110)
6.2. Laplace problem Consider the following boundary value problem defined over the previously described rectangular domain and with nonhomogeneous boundary conditions (see Fig. 6) : A,a(x, @) = 0
for
0< x, < L
cY(x,@) = 0
for
x, = 0,
a(~,@)
= T0Lx2
for
x,=L,
V,cl(w,@)n(x,@)
= T,,x,
for
x2 = 0,
V,a(x,Q)*n(x,Q)
= Tax,
for
x2=w.
and
0 < x2 < w,
(111)
Using separation of variables (Boyce and Diprima, 1986), the solution is readily determined a(~,@) = TOx,x2.
(112)
The referential response may be evaluated from eqns (28) and (112) : oi(X,@) = a(x(W),aq
= T,LwX,X2.
(113)
We may also verify eqn (66), using the above result for oi [i.e. eqn (113)] div(J*V&)
=E$+k$=O I
cqx,aq= 0 sAs30:9E
for
O
for
X, =O,
2
and
O
D. A. TORTORELLIand Z. WANG
1208
oi(X,@) &W,cP)
= ToLwX2
for
X,=1,
= TOLwX,
for
X2 = 0,
for
X2 = 1.
2
(114)
Here [email protected], = D,oi * El, i.e. the partial derivative of oi with respect to Xi of the (0 ; Ei} coordinate system. 6.3. Response function We choose to evaluate the net flux over surface S, (see Fig. 6), therefore the following response functional is defined
G(Q) =
s
Vu *n dax,
(115)
s1
where 9 = Va  n in eqn (72). In the reference domain, this functional is expressed as [see eqn (73)]
G(9)
=
s
Voi J*N dax,
S!
(116)
[email protected] = J‘Voi*JJTN/K= Voi*J*N/K. Using eqns (112), (113) and (101) we confirm that eqns (115) and (116) yield identical results :
G(@) =
s
*da,
S,axl
aoi s s,Lax,
Tow2
w
_ dax.
G2=
(117)
Here as/ax, = D ,a *ei, i.e. the partial derivative of c1with respect to Xiof the (CD ; e,} coordinate system. Next we consider the shape variation w + w + 6w and L + L + 6L and from eqn (112), we obtain 6&x,
m; cp) =
0.
(118)
We may also determine this result from eqn (76)
where, a direct calculation for eqn (113) gives b20i(X,(D; cp) = T~LX,X~CSW+TOWX,X~~L. 6.4. Direct diferentiation and aa’joint method We may verify the solution [email protected] in (119) by the direct differentiation d20isatisfies eqn (75). For our example, eqn (75) becomes
(119)
method, where
1209
A systematic approach to shape sensitivity analysis (_~+!%)!?+($_!E)$
for
O
syqx, 0)
and
O
= 0
[email protected](X,@) = T0LX26w+T0wX26L
for
XI = 0,
for
X, = 1,
V&Z* J*N +V,oi*&J*N
= 2T,,LX,&
for
X*=0,
V,820i.J*N+V,oi*62J*N
= 2T,Lx,6L
for
Xz = 1.
(120)
Where &J* is given in eqn (108). A straightforward calculation may be used to validate eqn (119) as the solution to the above problem. From eqn (74), we obtain (with Vig = V3g = 0)
6G=
[Vz~(V,62ai.J*N+V,ai.62J*NV,oi.J*N62K/K)+~621Y1dax.
(121)
s s1 The variation of G is also determined directly from eqn (117)
S,G(@;q) = TowSw. To verify the direct differentiation method, we substitute eqns (108), (110) and (119) into eqn (121), and obtain the identical result, i.e. 6G = TOwSw. We may also obtain the same result by the adjoint method. In this technique, we first evaluate the adjoint response 1, which from eqns (88) and (117), must satisfy div, (JJlJTV,;i)
JTV,l*JJTN
= 0
for
0 < Xi < 1 and
fi=O
for
XI =0,
fi=l
for
X,=1,
for
X2 =0
= 0
and
0 < XZ < 1,
Xz = 1.
(122)
The solution is determined by separation of variables as
LX,.
(123)
Using the above result, and eqns (99), (lOO), (lOl), (102), (103), (106), (107), (108), (110) and (113), we may evaluate the adjoint sensitivity via either eqn (89) or (91)
SGE+=
s Sl
d&K dax +
s s1
V24(V,oi* J*N&K/K)
JTV,&82~JTV,0iJs a4
dax + s dB
TV,@X12~JTV,&JJTN
(JTV,fi* JTV,oiSz~ dax
~(JTV,~P*~Z~K+82qPK+~P~2K)dax
+ s,(JTV,aiP.62~+62aP)J*V,L.N f = T,wSw.
dax (124)
1210
D. A. TORTORELLI and
Z. WANG
7. CONCLUSION
Shape sensitivities have been derived for the Laplace problem via the direct and adjoint approaches. The formulation is developed from the domain parameterization methodology, of which the material derivative method may be viewed as a special case. A generic derivation is developed and is readily adopted for both finite and boundary element applications. In the case of the boundary element approach, the need to differentiate the fundamental solution is circumvented ; and in the finite element method it is shown that the sensitivities may be expressed over the boundary for the adjoint approach, but not for the direct method. AcknowledgemenfsThe authors would like to thank Professor Donald E. Carlson and Dr Gunaseelan Krishnasamy of the Department of Theoretical and Applied Mechanics, University of Illinois at UrbanaChampaign, for their assistance. The first author would also like to thank Professors Bernard Rousselet and Denise Chenais of the Laboratoire de Mathematiques at the Universite de Nice for their helpful comments.
REFERENCES Aithal, R. and Saigal, S. (1990). Adjoint structure approach for shape sensitivity analysis using BEM. J. Engng Mech. 116( 12). 26632680. Arora, J. S., Lee, T. H. and Cardoso, J. B. (1991). Structure shape design sensitivity analysis: A unified viewpoint. AlAA/ASME/ASCE/ASC 32nd Structure, Structural Dynamics, and Materials Conference, Baltimore, MD, pp. 675683. Barone, M. R. and Yang, R. J. (1988). Boundary integral equations for recovery of design sensitivities in shape optimization. AIAA Jl26(5), 589594. Barone, M. R. and Yang, R. J. (1989). A boundary element approach for recovery of shape sensitivities in threedimensional elastic solids. Comp. Meth. Appl. Mech. Engng 74,6982. Becker, E. B., Carey, G. F. and Oden, J. T. (1981). Finite Elements: An Introduction. PrenticeHall, Englewood Cliffs, NJ. Boyce, W. E. and Diprima, R. C. (1986). Elementary Dtrerential Equations and Boundary Value Problems. Wiley, New York. Cardoso, J. B. and Arora, J. S. (1988). Variational method for design sensitivity analysis in nonlinear structural mechanics. AIAA Jl26(5), 595603. Carey, G. F. and Oden, J. T. (1983). Finite Elements: A Second Course. PrenticeHall, Englewood Cliffs, NJ. Carlson, D. E. (1991). Private communication. The University of Illinois, Department of Theoretical and Applied Mechanics. Carlson, D. E. and Hoger, A. (1986). On the derivatives of the principal invarients of a secondorder tensor. J. Elasticity 16, 221224. Cea, J. (1981a). Problems of Shape Optimal Design, Optimization of Distributed Parameter Structures (Edited by E. J. Haug and J. Cea), Vol. II, pp. 10051048. Sijthoff & Noordhoff, The Netherlands. Cea, J. (198lb). Numerical Methods in Shape Optimal Design, Optimization of Distributed Parameter Structures (Edited by E. J. Haug and J. Cea), Vol. II, pp. 10491987. Sijthotf& Noordhoff, The Netherlands. Chadwick, P. (1976). Continuum Mechanics: Concise Theory and Problems. Wiley, New York. Choi, K. K. and Huag, E. J. (1983). Shape design sensitivity analysis of elastic structures. J. Struct. Mech. 11(2), 231269. Choi, J. H. and Kwak, B. M. (1988a). Shape design sensitivity analysis of elliptic problems in boundary integral equation formulation. Mech. Struct. Mach. 16(2), 147165. Choi, J. H. and Kwak, B. M. (1988b). Boundary integral equation method for shape optimization of elastic structures. Int. J. Numer. Meth. Engng 26, 157991595. Ciarlet, P. (1988). Mathematical Elasticity. Elsevier Science Publishers, New York. Dems, K. and Mroz, Z. (1984). Variational approach by means of an adjoint system to structural optimization and sensitivity analysisII. Structural shape variation. Int. J. Solids Structures 20(6), 527552. DoCarmo, M. P. (1976). Differential Geometry of Curves and Surfaces. PrenticeHall, Englewood Cliffs, NJ. Gurtin, M. E. (1981). An Introduction to Continuum Mechanics. Academic Press, New York. Haber, R. B. (1987). A new variational approach to structural shape design sensitivity analysis. In Computer Aided Optimal Design : Structural and Mechanical Systems (Edited by C. A. Mota Soares). Springer, Berlin. Haug, E. J., Choi, K. K. and Komkov, V. (1986). Design Sensitivity Analysis of Structural Systems. Academic Press, New York. Merit, R. A. (1988). Shape optimization and identification of solid geometries considering discontinuities. J. Heat Transfer 110, 547550. Mota Soares, C. A., Rodrigues, H. C. and Choi, K. K. (1984). Shape optimal structural design using boundary elements and minimum compliance techniques. J. Mech. Tran. Auto. Des. 106,518523. Park, C. W. and Yoo, Y. M. (1988). Shape design sensitivity analysis of a twodimensional heat transfer system using the boundary element method. Comput. Struct. 28(4), 543560. Petryk, H. and Mroz, Z. (1986). Time derivatives of integral and functions defined on varying volume and surface domains. Arch. Mech. 38(5d), 6977724. Phelan, D. G. and Haber, R. B. (1989). Sensitivity analysis of linear elastic system using domain parameterization and mixed mutual energy principle. Comput. Meth. Appl. Mech. Engng 77,3159. Phelan, D. G., Vidal, C. and Haber, R. B. (1989). Explicit Sensitivity Analysis of Nonlinear Elastic Systems, Computer Aided Optimum Design : Recent Advances. Computational Mechanics Publications/Springer, Berlin. Rousselet, B. (1991). Introduction to shape sensitivity threedimensional and surface systems. Proceedings
A systematic
approach
to shape sensitivity
1211
analysis
NATOIDFG ASI Optimization oflarge Structural Systems, Berchtesgaden, Germany, 23 September4 October 1991. Tortorelli, D. A. and Haber, R. B. (1989). First order design sensitivities for transient conduction problems by an adjoint method. ht. J. Numer. Meth. Engng U(4), 733752. Tortorelli, D. A., Haber, R. B. and Lu, S. C.Y. (1990). Design sensitivites analysis for nonlinear transient thermal systems. Comput. Meth. Appl. Mech. Engng 75,6178. Tortorelli, D. A., Lu, S. C.Y. and Haber, R. B. (1990). Design sensitivities for elastodynamic systems. Mech. Struct. Mach. 18(l), 77106. Tortorelli, D. A., Haber, R. B. and Lu, S. C.Y. (1991). Adjoint sensitivity analysis for nonlinear dynamic thermoelastic systems. AIAA JI 29(2), 253263. Tortorelli, D. A., Subraamani, G. S., Lu, C. Y. and Haber, R. B. (1991). Sensitivity analysis for coupled thermoelastic systems. Znt. J. Solid Structures 27(12), 14771497. Zhana. 0. and Mukheriee, S. (1991). Design sensitivitv coefficients for linear elastic bodies with zones and corners by thederivative boundary element method. Int. J.Solids Structures 27(8), 983998. Zhao, Z. (1991). Shape Design Sensitivity Analysis and Optimization Using the Boundary Element Method. Springer, Berlin. Zolesio, J. P. (1981). The material derivative (or speed) method for shape optimization. In Optimization of Distributed Parameter Structures (Edited by E. J. Haug and J. Cea), Vol. II, pp. 1089I 151. Sijthoff & Noordhoff, The Netherlands.
APPENDIX A. 1. Alternative proof qf eqn (42) This proof is due to Carlson and Hodger (1986) and relies on the CayleyHamilton which states that the second order tensor Asatisfies its characteristic equation, i.e. A’n,(A)A’+n,(A)An,(A)1 where ?I,,
n2 and n, are the principal n,(A)
z,(A)
=
= 0,
1981)
(Al)
= tr(A), = $(a:(A)n,(A’)),
= det (A).
the trace of eqn (Al) and substitution n,(A)
(Gurtin,
invariants
B*(A) = $((tr(A))*tr(A’))
Taking
theorem
of eqn (A2), after some rearranging,
gives
:[~I(A~)~.,(A)~,(A*)+:(~:(A)HIA*))~,(A)I.
(A3)
Using eqns (24) (27) and (83) we may readily verify Dn,(A)U
Taking
the derivative
= 2n,(AU),
(44)
Dn,(A’)U
= 3n,(A%J).
(A5)
of eqn (A3) via the product Dn,(A)U
= x,(U),
Dn,(A’)U
rule [eqn (24)] and substituting
eqn (A4) we obtain
= f[3x,(A2U)n,(U)n,(A’)27c,(A)rt,(AU)
+~(~~~,(A)~,(U)~~,(AU)~,(A)~,(A’)XI(U))I = [A2rcn,(A)A+n,(A)I]T*U = [A‘n,(A)IT*U, where we used eqns (5) and (Al) to obtain immediately from the above and eqn (5).
(A6)
the second
and third equalities,
respectively.
A.2. Alternative proof af Piola’s identity This proof of Piola’s identity is supplied by Carlson (1991) and is a componentfree presented in Chadwick (1976) and Ciarlet (1988). Here the previous notation remains are constants. The following relations are required V, J(X, @)II = J(X, Q) tr (J‘(X,cP)D,J(X,@)u), which follows from eqns (14) and (45). Also required D:x([email protected])(u,v)
Equation
(43) follows
alternative to the proof in effect and a, u, YE vCN
(47)
is =
P,JPLW)u
= (0, JW, @)u)v,
(‘48)
1212
D. A. TORTORELLI and Z. WANG
which is obtained from eqns (30) and (12). The last preliminary argument is O,(J‘(Y&U++
= J‘(X,@)D,J(X,O)uJ‘(X,@)a = J‘(X,B)D,J(X,@)(J‘(X,cD)a)u,
(A9)
which follows from eqns (40) and (A8) where Y = J ‘(X, @)a, The proof is now presented : div, (J(X,@)JT(X,@))*a
= div, (J(X,@)J‘(X,@)a) =D,J(X,B).J‘(X,O)a+J(X,@)div,(J’(X,@)a) =.l(X,@)tr(J‘(X,@)D,J(X,Q))[J‘(X,@)a]) +J(X,cD)tr(D,J‘(X,(P)a) = J(X,Q)tr(J‘(X,@)D,J(X,@)[J‘(X,Q)a]) J(X,@)tr(J‘(X,@)D,J(X,Q)[J‘(X,@)a]) = 0,
(AlO)
where eqns (23,), (25,), (131) (with II = J ‘(X, @)a), and (133) (with the arbitrariness of u) are systematically applied.