Layout optimization of stiffeners in stiffened composite plates with thermal residual stresses

Available online at www.sciencedirect.com

Finite Elements in Analysis and Design 40 (2004) 1233 – 1257 www.elsevier.com/locate/!nel

Layout optimization of sti'eners in sti'ened composite plates with thermal residual stresses Xiangmin Wanga , J.S. Hansena , D.C.D. Oguamanamb;∗ a

Institute for Aerospace Studies, University of Toronto, 4925 Duerin St., Toronto, Ont., Canada M3H 5T6 b Department of Mechanical and Industrial Engineering, Ryerson University, 350 Victoria St., Toronto, Ont., Canada M5B 2K3 Received 5 January 2003; received in revised form 8 May 2003; accepted 10 June 2003

Abstract The in7uence of manufacturing induced thermal residual stresses in sti'ened, symmetrically laminated plates on the optimal locations of the sti'eners is investigated using !nite element analysis method. The sti'eners are either discrete or continuous. The objective is to determine the optimal locations of these sti'eners in order to maximize the !rst natural frequency of the sti'ened plate. The optimization problem is solved using the method of moving asymptotes (MMA). A two stage solution approach is adopted: the !rst solves a thermal problem and the results are then used in the second stage to solve a free vibration problem. The results of the numerical simulations indicate that thermal residual stresses in7uence the magnitude of the optimum !rst natural frequencies of only the plates sti'ened with continuous sti'eners. However, the optimum locations for both the continuous and the discrete sti'eners are invariant to thermal residual stresses. Both the arrangement of the sti'eners and the magnitude of the optimum !rst natural frequency are strongly in7uenced by the ply-angles of the basic laminates and those of the sti'eners. ? 2003 Elsevier B.V. All rights reserved. Keywords: Sti'eners; Optimization; Laminated plates; Thermal; Residual stress; Sti'ened plates

1. Introduction The presence of residual stresses in structures is usually a consequence of externally applied loads, environmentally induced loads due to thermal e'ects or moisture absorption, or a ∗

Corresponding author. Tel.: +1-416-979-5000x7490; fax: +1-416-979-5265. E-mail address: [email protected] (D.C.D. Oguamanam).

0168-874X/$ - see front matter ? 2003 Elsevier B.V. All rights reserved. doi:10.1016/j.!nel.2003.06.003

1234

X. Wang et al. / Finite Elements in Analysis and Design 40 (2004) 1233 – 1257

Nomenclature A Ab [B] [BG ] [D] [DG ] {e} {eG } G {f} I Ib Jb {k} [K] [K]e [K]b [KGR ] [KGRe ] L m [M ] [Me ] [Mb ] [N ] {M } N {N } [Q] s T uG u U vG v wG w W {} G {} {e }

area of the plate cross-sectional area of the sti'ener matrix relating linear strain components and nodal displacements matrix relating non-linear strain components and nodal displacements constitutive matrix between middle-surface linear strains and resultant forces constitutive matrix between middle-surface non-linear strains and resultant forces vector of linear strain components vector of non-linear strain components vector of edge tractions moment of inertia of the plate cross-sectional moment of inertia of the sti'ener cross-sectional polar moment of inertia of the sti'ener vector of middle-surface curvatures global sti'ness matrix element sti'ness matrix sti'ness matrix of the sti'ener global geometric sti'ness matrix due to thermal residual stresses element geometric sti'ness matrix due to thermal residual stresses Lagrangian mass per unit area global mass matrix element mass matrix mass matrix of the sti'ener matrix of interpolation functions vector of resultant moments number of elements in !nite element model vector of resultant forces matrix of material sti'nesses length direction of the sti'ener kinetic energy displacement in x-direction middle-surface displacement in x-direction strain energy displacement in y-direction middle-surface displacement in y-direction displacement in z-direction middle-surface displacement in z-direction work of external forces vector of thermal coeIcients of expansion vector of nodal displacements vector of nodal displacements for element

X. Wang et al. / Finite Elements in Analysis and Design 40 (2004) 1233 – 1257

JT {} G {GL } {GN } b    { G} b x y

1235

change in temperature vector of strain components vector of linear strain components vector of non-linear strain components cross-sectional twist of the sti'ener eigenvalue of the plate total potential energy mass density vector of stresses cross-sectional rotation of the sti'ener middle-surface rotation variable in x-direction middle-surface rotation variable in y-direction

combination of environmental e'ects and boundary constraints. These factors are often observed to act simultaneously. The investigation of systems where residual stresses result from a combination of environmental e'ects and boundary constraints is reported in the studies by Sai Ram and Sinha [1,2]. The studies by Foldager et al. [3], de Faria and Hansen [4], and de Almeida and Hansen [5,6] are some of the more recent attempts to examine the role of manufacturing induced thermal residual stresses on the elastic buckling and vibration behaviour of sti'ened composite plates and cylinders. In particular, de Almeida and Hansen [5] explore the idea of sti'ening plates along the perimeter with materials whose thermal coeIcient and elastic properties are di'erent from those of the base plate. The motivation is to exploit the di'erences in thermal and mechanical properties to induce favourable residual-stress state in the plate. They show enhanced buckling capability of over 100% in some cases. This investigation is extended to free vibration in [6] where improvements in natural frequency of close to 100% is observed in some scenario. It is worth mentioning that both the locations and the stacking sequences of the sti'eners used in these studies are predetermined. Foldager et al. [3] re-examines the concept presented by de Almeida and Hansen [5] with regard to sti'ened, laminated cylinders under buckling load. The sti'eners are arranged such that they extend from one end of the cylinder to the other and their locations are predetermined. However, the stacking sequences of the sti'eners and the base cylinder are optimally determined for maximum buckling load. The study by de Faria and Hansen [4] is also another attempt at inducing thermal residual stresses in structures by exploiting the di'erences in thermal and elastic properties. However, rather than introduce sti'eners as in Refs. [3–6], de Faria and Hansen [4] examine the notion of having sub-regions in a composite structure with completely di'erent stacking sequences, which invariably results in piecewise continuous material properties across the structure. This discontinuity in material properties is then optimally explored during manufacturing to induce thermal residual stresses that enchance overall structural buckling capability. The current study, like the study in [6], attempts to induced favourable thermal residual stresses in a sti'ened, laminated plate in order to enhance the free vibration properties of the structure.

1236

X. Wang et al. / Finite Elements in Analysis and Design 40 (2004) 1233 – 1257

