- Email: [email protected]
Sensitivity analysis in shape optimization of continuum structures
034s79491SJ 53.00 - .w PclgamlJo Press Ltd.
SENSITIVITY ANALYSIS IN SHAPE OPTIMIZATION OF CONTINUUM STRUCTURES SHU-YU
WANG.~
YASBING SUN: and R. H. GALLAGHER$
University of Arizona, Tucson, .4Z 85721, U.S.A. (Receiced 15 .Vocember 1983; received for publication 19 January 1983)
Abstract-This paper describes an efficient method for sensitivity analysis in shape optimum design. One feature is the use of limited number of master nodes to characterize the surfaces of a set of isoparametric finite elements, and the adoption of their coordinates as design variables of the shape optimization. Another is the derivation of analytical formulations of the gradients of both the stiffness terms and the load vectors with respect to the design variables. A finite element analysis code is adapted to the purposes of the method and numerical examples are performed and comparisons made with sensitivity analysis based on forward finite differences. 1.INTRODUCTION
Mathematical programming approaches used in structural optimization generally fall into two categories. One involves gradient-based methods which use the derivatives of the objective function and constraint functions, as well as the values of functions themselves. Alternatively, nongradient techniques only use the values of the functions. Since the first category utilizes more information from analysis at a design point than the second one, a gradient-based method is expected to be more efficient than the latter, and it demands much fewer cases of structural reanalysis. The derivatives of a function with respect to the design variables (the gradient) represent and characterize the trends of the function variation when design variables are changing. They provide important information for choosing a search direction to obtain an improved, feasible, new design point. The calculation of derivatives of the objective function and constraints with respect to the design variables, known as design sensitivity analysis, plays a central role in optimization algorithms. Design sensitivity analysis is now applied not only in most of the efficient optimization methods but also in developing explicit approximations of the constraint functions and in some reanalysis methods[l]. Numerous papers have been published on design sensitivity analysis. Generally, they are restricted in application to trusses or frameworks in which the stiffness matrix can be expressed or approximated as a linear function of the design variables (member sizes). Sensitivity analysis applications to continuum structures are relatively rare. The reason may be that the structural response parameters (such as stress, displacement, natural frequency, etc.) of these structures, caused by external loading, are generally implicit functions of the design variables; there are no linear or other obvious relations between responses and design variables. tvisiting Scholar: Associate Professor, Department of Civil Engineering, Zhejiang University, Hangzhou, Zhejiang, The People’s Republic of China. :Visiting Scholar: Engineer, Xian Aircraft Company, Xian. Shaanxi Province, The People’s Republic of China. BFormerly, Dean, College of Engineering. Presently, VicePresident, Worcester Poly. Inst., Worcester, MA, USA.
In the case of shape optimization, the design variables may involve configuration parameters, nodal locations, boundary shapes, and overall topology[Z]. This design optimization problem is more complicated than in the case of member size selection. Consequently, the conventional tool for calculation of the function derivatives is by means of the finite difference technique. Finite difference calculations of the derivatives with respect to design variables for a problem with M design variables requires performance of the analysis for M + I different stiffness matrices. Each function derivative evaluation requires a completely new analysis of the structure. Moreover, these calculations are sensitive to small errors in the initial data. The difficulties of sensitivity analysis are obstacles in the way of efficient shape optimization and to improving the efficiency of existing algorithms. Because a large part of the information which sensitivity analysis requires has been obtained in the process of structural analysis with the finite element method, one way to improve is to incorporate sensitivity analysis into a finite element analysis code and to derive analytical formulations. Perhaps the first analytical work in design sensitivity analysis was by Zienkiewicz and Campbell in 1973[3] and Ramakrishnan and Francavilla in 1974[4]. They present exact formulations for the evaluation of stress derivatives, but restricted it to a discussion fo a ?-dimensional problem. Also, this work was lacking in details, especially for the derivatives of the loading vector. The objectives of this paper are to introduce an efficient sensitivity formulation for shape optimization of continuum structures. It consists of two features. The first is the use of a limited number of master nodes to characterize the surfaces of a set of isoparametric finite elements, and the adoption of their coordinates as design variables of the shape optimization. The second is the derivation of analytical formulations of the gradients. Emphases are put on the calculation of the derivatives of load vectors with respect to the design variables. The sensitivity analysis is incorporated into a finite element analysis code. Numerical examples are performed and comparisons are made with sensitivity analysis based on forward finite differences.
$56
SfiU-w 2. PROBLEM
ST.ATE.WEhT
.A.VD ASSUMPTIONS
In order to provide a common terminology for discussion. the optimization problem will be stated mathematically as follows: Minimize Subject to
F(c) g,(G) 2 0
i=I...I
h,(C) = 0
j=l...J
a,L
(1)
[~I{~} = {R 1
(2)
where [K] is the stiffness matrix, formed by assembling all element stiffness matrices [K],; {u} is the global nodal displacement vector, and {R} is the nodal force vector caused by applied loads. Using the Choleski decomposition method, we obtain
= [~l’[W
(3)
Let [r/](crj = {R’). Then, eqn (2) becomes [(l]r{R’} = {Rj.
(4)
The elements of matrix [L/l, an upper triangular matrix, are computed from [K] by a simple recursion algorithm, and {u } can easily be solved by a forward and a backward substitution. Once the necessary nodal displacements are determined, all other structural responses can be sequentially calculated. Therefore, the displacement behavior constraints and their derivatives with respect to design variables, can be expressed in the following form
I-2
20
(5)
{$}= A{$}
(6)
g,(Z)=
In general, both [K] and {R ] are dependent on the design variables. Hence, the derivatives of eqn (2) with respect to any design variable a, will be
let
m =i...M
in which 2 (or {al) is a M-dimensional vector of generalized design variables. All of the design parameters which can uniquely determine an engineering design will be written as a function of a’ (such as polynomial, fourier, exponential, etc.), according to specified practical problems. F is the objective function or merit function of the problem; g,, h, are inequality and equality constraints, respectively. Both of them will construct a feasible region in Mdimensional design space. Inequality constraints can be divided into behavior constraints, which are generally implicit in form, and side constraints, which limit the ranges of the design variables. For simplicity, this paper will consider only behavior constraints on stress and displacement. Most practical structures are modeled, for purposes of analysis, by the finite element method. We therefore adopt this method in our treatment of design sensitivity analysis procedures. In terms of displacements, the relevant global equations are
[Kl
WANGer al
then
Because eqn (9) and eqn (2) have the same coefficient matrix [K], the decomposed form of eqn (3) is still available. Only the forward and backward substitution are needed to obtain {&/cYa,}[5]. According to eqn (8) and (9), it is obvious that the key to sensitivity analysis is how to determine {R^I: namely, how to calculate the derivatives of [K] and {R} with respect to the design variables. The relationship between [K] or {R} and the design variables is, in general, nonlinear. 3. SELECHON
OF THE DESIGN
24
Master Nodes b
0
where U denotes the allowable displacement point for which the constraint is written.
VARIABLES
In the case of shape optimization, the design variables involve not only member sizes, thickness and cross section area, but also the boundary shapes, configuration parameters, and nodal point locations. The choice of design variables will change the character of the problem by changing the degree of nonlinearity of the objective and constraint functions. For shape optimization, one choice of design variables might be the nodal coordinates of the boundary. The advantage of this choice is that one would be able to obtain a general curved boundary, consistent with the finite element model, in which the structure is allowed to assume whatever shape is necessary to obtain the optimum design. The major problem is the resulting large number of design variables. Another choice would be the amplitude of a prescribed function, or the summation of several functions, each of which relates the nodal points along the boundary[6, 71. Here, we select a limited number of master nodes (the open dots in Fig. 1) which control the configuration of the structure[8]. The coordinates of these nodes are the design variables. Connecting these master nodes with straight lines, or quadratic or other kinds of curves, forms a coarse mesh and a set of subregions. Each subregion is then automatically
at the Fig. 1.
Sensitivity analysis in shape optimization of continuum StructureS
subdivided into structural elements specified by the number of divisions desired in the local coordinate directions. After the coordinates of the master nodes are determined, the locations of all internal nodes (the black dots in Fig. I) in these subregions can be readily computed, using isoparametric mapping or other interpolation techniques. In accordance with isoparametric mapping
Here, X, Y and Z are the global coordinates of a nodal point; 5, r] and { are the local or isoparametric mapping coordinates of the point (- I 5 5, 7, 5 5 1); N, is a shape function; Xi, Vi and Zi are the coordinates of the point i specifying the subregion; and n, is the total number of master nodes in a subregion or relevant to the subregion. Using this technique, the number of design variables can be drastically reduced. Consider, for example, the region shown in Fig. I(b). We take 3 points (I,4 and 7) as the master nodes to describe an upper boundary with 7 points. After the other boundaries are defined, we specify 6 divisions in the < direction and 3 divisions in the 4 direction. We can then calculate the coordinates of all internal nodes including those boundary points (2, 3, 5 and 6) which have the local coordinates 5 = -0.667, -0.333, +0.333 and +0.667, by using the two-dimensional version of eqn (10). Thus, the structure is divided into 18 elements, each of which is a linear element with 4 corner nodes (see Fig. lb left). If we require the results of structural analysis to be more accurate, we can form 8-node quadratic elements by appending corresponding nodes both on boundaries and in the interior (see Fig. lb right). If one wishes to use 1Znode cubic elements, however, the boundary still has a parabolic shape because it is determined by only 3 master nodes and the isoparametric mapping equation remains unchanged. For a cubic representation of the boundary one can take 4 points (I, 3, 5 and 7) as the master nodes. It will be appropriate to divide the structure into several subregions when the number of master nodes on a bounary is greater than 4. The problem arising from this case is that the boundary curve connecting adjacent subregions will not be smooth. One can smooth it by using “spline functions” or other approaches after the optimization process.
[B]~D][B] dX d Y dZ
PI, =
where [B] is a matrix relating element strains to assumed displacement field parameters, and generally composed of derivatives of the shape functions N,; [D] an elasticity matrix which connects the stresses and strains; and IJJ is the determinant of the Jacobian matrix [J], which represents the mapping between global coordinates (X, Y and Z) and local curvilinear coordinates (5. q, 0. Suppose the mechanical properties of the material are prescribed and do not change during the optimization process; then, the derivatives of the element stiffness matrix with respect to urn can be written as
+ [IqTID]
!g
IJ 1
m
+
[B]‘[DI[B] fg m
>
dt dq 4.
(a) Derivatives of [K] with respect 10 Ihe design vur-
The derivative of a structural stiffness rfiatrix with respect to a design variable a, can be obtained by assembling all derivatives of element stiffness matrices. Suppose a n-node isoparametric brick element is adopted; then, the element stiffness matrix [/cl, is as follows (see Ref. [9])
(12)
Because the integrand of the above equation is too complicated to be integrated, numerical integration must be performed. Thus, dk
il I da,
NG NG NG
(
=ccc WjWkWlfgJblIBIIJI m I k 1
+ Plr’[~l (13)
where W,, W,, W, are the weighting factors of numerical (Gaussian) integration, NC is the number of the Gauss rule adopted, and ( )$.. is the integrand function of a Gauss sampling point with local coordinates c,, qk and [,. The calculation of [B] and 1J I are standard aspects of finite element analysis. Now, we must construct a[B] and FIJI aa, aa, In three dimensional
elasticity
[B] = [B,B,. . . BN]
r
Ni..r 0
4. FORMULATIONS
iable a, (nnmely, [JK/&I,I)
857
0
N.,
0 0
1 (14)
N,,X,N,., and Ni., are derivatives of Ni with respect to X, Y and Z, respectively.
SJiU-YU WANGet al
The derivatrve of [Bf with respect to a, will be
and each submatrix
Fig. 2.
is
?[BJ
--zz
?a,
For a 20-node brick element we have N, =
$(I+j&l +
4 &J)(1 + Co)CSo + tlo + io - 2)
$1- <*)(I + q&l
+ &))(I - S,*fq,*i,’
+ f ( 1+ 50x1- B % 1+ 6”N1 -
where50 = fd, TO =
WI
and
50
=
s,*)i,*s,*
Lf.
According to the chain differentiation rule for a composite function, we obtain
Obviously, (L;iV,.,/aa,), (a”NiJi3a,,,) and (~ivi,z/~a,) can be obtained by differentiating eqn (16). Notice that IV*.:,IV,, and I%‘~,; are functions only of the local coordinates of the Gauss sampling point (C, q, i). Because the number of the Gauss rule adopted in the element does not change during the optimization process, the location and coordinates of Gauss sampling points are fixed. That means that IV,.:, Ni,q and N,,: are independent of the design variables. Hence
where Ntr is the derivative of fVi with respect to <. Therefore, in view of eqn (IO) in which
ax:
z-
and similar relationships can be obtained for Y.: . I 2,;. in view of eqn (15), IV,,: can be written in the following explicit form
N,,: =
2
_T
=$<,(I
m
aY,t
aa ”
dZ.
_1
da,
(21)
+rlo)(l +io) and
+ 2 f,( 1 - q 3( 1 + &I)(1 - v,2)i*2r,’
+; 4,(1 -
_
(22) (lg)
where (~X~&z,,,)denoted the partial derivative of the X coordinate of the ith node with respect to a,,,, IX,,,,) denotes a vector LdX,/aa, . , ax,jaa,_? and !_N,,_I can be evaluated from eqn (IS).
859
Sensitivity analysis in shape optimization of continuum structures
As mentioned previously, we have already prescribed the relations between nodal coordinates and design variables by using the concept of master nodes and isoparametric mapping (see eqn 10). so that the derivatives of element nodal coordinates {X,,) can be analytically determined. Similarly for evaluating (I?Y,&?a,), (c?Z,~/C?U~,etc. Hence
caused by a concentrated force {PI at a point C; {Q jc = 1’j&N]‘{4 1 dS; the equivalent nodal forces caused by a distributed surface load {q} acting on surface S; {G)< = j [y[N]r{g) dV; the equivalent nodal forces caused by a body force {g}; and {Hje = j jv[B]‘[D]{cO,j dV; the equivalent nodal forces caused by an initial strain (~1, such as a change of temperature.
(aN,,,/aa,),
For simplicity of explanation, when varying the design variable a, we only consider the change of the equivalent forces at node i of the element. (1) Concenfrotedforce. If there is a concentrated force {P} = LPXPyPZJT acting at point C of the considered element, the equivalent nodal forces at the ith node are given as
After
obtaining
(aNi,,/aa,,,)
and
(&Vi&Ya,), (d/da,)[B]
The determinant
can be formed readily. of [J] is x7< y,: Z,(
IJI = X,, Y,, z,’
(24)
IX,, Y,! z, i The derivative of 1J 1 with respect to a, therefore equals
(61, = M.~K~F,~I~ =
(w,~w~bJT
(27)
where (N,), denotes the shape function Ni at point C on which {P} is acting. Suppose, as in most engineering problems, that the value of concentrated force and its location relative to the element nodes, namely the local coordinates remains unchanged. Then {aF,/aa,}, = 0. In general, the value of force does not change, but the local coordinates of the point at which it acts will vary with the change of a,. Thus
m
=3
da,
/_P,VP,P,Jr.
With the chain differentiation a (N,), -= (2% Each component has ken calculated previously in eqns (17) and (22). Substituting the results of eqns (19) and (25) for (a[f?]/aa,) and (a IJ I /au,) into (13) the derivative of the element stiffness matrix [ak/aa,]~ can be determined. (b) Derivatives of {I?} with respect to the design variable (I,,, (namely, {dltlda,})
The total nodal force vector {R} in eqn (2) can be formed by assembling the equivalent nodal forces {R}c of the respective elements
da,
N:$N.
‘.*da,
(28)
rule we obtain 3+N_
“‘I2a,
_%
‘.‘aam >c
(29)
where N,,c, N,,4 and Ni,( can be calculated from eqn (18). Note that <, qc and < herein are the local coordinates of the force location (t,, qc and I;,), instead of those of Gauss sampling points. Using the relations between the global Cartesian and the local curvilinear coordinate system (coordinate transformation, eqn 17) the following expressions exists
ax, a K = da, (LNJ,{X’,) = LNmJ@} + LNJc{X,)
(30) in which ran is the total number general, IRL=
of elenients.
{FIc+{Q}e+{GI<+
where {F}, = [NJCr{Pl, the equivalent
In
(f$ nodal forces
where X, denotes the global X coordinate of the force location; subscript c denotes the calculation of functions for the local coordinates of point c (5,. qc and I,); {A’,,,,} denotes a vector L(aX,/au,)(aX@u,) . . . (aX,/aa,)J’; and LNL is a row vector LN,, N2.. . N,] of the shape function at point c. Therefore
WANGel ai.
SWJ-w
Similarly (Y.:),Z
iz.,),;
m
(U.,),$
(Z,,)$ m
'm
+ (YJc;" = 2* - LNJ<,(K,j m m w
-I- (Z.;),$& = 2 - i_Nl,{Z,,). m m
Obviously, the coefficients of the simultaneous equations (31) and (32) are exactly the transportation of Jacobian matrix, [J]r, shown in eqn (16), but replacing <, q, i with <,, ‘I,, i,. According to en~nee~ng requirements, we can establish the function or the relation between the force global location X,, I’,, Z, and design variables ri or {u}, calculate (8X,/&,), (aY,/aa,) and (aZ,/&,,), and form the r.h.s. of eqns (31) and (32); then, (@/&I,), (&7/&z,,,)and (aC,/&z,,,)can be solved from this simuttaneous equation. Substituting (a
For example, when {q ] acts on the surface < = 1, we have dS = IS x TI,,, di dq where S = LX,, Y,,Z,, Jr and T = LX., Y,,.Z,,j’.
Thus, for example for the S-direction
Q,=
’ i‘l
’ (:V;q,r. -I
-I
/S xTI);=,d<
dg.
{Q,), = LQ,xQi~Q~z_le' (qiY S x V;,
=f
z.Z
(Gt), = LGi~rGrYGiZJeT=
I
I
I
-1
-I
*, dl dq
J”j
(qN,L( Y,,.z,
z,
If {g] is always constant as in the case of gravity, and if one still uses the same Gauss points in the numerical integration, the local coordinates <, PJand
force
If considering (q) normal to the curved surface { = + 1, then we can transform the first kind of surface integration (the integrand is a vector function) into the second kind of surface integration (the integrand is a scalar function). Thus
-
XNiLg,g,g,JrlJId5d'ldi.(33)
(36)
-
X,,Z,,)f
Y,:x>,)I:c 2, dS dv
(37)
similar relationships can be written for normal surfaceloadson<=kl andq=Jrl. If 4 does not change, such as water pressure acting on the upstream surface of a dam and the mesh generated along the elevation, we have
(38)
i will remain changed. So
unchanged.
That
means
Ni is un-
where 2X__Laz,, aa,
da,
are aiready given by eqn (22). (4) Temperature force. Considering temperature effects, with T the temperature change from the stress-free state, the equivalent nodal force will be in which (a/&z,,,)]J ( can be obtained from eqn (25). force a surface (3) Surface force. Suppose {q) = LqrqrqzJr acts on a certain surface S of the element. The equivalent nodaI force at the ith node
Xhl
Sensitivity analysis in shape optimization of continuum structures Hence
{Q&r,), should bed one only for the elements subjected to surface loads. Then, sum them to obtain (6) Form the r.h.s. {/?I of eqn (9). Call a solver subroutine to get {iiu/du,,,~. Also
where (dNL.r/&s,,,),etc. are given by eqn (20) and (8 I J I /hJ by ew (25). 5. .MAINSTEPSOF THE SENSlT1Vl-l-Y ANALYSIS As mentioned above, we adopt an automatic geometric generation procedure to generate the internal nodal locations of only the points in a subregion influenced by the specific master nodes, whose coordinates are taken as design variables. The derivatives of all nodal coordinates with respect to any design variable can be obtained analytically by differentiating eqn (IO), and they will be constant if the isoparametric mapping adopted does not change in the process. Therefore, they can he kept for use. The number of Gauss sampling points and their local coordinates in an element are constant because the adopted Gauss order rule generally remains unchanged in the whole procedure. Hence, neither the shape functions LNJ (see eqn 15) at a Gauss point nor their derivatives with respect to local coordinates LN,:f. LN,,J and l_N,J (eqn I$) will vary. They can be evatuated and stored at the outset and read when needed in subsequent operations. The main steps of a finite element-based sensitivity analysis can be summarized as follows: (1) In the typical (kth) iteration, generate the nodal coordinates of the structure in terms of the coordinates of master nodes, which can be determined according to the current design variables {u”). (2) Calculate the structural response. In the process determine X,< . . . Z,, (eqn 17). [J], [J]-’ and IJ I at each Gauss point using the current (kth) nodal coordinates, solve the Cartesian derivatives of the shape functions LN,xJ LN,rJ and LN,,J in terms of eqn (16). construct [S) and [ic&and assemble them to obtain [K]. Finally, the joint displacements (a) can be solved from eqn (2). The element stiffness matrix [k],, total nodal displacement {a}, as well as [B] and (o} at each Gauss point, must be stored for later use. (3) Read LN,t_l L&J LNd andl_.K,,_J 1 f’,,J l_Z,,J. Calculate (~~~u~)[~] and (~Ni,~~~u~), etc. (See eqns 20-23). Then, construct (~,@a,,#] and evaluate (a/&,)]J] according to eqns (19) and (25) respectively. (4) In terms of eqn (13), perform the multiplication of the relevant matrices at each Gauss point. Sum them to get [i?k/&,& Then, calculate the second term of the r.h.s. of eqn (8). This can be simphfied by considering the relation {o), = [D]jB]{u),. Thus
(5) At the same time, (aG/d~~}~ and {aH/&z,), can be obtained in terms of eqns (34) and (39), element by element, in the loop of calculating [dk/hz,&{u},.
(7) Construct the constraint eqn (6) etc. 6.
gradient according to
Ex&MPLES AND COMPARISON WITH THE FINITE DIFFFRWCE
As
APPROACH
mentioned above, a large part of the information which design sensitivity analysis requires has been obtained in the process of structural analysis and can be stored for use. The computational efforts of sensitivity analysis with analytical formulations are slightly less than the efforts required in a finite difference approach. For example, the computational efforts of [~~~~u~~~~~~,using the analytical method, require 72n* f 54n multiplications and 54n’ + 3On additions, instead of 72n* + 108n and 54n2 + 87~ respectively, for the forward difference method. In a system with n = 20, this means that the analytical method needs 29,880 multiplications and 22,200 additions, and the forward difference method requires 30,960 and 23,340, respectively. To demonstrate the applicability of the analytical formulation, a shape optimization program has been written to solve plane stress/strain problems, The program employs an g-node isoparametric element. The approach to optimization that is adopted is the method of successive or sequential linear programmingflO] with an improved move limit, which can make an infeasible design point become feasible by moving in the gradient direction of the most violated constraint [ I I]. Example I, shown in Figs. 3 and 4, is a horizontal cross-section of an arch dam. It consists of a concrete arch body and rock foundations on both banks of a symmetrical valley. The arch is subjected to uniformly distributed loads. The finite element idealization has 14 elements and 73 nodes. We choose the Y-coordinates of master nodes 31, 39, 47, 57, 65 and 73 as design variables. These variables then are scaled in terms of the prescribed ranges for the purpose of obtaining the same order of magnitude. After they are determined, all nodal coordinates can be evaluated by a geometric generation surboutine. The first three variables specify the location of the structure and the configuration of upper surface; the remainder specify the shape of lower surface, as well as the structure sizes or thicknesses. The volume of the arch body (composed of 8 elements, thickness = 1M) is taken as the objective function. Principal stresses at the six master nodes, limited to -t 100 and -8OOT/M*, are taken as the behavior constraints. Stresses are calculated at these nodes using stress smoothing and averaging techniques[l2]. The lower and upper bounds of the design variables (a,‘ and a,“), which are referred to as “side constraints”, equal 0 and 1, respectively, for all design variables, because these variables have been scaled in the optimization process. The total constraints number 24.
------r----T----r----l--7-500
--400
-300
-200
-100
-+
X Cm)
30
Fig. 3
Y Cm!
t
“i
Symnelry
------T--------T----?--600
-500
-400
-300
-200
--do- -
X(m)
Fig. 4.
The initial design point is based on engineering judgment and a simplified analysis. After six iterations, where each one includes geometric generation, analysis for displacements and stresses, and evaluation of design derivatives, the objective function decreases from 1025m’ to 195m’. It achieves this both by moving the whole structure to the lower bound of Y coordinates at nodes 31-and 47, prescribed in advance, and by reducing body thickness. The principal compressive stresses at nodes 39, 57 and 73 achieve their limits. The optimum shape and its corresponding data are shown in Figs. 3, 4 and Tables l-3.
In order to examine the efficiency of the analytical method in solving practical problems, we used both this method and the forward finite difference method to obtain derivatives of objective function, displacements and stress constraints. The CPU times on the CYBER 175 are almost identical. The major results of sensitivity by both methods are listed in Table 4. They are in close agreement. It is noteworthy that any perturbation of initial data (e.g. if the angle of the surface has a 0.01” deviation or the nodal coordinates are slightly in error), will cause these derivatives to have significant deviations when they are calculated via forward finite differ-
Sensitivity analysis in shape optimization of continuum structures Table 1. Master nodal coordinates for Example I. Total iterations: 6 Final
Initial Y CM)
x (MI 31
I
19.290 30.611
0.000
i
57.500
O.Or!O
y (MI
-22.9EO
36.970
-44.042
39
X (M)
47
44.042
36,970
22.980
19.290
57
-32.558
27.330
-ia.a72
15.642
o.oco
43.500
65 71
L
37
5%
77
0.000
i
27.611
aJ3-77
7-m
/
15
.!?A?
Table 2. Design variables (scaled) for Example I initial
Final 0.30000
0.68740
al a2
0.68750
0.01528
"3
0.68740
0.00000
“4
0.66625
a5
0.85714
0.02421 I
0.00000
0.66625
-66.
0.02421
Table 3. Principal stresses* (T/M’) for Example I. Tension (+), compression (-) Initial 31 39
33.16 -91 .Ol
Final 11.73
0.00 -799.99
54.62 -94.73
-611.39
33.16
11.73
54.62
0.00
-a7.04
-539.25
-114.66
-500.02
65
2.19
-283.88
-0.26
-752.25
73
-87.04
539.75
114.6fi
300.02
47 i 57
* Allowable stress [tr+]= 100 T/M*
, [o-] = -800
T/d
t
Y(m)
$_
r----r--
---
Symmetry
-X(m)
Fig. 5.
Sw-YU
WANG et al.
Table 4 Comparison between analytical formulation formulation
Analytical +a?,? 85.3635
85.4528
612.4a23
612.2667
266.scaa
266.7636
3G3.3745
303.9867
- 0.6535;.01
-
0.5699E-01
o.oa
85.3764
142.56i!
66.39
.
0.04
612.4753
669.1596
9.25
.
0.02
266.8017
323.:359
21.25
0.04
303.8674
360.551i
13.65
0.11
- 0.6725E-01
0.3374E-01
49.a3
0.5707E-01
0.14
0.5879E-01
0.9230E-01
57.00
O.l601E-01
0.13
O.l755E-01
0.5106E-01
90.94
0.6603E-01
-
-
- 0.9664E-02
0.7411E+02
0.7399E+02
-
0.16
0.5045E+O2
o.a077E+02
- o.a320E+02
- o.a299E+02
0.25
- O.lOllE+03
- 0.7078E*02
29.99
O.l722E+OZ
O.l31lE+OZ
176.13
0.04
o.s21aE+02
o.a251E+02
58.13
- 0.04
- 0.239aE+04
- O.l551E+04
35.32
0.4475E+C4
23.34
- O.l664E+02
- 0.2309E-rO4
546.69
0.238x-01
60.10
0.3490E+04
0.3492E+04
0.06
0.362aE+04
0.3618E+03
0.3620E+03
0.06
0.4445E+03
O.l292E+04
190.66
- 0.3569E+03
0.4904E+03
237.41
- 0.3134E+03
Sec.
-
0.06
0.4829E+02
- 0.230aE+04
- 0.3136E+03
0.06
7.159
7.034
l
Step lenath is 0.001
*
Because of symnetric structure. x
2
l
ences, but they will not influence the results in the analytical method. Example 2, shown in Fig. 5, examines the ability to optimize structure shape and overall location when the initial design is of different shape. The problem is the same as in Example I except that a straight beam configuration is selected as the initial design. The initial volume is 1073 m’. After nine iterations (see Fig. 6 and Tables S-7), it gradually becomes a curvilinear arch and moves to the lower bound of the
Initial
1
= $-
3
and ta4
=k
6
Y coordinates. In the final configuration the principal stresses become active at nodes 39 and 65 (the thickness of abutments are limited by lower bounds). The volume of the final configuration is 99.6m’. Examples 3 and 4 illustrate the capability of the analytical method to handle different load cases. The data of Example 3 are the same as in Example 2, but we add a - 5°C temperature load at nodes 14, 29 and 57-72. Results are shown in Fig. 7 and in Tables 5, 6 and 8. The initial objective function value is again
Table 5. IMaster nodal coordinates \
Node
for Examples
and 3 Final
-
X(H)
Y(H)
-:I.595
35.000
X(M)
Y(M)
\\
31
-22.980
j9
0.000
0.000
2a.o:3'
"1.595
35.000
22.980
19.290
5;
-21.4"3
la.000
-19.916
18.000
0.000
-.
x.000
19.290
A? 55
: LULL
0.10
0.4827E+02
R
,licht deviatit If T""?cI!,3 1
-
- n.i665E+02
rk
initial
data
- o.a96aE-02
L
CPU
Correct
finl:e differeace'
- 0.8959E-02
a0 /aa
a /aa t !
Fonard
#light deviatlo ,f rnntnur
ilitlal
lrrec?
and forward finite difference
3.5'30 _‘1 ,:1 3 l
la.000 28.506 for Example 3.
l
i6.718 27.187”
19.916 16.718 * 27.691 for Exa-ale 3.
1
Sensitivity analysis in shape optimization of continuum structures
865
Table 6. Design variables (scaled) for Example 2 Flrarl Example 3
Example 2
Initial a1
0.61081
0.00000
0.00000
a2
0.16667
0.00030
0.01204
a3
0.61081
0.00000
0.00000
a4
0.77384
0.00003
0.00000
as a6
0.67347
0.01333
0.01283
0.77314
0.0000Q
o_nnnoo
Table 7. Principal stresses* (T/M’) for Example 2. Tension (+), compression (-) Initial
Final
o1
0
0,
! \
0)
31
25.16
-5.19
49.89
-15.01
39
-19.27
-94.99
-28.67
-800.00
47
25.16
-5.19
49.09
-15.01
57
-32.47
-67.60
-40.13
-315.02
65 ~ 73
73.00 32.47
-0.94 67.60
9.60 3
-800.00 315.07
_
'Allowable stress [o*] = 100 l/M', [o-] = -800 T/M2
Table 8. Principal stresses* (T/M’) for Example 3. Tension (+), compression (-) Final
Initial
31
01 -_ 35.08
o2 -14 . 8-3-
Ol
o2
37.50
-1.48 -800.05
39
-23.05
-142.70
-29.95
47
35.08
-14.83
37.50
-1.48
57
14.06
-51.70
-14.17
-295.07
-1.91 -51.70
8.05 -14.17
-800.03 -295.07
65 73
119.53** 14.06
* Allowable stress [o'] = 100 T/M', [a-] = -800 T/N2 .* The constraint at this point has been violated in the initial Stage.
Y(m)
t $
Symmetry
lnitiol
C---
FinalDesign Iteration 3
--
Fig. 6.
Shope
X(m)
SIC’-YUWANGel al.
Fig. 7. Table 9. Master nodal coordinates for Example 4. Total iterations: 8 Initial
31
Final Y(M)
X(Ml
15.cKlo
-6.878
x(M)
y(M)
-5.537
10,000
39
0.000
15.100
0.300
10.919
47
+6.878
15.000
5.537
1o.ooa
57
-5.805
il.000
-5.218
9.034
65
0.000
11.100
0.030
10.319
11.000
5.278
71
+r &p
4
Table IO. Design variables (scaled) for Example 4 Initial
Final
a1
0.71429
0.00000
'2
0.72857
0.131?3
a3
c.71429
0.00000
a4
0.62830
0.00000
a5 a6
0.77273
0.00000
0 fi7R7n
n ml"""
Table 1I. Principal stresses* (T/M’) for Example 4. Tension (+), compression (--) Final
Initial Node 31 39
0
0
c
70.57
95.37
3e.lc
_ -313.73
10.85
-440.88
70.57
95.68
-320.96
-50.65
S 236k 9.65
47
236.82
57
-100.26
65
271.49
73
-100.26 *Allowable
-
0.54
-320.96 stress
[o']
: 100 T/M',
2.43 -50
57
[T-I = -800 T/M2
38.18 -353.14 -
2.86
.,
7c
Sensitivity analysis in shape optimization of continuum structures function
867
reduce to 3.80m’. The in the arch body appears more
value will ftrther
stress distribution uniform.
CONCLUSlON
This paper has derived analytical formulations
Grovitationol
load
9~9
‘lm )
60
for calculating design sensitivity analysis. Four exampoles show that it is possible to use this method for any gradient-based optimization algorithm. The results obtained are inspiring and desirable. That means the methodology is basically correct, feasible, and reasonable. The next step will be the devel-
opment of a program for 3-dimensional problems and the exploitation of it to solve other kinds of structure responses, in addition to stresses and displacements. Acknowledgemenrs-The authors wish to thank Dr. T. Triffet, Associate Dean, College of Engineering, University of Arizona. for offering computation fund, and Mrs. Helen Gallery for her careful typing of the manuscript. REFERENCES
I20
I00
60
60 X(m)
40
20
00
40
Fig. 8.
1073 m’ and the final value is 101.3 m’ after 11 iterations. It is noteworthy that the constraint at node 65 is violated initially; the program can make an infeasible design point become feasible. After 11 iterations the shape converges to a curvilinear arch. The fourth and final example studies a problem involving gravity and point loads. The data and results of this example, which is subjected to gravitational force as well as a concentrated load 5’at the node 39, are shown in Fig. 8 and in Tables 9-l 1. The initial value of objective function is 50.87 m’ and the final value is 8.19 m’ after 8 iterations. As in Example 3 the initial design is infeasible. Half of the principal stress constraints are violated. They become feasible, however, within the first 3 iterations. The principal stresses of the final design (see Table 11) do not achieve their limits because the prescribed lower bounds on thickness are too large (0.6m at the middle of arch and 1.Om near the abutments). If we decrease the lower bounds of thickness to 0.2 m and 0.6 m respectively, the cr, stresses at nodes 3 1 and 47 will achieve their allowable values (+ 100 T/M*) and or at node 65 will be -5lO.lOT/M* after two additional iterations. At the same time the objective
I. U. Kirsch and G. Toledano, Approximate reanalysis for modifications of structural geometry. Compur. Srruclures 16, 269-277 (1983). 2. G. N. Vanderplaats, Structual optimization-past, oresent and future. AIAA J. 20. 992-1000 (1981). 3. b. C. Zienkiewicz and J. S. Camphell, Shape opti-
mization and sequential linear programming. Opfimum Srrucrural Desipn (Edited bv R. H. Gallagher and 0. C. iienkiewicz). ‘%ley, New York (1973). 4. C. V. Ramakrishnan and A. Francavilla, Structural shape optimization using penalty functions. J. Sfrucfural Mech. 3, 74-75, 4OW22. 5. U. Kirsch, Optimum Sfrucfural Design: Concepts, Mefhods and Applicafions. McGraw-Hill, New York (1981). 6. M. E. Botkin, Shape optimization of plate and shell structures. A/AA /. 20, 268-273 (1981). 7. R. Sharpe. Optimum design of arch dams. Proc. Isfn. Cio. Engrs., Paper 72OOS,Suppl. Vol. 7389 (1969). 8. Y. K. Cheung and M. F. Yeo, A Pracfical Infroducfion fo Finife Elemenf Analysis. Pitman, London (1979). 9. R. H. Gallagher, Finife Elemenf Analysis Fundamentals. Prentice-Hall. Enelewood Cliffs. New Jersey (1975). IO. J. N. Siddall, bpfi;al Engineering Design. Dekker, New York (1982). II. S. S. Bhavikatti and C. V. Ramakrishnan, Computational efficiency of improved move limit method of sequential linear programming for structural optimization. Compuf. Sfrucfures 11, 191-196 (1980). 12. E. Hinton and J. S. Campbell, Local and globalsmoothing of discontinuous finite element functions using a least squares smoothing method. Inr. I. Iv’um. Mefh. Engng 8, 461-480 (1974).
SENSITIVITY ANALYSIS IN SHAPE OPTIMIZATION OF CONTINUUM STRUCTURES SHU-YU
WANG.~
YASBING SUN: and R. H. GALLAGHER$
University of Arizona, Tucson, .4Z 85721, U.S.A. (Receiced 15 .Vocember 1983; received for publication 19 January 1983)
Abstract-This paper describes an efficient method for sensitivity analysis in shape optimum design. One feature is the use of limited number of master nodes to characterize the surfaces of a set of isoparametric finite elements, and the adoption of their coordinates as design variables of the shape optimization. Another is the derivation of analytical formulations of the gradients of both the stiffness terms and the load vectors with respect to the design variables. A finite element analysis code is adapted to the purposes of the method and numerical examples are performed and comparisons made with sensitivity analysis based on forward finite differences. 1.INTRODUCTION
Mathematical programming approaches used in structural optimization generally fall into two categories. One involves gradient-based methods which use the derivatives of the objective function and constraint functions, as well as the values of functions themselves. Alternatively, nongradient techniques only use the values of the functions. Since the first category utilizes more information from analysis at a design point than the second one, a gradient-based method is expected to be more efficient than the latter, and it demands much fewer cases of structural reanalysis. The derivatives of a function with respect to the design variables (the gradient) represent and characterize the trends of the function variation when design variables are changing. They provide important information for choosing a search direction to obtain an improved, feasible, new design point. The calculation of derivatives of the objective function and constraints with respect to the design variables, known as design sensitivity analysis, plays a central role in optimization algorithms. Design sensitivity analysis is now applied not only in most of the efficient optimization methods but also in developing explicit approximations of the constraint functions and in some reanalysis methods[l]. Numerous papers have been published on design sensitivity analysis. Generally, they are restricted in application to trusses or frameworks in which the stiffness matrix can be expressed or approximated as a linear function of the design variables (member sizes). Sensitivity analysis applications to continuum structures are relatively rare. The reason may be that the structural response parameters (such as stress, displacement, natural frequency, etc.) of these structures, caused by external loading, are generally implicit functions of the design variables; there are no linear or other obvious relations between responses and design variables. tvisiting Scholar: Associate Professor, Department of Civil Engineering, Zhejiang University, Hangzhou, Zhejiang, The People’s Republic of China. :Visiting Scholar: Engineer, Xian Aircraft Company, Xian. Shaanxi Province, The People’s Republic of China. BFormerly, Dean, College of Engineering. Presently, VicePresident, Worcester Poly. Inst., Worcester, MA, USA.
In the case of shape optimization, the design variables may involve configuration parameters, nodal locations, boundary shapes, and overall topology[Z]. This design optimization problem is more complicated than in the case of member size selection. Consequently, the conventional tool for calculation of the function derivatives is by means of the finite difference technique. Finite difference calculations of the derivatives with respect to design variables for a problem with M design variables requires performance of the analysis for M + I different stiffness matrices. Each function derivative evaluation requires a completely new analysis of the structure. Moreover, these calculations are sensitive to small errors in the initial data. The difficulties of sensitivity analysis are obstacles in the way of efficient shape optimization and to improving the efficiency of existing algorithms. Because a large part of the information which sensitivity analysis requires has been obtained in the process of structural analysis with the finite element method, one way to improve is to incorporate sensitivity analysis into a finite element analysis code and to derive analytical formulations. Perhaps the first analytical work in design sensitivity analysis was by Zienkiewicz and Campbell in 1973[3] and Ramakrishnan and Francavilla in 1974[4]. They present exact formulations for the evaluation of stress derivatives, but restricted it to a discussion fo a ?-dimensional problem. Also, this work was lacking in details, especially for the derivatives of the loading vector. The objectives of this paper are to introduce an efficient sensitivity formulation for shape optimization of continuum structures. It consists of two features. The first is the use of a limited number of master nodes to characterize the surfaces of a set of isoparametric finite elements, and the adoption of their coordinates as design variables of the shape optimization. The second is the derivation of analytical formulations of the gradients. Emphases are put on the calculation of the derivatives of load vectors with respect to the design variables. The sensitivity analysis is incorporated into a finite element analysis code. Numerical examples are performed and comparisons are made with sensitivity analysis based on forward finite differences.
$56
SfiU-w 2. PROBLEM
ST.ATE.WEhT
.A.VD ASSUMPTIONS
In order to provide a common terminology for discussion. the optimization problem will be stated mathematically as follows: Minimize Subject to
F(c) g,(G) 2 0
i=I...I
h,(C) = 0
j=l...J
a,L
(1)
[~I{~} = {R 1
(2)
where [K] is the stiffness matrix, formed by assembling all element stiffness matrices [K],; {u} is the global nodal displacement vector, and {R} is the nodal force vector caused by applied loads. Using the Choleski decomposition method, we obtain
= [~l’[W
(3)
Let [r/](crj = {R’). Then, eqn (2) becomes [(l]r{R’} = {Rj.
(4)
The elements of matrix [L/l, an upper triangular matrix, are computed from [K] by a simple recursion algorithm, and {u } can easily be solved by a forward and a backward substitution. Once the necessary nodal displacements are determined, all other structural responses can be sequentially calculated. Therefore, the displacement behavior constraints and their derivatives with respect to design variables, can be expressed in the following form
I-2
20
(5)
{$}= A{$}
(6)
g,(Z)=
In general, both [K] and {R ] are dependent on the design variables. Hence, the derivatives of eqn (2) with respect to any design variable a, will be
let
m =i...M
in which 2 (or {al) is a M-dimensional vector of generalized design variables. All of the design parameters which can uniquely determine an engineering design will be written as a function of a’ (such as polynomial, fourier, exponential, etc.), according to specified practical problems. F is the objective function or merit function of the problem; g,, h, are inequality and equality constraints, respectively. Both of them will construct a feasible region in Mdimensional design space. Inequality constraints can be divided into behavior constraints, which are generally implicit in form, and side constraints, which limit the ranges of the design variables. For simplicity, this paper will consider only behavior constraints on stress and displacement. Most practical structures are modeled, for purposes of analysis, by the finite element method. We therefore adopt this method in our treatment of design sensitivity analysis procedures. In terms of displacements, the relevant global equations are
[Kl
WANGer al
then
Because eqn (9) and eqn (2) have the same coefficient matrix [K], the decomposed form of eqn (3) is still available. Only the forward and backward substitution are needed to obtain {&/cYa,}[5]. According to eqn (8) and (9), it is obvious that the key to sensitivity analysis is how to determine {R^I: namely, how to calculate the derivatives of [K] and {R} with respect to the design variables. The relationship between [K] or {R} and the design variables is, in general, nonlinear. 3. SELECHON
OF THE DESIGN
24
Master Nodes b
0
where U denotes the allowable displacement point for which the constraint is written.
VARIABLES
In the case of shape optimization, the design variables involve not only member sizes, thickness and cross section area, but also the boundary shapes, configuration parameters, and nodal point locations. The choice of design variables will change the character of the problem by changing the degree of nonlinearity of the objective and constraint functions. For shape optimization, one choice of design variables might be the nodal coordinates of the boundary. The advantage of this choice is that one would be able to obtain a general curved boundary, consistent with the finite element model, in which the structure is allowed to assume whatever shape is necessary to obtain the optimum design. The major problem is the resulting large number of design variables. Another choice would be the amplitude of a prescribed function, or the summation of several functions, each of which relates the nodal points along the boundary[6, 71. Here, we select a limited number of master nodes (the open dots in Fig. 1) which control the configuration of the structure[8]. The coordinates of these nodes are the design variables. Connecting these master nodes with straight lines, or quadratic or other kinds of curves, forms a coarse mesh and a set of subregions. Each subregion is then automatically
at the Fig. 1.
Sensitivity analysis in shape optimization of continuum StructureS
subdivided into structural elements specified by the number of divisions desired in the local coordinate directions. After the coordinates of the master nodes are determined, the locations of all internal nodes (the black dots in Fig. I) in these subregions can be readily computed, using isoparametric mapping or other interpolation techniques. In accordance with isoparametric mapping
Here, X, Y and Z are the global coordinates of a nodal point; 5, r] and { are the local or isoparametric mapping coordinates of the point (- I 5 5, 7, 5 5 1); N, is a shape function; Xi, Vi and Zi are the coordinates of the point i specifying the subregion; and n, is the total number of master nodes in a subregion or relevant to the subregion. Using this technique, the number of design variables can be drastically reduced. Consider, for example, the region shown in Fig. I(b). We take 3 points (I,4 and 7) as the master nodes to describe an upper boundary with 7 points. After the other boundaries are defined, we specify 6 divisions in the < direction and 3 divisions in the 4 direction. We can then calculate the coordinates of all internal nodes including those boundary points (2, 3, 5 and 6) which have the local coordinates 5 = -0.667, -0.333, +0.333 and +0.667, by using the two-dimensional version of eqn (10). Thus, the structure is divided into 18 elements, each of which is a linear element with 4 corner nodes (see Fig. lb left). If we require the results of structural analysis to be more accurate, we can form 8-node quadratic elements by appending corresponding nodes both on boundaries and in the interior (see Fig. lb right). If one wishes to use 1Znode cubic elements, however, the boundary still has a parabolic shape because it is determined by only 3 master nodes and the isoparametric mapping equation remains unchanged. For a cubic representation of the boundary one can take 4 points (I, 3, 5 and 7) as the master nodes. It will be appropriate to divide the structure into several subregions when the number of master nodes on a bounary is greater than 4. The problem arising from this case is that the boundary curve connecting adjacent subregions will not be smooth. One can smooth it by using “spline functions” or other approaches after the optimization process.
[B]~D][B] dX d Y dZ
PI, =
where [B] is a matrix relating element strains to assumed displacement field parameters, and generally composed of derivatives of the shape functions N,; [D] an elasticity matrix which connects the stresses and strains; and IJJ is the determinant of the Jacobian matrix [J], which represents the mapping between global coordinates (X, Y and Z) and local curvilinear coordinates (5. q, 0. Suppose the mechanical properties of the material are prescribed and do not change during the optimization process; then, the derivatives of the element stiffness matrix with respect to urn can be written as
+ [IqTID]
!g
IJ 1
m
+
[B]‘[DI[B] fg m
>
dt dq 4.
(a) Derivatives of [K] with respect 10 Ihe design vur-
The derivative of a structural stiffness rfiatrix with respect to a design variable a, can be obtained by assembling all derivatives of element stiffness matrices. Suppose a n-node isoparametric brick element is adopted; then, the element stiffness matrix [/cl, is as follows (see Ref. [9])
(12)
Because the integrand of the above equation is too complicated to be integrated, numerical integration must be performed. Thus, dk
il I da,
NG NG NG
(
=ccc WjWkWlfgJblIBIIJI m I k 1
+ Plr’[~l (13)
where W,, W,, W, are the weighting factors of numerical (Gaussian) integration, NC is the number of the Gauss rule adopted, and ( )$.. is the integrand function of a Gauss sampling point with local coordinates c,, qk and [,. The calculation of [B] and 1J I are standard aspects of finite element analysis. Now, we must construct a[B] and FIJI aa, aa, In three dimensional
elasticity
[B] = [B,B,. . . BN]
r
Ni..r 0
4. FORMULATIONS
iable a, (nnmely, [JK/&I,I)
857
0
N.,
0 0
1 (14)
N,,X,N,., and Ni., are derivatives of Ni with respect to X, Y and Z, respectively.
SJiU-YU WANGet al
The derivatrve of [Bf with respect to a, will be
and each submatrix
Fig. 2.
is
?[BJ
--zz
?a,
For a 20-node brick element we have N, =
$(I+j&l +
4 &J)(1 + Co)CSo + tlo + io - 2)
$1- <*)(I + q&l
+ &))(I - S,*fq,*i,’
+ f ( 1+ 50x1- B % 1+ 6”N1 -
where50 = fd, TO =
WI
and
50
=
s,*)i,*s,*
Lf.
According to the chain differentiation rule for a composite function, we obtain
Obviously, (L;iV,.,/aa,), (a”NiJi3a,,,) and (~ivi,z/~a,) can be obtained by differentiating eqn (16). Notice that IV*.:,IV,, and I%‘~,; are functions only of the local coordinates of the Gauss sampling point (C, q, i). Because the number of the Gauss rule adopted in the element does not change during the optimization process, the location and coordinates of Gauss sampling points are fixed. That means that IV,.:, Ni,q and N,,: are independent of the design variables. Hence
where Ntr is the derivative of fVi with respect to <. Therefore, in view of eqn (IO) in which
ax:
z-
and similar relationships can be obtained for Y.: . I 2,;. in view of eqn (15), IV,,: can be written in the following explicit form
N,,: =
2
_T
=$<,(I
m
aY,t
aa ”
dZ.
_1
da,
(21)
+rlo)(l +io) and
+ 2 f,( 1 - q 3( 1 + &I)(1 - v,2)i*2r,’
+; 4,(1 -
_
(22) (lg)
where (~X~&z,,,)denoted the partial derivative of the X coordinate of the ith node with respect to a,,,, IX,,,,) denotes a vector LdX,/aa, . , ax,jaa,_? and !_N,,_I can be evaluated from eqn (IS).
859
Sensitivity analysis in shape optimization of continuum structures
As mentioned previously, we have already prescribed the relations between nodal coordinates and design variables by using the concept of master nodes and isoparametric mapping (see eqn 10). so that the derivatives of element nodal coordinates {X,,) can be analytically determined. Similarly for evaluating (I?Y,&?a,), (c?Z,~/C?U~,etc. Hence
caused by a concentrated force {PI at a point C; {Q jc = 1’j&N]‘{4 1 dS; the equivalent nodal forces caused by a distributed surface load {q} acting on surface S; {G)< = j [y[N]r{g) dV; the equivalent nodal forces caused by a body force {g}; and {Hje = j jv[B]‘[D]{cO,j dV; the equivalent nodal forces caused by an initial strain (~1, such as a change of temperature.
(aN,,,/aa,),
For simplicity of explanation, when varying the design variable a, we only consider the change of the equivalent forces at node i of the element. (1) Concenfrotedforce. If there is a concentrated force {P} = LPXPyPZJT acting at point C of the considered element, the equivalent nodal forces at the ith node are given as
After
obtaining
(aNi,,/aa,,,)
and
(&Vi&Ya,), (d/da,)[B]
The determinant
can be formed readily. of [J] is x7< y,: Z,(
IJI = X,, Y,, z,’
(24)
IX,, Y,! z, i The derivative of 1J 1 with respect to a, therefore equals
(61, = M.~K~F,~I~ =
(w,~w~bJT
(27)
where (N,), denotes the shape function Ni at point C on which {P} is acting. Suppose, as in most engineering problems, that the value of concentrated force and its location relative to the element nodes, namely the local coordinates remains unchanged. Then {aF,/aa,}, = 0. In general, the value of force does not change, but the local coordinates of the point at which it acts will vary with the change of a,. Thus
m
=3
da,
/_P,VP,P,Jr.
With the chain differentiation a (N,), -= (2% Each component has ken calculated previously in eqns (17) and (22). Substituting the results of eqns (19) and (25) for (a[f?]/aa,) and (a IJ I /au,) into (13) the derivative of the element stiffness matrix [ak/aa,]~ can be determined. (b) Derivatives of {I?} with respect to the design variable (I,,, (namely, {dltlda,})
The total nodal force vector {R} in eqn (2) can be formed by assembling the equivalent nodal forces {R}c of the respective elements
da,
N:$N.
‘.*da,
(28)
rule we obtain 3+N_
“‘I2a,
_%
‘.‘aam >c
(29)
where N,,c, N,,4 and Ni,( can be calculated from eqn (18). Note that <, qc and < herein are the local coordinates of the force location (t,, qc and I;,), instead of those of Gauss sampling points. Using the relations between the global Cartesian and the local curvilinear coordinate system (coordinate transformation, eqn 17) the following expressions exists
ax, a K = da, (LNJ,{X’,) = LNmJ@} + LNJc{X,)
(30) in which ran is the total number general, IRL=
of elenients.
{FIc+{Q}e+{GI<+
where {F}, = [NJCr{Pl, the equivalent
In
(f$ nodal forces
where X, denotes the global X coordinate of the force location; subscript c denotes the calculation of functions for the local coordinates of point c (5,. qc and I,); {A’,,,,} denotes a vector L(aX,/au,)(aX@u,) . . . (aX,/aa,)J’; and LNL is a row vector LN,, N2.. . N,] of the shape function at point c. Therefore
WANGel ai.
SWJ-w
Similarly (Y.:),Z
iz.,),;
m
(U.,),$
(Z,,)$ m
'm
+ (YJc;" = 2* - LNJ<,(K,j m m w
-I- (Z.;),$& = 2 - i_Nl,{Z,,). m m
Obviously, the coefficients of the simultaneous equations (31) and (32) are exactly the transportation of Jacobian matrix, [J]r, shown in eqn (16), but replacing <, q, i with <,, ‘I,, i,. According to en~nee~ng requirements, we can establish the function or the relation between the force global location X,, I’,, Z, and design variables ri or {u}, calculate (8X,/&,), (aY,/aa,) and (aZ,/&,,), and form the r.h.s. of eqns (31) and (32); then, (@/&I,), (&7/&z,,,)and (aC,/&z,,,)can be solved from this simuttaneous equation. Substituting (a
For example, when {q ] acts on the surface < = 1, we have dS = IS x TI,,, di dq where S = LX,, Y,,Z,, Jr and T = LX., Y,,.Z,,j’.
Thus, for example for the S-direction
Q,=
’ i‘l
’ (:V;q,r. -I
-I
/S xTI);=,d<
dg.
{Q,), = LQ,xQi~Q~z_le' (qiY S x V;,
=f
z.Z
(Gt), = LGi~rGrYGiZJeT=
I
I
I
-1
-I
*, dl dq
J”j
(qN,L( Y,,.z,
z,
If {g] is always constant as in the case of gravity, and if one still uses the same Gauss points in the numerical integration, the local coordinates <, PJand
force
If considering (q) normal to the curved surface { = + 1, then we can transform the first kind of surface integration (the integrand is a vector function) into the second kind of surface integration (the integrand is a scalar function). Thus
-
XNiLg,g,g,JrlJId5d'ldi.(33)
(36)
-
X,,Z,,)f
Y,:x>,)I:c 2, dS dv
(37)
similar relationships can be written for normal surfaceloadson<=kl andq=Jrl. If 4 does not change, such as water pressure acting on the upstream surface of a dam and the mesh generated along the elevation, we have
(38)
i will remain changed. So
unchanged.
That
means
Ni is un-
where 2X__Laz,, aa,
da,
are aiready given by eqn (22). (4) Temperature force. Considering temperature effects, with T the temperature change from the stress-free state, the equivalent nodal force will be in which (a/&z,,,)]J ( can be obtained from eqn (25). force a surface (3) Surface force. Suppose {q) = LqrqrqzJr acts on a certain surface S of the element. The equivalent nodaI force at the ith node
Xhl
Sensitivity analysis in shape optimization of continuum structures Hence
{Q&r,), should bed one only for the elements subjected to surface loads. Then, sum them to obtain (6) Form the r.h.s. {/?I of eqn (9). Call a solver subroutine to get {iiu/du,,,~. Also
where (dNL.r/&s,,,),etc. are given by eqn (20) and (8 I J I /hJ by ew (25). 5. .MAINSTEPSOF THE SENSlT1Vl-l-Y ANALYSIS As mentioned above, we adopt an automatic geometric generation procedure to generate the internal nodal locations of only the points in a subregion influenced by the specific master nodes, whose coordinates are taken as design variables. The derivatives of all nodal coordinates with respect to any design variable can be obtained analytically by differentiating eqn (IO), and they will be constant if the isoparametric mapping adopted does not change in the process. Therefore, they can he kept for use. The number of Gauss sampling points and their local coordinates in an element are constant because the adopted Gauss order rule generally remains unchanged in the whole procedure. Hence, neither the shape functions LNJ (see eqn 15) at a Gauss point nor their derivatives with respect to local coordinates LN,:f. LN,,J and l_N,J (eqn I$) will vary. They can be evatuated and stored at the outset and read when needed in subsequent operations. The main steps of a finite element-based sensitivity analysis can be summarized as follows: (1) In the typical (kth) iteration, generate the nodal coordinates of the structure in terms of the coordinates of master nodes, which can be determined according to the current design variables {u”). (2) Calculate the structural response. In the process determine X,< . . . Z,, (eqn 17). [J], [J]-’ and IJ I at each Gauss point using the current (kth) nodal coordinates, solve the Cartesian derivatives of the shape functions LN,xJ LN,rJ and LN,,J in terms of eqn (16). construct [S) and [ic&and assemble them to obtain [K]. Finally, the joint displacements (a) can be solved from eqn (2). The element stiffness matrix [k],, total nodal displacement {a}, as well as [B] and (o} at each Gauss point, must be stored for later use. (3) Read LN,t_l L&J LNd andl_.K,,_J 1 f’,,J l_Z,,J. Calculate (~~~u~)[~] and (~Ni,~~~u~), etc. (See eqns 20-23). Then, construct (~,@a,,#] and evaluate (a/&,)]J] according to eqns (19) and (25) respectively. (4) In terms of eqn (13), perform the multiplication of the relevant matrices at each Gauss point. Sum them to get [i?k/&,& Then, calculate the second term of the r.h.s. of eqn (8). This can be simphfied by considering the relation {o), = [D]jB]{u),. Thus
(5) At the same time, (aG/d~~}~ and {aH/&z,), can be obtained in terms of eqns (34) and (39), element by element, in the loop of calculating [dk/hz,&{u},.
(7) Construct the constraint eqn (6) etc. 6.
gradient according to
Ex&MPLES AND COMPARISON WITH THE FINITE DIFFFRWCE
As
APPROACH
mentioned above, a large part of the information which design sensitivity analysis requires has been obtained in the process of structural analysis and can be stored for use. The computational efforts of sensitivity analysis with analytical formulations are slightly less than the efforts required in a finite difference approach. For example, the computational efforts of [~~~~u~~~~~~,using the analytical method, require 72n* f 54n multiplications and 54n’ + 3On additions, instead of 72n* + 108n and 54n2 + 87~ respectively, for the forward difference method. In a system with n = 20, this means that the analytical method needs 29,880 multiplications and 22,200 additions, and the forward difference method requires 30,960 and 23,340, respectively. To demonstrate the applicability of the analytical formulation, a shape optimization program has been written to solve plane stress/strain problems, The program employs an g-node isoparametric element. The approach to optimization that is adopted is the method of successive or sequential linear programmingflO] with an improved move limit, which can make an infeasible design point become feasible by moving in the gradient direction of the most violated constraint [ I I]. Example I, shown in Figs. 3 and 4, is a horizontal cross-section of an arch dam. It consists of a concrete arch body and rock foundations on both banks of a symmetrical valley. The arch is subjected to uniformly distributed loads. The finite element idealization has 14 elements and 73 nodes. We choose the Y-coordinates of master nodes 31, 39, 47, 57, 65 and 73 as design variables. These variables then are scaled in terms of the prescribed ranges for the purpose of obtaining the same order of magnitude. After they are determined, all nodal coordinates can be evaluated by a geometric generation surboutine. The first three variables specify the location of the structure and the configuration of upper surface; the remainder specify the shape of lower surface, as well as the structure sizes or thicknesses. The volume of the arch body (composed of 8 elements, thickness = 1M) is taken as the objective function. Principal stresses at the six master nodes, limited to -t 100 and -8OOT/M*, are taken as the behavior constraints. Stresses are calculated at these nodes using stress smoothing and averaging techniques[l2]. The lower and upper bounds of the design variables (a,‘ and a,“), which are referred to as “side constraints”, equal 0 and 1, respectively, for all design variables, because these variables have been scaled in the optimization process. The total constraints number 24.
------r----T----r----l--7-500
--400
-300
-200
-100
-+
X Cm)
30
Fig. 3
Y Cm!
t
“i
Symnelry
------T--------T----?--600
-500
-400
-300
-200
--do- -
X(m)
Fig. 4.
The initial design point is based on engineering judgment and a simplified analysis. After six iterations, where each one includes geometric generation, analysis for displacements and stresses, and evaluation of design derivatives, the objective function decreases from 1025m’ to 195m’. It achieves this both by moving the whole structure to the lower bound of Y coordinates at nodes 31-and 47, prescribed in advance, and by reducing body thickness. The principal compressive stresses at nodes 39, 57 and 73 achieve their limits. The optimum shape and its corresponding data are shown in Figs. 3, 4 and Tables l-3.
In order to examine the efficiency of the analytical method in solving practical problems, we used both this method and the forward finite difference method to obtain derivatives of objective function, displacements and stress constraints. The CPU times on the CYBER 175 are almost identical. The major results of sensitivity by both methods are listed in Table 4. They are in close agreement. It is noteworthy that any perturbation of initial data (e.g. if the angle of the surface has a 0.01” deviation or the nodal coordinates are slightly in error), will cause these derivatives to have significant deviations when they are calculated via forward finite differ-
Sensitivity analysis in shape optimization of continuum structures Table 1. Master nodal coordinates for Example I. Total iterations: 6 Final
Initial Y CM)
x (MI 31
I
19.290 30.611
0.000
i
57.500
O.Or!O
y (MI
-22.9EO
36.970
-44.042
39
X (M)
47
44.042
36,970
22.980
19.290
57
-32.558
27.330
-ia.a72
15.642
o.oco
43.500
65 71
L
37
5%
77
0.000
i
27.611
aJ3-77
7-m
/
15
.!?A?
Table 2. Design variables (scaled) for Example I initial
Final 0.30000
0.68740
al a2
0.68750
0.01528
"3
0.68740
0.00000
“4
0.66625
a5
0.85714
0.02421 I
0.00000
0.66625
-66.
0.02421
Table 3. Principal stresses* (T/M’) for Example I. Tension (+), compression (-) Initial 31 39
33.16 -91 .Ol
Final 11.73
0.00 -799.99
54.62 -94.73
-611.39
33.16
11.73
54.62
0.00
-a7.04
-539.25
-114.66
-500.02
65
2.19
-283.88
-0.26
-752.25
73
-87.04
539.75
114.6fi
300.02
47 i 57
* Allowable stress [tr+]= 100 T/M*
, [o-] = -800
T/d
t
Y(m)
$_
r----r--
---
Symmetry
-X(m)
Fig. 5.
Sw-YU
WANG et al.
Table 4 Comparison between analytical formulation formulation
Analytical +a?,? 85.3635
85.4528
612.4a23
612.2667
266.scaa
266.7636
3G3.3745
303.9867
- 0.6535;.01
-
0.5699E-01
o.oa
85.3764
142.56i!
66.39
.
0.04
612.4753
669.1596
9.25
.
0.02
266.8017
323.:359
21.25
0.04
303.8674
360.551i
13.65
0.11
- 0.6725E-01
0.3374E-01
49.a3
0.5707E-01
0.14
0.5879E-01
0.9230E-01
57.00
O.l601E-01
0.13
O.l755E-01
0.5106E-01
90.94
0.6603E-01
-
-
- 0.9664E-02
0.7411E+02
0.7399E+02
-
0.16
0.5045E+O2
o.a077E+02
- o.a320E+02
- o.a299E+02
0.25
- O.lOllE+03
- 0.7078E*02
29.99
O.l722E+OZ
O.l31lE+OZ
176.13
0.04
o.s21aE+02
o.a251E+02
58.13
- 0.04
- 0.239aE+04
- O.l551E+04
35.32
0.4475E+C4
23.34
- O.l664E+02
- 0.2309E-rO4
546.69
0.238x-01
60.10
0.3490E+04
0.3492E+04
0.06
0.362aE+04
0.3618E+03
0.3620E+03
0.06
0.4445E+03
O.l292E+04
190.66
- 0.3569E+03
0.4904E+03
237.41
- 0.3134E+03
Sec.
-
0.06
0.4829E+02
- 0.230aE+04
- 0.3136E+03
0.06
7.159
7.034
l
Step lenath is 0.001
*
Because of symnetric structure. x
2
l
ences, but they will not influence the results in the analytical method. Example 2, shown in Fig. 5, examines the ability to optimize structure shape and overall location when the initial design is of different shape. The problem is the same as in Example I except that a straight beam configuration is selected as the initial design. The initial volume is 1073 m’. After nine iterations (see Fig. 6 and Tables S-7), it gradually becomes a curvilinear arch and moves to the lower bound of the
Initial
1
= $-
3
and ta4
=k
6
Y coordinates. In the final configuration the principal stresses become active at nodes 39 and 65 (the thickness of abutments are limited by lower bounds). The volume of the final configuration is 99.6m’. Examples 3 and 4 illustrate the capability of the analytical method to handle different load cases. The data of Example 3 are the same as in Example 2, but we add a - 5°C temperature load at nodes 14, 29 and 57-72. Results are shown in Fig. 7 and in Tables 5, 6 and 8. The initial objective function value is again
Table 5. IMaster nodal coordinates \
Node
for Examples
and 3 Final
-
X(H)
Y(H)
-:I.595
35.000
X(M)
Y(M)
\\
31
-22.980
j9
0.000
0.000
2a.o:3'
"1.595
35.000
22.980
19.290
5;
-21.4"3
la.000
-19.916
18.000
0.000
-.
x.000
19.290
A? 55
: LULL
0.10
0.4827E+02
R
,licht deviatit If T""?cI!,3 1
-
- n.i665E+02
rk
initial
data
- o.a96aE-02
L
CPU
Correct
finl:e differeace'
- 0.8959E-02
a0 /aa
a /aa t !
Fonard
#light deviatlo ,f rnntnur
ilitlal
lrrec?
and forward finite difference
3.5'30 _‘1 ,:1 3 l
la.000 28.506 for Example 3.
l
i6.718 27.187”
19.916 16.718 * 27.691 for Exa-ale 3.
1
Sensitivity analysis in shape optimization of continuum structures
865
Table 6. Design variables (scaled) for Example 2 Flrarl Example 3
Example 2
Initial a1
0.61081
0.00000
0.00000
a2
0.16667
0.00030
0.01204
a3
0.61081
0.00000
0.00000
a4
0.77384
0.00003
0.00000
as a6
0.67347
0.01333
0.01283
0.77314
0.0000Q
o_nnnoo
Table 7. Principal stresses* (T/M’) for Example 2. Tension (+), compression (-) Initial
Final
o1
0
0,
! \
0)
31
25.16
-5.19
49.89
-15.01
39
-19.27
-94.99
-28.67
-800.00
47
25.16
-5.19
49.09
-15.01
57
-32.47
-67.60
-40.13
-315.02
65 ~ 73
73.00 32.47
-0.94 67.60
9.60 3
-800.00 315.07
_
'Allowable stress [o*] = 100 l/M', [o-] = -800 T/M2
Table 8. Principal stresses* (T/M’) for Example 3. Tension (+), compression (-) Final
Initial
31
01 -_ 35.08
o2 -14 . 8-3-
Ol
o2
37.50
-1.48 -800.05
39
-23.05
-142.70
-29.95
47
35.08
-14.83
37.50
-1.48
57
14.06
-51.70
-14.17
-295.07
-1.91 -51.70
8.05 -14.17
-800.03 -295.07
65 73
119.53** 14.06
* Allowable stress [o'] = 100 T/M', [a-] = -800 T/N2 .* The constraint at this point has been violated in the initial Stage.
Y(m)
t $
Symmetry
lnitiol
C---
FinalDesign Iteration 3
--
Fig. 6.
Shope
X(m)
SIC’-YUWANGel al.
Fig. 7. Table 9. Master nodal coordinates for Example 4. Total iterations: 8 Initial
31
Final Y(M)
X(Ml
15.cKlo
-6.878
x(M)
y(M)
-5.537
10,000
39
0.000
15.100
0.300
10.919
47
+6.878
15.000
5.537
1o.ooa
57
-5.805
il.000
-5.218
9.034
65
0.000
11.100
0.030
10.319
11.000
5.278
71
+r &p
4
Table IO. Design variables (scaled) for Example 4 Initial
Final
a1
0.71429
0.00000
'2
0.72857
0.131?3
a3
c.71429
0.00000
a4
0.62830
0.00000
a5 a6
0.77273
0.00000
0 fi7R7n
n ml"""
Table 1I. Principal stresses* (T/M’) for Example 4. Tension (+), compression (--) Final
Initial Node 31 39
0
0
c
70.57
95.37
3e.lc
_ -313.73
10.85
-440.88
70.57
95.68
-320.96
-50.65
S 236k 9.65
47
236.82
57
-100.26
65
271.49
73
-100.26 *Allowable
-
0.54
-320.96 stress
[o']
: 100 T/M',
2.43 -50
57
[T-I = -800 T/M2
38.18 -353.14 -
2.86
.,
7c
Sensitivity analysis in shape optimization of continuum structures function
867
reduce to 3.80m’. The in the arch body appears more
value will ftrther
stress distribution uniform.
CONCLUSlON
This paper has derived analytical formulations
Grovitationol
load
9~9
‘lm )
60
for calculating design sensitivity analysis. Four exampoles show that it is possible to use this method for any gradient-based optimization algorithm. The results obtained are inspiring and desirable. That means the methodology is basically correct, feasible, and reasonable. The next step will be the devel-
opment of a program for 3-dimensional problems and the exploitation of it to solve other kinds of structure responses, in addition to stresses and displacements. Acknowledgemenrs-The authors wish to thank Dr. T. Triffet, Associate Dean, College of Engineering, University of Arizona. for offering computation fund, and Mrs. Helen Gallery for her careful typing of the manuscript. REFERENCES
I20
I00
60
60 X(m)
40
20
00
40
Fig. 8.
1073 m’ and the final value is 101.3 m’ after 11 iterations. It is noteworthy that the constraint at node 65 is violated initially; the program can make an infeasible design point become feasible. After 11 iterations the shape converges to a curvilinear arch. The fourth and final example studies a problem involving gravity and point loads. The data and results of this example, which is subjected to gravitational force as well as a concentrated load 5’at the node 39, are shown in Fig. 8 and in Tables 9-l 1. The initial value of objective function is 50.87 m’ and the final value is 8.19 m’ after 8 iterations. As in Example 3 the initial design is infeasible. Half of the principal stress constraints are violated. They become feasible, however, within the first 3 iterations. The principal stresses of the final design (see Table 11) do not achieve their limits because the prescribed lower bounds on thickness are too large (0.6m at the middle of arch and 1.Om near the abutments). If we decrease the lower bounds of thickness to 0.2 m and 0.6 m respectively, the cr, stresses at nodes 3 1 and 47 will achieve their allowable values (+ 100 T/M*) and or at node 65 will be -5lO.lOT/M* after two additional iterations. At the same time the objective
I. U. Kirsch and G. Toledano, Approximate reanalysis for modifications of structural geometry. Compur. Srruclures 16, 269-277 (1983). 2. G. N. Vanderplaats, Structual optimization-past, oresent and future. AIAA J. 20. 992-1000 (1981). 3. b. C. Zienkiewicz and J. S. Camphell, Shape opti-
mization and sequential linear programming. Opfimum Srrucrural Desipn (Edited bv R. H. Gallagher and 0. C. iienkiewicz). ‘%ley, New York (1973). 4. C. V. Ramakrishnan and A. Francavilla, Structural shape optimization using penalty functions. J. Sfrucfural Mech. 3, 74-75, 4OW22. 5. U. Kirsch, Optimum Sfrucfural Design: Concepts, Mefhods and Applicafions. McGraw-Hill, New York (1981). 6. M. E. Botkin, Shape optimization of plate and shell structures. A/AA /. 20, 268-273 (1981). 7. R. Sharpe. Optimum design of arch dams. Proc. Isfn. Cio. Engrs., Paper 72OOS,Suppl. Vol. 7389 (1969). 8. Y. K. Cheung and M. F. Yeo, A Pracfical Infroducfion fo Finife Elemenf Analysis. Pitman, London (1979). 9. R. H. Gallagher, Finife Elemenf Analysis Fundamentals. Prentice-Hall. Enelewood Cliffs. New Jersey (1975). IO. J. N. Siddall, bpfi;al Engineering Design. Dekker, New York (1982). II. S. S. Bhavikatti and C. V. Ramakrishnan, Computational efficiency of improved move limit method of sequential linear programming for structural optimization. Compuf. Sfrucfures 11, 191-196 (1980). 12. E. Hinton and J. S. Campbell, Local and globalsmoothing of discontinuous finite element functions using a least squares smoothing method. Inr. I. Iv’um. Mefh. Engng 8, 461-480 (1974).