Unlike the study in [6], however, no prior assumption is made on the locations of the sti'eners. The only assumptions are on the number of sti'eners, the type of sti'eners (i.e. their discreteness or continuousness, and orientation), and the stacking sequences of both the sti'eners and the base plate. The investigations by Grall [7] and Jaunky and Knight [8] are examples of studies dedicated to the optimal design of sti'ened structures in regard to buckling and vibration via the use of the sti'ener distribution as design variables. In a di'erent approach to the problem, Kataoka-Filho [9] shows that the face-sheet distribution in optimal design of composite sandwich structures takes the form of sti'ened structures. The sti'eners considered in the present study are always located symmetrically about the midplane of the base laminated plate. The continuous sti'eners are unidirectional laminates because unidirectional laminates have high sti'ness and nearly zero thermal expansion coeIcient along the direction of !bres. For any given type of sti'eners, their locations are optimally determined based on the maximization of the !rst natural frequency of the resulting structure. The optimization process is implemented using the method of moving asymptotes (MMA) [10,11]. It is assumed that the structure is stress free at the elevated temperature (cure or consolidation temperature) and the thermal residual stresses are developed during the cooling phase to lower temperature (operating or room temperature). During the cooling process, the sti'ener contracts less than the basic laminate and hence resists the tendency of the basic laminate to shrink further. This e'ect induces tensile residual stresses in the basic laminate and compressive residual stresses in the sti'ener. The magnitude and distribution of the residual stresses depend on many factors such as material properties, temperature di'erence, layer-angles of the basic laminate and sti'eners, and the geometric arrangement of the sti'eners. The resulting distribution of residual stresses may be either bene!cial or detrimental to the mechanical behaviour of the composite sti'ened plate depending on the particular condition and the di'erent arrangements of the sti'eners. Based on the results of the numerical simulations implemented, the optimum locations of the sti'eners, either the discrete or continuous type, are observed to be invariant to thermal residual stresses but are more strongly in7uenced by their stacking sequences. Both thermal residual stresses and stacking sequences a'ect the magnitude of the optimum !rst natural frequency with continuous sti'eners. 2. Problem description and formulation For ease of presentation and brevity, the problem of interest is formally described in this section in conjunction with the non-!nite dimensional formulation. The objective of this study is to maximise the fundamental (or !rst) natural frequency of a sti'ened, symmetrically laminated composite plate with thermal residual stresses by determining the optimal locations of the sti'eners. The distribution of the thermal residual stresses is a consequence of the cure or consolidation process during the manufacturing of the plate. Further, the sti'ened plate is unconstrained hence the thermal residual stresses are solely a consequence of thermal coeIcient of expansion mismatches. These mismatches are either through the thickness of the plate or spatial in the plane of the plate. While thermal coeIcients of expansion are not temperature invariant [12], they are, however, assumed constant for the purposes of this study.

X. Wang et al. / Finite Elements in Analysis and Design 40 (2004) 1233 – 1257

1237

The Reissner-Mindlin plate theory [13] is assumed to be adequate to model the plate. The solution approach can be divided into three processes. In the !rst, the thermal residual stresses are obtained using linear strain–displacement relations. Thereafter, the full nonlinear strain–displacement is used to determine the reference state for the free vibration problem. This stage incorporates the thermal residual stresses determined in the preceding step. Finally, a formulation that is based on the stationarity of the Rayleigh quotient is implemented to determine the eigenvalue (or natural frequency) sensitivity of the plate with respect to the sti'ener locations. The major aspects of the solution approach as discussed above are now detailed in the following subsections. 2.1. Basic equation A displacement !eld that is appropriate for the Reissner-Mindlin plate hypothesis may be written as u(x; G y; z) = u(x; y) + z x (x; y); v(x; G y; z) = v(x; y) + z

y (x; y);

w(x; G y; z) = w(x; y);

(1)

where u(x; y), v(x; y) and w(x; y) are midplane displacements, respectively, de!ned along the x-, yand z-axes; the rotation of the cross section is denoted as x , when it is about the y-axis, and y , when it is about the x-axis. The convention adopted in this study in de!ning the variables is such that quantities with over-line are de!ned in the three-dimensional space (x, y and z) while those without the over-line are restricted to the midplane of the plate (i.e., a planar xy space). The strain vector G can be decomposed into a sum of linear and non-linear components. Thus {} G = {} G L + {} G N:

(2)

From the non-linear strain–displacement relations [14], the vector of the linear strain component can be expanded as GxL = u; x + zkxx ; GyL = v; y + zkyy ; L = u; y + v; x + zkxy ; #Gxy

#GxzL = w; x +

x;

L #Gyz = w; y +

y;

(3)

where the midplane curvatures of the plate which are denoted as kxx , kyy and kxy can be written as kxx =

x; x

;

kyy =

y; y ;

kxy =

x; y

+

y; x :

Note that ( )a; b represents the partial derivative of ( )a with respect to b.

(4)

1238

X. Wang et al. / Finite Elements in Analysis and Design 40 (2004) 1233 – 1257

The vector of non-linear strain component can be expanded and written as GxN = 12 [(u;2x + v;2x + w;2x ) + 2z(u; x

x; x

GyN = 12 [(u;2y + v;2y + w;2y ) + 2z(u; y

+ v; x

x; y

+ v; y

N = (u; x v; y + v; x v; y + w; x w; y ) + z(u; x #Gxy

+z 2 (

x; x x; y

+

y; x )

x; y

+ z2 (

y; y )

+ u; y

2 x; x

+ z2 ( x; x

+

2 x; y

+ v; x

2 y; x )]

+

2 y; y )]

y; y

+ v; y

y; x )

y; x y; y )

#GxzN = (u; x

x

+ v; x

y)

+ z(

x x; x

N #Gyz = (u; y

x

+ v; y

y)

+ z(

x x; y

+ +

y y; x ) y y; y ):

(5)

It is worthwhile to de!ned vectors for the in-plane linear strain component T , the midplane linear strain components {e}T , and the midplane non-linear strain components {eG }T . Thus {}T , [u; x v; y u; y + v; x ];

(6)

{e}T , [u; x v; y u; y + v; x kxx kyy kxy w; x + {eG }T , [u; x u; y v; x v; y w; x w; y

x; x

x w; y

+

x; y y; x y; y

y ];

(7)

x y ]:

(8)

2.2. Thermal problem The determination of the thermal residual stresses in the plate is the !rst of three steps in determining the optimal locations of the sti'eners. As mentioned earlier, the thermal residual stresses result from mismatches in thermal coeIcients of expansion which is either through the thickness of the plate or due to special preassigned variations. The thermal coeIcients of expansion are also assumed to be temperature invariant. At this stage of the solution process, the plate is free of any mechanically imposed boundary condition(s) but rigid-body motion is eliminated. There are no mechanical loads and linear strain–displacement relation is assumed. Further, the symmetric sequence of the laminates implies that there is no membrane–bending coupling. The stress–strain relation adopted involves thermal e'ects directly and is represented by G } { G } = [Q]({ G − JT ); G

(9)

G is the material sti'ness matrix in the global coordinate, {} where [Q] G is the vector of total strains, {} G is the vector of thermal coeIcients of expansion in the global coordinate, and change in temperature is denoted as JT . The corresponding total potential energy of the plate is now given as

1 T R = {e} [D]{e} dA − ({0 }T {N } + {k}T {M }) dA; (10) 2 A A

X. Wang et al. / Finite Elements in Analysis and Design 40 (2004) 1233 – 1257

1239

where A is the area of the plate. Further, following the accepted convention in the analysis of laminate structures [15], other variables are de!ned as follows: {e} , {0 k #0 }T ;   A B 0    [D] ,  B D 0 ; 0 0 As

(11)

(12)

{0 } , {x y #xy }T ;

(13)

{k} , {kxx kyy kxy + kyx }T ;

(14)

{#0 } , {#xz #yz }T ;

(15)

n


(Aij ; Bij ; Dij ) , Asij ,

n
 k=1

k=1 zk

zk − 1

zk

zk − 1

[QG ij ]k (1; z; z 2 ) d z

[QG sij ]k d z

i; j = 4; 5;

i; j = 1; 2; 6;

(16) (17)

where n is the total number of laminae, and zk −1 , zk , respectively, denote the coordinates, with respect to the plate midplane, of the bottom and the top of the kth layer. The thermal load and moment vectors are de!ned as n
zk  G } JT [Q]{ G d z; (18) {N } , k=1

{M } ,

zk − 1

n
 k=1

zk

zk − 1

G } JTz[Q]{ G d z:

(19)

Although this study is restricted to laminates without membrane–bending coupling, the membrane– bending coupling matrix [B] and the corresponding thermal moments {M } are retained in the above formulation for completeness. 2.3. Free vibration problem The free vibration problem is viewed as a linear perturbation about a non-linear reference state which results from the thermal calculation. This permits the incorporation of the e'ects of the thermal residual stresses. The governing equations of motion of the plate are derived using Hamilton’s principle. The Lagrangean L is given as L = T − ;

(20)

1240

X. Wang et al. / Finite Elements in Analysis and Design 40 (2004) 1233 – 1257

where T is the kinetic energy and  is the potential energy of the plate. The latter is further decomposed into the component due to elastic strain energy U and that due to the work done by the applied forces W. The stresses at any point in the plate { G } may be written as G G + GR }; { G } = [Q]{

(21)

where {} G is the total strain vector due to the vibrational motion and {GR } is thermal strain vector. G is the transformed reduced sti'ness matrix of the particular lamina of interest. The matrix [Q] The potential energy of the plate can be calculated from

1 G T {u} =U−W= {G + GR }T { G } dv − {f} G ds; (22) 2 v s G where v is the entire volume of the plate and s is its surface; {u} G are the displacements, and {f} is the vector of edge traction. The fact that there are no edge/surface tractions in a free vibration problem makes the work contribution W ≡ 0, and hence  = U. Eqs. (2) and (21) are now substituted into Eq. (22). The higher-order terms in the resulting expression are ignored and use is made of the fact that { G R } (i.e., the thermal residual stresses) correspond to a state of equilibrium to obtain the following:

1 G G L } dv + {G N }T { G R } dv: U= {G L }[Q]{ (23) 2 v v The above strain energy is integrated through the thickness (i.e., with respect to z) to obtain the following:

1 U= {e}T [D]{e} dA + {eG }T [DGR ]{eG }} dA: (24) 2 A A The second integral in the expression is a consequence of the coupling between linear and non-linear terms and it is denoted by the subscript G. By denoting the linear terms computed from the thermal analysis by the superscript R, the matrix [DGR ] can be written as  R  R R Nxx Nxy 0 0 0 0 MxxR Mxy 0 0 QxR 0  R  R R R 0 0 0 0 Mxy Myy 0 0 QyR 0   Nxy Nyy   R R  0 0 NxxR Nxy 0 0 0 0 MxxR Mxy 0 QxR      R R R R 0 Nxy Nyy 0 0 0 0 Mxy Myy 0 QyR   0   R  0 0 0 0 NxxR Nxy 0 0 0 0 0 0      R R  0 0 0 0 Nxy Nyy 0 0 0 0 0 0  R  ; (25) [DG ] =  R R 0 0 0 0 LRxx LRxy 0 0 TxR 0   Mxx Mxy   R  R  Mxy Myy 0 0 0 0 LRxy LRyy 0 0 TyR 0     0 R 0 0 0 0 LRxx LRxy 0 TxR  0 MxxR Mxy     R R R R R  0 0 Mxy Myy 0 0 0 0 Lxy Lyy 0 Ty     R  R R R Qy 0 0 0 0 Tx Ty 0 0 0 0   Qx QyR 0 0 0 0 TxR TyR 0 0 0 0 QxR

X. Wang et al. / Finite Elements in Analysis and Design 40 (2004) 1233 – 1257

where (Nijs ; Mijs ; Lsij )

,

n
 k=1

(Qis ; Tis )

n


,

k=1

zk

zk − 1

zk

zk − 1

G ij (1; z; z 2 ) d z

G iz (1; z) d z

for ij = xx; yy; xy;

1241

(26)

for i = x; y:

(27)

The higher-order moments about the midplane of the plate are denoted as LRij (ij = xx, yy or xy). While these parameters assume non-zero values for the case of nonsymmetric laminates and for laminates with non-planar loads, they also, in general, assume non-zero values for laminated (symmetric or otherwise) plates under in-plane loads. Further, the kinematic hypotheses imply that the stress resultants TiR are identically zero for symmetric laminates but may be non-zero for non-symmetric laminates. The kinetic energy of the plate T is  T   uG˙  uG˙      
     1 ˙ ˙ (28) T= vG vG  dv;    2 v      ˙  ˙ wG wG where  is the density. Since the mass distribution is uniform over the plate cross-section, the kinetic energy can be written in terms of the midplane displacements as  T    u˙  u˙   m 0 0 0 0                        v ˙ v ˙       0 m 0 0 0   
       1   w˙ T= (29)  0 0 m 0 0  w˙  dA;     2 A           ˙x  ˙x        0 0 0 I 0              ˙    ˙ 0 0 0 0 I y y where


(m; I ) =

zk

z0

(1; z 2 ) d z:

(30)

2.4. Sensitivity of eigenvalues Following the study by Liu et al. [16], the free vibration problem of the plate with a sti'ener, as depicted in Fig. 1, is formulated based on the stationarity of the appropriate Rayleigh quotient in conjunction with a kinematic compatibility constraint [17]. The free vibration problem ignores in-plane deformation. Hence an harmonic solution can be assumed for the relevant !eld variables such that w(x; y; t) = W (x; y)ei,t ;

x (x; y; t)

= -x (x; y)ei,t ;

and

y (x; y; t)

= -y (x; y)ei,t :

1242

X. Wang et al. / Finite Elements in Analysis and Design 40 (2004) 1233 – 1257 Y

S

L α (X0 ,Y0)

0

X

Fig. 1. Schematic of plate with a symmetric-layered sti'ener.

The above permits a rede!nition of the potential energy as U = 12 U ei,t where
 k T  D∗ 0   k  U [W; x ; y ] = dA; 0 A # # A 0 s 0   A B ∗ D = ; B D where k, #0 , A, B, D and As are the same as those in Eqs. (14)–(17). Similarly, the kinetic energy is rede!ned as T = 12 ,2 T ei2,t where

tp3 2 2 T [W; x ; y ] = tp W dA + ( x + y2 ) dA; A A 12 where tp is the plate thickness and  is the mass density of the plate. The Rayleigh quotient is given as U [W; x ; y ] + Ub [Wb ; b ; b ] ; R= T [W; x ; y ] + Tb [Wb ; b ; b ]

(31) (32)

(33)

(34)

where the strain energy of the plate is denoted as U [W; x ; y ].   The terms Ub [Wb ; b; b ] and Tb [Wb ; b ; b ] are from the strain energy 12 Ub ei2,t and kinetic energy of the sti'ener 12 Tb ei2,t , respectively. The former may be written as     
  d b 2 db 2 dWb 2 EIb ds; (35) Ub [Wb ; b ; b ] = + GAb + GJb b− ds ds ds s where Wb is the vertical displacement along the local axis of the sti'ener. The cross-sectional rotation and the twist angle of the cross-section of the sti'ener are, respectively, denoted as b and b . The modulus of elasticity is E, Ib is the second moment of the cross-sectional area, G is the shear modulus and Ab is the cross-sectional area. The polar moment of inertia of the cross-section of the sti'ener is denoted as Jb .

X. Wang et al. / Finite Elements in Analysis and Design 40 (2004) 1233 – 1257

The latter, that is the kinetic energy term, Tb [Wb ; b ; b ], may be expressed as
Tb [Wb ; b ; b ] = (b Ab Wb2 + b Ib b2 + b Jb 2b ) ds; s

1243

(36)

where b is the sti'ener mass density. The kinematic compatibility constraints require the equality of the deformation on the common points on the intersecting surface (a line to be speci!c) of the sti'ener and the plate. These constraints may be written as     W (x0 + s(cos ); y0 + s(sin )) Wb (s)      x (x0 + s(cos ); y0 + s(sin ))  − P  b (s)  = 0; (37)     b (s) y (x0 + s(cos ); y0 + s(sin )) where P is a coordinate transformation matrix which relates the local coordinate system of the sti'ener to the global system attached to the plate. The location of the sti'ener on the plate is completely speci!ed by the cartesian coordinates (x0 ; y0 ) and the angle  measured from the horizontal axis in a counterclockwise direction. The stationarity of the Rayleigh quotient R is subject to the kinematic compatibility constraints as de!ned in Eq. (37). Hence the variables W , x , y , Wb , b , and b are not independent variables. In order to facilitate the derivation of the eigenvalue sensitivity, the stationary problem of the generalized Rayleigh quotient is introduced in the form U [W; x ; y ] + Ub [Wb ; b ; b ] − 2L[ ; W; x ; y ; Wb ; b ; b ] ; (38) RG = T [W; x ; y ] + Tb [Wb ; b ; b ] where L is de!ned as     W (x0 + s(cos ); y0 + s(sin )) Wb (s)
    T  x (x0 + s(cos ); y0 + s(sin ))  − P  b (s)  ds L, (s)      s

y (x0

+ s(cos ); y0 + s(sin ))

(39)

b (s)

and the vector of Lagrange multipliers is denoted as (s). The generalized Rayleigh quotient RG has seven independent variables ; W; x ; y ; Wb ; b ; b . It is known from variational principle states that the entire set of Euler equations for generalized Rayleigh quotient RG is the same as the Euler equations of the ordinary Rayleigh quotient R in conjunction with the kinematic compatibility constant Eq. (37). The implication, therefore, is that the stationarity of RG and the stationarity of R are equivalent in that they give the same set of eigenfunctions. However, they are di'erent with respect to the constraint condition and the number of independent arguments. Another view to understanding the two quotients is to consider the stationarity of RG as an unconstrained problem with seven variables while the stationarity of R is a constrained problem. By denoting the vector of Lagrange multiplier that makes the generalized Rayleigh quotient stationary as ˆ and the corresponding eigenfunctions as Wˆ , ˆ x , ˆ y , Wˆ b , ˆ b , and ˆ b , then U [Wˆ ; ˆ x ; ˆ y ] + Ub [Wˆb ; ˆ b ; ˆ b ] − 2L[ ˆ; Wˆ ; ˆ x ; ˆ y ; Wˆ b ; ˆ b ; ˆ b ] = : (40) T [Wˆ ; ˆ x ; ˆ y ] + Tb [Wˆ b ; ˆ b ; ˆ b ]

1244

X. Wang et al. / Finite Elements in Analysis and Design 40 (2004) 1233 – 1257

The relationship between  and the position of the sti'ener, which is determined by (x0 ; y0 ; ) is provided by Eqs. (37) and (40). Two important deductions can be inferred. The !rst is that ˆ, Wˆ , ˆ x , ˆ y , Wˆ b , ˆ b , ˆ b change with any change in the coordinates of the sti'ener (i.e., any change to either or a combination of x0 ; y0 and ), but variations in 2 ˆ, 2Wˆ , 2 ˆ x , 2 ˆ y , 2Wˆ b , 2 ˆ b , 2ˆ b will not result in any change to  to the !rst order because of the stationarity of RG . Hence, a di'erentiation of Eq. (40) with respect to the location of the sti'ener (i.e., (x0 ; y0 ; )) will not involve di'erentiating ˆ, Wˆ , ˆ x , ˆ y , Wˆ b , ˆ b , ˆ b since they are independent of the sti'ener location. The second deduction is that though a perturbation of the sti'ener changes quotient L as de!ned by Eq. (39), the perturbation does not result in any !rst order change in U , T , Ub , Tb . Based on the aforementioned observations, the di'erentiation of Eq. (40) with respect to any coordinate of the sti'ener location (x0 , y0 or ) yields −2(9L=9p) 9 = ˆ ˆ 9p T [W ; x ; ˆ y ] + Tb [Wˆ b ; ˆ b ; ˆ b ]

(p = x0 ; y0 ; ):

(41)

By normalizing the eigenfunctions such that T [Wˆ ; ˆ x ; ˆ y ] + Tb [Wˆ b ; ˆ b ; ˆ b ] = 1

(42)

a simpli!ed form of the above di'erentiation, i.e., Eq. (41), may be written as 9 9L = −2 : 9p 9p

(43)

Following some basic algebraic manipulations, the eigenvalue sensitivity with respect to the location of the sti'ener (as de!ned by the coordinates x0 , y0 and ) may be expressed as:  ˆ 


d 9L = −2 = −2 d x0 9x0

s

  ˆT (s)   

d 9L = −2 = −2 dy0 9y0


s

  ˆ (s)    T

 d 9L = −2 = −2 d d


s

  ˆT (s)  

9W 9x

9 ˆx 9x 9 ˆy 9x 9Wˆ 9y 9 ˆx 9y 9 ˆy 9y 9Wˆ 9n 9 ˆx 9n 9 ˆy 9n

   ds; 

(44)

    ds;   

(45)

  dP     ˆ b (s)  ds: −   d   ˆ b (s) 

Wˆ b (s)

(46)

The integrals are evaluated along the sti'ener, and [9Wˆ =9n; 9 ˆ x =9n; 9 ˆ y =9n]T are the directional derivative of the plate eigenfunctions at the nodes along the interface of the sti'ener and plate in a direction normal to the sti'ener. The normal is positive in the direction of increasing angle .

X. Wang et al. / Finite Elements in Analysis and Design 40 (2004) 1233 – 1257

1245

The vector of Lagrange multipliers ˆ(s) can be determined by equating the !rst variation of the right-hand side of Eq. (40) to zero and then selecting ˆ(s) such that the coeIcients of the !rst variation of Wˆ b , ˆ b , ˆ b vanish. The result is      b ˆ − dsd GAb ˆ b − dW + A W b b ds         ˆ   d d ˆb d W b ˆ ˆ (47) ˆ(s) = P  ds EIb ds − GAb b − ds + Ib b  :       ˆ b d GJb dds + Jb ˆ b ds The components of ˆ(s) can be interpreted as the internal reaction forces acting at the interface of the plate and sti'ener. This leads to the conclusion that the eigenvalue sensitivity is proportional to these forces and to the derivative of eigenfunctions.

3. Finite element method formulation The derivations of the preceding section are now reformulated in a form that is amenable to the application of !nite element analysis. 3.1. Basic problem The !nite element formulation to solve the thermal, free vibration and eigenvalue sensitivity problems is based on a 16-node isoparametric, bi-cubic Lagrangean element. The geometry of the element is shown in Fig. 2. Each node of the element has !ve degrees of freedom (u, v, w, x , and y ). The displacements at an arbitrary point within the domain of the eth element is given by {} = [u v w

T x y]

= [N ]5×80 {e }1×80 ;

(48)

where [N ] is the matrix of interpolation functions. The vector of nodal displacements is denoted as {e } and it is of the form {e }T = [u1 v1 w1

x1 y1

· · · u16 v16 w16

x16

y16 ]:

(49)

The linear strain vector {e} for an element can be expressed as {e} = [B]8×80 {e };

(50)

where [B] is the linear strain–displacement transformation matrix. Also, the non-linear strain vector {eG } for an element can be expressed as {eG } = [BG ]12×80 {e }; where matrix [BG ] is the non-linear strain–displacement transformation matrix.

(51)

1246

X. Wang et al. / Finite Elements in Analysis and Design 40 (2004) 1233 – 1257

η 13

14

9

10

15

(-1, 1)

11

16 (1, 1)

12

ξ 5

1 (-1,-1)

6

7

8

2

3

4 (1, -1)

Fig. 2. Bi-cubic Lagrangian isoparametric element.

The matrix of element sti'ness and of geometric sti'ness due to thermal residual stresses are, respectively, given as
(52) [Ke ] = [B]T [D][B] dA; A

and R [KGe ]=


A

[BG ]T [DG ][BG ] dA:

(53)

R The latter, i.e., [KGe ], is computed by assuming a thermal loading that corresponds to the temperature di'erence.

3.2. The thermal problem The thermal problem is solved by minimising the total potential energy, i.e., Eq. (10), and computing the displacement !elds. Thus    N 

   2R = {2e}T [D]{e} dA − {2e}T M dA = 0: (54)   A A   0   The use of the linear strain vector from Eq. (50) yields an equation of the form [Ke ]{e } − {fe } = 0; where the sti'ness matrix of the element
[Ke ] , [B]T [D][B] dA A

(55) (56)

X. Wang et al. / Finite Elements in Analysis and Design 40 (2004) 1233 – 1257

and the force vector of the element    N  
  T M {fe } , [B] dA:   A   0  

1247

(57)

The global governing equations are assembled in the usual !nite element method convention [18] and the global nodal displacements are computed. A postprocess operation is then performed to determined the elemental strains and stresses. Thereafter, the eight linear stress resultants (NxR ; NyR , R R Nxy ; MxR , MyR ; Mxy , QxR ; QyR ) and !ve non-linear stress resultants (LRx ; LRy ; LRxy ; TxR ; TyR ) of each element are evaluated at the Gauss integration points. These are then used to develop the global non-linear strain–displacement transformation matrix [BG ]. 3.3. The free vibration problem The element mass matrix [Me ] is derived from Eqs. (28) and (48) and it is expressed as
[Me ] = [N ]T [J ][N ] dA; A

where the inertia  m  0   [J ] =  0  0  0

matrix [J ] is given as  0 0 0 0  m 0 0 0   0 m 0 0:  0 0 I 0  0

0

0

(58)

(59)

I

Thus the plate (or global) kinetic energy can be written as N 1 T  [Me ]e ; TN = 2 e=1 e where N is the total number of the elements in the model. Similarly, the expression for the strain energy may be written as N N 1 T 1 T R e [Ke ]e +  [K ]e ; UN = 2 e=1 2 e=1 e Ge

(60)

(61)

R where [KGe ] is due to thermal residual stresses. The above energy expressions are used in the Hamilton’s principle and variations are observed over the free variables to yield the system governing equation. By assuming a harmonic expansion of the free variables the resultant eigenvalue problem becomes

[[K] + [KGR ] − [M ]]{} = 0;

(62)

where [M ], [K], [KGR ] are the global matrices for the plate. The natural frequencies and corresponding vibration modes are solved using subspace iteration algorithm [18].

1248

X. Wang et al. / Finite Elements in Analysis and Design 40 (2004) 1233 – 1257

3.4. Sensitivity of eigenvalues The vector of Lagrange multipliers, Eq. (47), can be evaluated using an equivalent !nite element version of the form G b ]UG b ; G = [ − Kb + M

(63)

G UG b are the discretized versions of the Lagrange multipliers, the eigenvalue and the where G , , eigenfunctions, and Kb , Mb are the corresponding sti'ness and mass matrices of the sti'ener. The above equation is now used to write the corresponding !nite element formulation equivalence of Eqs. (44)–(46), which, respectively, are G d G G b ) dU ; = −2UG Tb (−Kb + M d x0 dx

(64)

where d UG =d x is the derivative of the plate eigenvector at the interface nodes between the plate and the sti'ener, G d G G b ) dU ; = −2UG Tb (−Kb + M dy0 dy

(65)

where d UG =dy is the derivative of the plate eigenvector at the interface nodes between the plate and the sti'ener,      G d G G b dP G b )Is 9U + 2UG T −Kb dP + M Is UG b ; (66) = −2UG Tb (−Kb + M b d 9n d d where Is is a diagonal matrix with diagonal entries [si ] which denote the s-coordinate of the nodes along the sti'ener, and dp=d is a block-diagonal matrix with its block as [dP=d] which is a 5 by 5 matrix, and 9UG =9n is the directional derivative of the plate eigenvectors at the nodes along the interface in a direction normal to the sti'ener. The normal is assumed to be positive in the direction of increasing angle . The [si ] and [dP=d] matrices are, respectively, given as   xi 0 0 0 0    0 yi 0 0 0      0 0 1 0 0 [si ] =  (67) ;    0 0 0 1 0   0 





dP  = d

0

0

0

−sin  −cos  cos  −sin 

03 × 2

1 

 02×3  : 03× 3

(68)

The directional derivative of the plate eigenvectors 9UG =9n can be computed in terms of derivatives of the plate eigenvector at the interface nodes between the plate and the sti'ener with respect to the

X. Wang et al. / Finite Elements in Analysis and Design 40 (2004) 1233 – 1257

coordinates x and y. Speci!cally, d UG d UG 9UG =− sin  + cos : 9n dx dy

1249

(69)

4. Optimisation of natural frequencies It is considered prudent at this juncture to elaborate on the process of optimising the natural frequencies of the plate before presenting the results of the numerical simulations in Section 5. Four di'erent types of sti'ened composite plates are considered in order to illustrate the e'ect of the di'erent geometric arrangements of the sti'eners on the thermal residual stress resultants and natural frequencies of the plate. These are (1) the plate with a single continuous sti'ener, (2) the plate with continuous vertical and horizontal sti'eners, (3) the plate with continuous diagonal sti'eners, and (4) the plate with discrete sti'eners. The objective function is the maximisation of the !rst natural frequency of the plate and the design variables are the location parameters of the sti'eners. The location of each sti'ener is identi!ed by its cartesian coordinates x0 , y0 , and, where applicable, by angle . In particular, the plate with a single continuous sti'ener permits the arbitrarily location of the sti'ener in the domain of the plate. Thus, the location of a sti'ener is determined by three positional parameters x0 , y0 and . The location of the sti'eners in a plate with both vertical and horizontal sti'eners is determined from the respective x0 or y0 coordinate. Each of the sti'eners can move in the x or the y direction of the plate, depending on whether it is originally vertical or horizontal. The location of each sti'ener that comprise the plate with diagonal sti'eners is identi!ed by the x0 and y0 coordinates. For the plate with discrete sti'eners, the location of each sti'ener can be changed arbitrarily along x and y directions in the domain of the plate. The present optimisation problem falls into the category of unconstrained function optimisation [19]. It is solved using MMA [10,11], a globally convergent iterative optimisation technique— a convex subproblem that approximates the original problem is generated and solved during each iteration. The MMA technique requires the derivatives of the objective function with respect to the design variables. These are determined from the expressions derived in Sections 2.4 and 3.4. 5. Numerical simulation The 20 ply graphite/epoxy T300/5208 laminated plates used in the numerical simulations are each square shaped and are either [0; 90]s or [ ± 45]s . Each side of the plate is 360 mm in length and the sti'eners are symmetrically located on the upper and lower surfaces of the base laminated plate. The material properties are assumed to be temperature invariant and are tabulated in Table 1. Unless otherwise stated, the width of sti'eners used in the simulations is 6 mm. The boundary conditions of the thermal problem are only such as to eliminate rigid-body motion of the plate. This is because the thermal residual stress resultants are strictly a consequence of the mismatch in the thermal coeIcient of expansions of the sti'eners and the base laminated plate. The free vibration problem, on the other hand, is based on simply supported conditions on the edges of the plate. The application of these boundary conditions does not introduce additional residual stresses.

1250

X. Wang et al. / Finite Elements in Analysis and Design 40 (2004) 1233 – 1257

Table 1 T300/5208 Graphite-Epoxy material properties Properties

Values

Longitudinal modulus of elasticity, E1 Transversal modulus of elasticity, E2 In-plane Poisson’s ratio, V 12 In-plane shear modulus, G12 Transversal shear modulus, G13 Transversal shear modulus, G23 Longitudinal thermal expansion coeIcient, 1 Transversal thermal expansion coeIcient, 2 Layer thickness, t Mass density, 

154:5 GPa 11:13 GPa 0.304 6:98 GPa 6:98 GPa 3:36 GPa −0:17E − 06◦ C−1 2:31E − 05◦ C−1 0:15 mm 1:56 g=cm3

While the temperature di'erence used in the thermal problem is given as a positive number it is, by de!nition, the di'erence between the lower operating or room temperature and the elevated (cure or consolidated) temperature. The simulations using each of the four types of sti'ened laminated plates are individually presented in the following subsections. Perhaps it is worth mentioning that repeated eigenvalues are not encountered in the these simulations. The sensitivity of the eigenvalues would be determined from subeigenproblems [20] in the event of repeated eigenvalues. 5.1. Type A: single continuous stiener The sti'ener used in this analysis has a length of 120 mm. The symmetric arrangement of the sti'ener on the base laminate is ensured by placing ten plies at the top of the base laminate and aIxing ten plies directly below at the bottom of the base laminate. The basic sti'ener size remains constant but the location, as identi!ed by the parameters x0 , y0 and , varies over the domain of the base laminate. Further, taking into cognizance the periodicity of trigonometric sine and cosine functions, the directional angle  is bounded between −90◦ and 90◦ . The number of elements used in the analysis is 38. Also worthy of mention is the fact that the !bres of the laminates of the sti'ener are aligned with the local sti'ener x-axis (i.e., longitudinal axis). The optimum location of the sti'ener is observed to be a corner of the plate and it is at an angle of −45◦ with the x-axis of the base laminate. This is illustrated in Fig. 3. The variation of the optimum !rst natural frequencies of the plate with respect to a change in temperature is depicted in Fig. 4. The two base laminate plates considered are [0; 90]s and [ ± 45]s and they are, henceforth, identi!ed in this and identical !gures as C1 and C2, respectively. The optimum !rst natural frequency from the former when there is no change in temperature is 24.467 and 30:619 Hz for the latter. These values increase slightly with increasing change in temperature. For example, a 200◦ C change in temperature yields optimum !rst natural frequencies of 24.680 and 30:642 Hz, respectively. The [0; 90]s ply base laminate shows the faster rate of change in the natural frequency.

X. Wang et al. / Finite Elements in Analysis and Design 40 (2004) 1233 – 1257

1251

350

300

250

Y

200

150

100

50

0 0

100

200

300

X

Fig. 3. Optimum sti'ener location (a single continuous sti'ener scenario).

36 C1 C2

35 34

Natural Frequency (Hz)

33 32 31 30 29 28 27 26 25 24 23 22 21 0

50

100

150

200

Temperature Difference (°C)

Fig. 4. Variation of the 1st natural frequencies with respect to temperature change (for a single continuous sti'ener).

5.2. Type B: base plate with vertically and horizontally located continuous stieners Each sti'ener in the Type B scenario runs from one side of the plate to the opposite side. The arrangement of the sti'eners on the base laminate is such that symmetry is maintained about the midplane of the base laminate. Unlike Type A arrangement, the positions of the sti'eners in this scenario are completely de!ned by the two cartesian coordinates x0 and y0 . Figs. 5 and 6, respectively, illustrate the optimum sti'ener arrangement and the variation of the 1st natural frequency with change in temperature. The sti'eners are observed to cluster around the middle of the base plate. This increases the plate sti'ness which, in turn, increases the !rst natural

1252

X. Wang et al. / Finite Elements in Analysis and Design 40 (2004) 1233 – 1257 350

300

250

Y

200

150

100

50

0

0

100

200

300

X

Fig. 5. Optimum sti'ener location (case of vertically and horizontally located continuous sti'eners).

38 37 C1 C2

36

Natural Frequency (Hz)

35 34 33 32 31 30 29 28 27 26 25 24 23 0

50

100

150

200

Temperature Difference (°C)

Fig. 6. Variation of the 1st natural frequencies with respect to temperature change (for vertically and horizontally located continuous sti'eners).

frequency. The !rst natural frequencies decrease with increasing change in temperature for both base laminate ply sequences. The above is a further indication that while the arrangement of the sti'eners can be tailored to enhance vibrational performance of the entire system/structure via the introduction of residual stresses in the base plate, it is possible, however, to select arrangements that, in fact, deteriorate the vibrational performance.

X. Wang et al. / Finite Elements in Analysis and Design 40 (2004) 1233 – 1257

1253

350

300

250

Y

200

150

100

50

0

0

100

200

300

X

Fig. 7. Optimum sti'ener location (case of diagonal located continuous sti'eners).

39 C1 C2

Natural Frequency (Hz)

38 37 36 35 34 33 32 31 30

0

50

100

150

200

Temperature Difference(°C)

Fig. 8. Variation of the 1st natural frequencies with respect to temperature change(case of diagonal located continuous sti'eners).

5.3. Type C: base plate with diagonally located continuous stieners This scenario uses 32 elements and considers four diagonally located continuous sti'eners. Given that the sti'eners are de!ned in the domain of the base plate, the length of each is dependent upon its position on the base plate. The optimum arrangement of the sti'eners, like the observation in the preceding scenario, is such that the sti'eners cluster around the centre (see Fig. 7). For this case study, however, increments to the temperature change act to increase the !rst natural frequency as can be inferred from Fig. 8,

1254

X. Wang et al. / Finite Elements in Analysis and Design 40 (2004) 1233 – 1257

Fig. 9. Optimum sti'ener locations {[0; 90]ps ; [0; 90]ss }.

with the greater in7uence observed with the [ ± 45]s laminate. This implies that the residual tensile thermal stress in the [ ± 45]s laminate is higher stress than that in the [ ± 90]s laminate. 5.4. Type D: base plate with discrete stieners Here, while the sti'eners each have 20 plies, their width and length are each 18 mm. The sti'eners are vertically or horizontally positioned on the base laminate. Thus, and in light of the square shape of the sti'eners, each sti'ener is identi!ed by the cartesian coordinate of its centroid, i.e., (x0 ; y0 ). The determining factor is the relationship between the ply-angles of the base laminate and those of the sti'eners. This is demonstrated via the use of four sets of ply-angles combinations: {[0; 90]ps ; [0; 90]ss }, {[ ± 45]ps ; [0; 90]ss }, {[0; 90]ps ; [ ± 45]ss }, {[ ± 45]ps ; [ ± 45]ss }, where the superscript ‘p’ refers to base laminate or plate and ‘s’ refers to sti'ener. The corresponding optimal sti'ener locations are respectively depicted in Figs. 9–12. The optimum 1st natural frequencies are 29.518, 32.589, 32.135 and 37:060 Hz respectively. The discrete sti'ener results are reminiscent of the type of results generated by structural topology optimisation [21]. 6. Summary The studies in Refs. [3–6] show that thermal residual stresses can be used to improve both buckling and vibrational performance of (sti'ened) composite structures. The current e'ort extends these studies by exploring the problem of the systematic arrangement of the sti'eners in order to obtain the best possible enhancement in performance. The thermal residual stresses are solely a consequence of thermal coeIcient of expansion mismatches.

X. Wang et al. / Finite Elements in Analysis and Design 40 (2004) 1233 – 1257

1255

Fig. 10. Optimum sti'ener locations {[45; −45]ps ; [0; 90]ss }.

Fig. 11. Optimum sti'ener locations {[0; 90]ps ; [45; −45]ss }.

The factors that determine whether thermal residual stresses are bene!cial or detrimental to the performance of a sti'ened composite plate include the geometric arrangement of the sti'eners, the ply-angles of the basic laminate and the sti'eners, material properties, and the boundary conditions. The !rst two factors are the primary concern of the present study where the objective is the maximisation of the !rst natural frequency of the plate.

1256

X. Wang et al. / Finite Elements in Analysis and Design 40 (2004) 1233 – 1257

Fig. 12. Optimum sti'ener locations {[45; −45]ps ; [45; −45]ss }.

The sti'eners considered are classi!ed as either continuous or discrete. The results of the numerical simulations indicate that thermal residual stresses in7uence the magnitude of the optimum !rst natural frequencies of only the plates with continuous sti'eners. However, the optimum locations of the sti'eners in both plate types are invariant to thermal residual stresses. Both the arrangement of the sti'eners and the magnitude of the optimum !rst natural frequency are strongly in7uenced by the ply-angles of the basic laminates and those of the sti'eners. For example, the use of [ ± 45]s basic laminate or sti'eners always results in higher natural frequency when compared to [0; 90]s . Finally, the idealistic nature of the sti'ener arrangements examined in this study is acknowledged. However, they are so selected to highlight the e'ects of thermal residual stresses on the vibration performance of a composite sti'ened plate. A more realistic design process requires a comprehensive analysis of the overall structural mechanical behaviour. Thus the process of enhancing both buckling and vibration performance would be done in conjunction with other structural performance considerations such as static strength and damage tolerance. Acknowledgements This work was partially supported by the Natural Sciences and Engineering Research Council of Canada (NSERC). References [1] K.S. Sai Ram, P.K. Sinha, Hygrothermal e'ects on the buckling of composite laminated plates, Composite Struct. 21 (1992) 233–247.

X. Wang et al. / Finite Elements in Analysis and Design 40 (2004) 1233 – 1257

1257

[2] K.S. Sai Ram, P.K. Sinha, Hygrothermal e'ects on the free vibration of laminated composite plates, ASME J. Appl. Mech. 18 (1992) 31–38. [3] J.P. Foldager, J.S. Hansen, N. Olho', Optimization of the buckling load for composite structures taking thermal e'ects into account, Struct. Multdisciplinary Optim. 21 (2001) 14–31. [4] A.R. de Faria, J.S. Hansen, Optimal buckling loads of nonuniform composite plates with thermal residual stresses, ASME J. Appl. Mech. 66 (1999) 388–395. [5] S.F.M. de Almeida, J.S. Hansen, Enhanced elastic buckling loads of composite plates with tailored residual Stresses, ASME J. Appl. Mech. 64 (1997) 772–780. [6] S.F.M. de Almeida, J.S. Hansen, Natural frequencies of composite plates with tailored thermal residual-stresses, Int. J. Solids Struct. 36 (1999) 3517–3539. [7] B. Grall, Z. Gurdal, Structural analysis and design of geodesically sti'ened reference: panels with variable sti'ener distribution (NASA-CR-190608 1992). [8] N. Jaunky, N.F. Knight, Jr., Optimal design of grid-sti'ened composite panels with global and local buckling analysis, Proceedings of the 37th AIAA/ASME/ASCE/ASC Structures, Structural Dynamics and Materials Conference, Salt Lake City, Part 4, 1976, pp. 2315 –2325. [9] M. Kataoka-Filho, Optimization of nonhomogeneous facesheets in composite sandwich plates, Ph.D. Dissertation, University of Toronto, 1997. [10] K. Svanberg, The Method of Moving Asymptotes (MMA) with some extensions, Optimization Large Struct. Systems 1 (1991) 555–566. [11] K. Svanberg, A globally convergent version of MMA without linesearch, in: N. Olho', G. Rozvany (Eds.), Proceedings of the 1st World Congress of Structural and Multidisciplinary Optimization, Pergamon Press, Elmsford, NY, 1995, pp. 9 –16. [12] D.S. Adams, D.E. Bowles, C.T. Herakovich, Thermally induced transverse cracking in graphite-epoxy cross-ply laminates, in: G.S. Springer (Ed.), Environmental E'ects on Composite Materials, Vol. 3, Technomic Publ. Co. Inc., Lancaster, PA, 1988, pp. 247–274. [13] R.D. Mindlin, In7uence of rotary inertia and shear on 7exural motions of isotropic elastic plates, ASME J. Appl. Mech. 18 (1951) 31–38. [14] V.V. Novozhilov, Foundations of the Non-linear Theory of Elasticity, Graylock, Rochester, New York, 1953. [15] R.M. Jones, Mechanics of Composite Materials, Hemisphere Publishing Co., Washington DC, 1975. [16] Z.S. Liu, J.S. Hansen, D.C.D. Oguamanam, Eigenvalue sensitivity analysis of sti'ened plates with respect to the location of sti'eners, Struct. Optim. 16 (1998) 155–161. [17] H.C. Hu, Variational Principles of Theory of Elasticity with Applications, Science Press, Beijing, P. R. China, 1984. [18] K.-J. Bathe, Finite Element Procedures, Prentice-Hall, Englewood Cli's, NJ, 1996. [19] G.N. Vanderplatts, Numerical Optimization Techniques for Engineering Design with Applications, McGraw-Hill Inc., New York, 1984. [20] D.C.D. Oguamanam, Z.S. Liu, J.S. Hansen, Natural frequency sensitivity analysis with respect to lumped mass location, AIAA J. 37 (1999) 928–932. [21] M.P. Bendsoe, N. Kikuchi, Generating optimal topologies in structural design using a homogenization motion, Comput. Methods Appl. Mech. Eng. 71 (1988) 197–224.