- Email: [email protected]
Variational asymptotic modeling of the multilayer functionally graded cylindrical shells
Composite Structures 94 (2012) 966–976
Contents lists available at SciVerse ScienceDirect
Composite Structures journal homepage: www.elsevier.com/locate/compstruct
Variational asymptotic modeling of the multilayer functionally graded cylindrical shells q Zhong Yifeng a,b,⇑, Chen Lei a,b, Wenbin Yu c a
College of Civil Engineering, Chongqing University, Chongqing, PR China Key Laboratory of New Technology of Construction of Cities in Mountain Area, Chongqing University, Ministry of Education, Chongqing 400045, PR China c Utah State University, Logan, UT 84322-4130, USA b
a r t i c l e
i n f o
Article history: Available online 13 October 2011 Keywords: Asymptotic series Functionally graded material Cylindrical shell Variational asymptotic method Simplified model
a b s t r a c t An efficient high-fidelity shell model is developed for heterogeneous multilayer cylindrical shells made of functionally graded material by using the variational asymptotic method (VAM). Taking advantage of the smallness parameters inherent in the shell structure, the VAM is applied to rigorously decouple the 3-D, anisotropic elasticity problem into a 1-D through-the-thickness analysis and a 2-D shell analysis. The through-the-thickness analysis servers as a link between the original 3-D analysis and the shell analysis by providing a constitutive model for the shell analysis and recovering the 3-D field variables in terms of global responses calculated by the shell analysis. The present model is valid for large displacements and global rotations and can capture all the geometric nonlinearity of a shell when the strains are small. A cylindrical bending example of a homogeneous substrate with a thin SiC–Al functionally graded coating under sinusoidal pressure on the top surface is used to validate this model. Ó 2011 Elsevier Ltd. All rights reserved.
1. Introduction Functionally graded materials (FGM) are microscopically inhomogeneous composites usually made from a mixture of metals and ceramics. By gradually varying the volume fraction of consistent material, its material properties exhibit a smooth and continuous change from one surface to another. The surface with high ceramic constituents can provide superior thermal-resistance for high temperature environments while the surface with high metal constituents offers strong mechanical performance. This kind of material arrangement reduces the risk of catastrophic fracture and mitigates thermal stress concentration under extreme environments. The advantage of FGM has been used to address a large variety of application fields, such as graded thermoelectrics and dielectrics, piezoelectrically graded materials for ultrasonic transducers, and tungsten–copper composites for high current connectors and diverter plates, to name but a few [1,2]. Recently, the concept of FGM is actively explored in the multilayered design of thermal coatings as well as sandwich panels to overcome the mismatch of the thermomechanical properties between the coating and the substrate or between the surface panels and the core [3–5].
q
This paper is contributed by all authors.
⇑ Corresponding author at: College of Civil Engineering, Chongqing University, Chongqing, PR China. Tel.: +86 023 65123292. E-mail address: [email protected] (Z. Yifeng). 0263-8223/$ - see front matter Ó 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.compstruct.2011.10.005
To use FGM effectively, an efficient and accurate model for structures made of such materials needs to be developed. These structures are characterized by one of their dimensions (the thickness) being much smaller than the other two. Although all structures made of FGM can be described using 3-D continuum mechanics, exact solutions exist only for a few specific problems with very idealized material types, geometry, and boundary conditions [6,7]. For more realistic cases, 3-D numerical simulation tools such as ANSYS and ABAQUS are often used to find approximate solutions. However, this approach is computation intensive and usually used in the detailed analysis due to its prohibitive computational cost. In view of the fact that the thickness is small, analysis of such structures can be simplified using 2-D models. Although many 2-D models have been developed to analyze FGM shells treating different topics, most of them rely on some apriori kinematic assumptions [8]. Examples may be found in the application of classical lamination theory (CLT) in the problem of maximization of buckling stability limits of Functional Graded (FG) rings/long cylindrical shells subjected to external pressure [9], the use of the first-order sheardeformation theory (FOSDT) in studying the stability problem of a circular cylindrical shell composed of FGM with elasticity modulus varying continuously in the thickness direction under combined external pressure and axial compression loads [10], the implementation of the high-order shell theory in dynamic buckling of imperfect FGM cylindrical shells with integrated surface-bonded sensor and actuator layers subjected to some complex combinations
Z. Yifeng et al. / Composite Structures 94 (2012) 966–976
967
Nomenclature
v
a column matrix containing 2-D generalized strains v ¼ ½T jT T
dW; dW 2D
, j Cij
c l l⁄, K⁄ X
x oX eab, Kab ck Cij D eijk fi, ma Fij Gk gk
3-D and 2-D virtual work, respectively column matrix of 2-D generalized strains three-dimensional strain tensor column matrix of transverse strain measures characteristic magnitude of the elastic constants effective shear and bulk modulus estimated by Mori– Tanaka scheme, respectively the reference surface inertial angular velocity of triad Bi boundary of the reference surface 2-D generalized shell strains through-the-thickness average of in-plane warping functions matrix of direction cosines 3-D material matrix the permutation symbol generalized forces and moments, respectively mixed-basis component of the deformation gradient tensor covariant basis vector of the deformed shell contravariant base vector of the undeformed shell
of thermo–electro-mechanical loads [11], and the use of Layerwise theories in analysis of a FGM cylindrical shell under dynamic loads [12]. CLT ignores transverse shear effects and provides reasonable results only for very thin shells. Moreover, in CLT, both plane strains and plane stresses are assumed which we know will not be true at the same time for materials having nonzero Poisson’s ratios. A number of shear deformation theories have been developed to overcome some drawbacks of CLT, with the simplest of which being FOSDT (equivalent to Reissner–Mindlin theory for shells made of isotropic homogeneous materials), where a constant distribution of shear strain through the thickness is assumed and a shear correction factor is required to account for the deviation of the real shear strain from the assumed constant one. The dependence of the shear correction factor on the geometry and material of the shell makes it difficult to guarantee the accuracy of FOSDT. By expanding the displacement field of the shell using higher-order polynomials, higher-order shear deformation shell theories are developed, which can account for both transverse normal and shear deformations without relying on shear correction factors. However, as indicated by Bian et al. [13] and Das et al. [5] for laminated shells, models based on higher-order theories cannot capture the discontinuous slope of in-plane and transverse displacement components in the thickness direction. By extending a generalized refined theory (referred as Soldatos shell theory [14]) which incorporates shape functions to guarantee continuity of transverse shear stresses at interfaces, these authors provide analytical solutions for simply- and multiply- spanned FG shells under cylindrical bending. Layerwise theories can produce reasonable results at the cost of complex models and expensive computation. Despite of the popularity of the aforementioned methods in analyzing many FG shells and FG laminated shells, these approaches have two major disadvantages: (1) the apriori assumptions which are naturally extended from the analysis of isotropic homogeneous structures cannot be easily justified for heterogeneous and anisotropic structures, such as FGMs; (2) it is difficult
h, R1 K, K2D l La
thickness and radius of the shell, respectively 3-D and 2-D kinetic energy, respectively characteristic wavelength of the shell deformation integration constants for determining the relationship between v and ck M, N shell moment resultants and force resultants, respectively n a small parameter used to denote the order of strains pi applied body force Qi tractions applied on lateral surfaces R strain measures of Reissner–Mindlin shell U free energy per unit area wi warping functions for the normal-line element xi Cartesian coordinate base vector of the coordinate system before and after bi, Bi deformation, respectively b ^r; R position vector of a material point before and after deformation, respectively triad vector for deformed Reissner–Mindlin shell Bi b R through-the-thickness average of the position vector R kk, k3, K3 Lagrange multipliers for enforcing in-plane and outplane warping constraints
for an analyst to determine the accuracy of the result and which assumption should be chosen for efficient and accurate analysis for a particular shell. In this paper, an efficient high-fidelity shell model for multilayer FG cylindrical shell is constructed by the variational asymptotic method (VAM). The original 3-D nonlinear elasticity problem is casted in an intrinsic form so that the theory is applicable for arbitrarily displacement and global rotation with the restriction that to the strain field is small. VAM is then used to rigorously split the 3-D elasticity problem into two problems: a nonlinear, 2-D, shell analysis over the reference surface to obtain the global deformation and a linear, 1-D, through-the-thickness analysis to provide the 2-D generalized constitutive law and the recovering relations to approximate the original 3-D results. The nonuniqueness of asymptotic theory correct up to a certain order is used to cast the obtained asymptotically correct second-order free energy into a Reissner– Mindlin type model to account for transverse shear deformation. The present theory extends a simple and accurate model developed recently for composite laminates by Yu et al. [15,16] so that the multilayer FG shell can be treated in the same framework. Since the procedure is quite similar, the authors have chosen to repeat some formulae and texts from their previous publications in order to make the present paper more self-contained. 2. 3-D Formulation The elastodynamic behavior of a solid is governed by the extended Hamilton principle:
Z
t2
½dðK UÞ þ dWdt ¼ 0
ð1Þ
t1
where t1 and t2 are arbitrary fixed times, K is the kinetic energy, U is a term related with the internal energy and dW is the virtual work of applied loads. The overbar is used to indicate that the virtual work does not necessarily represent variation of functionals.
968
Z. Yifeng et al. / Composite Structures 94 (2012) 966–976
ga ¼ aa þ x3 b3;a ;
g3 ¼ b3
ð7Þ
From the differential geometry of surface and following Ref. [17], the derivatives of unit vectors bi,a can be expressed as
bi;a ¼ Aa ðka2 b1 þ ka1 b2 þ ka3 b3 Þ bi
ð8Þ
where Kab refers to the usual out-of-plane curvatures and k12 = k21 = 0 since the coordinates are the lines of curvatures. The geodesic curvatures ka3 can be expressed in terms of the Lamé parameters as
k13 ¼
A1;2 ; A1 A2
k23 ¼
A2;1 A1 A2
ð9Þ
One can obtain the explicit expressions for the 3-D covariant and contravariant base vectors by using Eq. (8):
ga ¼ Aa ð1 þ x3 kaa Þba ; ga ¼ Fig. 1. Schematic of shell deformation.
A cylindrical shell is a 3-D body with a relatively small thickness h and a smooth reference surface (see Fig. 1). The geometry of the reference surface can be described by a set of arbitrary curvilinear coordinates, xa. (Here and throughout the paper, Greek indices assume values 1 and 2 while Latin indices assume 1, 2, and 3; repeated indices are summed over their range except where explicitly indicated). Without loss of generality, the lines of curvatures are chosen to be the curvilinear coordinates to simplify the formulation. To uniquely represent an arbitrary point in a 3-D medium, another coordinate off the reference surface is needed. It is a natural and convenient choice to take the third coordinate x3 as the transverse normal coordinate. Letting bi denote the unit vector along xi for the undeformed shell, the position of any material point in the undeformed configuration can then be described by its position vector ^r from a fixed point O as
^rðx1 ; x2 ; x3 ; tÞ ¼ rðx1 ; x2 ; tÞ þ x3 b3
ð3Þ
where the angle brackets denote the definite integral through the thickness of the shell and will be used throughout the rest of the development. The 2-D base vectors are defined conventionally as
aa ðx1 ; x2 Þ ¼ r;a
ð4Þ
@ðÞ where ðÞ;a ¼ @x . a The so-called Lamé parameters can be obtained from this definition as
Aa ðx1 ; x2 Þ ¼
pffiffiffiffiffiffiffiffiffiffiffiffiffiffi aa aa
ð5Þ
Note here the summation convention is not applied here because a is not a dummy index; this rule will apply in similar situations throughout the rest of development. Then the unit vectors along coordinates xa can be obtained as:
ba ðx1 ; x2 Þ ¼
aa Aa
ð10Þ
ba Aa ð1 þ x3 kaa Þ
When the shell deforms, the particle that had position vector ^r b in the deformed in the undeformed state now has position vector R shell. The latter can be uniquely determined by the deformation of the 3-D body. Similarly, another triad Bi is introduced for the deformed configuration. Note that the Bi unit vectors are just tools to express vectors and tensors in their component form during the derivation. They are not necessarily tangent to the coordinates of the deformed shell. The relation between Bi and bi can be specified by an arbitrarily large rotation specified in terms of the matrix of direction cosines C(x1, x2) so that
Bi ¼ C ij bj ;
C ij ¼ Bi bj
ð11Þ
subjecting to the requirement that Bi is coincident with bi when the structure is undeformed. b can be represented as Now the position vector R
b i ðx1 ; x2 ; x3 ; tÞ ¼ Rðx1 ; x2 ; tÞ þ x3 B3 ðx1 ; x2 ; tÞ R þ wi ðx1 ; x2 ; x3 ; tÞBi ðx1 ; x2 ; tÞ
ð12Þ
ð2Þ
where r is the position vector from O to the point located by xa on the reference surface at a specific time t. When the origin of the coordinate system is located on the middle surface, it naturally follows that
h^rðx1 ; x2 ; x3 ; tÞi ¼ hrðx1 ; x2 ; tÞ
g3 ¼ g3 ¼ b 3 ;
ð6Þ
It is obvious that the unit vectors bi form an orthogonal triad. The 3-D covariant base vectors gi ¼ ^r;i for the chosen coordinate system can be obtained by using Eq. (4) as
where wi are the warping functions introduced to accommodate all possible deformations. Eq. (12) is six times redundant because of the way warping introduced. Proper definitions of R and Bi are needed to introduce six constraints to ensure a one-to-one mapping between R, Bi and wi. R can be defined similarly as Eq. (3) to be the average position through the thickness, from which it follows that the warping functions must satisfy the three constraints
hwi ðx1 ; x2 ; x3 ; tÞi ¼
ck ; 0
ck ¼
c1 c2
ð13Þ
where c1 and c2 are functions of the in-plane coordinates xa and time t, introduced for providing free variables for the construction of an optimal Reissner–Mindlin model. Two other constraints can be specified by taking B3 as the normal to the reference surface of the deformed shell. It should be noted that this choice has nothing to do with the well-known Kirchhoff hypothesis. In the Kirchhoff assumption, no local deformation of the transverse normal is allowed. However, all possible deformation are allowed by using the warping functions in present derivation. Because Ba can freely rotate around B3, the last constraint can be introduced as
B1 R;2 ¼ B2 R;1
ð14Þ
Based on the concept of decomposition of rotation tensor [18], the Jauman–Biot–Cauchy strain components for small local rotation are given by
969
Z. Yifeng et al. / Composite Structures 94 (2012) 966–976
1 2
Cij ¼ ðF ij þ F ji Þ dij
ð15Þ U¼
X
where dij is the Kronecker symbol, and Fij is the mixed-basis component of the deformation gradient tensor such that
F ij ¼ Bi Gk gk bj
ð16Þ
@b R where Gk ¼ @x is the covariant basis vector of the deformed configuk ration, which can be expressed in terms of the 2-D generalized strains following Ref. [19] as:
R;a ¼ Aa ðBa þ eab Bb Þ
ð17Þ
Bi;a ¼ Aa ðK a2 B1 þ K a1 B2 þ K a3 B3 Þ Bi
ð18Þ
where eab are the 2-D in-plane strains and of which the order is denoted by n, and Kij are the curvatures of the deformed surface which are the summation of curvatures of undeformed geometry kij and curvatures introduced by the deformation jij of which the order is denoted by n/h. Using this definition, one can show that Eq. (14) implies e12 = e21, i.e.
B1 R;2 B2 R;1 ¼ A2 A1
ð19Þ
which can serve as another constraint to specify the global rotation of the triad Bi. For geometrically nonlinear analysis, the strains are small compared to unity and warpings are of the order of strain or smaller, which have the effect of removing all the terms that are products of the warping and the generalized strains. With the help of Eqs. (12), (15), (16), (17) and (18), we can express the 3-D strain field as
Ce ¼ þ x3 j þ Ia wk;a ; @ðÞ where ðÞ;i ¼ @x i
1 Aa
2Cs ¼ wk;3 þ ea w3;a ;
; ðÞk ¼ ½ ðÞ1
Ct ¼ w3;3
ð20Þ
2C23 T ;
Ct ¼ C33 ;
¼ ½ e11 2e12 e22 T ; j ¼ ½ K 11 K 12 þ K 21 K 22 T 2
3
1 0 6 7 I 1 ¼ 4 0 1 5; 0 0 T
e1 ¼ ½ 1 0 ;
2
ð21Þ
ð24Þ where De, Des, Det, Ds, Dst, Dt are the corresponding partition matrices of the 6 6 elastic material matrix D at constant electric field, X is the surface stretched by the undeformed reference surface and / is geometric correction as 2
g g2 g3 h /¼ 1 ¼ 1 þ x3 ðk11 þ k22 Þ þ O 2 ja1 a2 j l
! ð25Þ
To calculate the kinetic energy, the absolute velocity of a generic point in the structure is obtained by taking a time derivative of Eq. (12) as
b ¼ V þ xðn þ wÞ þ w _ V
ð26Þ
where ( ) is the partial derivative with respect to time, V is the absolute velocity of a point in the deformed reference surface, x is the inertial angular velocity of Bi bases, and the notation () forms an antisymmetric matrix from a vector according to ð Þ ¼ eijk ðÞk using b ; V; x; w denote column the permutation symbol eijk. The symbols V matrices containing the components of corresponding vectors in Bi bases, and n ¼ ½ 0 0 x3 T . The kinetic energy of the shell structure can be obtained by
K¼
1 2
Z
q Vb T Vb dV ¼ K 2D þ K
ð27Þ
V
where q is the mass density and
Z 1 fnV þ xT jxÞdX V T V þ 2xT l ðl 2 X Z 1 ~ w þ wÞ ~ w þ wÞ ~ nÞT ðx ~ w þ wÞdV _ T ðx _ þ 2ðV þ x _ K ¼ q½ðx 2 V K 2D ¼
; l n; j are inertial constants commonly used in shell dynamwhere l ics, which can be trivially obtained through simple integrals over the thickness as:
0 ð22Þ
0
0
b þ s dR b þ þ b dR b ÞdX þ ðhP d Ri
X
1 þ x3 k11
e12 þ x3 j21 þ w1;2 w2 k23
e12 þ x3 j12 þ w2;1 þ w1 k13 þ 1 þ x3 k22 1 þ x3 k11 e22 þ x3 j22 þ w2;2 þ w3 k22 þ w1 k23 C22 ¼ 1 þ x3 k22 w3;1 w1 k11 2C13 ¼ w1;3 þ 1 þ x3 k11 w3;2 w2 k22 2C23 ¼ w2;3 þ 1 þ x3 k22 C33 ¼ w3;3 ð23Þ With the knowledge of mechanical strains, the energy enthalpy can be expressed as:
ð29Þ
The virtual work of the structure can be calculated as Z Z
dW ¼
T
e11 þ x3 j11 þ w1;1 þ w3 k11 w2 k13
2C12 ¼
DTst
l ¼ hqi; ln ¼ ½ 0 0 hx3 qi T ; 2 2 3 x3 q 0 0 6 7 j¼4 0 x23 q 0 5
The explicit expressions of Cij can be written as
C11 ¼
Ds
9 3T 8 + > Z < Ce > = 7 Dst 5 2Cs / dX ¼ U A /dX > > X : ; Ct Dt
Det
ð28Þ
3
0 0 6 7 I 2 ¼ 4 1 0 5; 0 1 e2 ¼ ½ 0 1
Des
ðÞ2 T , and
Ce ¼ ½ C11 2C12 C22 T ; 2Cs ¼ ½ 2C13
8 9 2 *> Ce >T De = 1 < 6 DT 2 Cs 4 es > 2 > : ; Ct DTes
Z
b hQ d Rids
ð30Þ
@X
where oX denotes the boundary of the reference surface, and ðÞ ¼ ðÞjx3 ¼h=2 . P = PiBi is the applied body force, s = siBi, b = biBi are tractions applied on the top and bottom surfaces, respectively, and Q = QiBi is the applied tractions along the lateral surfaces. The Lagrangian variation of the displacement field which can be expressed as
b ¼ dq i Bi þ x3 dB3 þ dwi Bi þ wj dBj dR
ð31Þ
in which the virtual displacement and rotation are defined by
2 B1 þ dw 1 B2 þ dw 3 B3 Þ Bi i ¼ dR Bi ; dBi ¼ ðdw dq
ð32Þ
i contain the components of the virtual displacei and dw where dq ment and rotation in the Bi system, respectively. Since the warping functions are small, one may safely ignore b products of the warping functions and the virtual rotations in d R and obtain the virtual work due to applied loads as
970
Z. Yifeng et al. / Composite Structures 94 (2012) 966–976
dW ¼ dW 2D þ dW
ð33Þ
where
dW 2D ¼
Z
a ÞdX þ i þ ma dw ðfi dq
X
dW ¼
Z X
Z
a Þds i þ hx3 Q a idw ðhQ 1 idq
ð34Þ
where UA0 can be obtained from Eq. (24) by dropping the derivatives with respect to xi in Eq. (20), resulting in
D E 2U A0 ¼ 2 ð þ x3 jÞT ðDes wk;3 þ Det w3;3 Þ þ wTk;3 Dst w3;3 þ hð þ x3 jÞT De ð þ x3 jÞi þ wTk;3 Ds wk;3 þ wT3;3 Dt w3;3
@X
hpi dwi i þ si dwþi þ bi dwi dX þ
Z
hQ i dwi ids
þ oððh=lÞ0 n2 Þ
ð35Þ
ð40Þ
X
with the generalized forces fi and moments ma defined as
h fi ¼ h/i i þ si þ bi ; ma ¼ hx3 /a i þ ðsa ba Þ 2
ð36Þ
T where ðÞk ¼ ½ ðÞ1 ðÞ2 ; ðÞ;3 ¼ @ðÞ=@x3 . It is clear that the warping functions wi only appear in UA0, which means wi can be solved from the following much simpler variational statement:
The second integration in Eq. (35) accounts for virtual work done through warping functions along the lateral boundaries of the shell. This term is necessary for the edge-zone problem, which is an important subject in its own right and beyond the scope of the present paper. Therefore, we will drop this term hereafter. With the knowledge of Eqs. (33), (24) and (27), the Hamiltons principle in Eq. (1) becomes
along with the constrains in Eq. (13). This variational statement has the following Euler–Lagrange equations:
Z
t2
½dðK 2D þ K UÞ þ dW 2D þ dW dt ¼ 0
ð37Þ
t1
3. Dimensional reduction The dimensional reduction from the original 3-D formulation to a 2-D shell model can only be done approximately. One way to do it is to take advantage of the small parameters in the formulation to construct the 2-D formulation so that the reduced model can achieve the minimum loss of accuracy in comparison to the original 3-D formulation. In order to apply VAM, we first need to assess the order of quantities in terms of small parameters. The quantities of interest have the following orders:
hp3 s3 b3 lðh=lÞ2 n; hpa sa ba lðh=lÞn; Q a ln; Q 3 lðh=lÞn
ð38Þ
where l represents the order of elastic material constants.
;3
ð þ x3 jÞ
T
Dkes
T
T þ wkk;3 Dks þ wk3;3 Dkst ¼ kk
ð42Þ
Z d K 2D U A0 /dX þ dW 2D dt ¼ 0 X
where k = 1, 2, 3, . . . , N with N denoting the total number of layers, kk and k3 are Lagrange multipliers corresponding to the in-plane and out-of-plane constraint equations expressed in Eq. (13); ()k denotes functions () for the kth layer. The boundary conditions as well as the interlaminar continuous conditions are:
½ð þ x3 jÞT Det þ ðwk;3 ÞT Dst þ w3;3 Dt þ ¼ 0; h i ð þ x3 jÞT Des þ ðwk;3 ÞT Ds þ w3;3 DTst ¼ 0
ð43Þ
and
h i ½w3 ¼ 0; ð þ x3 jÞT Det þ ðwk;3 ÞT Dst þ w3;3 DTt ¼ 0; X h i i ½wk ¼ 0; ð þ x3 jÞT Des þ ðwk;3 ÞT Ds þ w3;3 DTst ¼ 0
ð44Þ
Xi
where Xi denotes the interfaces between the ith layer and i + 1th layer; i = 1, . . . , N 1; the bracket [] denotes the jump of the enclosed argument on the interface. From these conditions, we can solve wi,3 from Eq. (42):
T
1 bk D bk wkk;3 ¼ ð þ x3 jÞT D ; es s
ð45Þ
b k =D bk wk3;3 ¼ ð þ x3 jÞT D et t with the hatted quantities being expressed as:
T b k ¼ Dk D b k Dk =Dk ; D es et es st t
1 k k k k b D et ¼ Det Des Ds =Dkst ;
T 1 b k ¼ Dk Dk bk D D Dkst t s t st
To clearly illustrate the application of VAM for FGM shells, a classical FGM shell model is first constructed. By applying VAM, the zeroth-order approximation of the variational statement in Eq. (37) can be obtained as
t1
T ð þ x3 jÞT Dket þ wkk;3 Dkst þ wk3;3 Dkt ¼ k3 ;
ð46Þ
Substituting Eq. (45) into Eq. (40), the zeroth-order strain energy can be expressed as
3.1. Zeroth-order approximation
t2
ð41Þ
;3
So far, a 3-D displacement-based formulation for FG shells has been developed in terms of 2-D displacements (represented by R r), rotations (represented by bi and Bi), and 3-D warping functions (wi). If we attempt to solve this problem directly, we will meet the same difficulty as solving any full 3-D problem with the additional difficulty coming from the anisotropy and heterogeneity of functional graded materials. The main complexity comes from the unknown 3-D warping functions wi. The common practice in the literature is to assume wi in terms of 2-D displacements and rotations to straightforwardly reduce the original 3-D continuum model into a 2-D shell model. However, for shells made with general anisotropic and heterogeneous materials such as FGMs, the imposition of such assumptions may introduce significant errors. Fortunately, VAM provides a useful technique to obtain wi through an asymptotical analysis of the variational statement in Eq. (37) in terms of small parameters inherent in the problem.
Z
dU A0 ¼ 0
ð39Þ
2U A0 ¼
j
T " b A bT B
b B b D
#
j
þo
0 ! h n2 l
ð47Þ
with
D E b ¼ hD b ¼ hx3 D b ¼ x2 D b e i; D b b e i; B A 3 e ; b e ¼ De Des D1 D bt b T Det D b T =D D e es et
ð48Þ
971
Z. Yifeng et al. / Composite Structures 94 (2012) 966–976
With UA0 expressed in Eq. (47), the original 3-D problem in Eq. (1) has been rigorously reduced to a 2-D formulation in Eq. (39) which approximates the original problem asymptotically correct to the order of (h/l)0. If we define the force resultants N and moment resultants M conjugate to and j, such that
N¼
@U A0 ; @
@U A0 @j
M¼
ð49Þ
then we obtain a 2-D constitutive model for the classical shell analysis of FGM shells as
N M
"
¼
b A bT B
b B b D
#
j
ð50Þ
It is clear that although the shell is made of FGM, the 2-D shell model of the zeroth-order remains the same. Despite the similarity with CLT, the present model is asymptotically correct and has the following features in contrast with CLT: 1. The normal line of undeformed shell does not remain straight and normal to the deformed shell; rather, it deforms in both the normal and in-plane directions in response to shell deformation ( and j). 2. This model can handle general FGMs with full anisotropy. 3. It can be easily observed that neither the normal strain nor the transverse shear strains vanish. The transverse normal and shear stresses can be shown to vanish, which are not assumed a priori but direct consequences of the model. For the zeroth-order approximation, the 3-D strain field can be recovered using Eq. (20) by neglecting those terms of order higher than (h/l)0: 0 e
C ¼ þ x3 j;
0 s
2C ¼ wk;3 ;
0 t
C ¼ w3;3
ð51Þ
The 3-D stress field can be computed using the original 3-D constitutive relations, which can be derived from Eq. (24) as
8 9 < re =
2
De T D rs ¼ 6 4 es : ; rt DTet
Des Ds DTst
38 0 9 Det > < Ce > = Dst 7 5 2C0s > > D : 0 ; t
ð53Þ
Up to this point, the zeroth-order solution wi in Eq. (45) as well as the 2-D strain energy in Eq. (47) are valid for FGM shells with fully populated 6 6 material matrix D. As most real materials have at least a monoclinic symmetry about their own mid-plane, hereafter, monoclinic material matrix characterized by 13 independent material properties, implying Des = 0 and Dst = 0, will be used for the rest of derivation. This leads to much simpler expressions for the zeroth-order approximation of warping function:
where
wk3 ¼ Dk? v
T Dk? ¼ Dkes =Dkt
To obtain the first-order approximation, we simply perturb the zeroth-order warping functions as
hn ; wk ¼ v k þ o l
hn w3 ¼ v 3 þ D? v þ o l
where Dk ¼ De Det DTet =Dt and v3 is decoupled from vk. Considering the warping constraint in Eq. (13), v3 only has a trivial solution. The stationary conditions of the functional given in Eq. (57) are
Dks v k;3 þ Dks ea Dk? v;a ¼ Dka;3 v;a þ g k;3 þ kk ;3
Ds v k;3 þ Ds ea D? v;a ¼ sk x ¼þh=2
3 Ds v k;3 þ Ds ea D? v;a ¼ bk
ð58Þ
x3 ¼h=2
where
h Dka;3 ¼ ITa Dkk
i x3 Dkk ;
g k;3 ¼ /kk
½v k ¼ 0; ½Ds ðv k;3 þ ea D? v;a Þ
¼0
ð59Þ
Dks;3 v k þ Dks ea Dk? v;a ¼ Dka v;a þ g k þ kk x3 þ constkk
ð61Þ
and the interface continuity condition in the second equation of Eq. (60) becomes
ð55Þ
v ¼ ½ j T ; hD? i ¼ 0 Note that inter-lamina continuity of Dk? must be maintained due to the continuity of warping functions to produce a continuous displacement field. It should also be pointed out that we have constrained in-plane warpings to vanish (hwki = 0) as the free constants in Eq. (13) will be absorbed in the first-order approximation.
ð60Þ
The inter-lamina continuity on Dka and gk are maintained by the second condition in Eq. (60). It should be mentioned that since the goal is to obtain an interior solution for the shell without considering edge effects, integration by parts with respect to the in-plane coordinates is used hereafter and throughout the rest of the paper, whenever it is convenient for the derivation. Integrating the first equation in Eq. (58), one obtains
ð54Þ
T x3 Dket =Dkt ;
ð56Þ
Substituting Eq. (56) back into Eq. (20), and using Eqs. (29), (34) and (37), one can obtain the leading terms for the first-order approximation of the variational statement in Eq. (37) as D T dP1 ¼ v Tk;3 Ds dv k;3 þ v 3;3 Dt dv 3;3 þ ðe þ x3 jÞ Dk Ia dv k;a E þðea D? v;a ÞT Ds dv k;3 pTk dv k sTk dv þk bTk dv k þ oððh=lÞ2 n2 Þ ð57Þ
Xi
with
wkk ¼ 0;
3.2. First-order approximation
The continuity conditions on the interfaces can be derived as
ð52Þ
Ct
re ¼ ½ r11 r12 r22 T ; rs ¼ ½ r13 r23 T ; rt ¼ r33
In view of Eq. (42), rs and rt are identically zero regardless the material properties, which means the classical model can not predict out-of-plane stresses. In order to predict the out-of-plane stresses, we need to carry out the asymptotic expansion of the 3D variational statement in Eq. (1) into higher order, which is done in next section.
Fig. 2. Schematic figure for laminated cylindrical shell.
972
Z. Yifeng et al. / Composite Structures 94 (2012) 966–976
constkþ1 constkk ¼ k
Z
zkþ1
zk
h ITa Dkk
i x3 Dkk dx3
ð62Þ
constkk
where is the integration constant for kth layer, zk is the x3 coordinate of the bottom of the kth layer (see Fig. 2). Dka ; g ka can be defined as:
Dka ðx3 Þ ¼
k
g a ðx3 Þ ¼
Z
x3
zk
Z
x3
zk
Dka;3 dz ¼
k
g a;3 dz ¼
Z
x3
zk
Z
x3
zk
h
ITa Dkk
i
zDkk dz;
ð63Þ
/kk dz
ð64Þ
Since Bi is uniquely determined by Bi and c, the kinematic identity between the strains measures R of Reissner–Mindlin shell and can be derived as
¼ R Da c;a where
D1 ¼ D2 ¼
0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 0 1
R ¼ ½ e11
Solving Eqs. (61) and (62), the two boundary equations in Eq. (58) as well as Eq. (13), one obtains the following Lagrange multipliers and warping functions:
1 ðsk þ bk hDa;3 iva hg ;3 iÞ h k ¼ Da þ La v;a þ gk
ð72Þ
2e12
T ; T
ð73Þ
;
e22 K 11 K 12 þ K 21 K 22 T
Now one can rewrite the strain energy expressed in Eq. (68), correct to the orders of interest according to Eq. (72), in terms of strains of the Reissner–Mindlin model as
kk ¼
ð65Þ
2PR ¼ RT AR 2RT ADa c;a þ RT;a Ba R;a þ 2RT;1 CR;2 2RT F
v
ð66Þ
The generalized Reissner–Mindlin model of practical use can be of the form
k k
with
Dka;3
1 bk; D ¼ Dks a
hDa i ¼ 0;
2PR ¼ RT AR R þ cT Gc þ 2RT F R þ 2cT F c
1 gk;3 ¼ Dks g^k ;
ð67Þ
La v;a ¼ ck =h
hgi ¼ 0;
b k ; g^k can be obtained from Eqs. (61), (62), (65) as well as the where D a boundary conditions in Eq. (58). Now we are ready to obtain an expression for the total energy that is asymptotically correct through the order of l(h/l)2n2, viz., T
2P1 ¼ v Av þ v
T ;a Ba
T ;1 C
T
v;a þ 2v v;2 2v F
To find an equivalent Reissner–Mindlin model Eq. (75) for Eq. (74), one has to eliminate all partial derivatives of the classical 2-D strain. The equilibrium equations in Ref. [19] are used to achieve this purpose. From the two equilibrium equations balancing bending moments with applied moments ma, one can obtain the following formula:
Gc F c ¼ DTa AR;a þ ½ m1
T
m2
ð76Þ
where FR,a is dropped because they are high order terms. Substituting Eq. (76) into Eq. (74), one can show that FR = F and Fc = 0. Finally one can rewrite Eq. (74) as
# hx3 Dk i ; A¼ hx3 Dk i x23 Dk D E b a þ LT hDa;3 i þ hDT iLa ; b T D1 D Ba ¼ Dsaa DT? D? D a s a a;3 D E D E T 1 T T b 2 þ L D2;3 þ D bTD D C ¼ Ds12 D? D? D 1 s 1 1;3 L2 ; D E D E T T b T D1 g^;a F ¼ Dþ? s3 þ D? b3 þ DT? /3 D a s
LTa sk þ bk þ hpk i hDk i
2P1 ¼ 2PR þ U
ð69Þ
Ba ¼ Ba þ ADa G1 DTa A;
ð78Þ
C ¼ C þ AD1 G1 DT2 A
ð79Þ
⁄
4. Transforming into the Reissner–Mindlin model Although Eq. (68) is asymptotically correct through the second order and straightforward use of this strain energy is possible, it involves more complex boundary conditions than necessary since it contains derivatives of the generalized strain measures. To obtain an energy functional that is of practical use, one can transform Eq. (68) into the Reissner–Mindlin model. In the Reissner–Mindlin model, there are two additional degrees of freedom, which are the transverse shear strains incorporated into the rotation of transverse normal. Another triad Bi for the deformed shell is introduced to define the 2-D strains as
ð70Þ ð71Þ 2c23 T .
U ¼ RT;a Ba R;a þ 2RT;1 CR;2 and
Eq. (68) is an energy functional expressed in terms of 2-D variables which can asymptotically approximate the original 3-D energy. It is noted that quadratic terms of the applied loads are neglected as they will not affect the 2-D model.
R;a ¼ Aa Ba þ eab Bb þ 2ca3 B3 Bi;a ¼ Aa K a2 B1 þ K a1 B2 þ K a3 B3 Bi
ð77Þ
where
;a
where the transverse shear strains are c ¼ ½ 2c13
ð75Þ
ð68Þ
where
"
ð74Þ
If we can drive U to be zero for any R, then we have found an asymptotically correct Reissner–Mindlin shell model. For general anisotropic shells, this term will not be zero; but we can minimize the error to obtain a Reissner–Mindlin model that is as close to asymptotical correctness as possible. The accuracy of the Reissner–Mindlin model depends on how close to zero one can drive this term of the energy. In other words, one needs to seek an optimal set of the 27 unknowns (3 unknowns for G and 24 unknowns for La) so that the value of the quadratic form in Eq. (78) is as close to zero as possible for arbitrary generalized strain measures. We let the distinct 78 terms in the symmetric 12 12 coefficient matrix equal to zeros to formulate 78 equations. It is a linear system with 27 unknowns. Then we use a least square technique to solve the overdetermined system for the constants. Mathematically, the overdetermined system (78 equations with 27 unknowns, indicated by MX = b) may be singular for some material properties. For example, the rank of MTM is only 26 for single-layer isotropic and orthotropic plates. In this situation, singular value decomposing technique can be applied to solve this least square problem. After minimizing U⁄, the Reissner–Mindlin model to be used for 2-D shell problem can be expressed as:
2PR ¼ RT AR þ cT Gc þ 2RT F
ð80Þ
973
Z. Yifeng et al. / Composite Structures 94 (2012) 966–976
where A, G, F capture the necessary material and geometric information obtained from the dimensional reduction process. It is worthy to emphasize that although the 2-D constitutive model is constructed in a way dramatically different from traditional Reissner–Mindlin models, the shell analysis remains the same, with no changes in the governing equations and boundary conditions except that the strain measures are now defined using Eq. (70).
T
Dkt yk3;3 þ Dket Ia v kk;a þ eTb Dks v kk;3 þ ea Dk? v;a þ /3 ¼ K3 ;
Dt y3;3 þ
DTet Ia
v k;a
þ
;b
;3
¼ s3 ;
Dt y3;3 þ DTet Ia v k;a ¼ b3 ; h i ½y3 ¼ 0; Dt y3;3 þ DTet Ia v k;a
Xi
¼0 ð85Þ
5. Recovery relations Thus far, a generalized Reissner–Mindlin model based on the asymptotically correct second-order energy for FG cylindrical shells has been constructed. Then, the original problem governed by the 3-D Hamiltons principle in Eq. (37) can be approximated by the following 2-D variational statement:
Z
t2
½dðK 2D PR Þ þ dW 2D dt ¼ 0
ð81Þ
t1
with K2D obtained from Eq. (28), PR from Eq. (80) and dW 2D from Eq. (34). It is noted that the variational statement in Eq. (80) is not asymptotically correct because some optimization technique has been used in the transform the asymptotically correct energy functional in Eq. (68) to the generalized Reissner–Mindlin model in Eq. (75), which is a price one has to pay for the simplicity of the 2-D shell analysis using the Reissner–Mindlin model. The 2-D analysis based on the variational statement of Eq. (76) can be routinely solved using shell elements in commercial finite element packages. This model can be used for various analyses of FGM shells, spanning from static, dynamic, buckling, to aeroelastic analyses. In many applications, however, the capability of predicting accurate 2-D displacement fields of FGM shells is not sufficient. Ultimately, the fidelity of a reduced-order model should be evaluated based on how well it can predict the 3-D displacement/strain/stress fields for the original 3-D problem. Therefore, it is necessary to provide recovery relations to express the 3-D displacement/strain/stress fields in terms of 2-D quantities and x3. Using Eqs. (2), (11) and (12), one can recover the 3-D displacement field through the first order as
U i ¼ ui þ x3 ðC 3i d3i Þ þ C ij wj
ð82Þ
where wa = va, w3 = C\v. From Eq. (20), the 3-D strain field can be recovered up to the first order as
Ce ¼ e þ x3 j;
2Cs ¼ v k;3 þ ea D? v;a ;
Ct ¼ D?;3 v
ð83Þ
Consequently, 3-D stresses rij can be obtained by applying the 3-D constitutive relations. Since we have obtained an optimal estimation of the shear stiffness matrix G, the recovered 3-D results up to the first order are better than CLT and FOSDT. However, the transverse normal stress (r33) is a second-order quantity and cannot be estimated within the first-order approximation. Despite that it is usually much smaller than other stress components, r33 is critical for predicting some structural failure mechanisms such as delamination. In order to obtain a reasonable recovery for the transverse normal stress, VAM procedure is applied once more to find the warping functions of second-order accuracy. Similarly, we perturb the warping functions as
wk ¼ v k þ yk þ oððh=lÞ2 nÞ; w3 ¼ D? v þ y3 þ ððh=lÞ2 nÞ
where K3 is the Lagrange multiplier to enforce the constraint hy3i = 0. The Euler–Lagrange equation and the inter-surface continuity equation in Eq. (85) can be expressed as:
T Dkt yk3;3 þ Dket Ia v kk;a ¼ Ekab vab þ Sk þ K3 x3 þ constky3
ð86Þ
k k k constkþ1 y3 consty3 ¼ Eab vab ðzkþ1 Þ þ S ðzkþ1 Þ
ð87Þ
where
Ekab ðx3 Þ ¼ Sk ðx3 Þ ¼
Z
Z
x3
zk x3
zk
k k eTa Dk b þ Dsab D? dz
k eTa g k ;a þ /3 dz
ð88Þ ð89Þ
The Lagrange multiplier K3 and the solution of y3 are given by
1 s3 þ b3 hEab;3 iv;ab hS;3 i h yk3 ¼ Ek;ab v;ab þ Sk
K3 ¼
ð90Þ ð91Þ
with
1 D E k Ekab;3 ¼ Dkt Ek ab ; Eab ¼ 0;
1 Sk;3 ¼ Dkt Sk ; hSi ¼ 0
ð92Þ
k⁄ where Ek ab and S can be obtained from Eqs. (86), (87), (90) as well as the interlaminar continuous conditions in the last equation of Eq. (85). Although y3 can help us obtain an energy expression asymptotically corrected up to the order of (h/l)4n2, such an energy expression is too complex for practical use. We will still use the Reissner– Mindlin model to carry out the 2-D shell analysis and use y3 for the second-order prediction of the 3-D displacement/strain/stress field. As will be shown latter, this approach achieves excellent predictions even though only the Reissner–Mindlin shell model is used for the 2-D shell analysis. Finally, we can recover the 3-D displacement field up to the second order as
U i ¼ ui þ x3 ðC 3i d3i Þ þ C ji wj þ d3i C 3i y3
ð93Þ
and the strains up to the second order as
Ce ¼ þ x3 j þ Ia v k;a ;
2Cs ¼ v k;3 þ ea D? v;a ;
Ct ¼ D?;3 v þ y3;3
ð94Þ
Finally, we can recover the 3-D stress field up to the second order as
½ r13
r12 r22 T ¼ Dk ð þ x3 jÞ þ Det y3;3 þ De Ia v k;a ; r23 T ¼ Ds ðv k;3 þ ea D? v;a Þ;
r33 ¼
DTet Ia
½ r11
ð95Þ
v k;a þ Dt y3;3
ð84Þ
where yk and y3 are the second-order warping functions. It can be shown that the in-plane components yk vanish and the equations governing y3 are
6. Numerical example A general formulation to treat multilayer FG cylindrical shells with properties as functions of x3 has been derived. In the
974
Z. Yifeng et al. / Composite Structures 94 (2012) 966–976
following, a cylindrical bending example of a homogeneous substrate with a thin SiC–Al functionally graded coating under sinusoidal pressure on the top surface will be used to demonstrate the performance of the present theory. The example is a double-layer coating/substrate system as discussed in Ref. [4] including a homogeneous substrate with a thin functionally graded coating (see Fig. 3). The coordinates are chosen in such a way that x1 2 [0, u], x2 2 [0, 1], and x3 2 [h/2, h/2]. In this system, the Poisson’s ratios for the substrate and the coating ð1Þ ð2Þ ð1Þ ð2Þ are assumed to be constants and equal, mH ¼ mH ¼ mF ¼ mF ¼ 0:3. The shear modulus of the substrate and the coating takes the form of (see Fig. 4): ð1Þ
ð1Þ
GF ¼ GH ;
ð2Þ
ð2Þ
GF ðx3 Þ ¼ nGH ehð1þx3 =hÞ ð2Þ
ð96Þ
ð1Þ
where n ¼ GH =GH is the ratio of the stiffness on the top of the FG coating to that of the homogeneous substrate; h = h/h(2)lnn is the inhomogeneous parameter of the FG coating; h(2) is the coating thickness. In the present analysis, n and h/h(2) are chosen to be 10 and 5, respectively. The top surface of the FGM coating is ceramic rich and the bottom surface is metal rich. The region between the two surfaces is made of the mixture of ceramic–metal materials with continually varying of the volume fractions of the ceramic and metal. The volume fraction of the ceramics Vc is assumed to vary according to a power law as
1 x3 p þ V c ¼ V c þ V þc V c 2 h
Fig. 4. Through-the-thickness variation of shear moduli G in two coating/substrate system.
ð97Þ
where V þ c and V c are the volume fractions of the ceramic on the top and the bottom surfaces, respectively; p P 0 is the volume fraction index, and h the thickness of the shell. The effective elastic moduli of the functionally graded metal– ceramic material are estimated by the Mori–Tanaka scheme as
K Km Vc ¼ ; c K m K c K m 1 þ ð1 V c Þ K Kþ4 l =3 m m
ð98Þ
l lm Vc ¼ lc lm 1 þ ð1 V c Þ llcmþslm1
Fig. 5. Through-the-thickness variation of effective Young’s modulus of FG metal– ceramic material estimated by Mori–Tanaka scheme for different values of p.
where s1 = lm(9Km + 8lm)/(6Km + 12lm), Km and lm represent the bulk and shear modulus of the metal and Kc and lc represent the bulk and shear modulus of the ceramic material [7,20]. The Young’s modulus, E(x3), and Poisson’s ratio, m(x3), are related to effective bulk and shear moduli by
K ðx3 Þ ¼
Eðx3 Þ ; 3½1 2mðx3 Þ
l ðx3 Þ ¼
E ð x3 Þ 2½1 þ mðx3 Þ
ð99Þ
For the purpose of illustration, we choose the properties of the constituent materials as Em = 70 GPa, vm = 0.3 for Al and Em = 427 GPa, vm = 0.17 for SiC. Figs. 5 and 6 plots through-the-thickness var-
Fig. 6. Through-the-thickness variation of Poisson’s ratio of FG metal–ceramic material estimated by Mori–Tanaka scheme for different values of p.
iation of effective Young’s modulus and Poisson’s ratio for various material index p. The shell is assumed to be simply supported and subjected to the following boundary conditions: Fig. 3. Sketch of a two double-layer coating/substrate cylindrical shell.
x1 ¼ 0;
p=3 : r11 ¼ 0; U 2 ¼ U 3 ¼ 0
ð100Þ
975
Z. Yifeng et al. / Composite Structures 94 (2012) 966–976
Fig. 7. Distribution of the longitudinal deflection U 1 vs the thickness coordinate at x1 = 0.
Fig. 8. Distribution of the transverse deflection U 3 vs the thickness coordinate at x1 = p/6.
Fig. 9. Distribution of the longitudinal normal stress r11 vs the thickness coordinate at x1 = p/6.
The top surface of the coating at x3 = h/2
issubjected to a sinusoidally distributed pressure s3 ¼ p20 sin nRp1x/1 ðR1 ¼ 10; n ¼ 2; u ¼ p=3Þ. There is no body force and the bottom surface is traction free. The stress and displacement components are nondimensionlized as r ij ¼ rij =p0 ; U i ¼ U i GHð2Þ =ðp0 hÞ. The through-thickness variations of various displacement and stress components for a thick (R1u/h = 3) double-layer FG coating/substrate system are plotted in Figs. 7–12. It can be observed
Fig. 10. Distribution of the in-plane shear stress x1 = 0.
r12 vs the thickness coordinate at
Fig. 11. Distribution of the transverse shear stress at x1 = 0.
r13 vs the thickness coordinate
Fig. 12. Distribution of the transverse normal stress r33 vs the thickness coordinate at x1 = p/6.
from these figures that all displacement and stress components except for U 3 match well with the 3-D exact solutions. It seems that there exist a constant shift (about 5% of the maximum value) between the results of the present model and the 3-D solution for the transfers displacement U 3 . This may be attributed to that the shell model is reduced from the original 3-D model and some information which cannot be captured by a 2-D model are lost during
976
Z. Yifeng et al. / Composite Structures 94 (2012) 966–976
the dimensional reduction process. Nevertheless, both results show similar variation trend for U 3 . It is important to note that by adopting FG coating the shear moduli G is continuous at the interface, the in-plane stress components such as r11 and r12 , therefore, are smooth throughout of the total thickness. In addition, no sharp peak appears for transverse shear stress r13 . Results presented in Figs. 7–12 clearly show that using FG coating will eliminate the mismatch of interface stress components to reduce the risk of cracking and rebounding of the coating. 7. Conclusions Using the variational asymptotic method (VAM), an efficient high-fidelity shell model for multilayer FG cylindrical shell has been developed. The VAM is used to systematically reduce the original, nonlinear 3-D model to a series of 2-D models in terms of the small parameter. No apriori assumptions have been adopted during the derivation. The theory is applicable to FG shells with material properties either being constant or changing continuously in each layer. A cylindrical bending example of a homogeneous substrate with a thin SiC–Al functionally graded coating under sinusoidal pressure on the top surface has been used to validate the present model. Excellent agreement between the present model and the 3-D exact solution have been observed for the 3-D distribution of displacements and stresses through the shell thickness. Moreover, the present model is valid for large displacements and global rotations and can capture all the geometric nonlinearity of a shell when the strains are small. The present paper has built on the second author’s previous work [15,21] with the following new contributions: 1. The present work has the capability of analyzing multilayer FGM cylindrical shell with material properties as functions of transverse locations or constants while in previous work these properties are treated as constants for each layer. 2. Interface continuity conditions are explicitly derived and solved to obtain the multilayer solutions. 3. Simplifications have been made in deriving B, C, D matrices. 4. Explicit analytical solutions for the second-order approximation of warping functions have been provided.
Acknowledgments This research is supported by the Fundamental Research Funds for the Central Universities, China, under Grant No. 0218005202017 and China Postdoctoral Science Foundation, under Grant No. 2011M501386. The views and conclusions contained herein are those of the writers and should not be interpreted as
necessarily representing the official policies or endorsement, either expressed or implied, of NSFC. References [1] Muller E, Drasar C, Schilz J, Kaysser WA. Functionally graded materials for sensor and energy applications. Mater Sci Eng A – Struct Mater Prop Microstruct Process 2003;362(1–2):17–39. [2] Mishnaevsky LL. Functionally gradient metal matrix composites: numerical analysis of the microstructure-strength relationships. Compos Sci Technol 2006;66(11–12):1873–87. [3] Shaw LL. Thermal residual stresses in plates and coatings composed of multilayered and functionally graded materials. Composites Part B – Eng 1998;29(3):199–210. [4] Kashtalyan M, Menshykova M. Three-dimensional elastic deformation of a functionally graded coating/substrate system. Int J Solids Struct 2007;44(16):5272–88. [5] Das M, Barut A, Madenci E, Ambur DR. A triangular plate element for thermoelastic analysis of sandwich panels with a functionally graded core. Int J Numer Meth Eng 2006;68(9):940–66. [6] Pan E, Heyliger PR. Exact solutions for magneto–electro-elastic laminates in cylindrical bending. Int J Solids Struct 2003;40(24):6859–76. [7] Vel SS, Batra RC. Exact solution for thermoelastic deformations of functionally graded thick rectangular plates. AIAA J 2002;40(7):1421–33. [8] Qatu MS, Sullivan RW, Wang W. Recent research advances on the dynamic analysis of composite shells: 2000–2009. Compos Struct 2010;93(1):14–31. [9] Maalawi KY. A generalized formulation for radially graded composite cylinders/ rings with maximized buckling stability limits. J Eng Appl Sci 2010;23(5):129–48. [10] Khazaeinejad P, Najafizadeh MM, Jenabi J, Isvandzibaei MR. On the buckling of functionally graded cylindrical shells under combined external pressure and axial compression. J Press Ves Technol–Trans ASME 2010;132(6):1551–7. [11] Shariyat M. Dynamic buckling of suddenly loaded imperfect hybrid FGM cylindrical shells with temperature-dependent material properties under thermo–electro-mechanical loads. Int J Mech Sci 2008;50(12):1561–71. [12] Yas MH, Shakeri M, Heshmati M, Mohammadi S. Layer-wise finite element analysis of functionally graded cylindrical shell under dynamic load. J Mech Sci Technol 2011;25(3):597–604. [13] Bian ZG, Chen WQ, Lim CW, Zhang N. Analytical solutions for single- and multi-span functionally graded plates in cylindrical bending. Int J Solids Struct 2005;42(24–25):6433–56. [14] Soldatos KP, Watson P. A method for improving the stress analysis performance of one- and two-dimensional theories for laminated composites. Acta Mech 1997;123(1–4):163–86. [15] Yu WB, Hodges DH, Volovoi VV. Asymptotic generalization of Reissner– Mindlin theory: accurate three-dimensional recovery for composite shells. Comput Meth Appl Mech Eng 2002;191(44):5087–109. [16] Yu WB, Hodges DH, Volovoi VV. Asymptotically accurate 3-D recovery from Reissner-like composite plate finite elements. Comput Struct 2003;81(7):439–54. [17] Hodges DH, Atilgan AR, Danielson DA. A geometrically nonlinear-theory of elastic plates. J Appl Mech–Trans ASME 1993;60(1):109–16. [18] Danielson DA, Hodges DH. Nonlinear beam kinematics by decomposition of the rotation tensor. J Appl Mech–Trans ASME 1987;54(2):258–62. [19] Atilgan AR, Hodges DH. On the strain-energy of laminated composite plates. Int J Solids Struct 1992;29(20):2527–43. [20] Gilhooley DF, Batra RC, Xiao JR, McCarthy MA, Gillespie JW. Analysis of thick functionally graded plates by using higher-order shear and normal deformable plate theory and MLPG method with radial basis functions. Compos Struct 2007;80(4):539–52. [21] Yu WB. Mathematical construction of a Reissner–Mindlin plate theory for composite laminates. Int J Solids Struct 2005;42(26):6680–99.
Contents lists available at SciVerse ScienceDirect
Composite Structures journal homepage: www.elsevier.com/locate/compstruct
Variational asymptotic modeling of the multilayer functionally graded cylindrical shells q Zhong Yifeng a,b,⇑, Chen Lei a,b, Wenbin Yu c a
College of Civil Engineering, Chongqing University, Chongqing, PR China Key Laboratory of New Technology of Construction of Cities in Mountain Area, Chongqing University, Ministry of Education, Chongqing 400045, PR China c Utah State University, Logan, UT 84322-4130, USA b
a r t i c l e
i n f o
Article history: Available online 13 October 2011 Keywords: Asymptotic series Functionally graded material Cylindrical shell Variational asymptotic method Simplified model
a b s t r a c t An efficient high-fidelity shell model is developed for heterogeneous multilayer cylindrical shells made of functionally graded material by using the variational asymptotic method (VAM). Taking advantage of the smallness parameters inherent in the shell structure, the VAM is applied to rigorously decouple the 3-D, anisotropic elasticity problem into a 1-D through-the-thickness analysis and a 2-D shell analysis. The through-the-thickness analysis servers as a link between the original 3-D analysis and the shell analysis by providing a constitutive model for the shell analysis and recovering the 3-D field variables in terms of global responses calculated by the shell analysis. The present model is valid for large displacements and global rotations and can capture all the geometric nonlinearity of a shell when the strains are small. A cylindrical bending example of a homogeneous substrate with a thin SiC–Al functionally graded coating under sinusoidal pressure on the top surface is used to validate this model. Ó 2011 Elsevier Ltd. All rights reserved.
1. Introduction Functionally graded materials (FGM) are microscopically inhomogeneous composites usually made from a mixture of metals and ceramics. By gradually varying the volume fraction of consistent material, its material properties exhibit a smooth and continuous change from one surface to another. The surface with high ceramic constituents can provide superior thermal-resistance for high temperature environments while the surface with high metal constituents offers strong mechanical performance. This kind of material arrangement reduces the risk of catastrophic fracture and mitigates thermal stress concentration under extreme environments. The advantage of FGM has been used to address a large variety of application fields, such as graded thermoelectrics and dielectrics, piezoelectrically graded materials for ultrasonic transducers, and tungsten–copper composites for high current connectors and diverter plates, to name but a few [1,2]. Recently, the concept of FGM is actively explored in the multilayered design of thermal coatings as well as sandwich panels to overcome the mismatch of the thermomechanical properties between the coating and the substrate or between the surface panels and the core [3–5].
q
This paper is contributed by all authors.
⇑ Corresponding author at: College of Civil Engineering, Chongqing University, Chongqing, PR China. Tel.: +86 023 65123292. E-mail address: [email protected] (Z. Yifeng). 0263-8223/$ - see front matter Ó 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.compstruct.2011.10.005
To use FGM effectively, an efficient and accurate model for structures made of such materials needs to be developed. These structures are characterized by one of their dimensions (the thickness) being much smaller than the other two. Although all structures made of FGM can be described using 3-D continuum mechanics, exact solutions exist only for a few specific problems with very idealized material types, geometry, and boundary conditions [6,7]. For more realistic cases, 3-D numerical simulation tools such as ANSYS and ABAQUS are often used to find approximate solutions. However, this approach is computation intensive and usually used in the detailed analysis due to its prohibitive computational cost. In view of the fact that the thickness is small, analysis of such structures can be simplified using 2-D models. Although many 2-D models have been developed to analyze FGM shells treating different topics, most of them rely on some apriori kinematic assumptions [8]. Examples may be found in the application of classical lamination theory (CLT) in the problem of maximization of buckling stability limits of Functional Graded (FG) rings/long cylindrical shells subjected to external pressure [9], the use of the first-order sheardeformation theory (FOSDT) in studying the stability problem of a circular cylindrical shell composed of FGM with elasticity modulus varying continuously in the thickness direction under combined external pressure and axial compression loads [10], the implementation of the high-order shell theory in dynamic buckling of imperfect FGM cylindrical shells with integrated surface-bonded sensor and actuator layers subjected to some complex combinations
Z. Yifeng et al. / Composite Structures 94 (2012) 966–976
967
Nomenclature
v
a column matrix containing 2-D generalized strains v ¼ ½T jT T
dW; dW 2D
, j Cij
c l l⁄, K⁄ X
x oX eab, Kab ck Cij D eijk fi, ma Fij Gk gk
3-D and 2-D virtual work, respectively column matrix of 2-D generalized strains three-dimensional strain tensor column matrix of transverse strain measures characteristic magnitude of the elastic constants effective shear and bulk modulus estimated by Mori– Tanaka scheme, respectively the reference surface inertial angular velocity of triad Bi boundary of the reference surface 2-D generalized shell strains through-the-thickness average of in-plane warping functions matrix of direction cosines 3-D material matrix the permutation symbol generalized forces and moments, respectively mixed-basis component of the deformation gradient tensor covariant basis vector of the deformed shell contravariant base vector of the undeformed shell
of thermo–electro-mechanical loads [11], and the use of Layerwise theories in analysis of a FGM cylindrical shell under dynamic loads [12]. CLT ignores transverse shear effects and provides reasonable results only for very thin shells. Moreover, in CLT, both plane strains and plane stresses are assumed which we know will not be true at the same time for materials having nonzero Poisson’s ratios. A number of shear deformation theories have been developed to overcome some drawbacks of CLT, with the simplest of which being FOSDT (equivalent to Reissner–Mindlin theory for shells made of isotropic homogeneous materials), where a constant distribution of shear strain through the thickness is assumed and a shear correction factor is required to account for the deviation of the real shear strain from the assumed constant one. The dependence of the shear correction factor on the geometry and material of the shell makes it difficult to guarantee the accuracy of FOSDT. By expanding the displacement field of the shell using higher-order polynomials, higher-order shear deformation shell theories are developed, which can account for both transverse normal and shear deformations without relying on shear correction factors. However, as indicated by Bian et al. [13] and Das et al. [5] for laminated shells, models based on higher-order theories cannot capture the discontinuous slope of in-plane and transverse displacement components in the thickness direction. By extending a generalized refined theory (referred as Soldatos shell theory [14]) which incorporates shape functions to guarantee continuity of transverse shear stresses at interfaces, these authors provide analytical solutions for simply- and multiply- spanned FG shells under cylindrical bending. Layerwise theories can produce reasonable results at the cost of complex models and expensive computation. Despite of the popularity of the aforementioned methods in analyzing many FG shells and FG laminated shells, these approaches have two major disadvantages: (1) the apriori assumptions which are naturally extended from the analysis of isotropic homogeneous structures cannot be easily justified for heterogeneous and anisotropic structures, such as FGMs; (2) it is difficult
h, R1 K, K2D l La
thickness and radius of the shell, respectively 3-D and 2-D kinetic energy, respectively characteristic wavelength of the shell deformation integration constants for determining the relationship between v and ck M, N shell moment resultants and force resultants, respectively n a small parameter used to denote the order of strains pi applied body force Qi tractions applied on lateral surfaces R strain measures of Reissner–Mindlin shell U free energy per unit area wi warping functions for the normal-line element xi Cartesian coordinate base vector of the coordinate system before and after bi, Bi deformation, respectively b ^r; R position vector of a material point before and after deformation, respectively triad vector for deformed Reissner–Mindlin shell Bi b R through-the-thickness average of the position vector R kk, k3, K3 Lagrange multipliers for enforcing in-plane and outplane warping constraints
for an analyst to determine the accuracy of the result and which assumption should be chosen for efficient and accurate analysis for a particular shell. In this paper, an efficient high-fidelity shell model for multilayer FG cylindrical shell is constructed by the variational asymptotic method (VAM). The original 3-D nonlinear elasticity problem is casted in an intrinsic form so that the theory is applicable for arbitrarily displacement and global rotation with the restriction that to the strain field is small. VAM is then used to rigorously split the 3-D elasticity problem into two problems: a nonlinear, 2-D, shell analysis over the reference surface to obtain the global deformation and a linear, 1-D, through-the-thickness analysis to provide the 2-D generalized constitutive law and the recovering relations to approximate the original 3-D results. The nonuniqueness of asymptotic theory correct up to a certain order is used to cast the obtained asymptotically correct second-order free energy into a Reissner– Mindlin type model to account for transverse shear deformation. The present theory extends a simple and accurate model developed recently for composite laminates by Yu et al. [15,16] so that the multilayer FG shell can be treated in the same framework. Since the procedure is quite similar, the authors have chosen to repeat some formulae and texts from their previous publications in order to make the present paper more self-contained. 2. 3-D Formulation The elastodynamic behavior of a solid is governed by the extended Hamilton principle:
Z
t2
½dðK UÞ þ dWdt ¼ 0
ð1Þ
t1
where t1 and t2 are arbitrary fixed times, K is the kinetic energy, U is a term related with the internal energy and dW is the virtual work of applied loads. The overbar is used to indicate that the virtual work does not necessarily represent variation of functionals.
968
Z. Yifeng et al. / Composite Structures 94 (2012) 966–976
ga ¼ aa þ x3 b3;a ;
g3 ¼ b3
ð7Þ
From the differential geometry of surface and following Ref. [17], the derivatives of unit vectors bi,a can be expressed as
bi;a ¼ Aa ðka2 b1 þ ka1 b2 þ ka3 b3 Þ bi
ð8Þ
where Kab refers to the usual out-of-plane curvatures and k12 = k21 = 0 since the coordinates are the lines of curvatures. The geodesic curvatures ka3 can be expressed in terms of the Lamé parameters as
k13 ¼
A1;2 ; A1 A2
k23 ¼
A2;1 A1 A2
ð9Þ
One can obtain the explicit expressions for the 3-D covariant and contravariant base vectors by using Eq. (8):
ga ¼ Aa ð1 þ x3 kaa Þba ; ga ¼ Fig. 1. Schematic of shell deformation.
A cylindrical shell is a 3-D body with a relatively small thickness h and a smooth reference surface (see Fig. 1). The geometry of the reference surface can be described by a set of arbitrary curvilinear coordinates, xa. (Here and throughout the paper, Greek indices assume values 1 and 2 while Latin indices assume 1, 2, and 3; repeated indices are summed over their range except where explicitly indicated). Without loss of generality, the lines of curvatures are chosen to be the curvilinear coordinates to simplify the formulation. To uniquely represent an arbitrary point in a 3-D medium, another coordinate off the reference surface is needed. It is a natural and convenient choice to take the third coordinate x3 as the transverse normal coordinate. Letting bi denote the unit vector along xi for the undeformed shell, the position of any material point in the undeformed configuration can then be described by its position vector ^r from a fixed point O as
^rðx1 ; x2 ; x3 ; tÞ ¼ rðx1 ; x2 ; tÞ þ x3 b3
ð3Þ
where the angle brackets denote the definite integral through the thickness of the shell and will be used throughout the rest of the development. The 2-D base vectors are defined conventionally as
aa ðx1 ; x2 Þ ¼ r;a
ð4Þ
@ðÞ where ðÞ;a ¼ @x . a The so-called Lamé parameters can be obtained from this definition as
Aa ðx1 ; x2 Þ ¼
pffiffiffiffiffiffiffiffiffiffiffiffiffiffi aa aa
ð5Þ
Note here the summation convention is not applied here because a is not a dummy index; this rule will apply in similar situations throughout the rest of development. Then the unit vectors along coordinates xa can be obtained as:
ba ðx1 ; x2 Þ ¼
aa Aa
ð10Þ
ba Aa ð1 þ x3 kaa Þ
When the shell deforms, the particle that had position vector ^r b in the deformed in the undeformed state now has position vector R shell. The latter can be uniquely determined by the deformation of the 3-D body. Similarly, another triad Bi is introduced for the deformed configuration. Note that the Bi unit vectors are just tools to express vectors and tensors in their component form during the derivation. They are not necessarily tangent to the coordinates of the deformed shell. The relation between Bi and bi can be specified by an arbitrarily large rotation specified in terms of the matrix of direction cosines C(x1, x2) so that
Bi ¼ C ij bj ;
C ij ¼ Bi bj
ð11Þ
subjecting to the requirement that Bi is coincident with bi when the structure is undeformed. b can be represented as Now the position vector R
b i ðx1 ; x2 ; x3 ; tÞ ¼ Rðx1 ; x2 ; tÞ þ x3 B3 ðx1 ; x2 ; tÞ R þ wi ðx1 ; x2 ; x3 ; tÞBi ðx1 ; x2 ; tÞ
ð12Þ
ð2Þ
where r is the position vector from O to the point located by xa on the reference surface at a specific time t. When the origin of the coordinate system is located on the middle surface, it naturally follows that
h^rðx1 ; x2 ; x3 ; tÞi ¼ hrðx1 ; x2 ; tÞ
g3 ¼ g3 ¼ b 3 ;
ð6Þ
It is obvious that the unit vectors bi form an orthogonal triad. The 3-D covariant base vectors gi ¼ ^r;i for the chosen coordinate system can be obtained by using Eq. (4) as
where wi are the warping functions introduced to accommodate all possible deformations. Eq. (12) is six times redundant because of the way warping introduced. Proper definitions of R and Bi are needed to introduce six constraints to ensure a one-to-one mapping between R, Bi and wi. R can be defined similarly as Eq. (3) to be the average position through the thickness, from which it follows that the warping functions must satisfy the three constraints
hwi ðx1 ; x2 ; x3 ; tÞi ¼
ck ; 0
ck ¼
c1 c2
ð13Þ
where c1 and c2 are functions of the in-plane coordinates xa and time t, introduced for providing free variables for the construction of an optimal Reissner–Mindlin model. Two other constraints can be specified by taking B3 as the normal to the reference surface of the deformed shell. It should be noted that this choice has nothing to do with the well-known Kirchhoff hypothesis. In the Kirchhoff assumption, no local deformation of the transverse normal is allowed. However, all possible deformation are allowed by using the warping functions in present derivation. Because Ba can freely rotate around B3, the last constraint can be introduced as
B1 R;2 ¼ B2 R;1
ð14Þ
Based on the concept of decomposition of rotation tensor [18], the Jauman–Biot–Cauchy strain components for small local rotation are given by
969
Z. Yifeng et al. / Composite Structures 94 (2012) 966–976
1 2
Cij ¼ ðF ij þ F ji Þ dij
ð15Þ U¼
X
where dij is the Kronecker symbol, and Fij is the mixed-basis component of the deformation gradient tensor such that
F ij ¼ Bi Gk gk bj
ð16Þ
@b R where Gk ¼ @x is the covariant basis vector of the deformed configuk ration, which can be expressed in terms of the 2-D generalized strains following Ref. [19] as:
R;a ¼ Aa ðBa þ eab Bb Þ
ð17Þ
Bi;a ¼ Aa ðK a2 B1 þ K a1 B2 þ K a3 B3 Þ Bi
ð18Þ
where eab are the 2-D in-plane strains and of which the order is denoted by n, and Kij are the curvatures of the deformed surface which are the summation of curvatures of undeformed geometry kij and curvatures introduced by the deformation jij of which the order is denoted by n/h. Using this definition, one can show that Eq. (14) implies e12 = e21, i.e.
B1 R;2 B2 R;1 ¼ A2 A1
ð19Þ
which can serve as another constraint to specify the global rotation of the triad Bi. For geometrically nonlinear analysis, the strains are small compared to unity and warpings are of the order of strain or smaller, which have the effect of removing all the terms that are products of the warping and the generalized strains. With the help of Eqs. (12), (15), (16), (17) and (18), we can express the 3-D strain field as
Ce ¼ þ x3 j þ Ia wk;a ; @ðÞ where ðÞ;i ¼ @x i
1 Aa
2Cs ¼ wk;3 þ ea w3;a ;
; ðÞk ¼ ½ ðÞ1
Ct ¼ w3;3
ð20Þ
2C23 T ;
Ct ¼ C33 ;
¼ ½ e11 2e12 e22 T ; j ¼ ½ K 11 K 12 þ K 21 K 22 T 2
3
1 0 6 7 I 1 ¼ 4 0 1 5; 0 0 T
e1 ¼ ½ 1 0 ;
2
ð21Þ
ð24Þ where De, Des, Det, Ds, Dst, Dt are the corresponding partition matrices of the 6 6 elastic material matrix D at constant electric field, X is the surface stretched by the undeformed reference surface and / is geometric correction as 2
g g2 g3 h /¼ 1 ¼ 1 þ x3 ðk11 þ k22 Þ þ O 2 ja1 a2 j l
! ð25Þ
To calculate the kinetic energy, the absolute velocity of a generic point in the structure is obtained by taking a time derivative of Eq. (12) as
b ¼ V þ xðn þ wÞ þ w _ V
ð26Þ
where ( ) is the partial derivative with respect to time, V is the absolute velocity of a point in the deformed reference surface, x is the inertial angular velocity of Bi bases, and the notation () forms an antisymmetric matrix from a vector according to ð Þ ¼ eijk ðÞk using b ; V; x; w denote column the permutation symbol eijk. The symbols V matrices containing the components of corresponding vectors in Bi bases, and n ¼ ½ 0 0 x3 T . The kinetic energy of the shell structure can be obtained by
K¼
1 2
Z
q Vb T Vb dV ¼ K 2D þ K
ð27Þ
V
where q is the mass density and
Z 1 fnV þ xT jxÞdX V T V þ 2xT l ðl 2 X Z 1 ~ w þ wÞ ~ w þ wÞ ~ nÞT ðx ~ w þ wÞdV _ T ðx _ þ 2ðV þ x _ K ¼ q½ðx 2 V K 2D ¼
; l n; j are inertial constants commonly used in shell dynamwhere l ics, which can be trivially obtained through simple integrals over the thickness as:
0 ð22Þ
0
0
b þ s dR b þ þ b dR b ÞdX þ ðhP d Ri
X
1 þ x3 k11
e12 þ x3 j21 þ w1;2 w2 k23
e12 þ x3 j12 þ w2;1 þ w1 k13 þ 1 þ x3 k22 1 þ x3 k11 e22 þ x3 j22 þ w2;2 þ w3 k22 þ w1 k23 C22 ¼ 1 þ x3 k22 w3;1 w1 k11 2C13 ¼ w1;3 þ 1 þ x3 k11 w3;2 w2 k22 2C23 ¼ w2;3 þ 1 þ x3 k22 C33 ¼ w3;3 ð23Þ With the knowledge of mechanical strains, the energy enthalpy can be expressed as:
ð29Þ
The virtual work of the structure can be calculated as Z Z
dW ¼
T
e11 þ x3 j11 þ w1;1 þ w3 k11 w2 k13
2C12 ¼
DTst
l ¼ hqi; ln ¼ ½ 0 0 hx3 qi T ; 2 2 3 x3 q 0 0 6 7 j¼4 0 x23 q 0 5
The explicit expressions of Cij can be written as
C11 ¼
Ds
9 3T 8 + > Z < Ce > = 7 Dst 5 2Cs / dX ¼ U A /dX > > X : ; Ct Dt
Det
ð28Þ
3
0 0 6 7 I 2 ¼ 4 1 0 5; 0 1 e2 ¼ ½ 0 1
Des
ðÞ2 T , and
Ce ¼ ½ C11 2C12 C22 T ; 2Cs ¼ ½ 2C13
8 9 2 *> Ce >T De = 1 < 6 DT 2 Cs 4 es > 2 > : ; Ct DTes
Z
b hQ d Rids
ð30Þ
@X
where oX denotes the boundary of the reference surface, and ðÞ ¼ ðÞjx3 ¼h=2 . P = PiBi is the applied body force, s = siBi, b = biBi are tractions applied on the top and bottom surfaces, respectively, and Q = QiBi is the applied tractions along the lateral surfaces. The Lagrangian variation of the displacement field which can be expressed as
b ¼ dq i Bi þ x3 dB3 þ dwi Bi þ wj dBj dR
ð31Þ
in which the virtual displacement and rotation are defined by
2 B1 þ dw 1 B2 þ dw 3 B3 Þ Bi i ¼ dR Bi ; dBi ¼ ðdw dq
ð32Þ
i contain the components of the virtual displacei and dw where dq ment and rotation in the Bi system, respectively. Since the warping functions are small, one may safely ignore b products of the warping functions and the virtual rotations in d R and obtain the virtual work due to applied loads as
970
Z. Yifeng et al. / Composite Structures 94 (2012) 966–976
dW ¼ dW 2D þ dW
ð33Þ
where
dW 2D ¼
Z
a ÞdX þ i þ ma dw ðfi dq
X
dW ¼
Z X
Z
a Þds i þ hx3 Q a idw ðhQ 1 idq
ð34Þ
where UA0 can be obtained from Eq. (24) by dropping the derivatives with respect to xi in Eq. (20), resulting in
D E 2U A0 ¼ 2 ð þ x3 jÞT ðDes wk;3 þ Det w3;3 Þ þ wTk;3 Dst w3;3 þ hð þ x3 jÞT De ð þ x3 jÞi þ wTk;3 Ds wk;3 þ wT3;3 Dt w3;3
@X
hpi dwi i þ si dwþi þ bi dwi dX þ
Z
hQ i dwi ids
þ oððh=lÞ0 n2 Þ
ð35Þ
ð40Þ
X
with the generalized forces fi and moments ma defined as
h fi ¼ h/i i þ si þ bi ; ma ¼ hx3 /a i þ ðsa ba Þ 2
ð36Þ
T where ðÞk ¼ ½ ðÞ1 ðÞ2 ; ðÞ;3 ¼ @ðÞ=@x3 . It is clear that the warping functions wi only appear in UA0, which means wi can be solved from the following much simpler variational statement:
The second integration in Eq. (35) accounts for virtual work done through warping functions along the lateral boundaries of the shell. This term is necessary for the edge-zone problem, which is an important subject in its own right and beyond the scope of the present paper. Therefore, we will drop this term hereafter. With the knowledge of Eqs. (33), (24) and (27), the Hamiltons principle in Eq. (1) becomes
along with the constrains in Eq. (13). This variational statement has the following Euler–Lagrange equations:
Z
t2
½dðK 2D þ K UÞ þ dW 2D þ dW dt ¼ 0
ð37Þ
t1
3. Dimensional reduction The dimensional reduction from the original 3-D formulation to a 2-D shell model can only be done approximately. One way to do it is to take advantage of the small parameters in the formulation to construct the 2-D formulation so that the reduced model can achieve the minimum loss of accuracy in comparison to the original 3-D formulation. In order to apply VAM, we first need to assess the order of quantities in terms of small parameters. The quantities of interest have the following orders:
hp3 s3 b3 lðh=lÞ2 n; hpa sa ba lðh=lÞn; Q a ln; Q 3 lðh=lÞn
ð38Þ
where l represents the order of elastic material constants.
;3
ð þ x3 jÞ
T
Dkes
T
T þ wkk;3 Dks þ wk3;3 Dkst ¼ kk
ð42Þ
Z d K 2D U A0 /dX þ dW 2D dt ¼ 0 X
where k = 1, 2, 3, . . . , N with N denoting the total number of layers, kk and k3 are Lagrange multipliers corresponding to the in-plane and out-of-plane constraint equations expressed in Eq. (13); ()k denotes functions () for the kth layer. The boundary conditions as well as the interlaminar continuous conditions are:
½ð þ x3 jÞT Det þ ðwk;3 ÞT Dst þ w3;3 Dt þ ¼ 0; h i ð þ x3 jÞT Des þ ðwk;3 ÞT Ds þ w3;3 DTst ¼ 0
ð43Þ
and
h i ½w3 ¼ 0; ð þ x3 jÞT Det þ ðwk;3 ÞT Dst þ w3;3 DTt ¼ 0; X h i i ½wk ¼ 0; ð þ x3 jÞT Des þ ðwk;3 ÞT Ds þ w3;3 DTst ¼ 0
ð44Þ
Xi
where Xi denotes the interfaces between the ith layer and i + 1th layer; i = 1, . . . , N 1; the bracket [] denotes the jump of the enclosed argument on the interface. From these conditions, we can solve wi,3 from Eq. (42):
T
1 bk D bk wkk;3 ¼ ð þ x3 jÞT D ; es s
ð45Þ
b k =D bk wk3;3 ¼ ð þ x3 jÞT D et t with the hatted quantities being expressed as:
T b k ¼ Dk D b k Dk =Dk ; D es et es st t
1 k k k k b D et ¼ Det Des Ds =Dkst ;
T 1 b k ¼ Dk Dk bk D D Dkst t s t st
To clearly illustrate the application of VAM for FGM shells, a classical FGM shell model is first constructed. By applying VAM, the zeroth-order approximation of the variational statement in Eq. (37) can be obtained as
t1
T ð þ x3 jÞT Dket þ wkk;3 Dkst þ wk3;3 Dkt ¼ k3 ;
ð46Þ
Substituting Eq. (45) into Eq. (40), the zeroth-order strain energy can be expressed as
3.1. Zeroth-order approximation
t2
ð41Þ
;3
So far, a 3-D displacement-based formulation for FG shells has been developed in terms of 2-D displacements (represented by R r), rotations (represented by bi and Bi), and 3-D warping functions (wi). If we attempt to solve this problem directly, we will meet the same difficulty as solving any full 3-D problem with the additional difficulty coming from the anisotropy and heterogeneity of functional graded materials. The main complexity comes from the unknown 3-D warping functions wi. The common practice in the literature is to assume wi in terms of 2-D displacements and rotations to straightforwardly reduce the original 3-D continuum model into a 2-D shell model. However, for shells made with general anisotropic and heterogeneous materials such as FGMs, the imposition of such assumptions may introduce significant errors. Fortunately, VAM provides a useful technique to obtain wi through an asymptotical analysis of the variational statement in Eq. (37) in terms of small parameters inherent in the problem.
Z
dU A0 ¼ 0
ð39Þ
2U A0 ¼
j
T " b A bT B
b B b D
#
j
þo
0 ! h n2 l
ð47Þ
with
D E b ¼ hD b ¼ hx3 D b ¼ x2 D b e i; D b b e i; B A 3 e ; b e ¼ De Des D1 D bt b T Det D b T =D D e es et
ð48Þ
971
Z. Yifeng et al. / Composite Structures 94 (2012) 966–976
With UA0 expressed in Eq. (47), the original 3-D problem in Eq. (1) has been rigorously reduced to a 2-D formulation in Eq. (39) which approximates the original problem asymptotically correct to the order of (h/l)0. If we define the force resultants N and moment resultants M conjugate to and j, such that
N¼
@U A0 ; @
@U A0 @j
M¼
ð49Þ
then we obtain a 2-D constitutive model for the classical shell analysis of FGM shells as
N M
"
¼
b A bT B
b B b D
#
j
ð50Þ
It is clear that although the shell is made of FGM, the 2-D shell model of the zeroth-order remains the same. Despite the similarity with CLT, the present model is asymptotically correct and has the following features in contrast with CLT: 1. The normal line of undeformed shell does not remain straight and normal to the deformed shell; rather, it deforms in both the normal and in-plane directions in response to shell deformation ( and j). 2. This model can handle general FGMs with full anisotropy. 3. It can be easily observed that neither the normal strain nor the transverse shear strains vanish. The transverse normal and shear stresses can be shown to vanish, which are not assumed a priori but direct consequences of the model. For the zeroth-order approximation, the 3-D strain field can be recovered using Eq. (20) by neglecting those terms of order higher than (h/l)0: 0 e
C ¼ þ x3 j;
0 s
2C ¼ wk;3 ;
0 t
C ¼ w3;3
ð51Þ
The 3-D stress field can be computed using the original 3-D constitutive relations, which can be derived from Eq. (24) as
8 9 < re =
2
De T D rs ¼ 6 4 es : ; rt DTet
Des Ds DTst
38 0 9 Det > < Ce > = Dst 7 5 2C0s > > D : 0 ; t
ð53Þ
Up to this point, the zeroth-order solution wi in Eq. (45) as well as the 2-D strain energy in Eq. (47) are valid for FGM shells with fully populated 6 6 material matrix D. As most real materials have at least a monoclinic symmetry about their own mid-plane, hereafter, monoclinic material matrix characterized by 13 independent material properties, implying Des = 0 and Dst = 0, will be used for the rest of derivation. This leads to much simpler expressions for the zeroth-order approximation of warping function:
where
wk3 ¼ Dk? v
T Dk? ¼ Dkes =Dkt
To obtain the first-order approximation, we simply perturb the zeroth-order warping functions as
hn ; wk ¼ v k þ o l
hn w3 ¼ v 3 þ D? v þ o l
where Dk ¼ De Det DTet =Dt and v3 is decoupled from vk. Considering the warping constraint in Eq. (13), v3 only has a trivial solution. The stationary conditions of the functional given in Eq. (57) are
Dks v k;3 þ Dks ea Dk? v;a ¼ Dka;3 v;a þ g k;3 þ kk ;3
Ds v k;3 þ Ds ea D? v;a ¼ sk x ¼þh=2
3 Ds v k;3 þ Ds ea D? v;a ¼ bk
ð58Þ
x3 ¼h=2
where
h Dka;3 ¼ ITa Dkk
i x3 Dkk ;
g k;3 ¼ /kk
½v k ¼ 0; ½Ds ðv k;3 þ ea D? v;a Þ
¼0
ð59Þ
Dks;3 v k þ Dks ea Dk? v;a ¼ Dka v;a þ g k þ kk x3 þ constkk
ð61Þ
and the interface continuity condition in the second equation of Eq. (60) becomes
ð55Þ
v ¼ ½ j T ; hD? i ¼ 0 Note that inter-lamina continuity of Dk? must be maintained due to the continuity of warping functions to produce a continuous displacement field. It should also be pointed out that we have constrained in-plane warpings to vanish (hwki = 0) as the free constants in Eq. (13) will be absorbed in the first-order approximation.
ð60Þ
The inter-lamina continuity on Dka and gk are maintained by the second condition in Eq. (60). It should be mentioned that since the goal is to obtain an interior solution for the shell without considering edge effects, integration by parts with respect to the in-plane coordinates is used hereafter and throughout the rest of the paper, whenever it is convenient for the derivation. Integrating the first equation in Eq. (58), one obtains
ð54Þ
T x3 Dket =Dkt ;
ð56Þ
Substituting Eq. (56) back into Eq. (20), and using Eqs. (29), (34) and (37), one can obtain the leading terms for the first-order approximation of the variational statement in Eq. (37) as D T dP1 ¼ v Tk;3 Ds dv k;3 þ v 3;3 Dt dv 3;3 þ ðe þ x3 jÞ Dk Ia dv k;a E þðea D? v;a ÞT Ds dv k;3 pTk dv k sTk dv þk bTk dv k þ oððh=lÞ2 n2 Þ ð57Þ
Xi
with
wkk ¼ 0;
3.2. First-order approximation
The continuity conditions on the interfaces can be derived as
ð52Þ
Ct
re ¼ ½ r11 r12 r22 T ; rs ¼ ½ r13 r23 T ; rt ¼ r33
In view of Eq. (42), rs and rt are identically zero regardless the material properties, which means the classical model can not predict out-of-plane stresses. In order to predict the out-of-plane stresses, we need to carry out the asymptotic expansion of the 3D variational statement in Eq. (1) into higher order, which is done in next section.
Fig. 2. Schematic figure for laminated cylindrical shell.
972
Z. Yifeng et al. / Composite Structures 94 (2012) 966–976
constkþ1 constkk ¼ k
Z
zkþ1
zk
h ITa Dkk
i x3 Dkk dx3
ð62Þ
constkk
where is the integration constant for kth layer, zk is the x3 coordinate of the bottom of the kth layer (see Fig. 2). Dka ; g ka can be defined as:
Dka ðx3 Þ ¼
k
g a ðx3 Þ ¼
Z
x3
zk
Z
x3
zk
Dka;3 dz ¼
k
g a;3 dz ¼
Z
x3
zk
Z
x3
zk
h
ITa Dkk
i
zDkk dz;
ð63Þ
/kk dz
ð64Þ
Since Bi is uniquely determined by Bi and c, the kinematic identity between the strains measures R of Reissner–Mindlin shell and can be derived as
¼ R Da c;a where
D1 ¼ D2 ¼
0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 0 1
R ¼ ½ e11
Solving Eqs. (61) and (62), the two boundary equations in Eq. (58) as well as Eq. (13), one obtains the following Lagrange multipliers and warping functions:
1 ðsk þ bk hDa;3 iva hg ;3 iÞ h k ¼ Da þ La v;a þ gk
ð72Þ
2e12
T ; T
ð73Þ
;
e22 K 11 K 12 þ K 21 K 22 T
Now one can rewrite the strain energy expressed in Eq. (68), correct to the orders of interest according to Eq. (72), in terms of strains of the Reissner–Mindlin model as
kk ¼
ð65Þ
2PR ¼ RT AR 2RT ADa c;a þ RT;a Ba R;a þ 2RT;1 CR;2 2RT F
v
ð66Þ
The generalized Reissner–Mindlin model of practical use can be of the form
k k
with
Dka;3
1 bk; D ¼ Dks a
hDa i ¼ 0;
2PR ¼ RT AR R þ cT Gc þ 2RT F R þ 2cT F c
1 gk;3 ¼ Dks g^k ;
ð67Þ
La v;a ¼ ck =h
hgi ¼ 0;
b k ; g^k can be obtained from Eqs. (61), (62), (65) as well as the where D a boundary conditions in Eq. (58). Now we are ready to obtain an expression for the total energy that is asymptotically correct through the order of l(h/l)2n2, viz., T
2P1 ¼ v Av þ v
T ;a Ba
T ;1 C
T
v;a þ 2v v;2 2v F
To find an equivalent Reissner–Mindlin model Eq. (75) for Eq. (74), one has to eliminate all partial derivatives of the classical 2-D strain. The equilibrium equations in Ref. [19] are used to achieve this purpose. From the two equilibrium equations balancing bending moments with applied moments ma, one can obtain the following formula:
Gc F c ¼ DTa AR;a þ ½ m1
T
m2
ð76Þ
where FR,a is dropped because they are high order terms. Substituting Eq. (76) into Eq. (74), one can show that FR = F and Fc = 0. Finally one can rewrite Eq. (74) as
# hx3 Dk i ; A¼ hx3 Dk i x23 Dk D E b a þ LT hDa;3 i þ hDT iLa ; b T D1 D Ba ¼ Dsaa DT? D? D a s a a;3 D E D E T 1 T T b 2 þ L D2;3 þ D bTD D C ¼ Ds12 D? D? D 1 s 1 1;3 L2 ; D E D E T T b T D1 g^;a F ¼ Dþ? s3 þ D? b3 þ DT? /3 D a s
LTa sk þ bk þ hpk i hDk i
2P1 ¼ 2PR þ U
ð69Þ
Ba ¼ Ba þ ADa G1 DTa A;
ð78Þ
C ¼ C þ AD1 G1 DT2 A
ð79Þ
⁄
4. Transforming into the Reissner–Mindlin model Although Eq. (68) is asymptotically correct through the second order and straightforward use of this strain energy is possible, it involves more complex boundary conditions than necessary since it contains derivatives of the generalized strain measures. To obtain an energy functional that is of practical use, one can transform Eq. (68) into the Reissner–Mindlin model. In the Reissner–Mindlin model, there are two additional degrees of freedom, which are the transverse shear strains incorporated into the rotation of transverse normal. Another triad Bi for the deformed shell is introduced to define the 2-D strains as
ð70Þ ð71Þ 2c23 T .
U ¼ RT;a Ba R;a þ 2RT;1 CR;2 and
Eq. (68) is an energy functional expressed in terms of 2-D variables which can asymptotically approximate the original 3-D energy. It is noted that quadratic terms of the applied loads are neglected as they will not affect the 2-D model.
R;a ¼ Aa Ba þ eab Bb þ 2ca3 B3 Bi;a ¼ Aa K a2 B1 þ K a1 B2 þ K a3 B3 Bi
ð77Þ
where
;a
where the transverse shear strains are c ¼ ½ 2c13
ð75Þ
ð68Þ
where
"
ð74Þ
If we can drive U to be zero for any R, then we have found an asymptotically correct Reissner–Mindlin shell model. For general anisotropic shells, this term will not be zero; but we can minimize the error to obtain a Reissner–Mindlin model that is as close to asymptotical correctness as possible. The accuracy of the Reissner–Mindlin model depends on how close to zero one can drive this term of the energy. In other words, one needs to seek an optimal set of the 27 unknowns (3 unknowns for G and 24 unknowns for La) so that the value of the quadratic form in Eq. (78) is as close to zero as possible for arbitrary generalized strain measures. We let the distinct 78 terms in the symmetric 12 12 coefficient matrix equal to zeros to formulate 78 equations. It is a linear system with 27 unknowns. Then we use a least square technique to solve the overdetermined system for the constants. Mathematically, the overdetermined system (78 equations with 27 unknowns, indicated by MX = b) may be singular for some material properties. For example, the rank of MTM is only 26 for single-layer isotropic and orthotropic plates. In this situation, singular value decomposing technique can be applied to solve this least square problem. After minimizing U⁄, the Reissner–Mindlin model to be used for 2-D shell problem can be expressed as:
2PR ¼ RT AR þ cT Gc þ 2RT F
ð80Þ
973
Z. Yifeng et al. / Composite Structures 94 (2012) 966–976
where A, G, F capture the necessary material and geometric information obtained from the dimensional reduction process. It is worthy to emphasize that although the 2-D constitutive model is constructed in a way dramatically different from traditional Reissner–Mindlin models, the shell analysis remains the same, with no changes in the governing equations and boundary conditions except that the strain measures are now defined using Eq. (70).
T
Dkt yk3;3 þ Dket Ia v kk;a þ eTb Dks v kk;3 þ ea Dk? v;a þ /3 ¼ K3 ;
Dt y3;3 þ
DTet Ia
v k;a
þ
;b
;3
¼ s3 ;
Dt y3;3 þ DTet Ia v k;a ¼ b3 ; h i ½y3 ¼ 0; Dt y3;3 þ DTet Ia v k;a
Xi
¼0 ð85Þ
5. Recovery relations Thus far, a generalized Reissner–Mindlin model based on the asymptotically correct second-order energy for FG cylindrical shells has been constructed. Then, the original problem governed by the 3-D Hamiltons principle in Eq. (37) can be approximated by the following 2-D variational statement:
Z
t2
½dðK 2D PR Þ þ dW 2D dt ¼ 0
ð81Þ
t1
with K2D obtained from Eq. (28), PR from Eq. (80) and dW 2D from Eq. (34). It is noted that the variational statement in Eq. (80) is not asymptotically correct because some optimization technique has been used in the transform the asymptotically correct energy functional in Eq. (68) to the generalized Reissner–Mindlin model in Eq. (75), which is a price one has to pay for the simplicity of the 2-D shell analysis using the Reissner–Mindlin model. The 2-D analysis based on the variational statement of Eq. (76) can be routinely solved using shell elements in commercial finite element packages. This model can be used for various analyses of FGM shells, spanning from static, dynamic, buckling, to aeroelastic analyses. In many applications, however, the capability of predicting accurate 2-D displacement fields of FGM shells is not sufficient. Ultimately, the fidelity of a reduced-order model should be evaluated based on how well it can predict the 3-D displacement/strain/stress fields for the original 3-D problem. Therefore, it is necessary to provide recovery relations to express the 3-D displacement/strain/stress fields in terms of 2-D quantities and x3. Using Eqs. (2), (11) and (12), one can recover the 3-D displacement field through the first order as
U i ¼ ui þ x3 ðC 3i d3i Þ þ C ij wj
ð82Þ
where wa = va, w3 = C\v. From Eq. (20), the 3-D strain field can be recovered up to the first order as
Ce ¼ e þ x3 j;
2Cs ¼ v k;3 þ ea D? v;a ;
Ct ¼ D?;3 v
ð83Þ
Consequently, 3-D stresses rij can be obtained by applying the 3-D constitutive relations. Since we have obtained an optimal estimation of the shear stiffness matrix G, the recovered 3-D results up to the first order are better than CLT and FOSDT. However, the transverse normal stress (r33) is a second-order quantity and cannot be estimated within the first-order approximation. Despite that it is usually much smaller than other stress components, r33 is critical for predicting some structural failure mechanisms such as delamination. In order to obtain a reasonable recovery for the transverse normal stress, VAM procedure is applied once more to find the warping functions of second-order accuracy. Similarly, we perturb the warping functions as
wk ¼ v k þ yk þ oððh=lÞ2 nÞ; w3 ¼ D? v þ y3 þ ððh=lÞ2 nÞ
where K3 is the Lagrange multiplier to enforce the constraint hy3i = 0. The Euler–Lagrange equation and the inter-surface continuity equation in Eq. (85) can be expressed as:
T Dkt yk3;3 þ Dket Ia v kk;a ¼ Ekab vab þ Sk þ K3 x3 þ constky3
ð86Þ
k k k constkþ1 y3 consty3 ¼ Eab vab ðzkþ1 Þ þ S ðzkþ1 Þ
ð87Þ
where
Ekab ðx3 Þ ¼ Sk ðx3 Þ ¼
Z
Z
x3
zk x3
zk
k k eTa Dk b þ Dsab D? dz
k eTa g k ;a þ /3 dz
ð88Þ ð89Þ
The Lagrange multiplier K3 and the solution of y3 are given by
1 s3 þ b3 hEab;3 iv;ab hS;3 i h yk3 ¼ Ek;ab v;ab þ Sk
K3 ¼
ð90Þ ð91Þ
with
1 D E k Ekab;3 ¼ Dkt Ek ab ; Eab ¼ 0;
1 Sk;3 ¼ Dkt Sk ; hSi ¼ 0
ð92Þ
k⁄ where Ek ab and S can be obtained from Eqs. (86), (87), (90) as well as the interlaminar continuous conditions in the last equation of Eq. (85). Although y3 can help us obtain an energy expression asymptotically corrected up to the order of (h/l)4n2, such an energy expression is too complex for practical use. We will still use the Reissner– Mindlin model to carry out the 2-D shell analysis and use y3 for the second-order prediction of the 3-D displacement/strain/stress field. As will be shown latter, this approach achieves excellent predictions even though only the Reissner–Mindlin shell model is used for the 2-D shell analysis. Finally, we can recover the 3-D displacement field up to the second order as
U i ¼ ui þ x3 ðC 3i d3i Þ þ C ji wj þ d3i C 3i y3
ð93Þ
and the strains up to the second order as
Ce ¼ þ x3 j þ Ia v k;a ;
2Cs ¼ v k;3 þ ea D? v;a ;
Ct ¼ D?;3 v þ y3;3
ð94Þ
Finally, we can recover the 3-D stress field up to the second order as
½ r13
r12 r22 T ¼ Dk ð þ x3 jÞ þ Det y3;3 þ De Ia v k;a ; r23 T ¼ Ds ðv k;3 þ ea D? v;a Þ;
r33 ¼
DTet Ia
½ r11
ð95Þ
v k;a þ Dt y3;3
ð84Þ
where yk and y3 are the second-order warping functions. It can be shown that the in-plane components yk vanish and the equations governing y3 are
6. Numerical example A general formulation to treat multilayer FG cylindrical shells with properties as functions of x3 has been derived. In the
974
Z. Yifeng et al. / Composite Structures 94 (2012) 966–976
following, a cylindrical bending example of a homogeneous substrate with a thin SiC–Al functionally graded coating under sinusoidal pressure on the top surface will be used to demonstrate the performance of the present theory. The example is a double-layer coating/substrate system as discussed in Ref. [4] including a homogeneous substrate with a thin functionally graded coating (see Fig. 3). The coordinates are chosen in such a way that x1 2 [0, u], x2 2 [0, 1], and x3 2 [h/2, h/2]. In this system, the Poisson’s ratios for the substrate and the coating ð1Þ ð2Þ ð1Þ ð2Þ are assumed to be constants and equal, mH ¼ mH ¼ mF ¼ mF ¼ 0:3. The shear modulus of the substrate and the coating takes the form of (see Fig. 4): ð1Þ
ð1Þ
GF ¼ GH ;
ð2Þ
ð2Þ
GF ðx3 Þ ¼ nGH ehð1þx3 =hÞ ð2Þ
ð96Þ
ð1Þ
where n ¼ GH =GH is the ratio of the stiffness on the top of the FG coating to that of the homogeneous substrate; h = h/h(2)lnn is the inhomogeneous parameter of the FG coating; h(2) is the coating thickness. In the present analysis, n and h/h(2) are chosen to be 10 and 5, respectively. The top surface of the FGM coating is ceramic rich and the bottom surface is metal rich. The region between the two surfaces is made of the mixture of ceramic–metal materials with continually varying of the volume fractions of the ceramic and metal. The volume fraction of the ceramics Vc is assumed to vary according to a power law as
1 x3 p þ V c ¼ V c þ V þc V c 2 h
Fig. 4. Through-the-thickness variation of shear moduli G in two coating/substrate system.
ð97Þ
where V þ c and V c are the volume fractions of the ceramic on the top and the bottom surfaces, respectively; p P 0 is the volume fraction index, and h the thickness of the shell. The effective elastic moduli of the functionally graded metal– ceramic material are estimated by the Mori–Tanaka scheme as
K Km Vc ¼ ; c K m K c K m 1 þ ð1 V c Þ K Kþ4 l =3 m m
ð98Þ
l lm Vc ¼ lc lm 1 þ ð1 V c Þ llcmþslm1
Fig. 5. Through-the-thickness variation of effective Young’s modulus of FG metal– ceramic material estimated by Mori–Tanaka scheme for different values of p.
where s1 = lm(9Km + 8lm)/(6Km + 12lm), Km and lm represent the bulk and shear modulus of the metal and Kc and lc represent the bulk and shear modulus of the ceramic material [7,20]. The Young’s modulus, E(x3), and Poisson’s ratio, m(x3), are related to effective bulk and shear moduli by
K ðx3 Þ ¼
Eðx3 Þ ; 3½1 2mðx3 Þ
l ðx3 Þ ¼
E ð x3 Þ 2½1 þ mðx3 Þ
ð99Þ
For the purpose of illustration, we choose the properties of the constituent materials as Em = 70 GPa, vm = 0.3 for Al and Em = 427 GPa, vm = 0.17 for SiC. Figs. 5 and 6 plots through-the-thickness var-
Fig. 6. Through-the-thickness variation of Poisson’s ratio of FG metal–ceramic material estimated by Mori–Tanaka scheme for different values of p.
iation of effective Young’s modulus and Poisson’s ratio for various material index p. The shell is assumed to be simply supported and subjected to the following boundary conditions: Fig. 3. Sketch of a two double-layer coating/substrate cylindrical shell.
x1 ¼ 0;
p=3 : r11 ¼ 0; U 2 ¼ U 3 ¼ 0
ð100Þ
975
Z. Yifeng et al. / Composite Structures 94 (2012) 966–976
Fig. 7. Distribution of the longitudinal deflection U 1 vs the thickness coordinate at x1 = 0.
Fig. 8. Distribution of the transverse deflection U 3 vs the thickness coordinate at x1 = p/6.
Fig. 9. Distribution of the longitudinal normal stress r11 vs the thickness coordinate at x1 = p/6.
The top surface of the coating at x3 = h/2
issubjected to a sinusoidally distributed pressure s3 ¼ p20 sin nRp1x/1 ðR1 ¼ 10; n ¼ 2; u ¼ p=3Þ. There is no body force and the bottom surface is traction free. The stress and displacement components are nondimensionlized as r ij ¼ rij =p0 ; U i ¼ U i GHð2Þ =ðp0 hÞ. The through-thickness variations of various displacement and stress components for a thick (R1u/h = 3) double-layer FG coating/substrate system are plotted in Figs. 7–12. It can be observed
Fig. 10. Distribution of the in-plane shear stress x1 = 0.
r12 vs the thickness coordinate at
Fig. 11. Distribution of the transverse shear stress at x1 = 0.
r13 vs the thickness coordinate
Fig. 12. Distribution of the transverse normal stress r33 vs the thickness coordinate at x1 = p/6.
from these figures that all displacement and stress components except for U 3 match well with the 3-D exact solutions. It seems that there exist a constant shift (about 5% of the maximum value) between the results of the present model and the 3-D solution for the transfers displacement U 3 . This may be attributed to that the shell model is reduced from the original 3-D model and some information which cannot be captured by a 2-D model are lost during
976
Z. Yifeng et al. / Composite Structures 94 (2012) 966–976
the dimensional reduction process. Nevertheless, both results show similar variation trend for U 3 . It is important to note that by adopting FG coating the shear moduli G is continuous at the interface, the in-plane stress components such as r11 and r12 , therefore, are smooth throughout of the total thickness. In addition, no sharp peak appears for transverse shear stress r13 . Results presented in Figs. 7–12 clearly show that using FG coating will eliminate the mismatch of interface stress components to reduce the risk of cracking and rebounding of the coating. 7. Conclusions Using the variational asymptotic method (VAM), an efficient high-fidelity shell model for multilayer FG cylindrical shell has been developed. The VAM is used to systematically reduce the original, nonlinear 3-D model to a series of 2-D models in terms of the small parameter. No apriori assumptions have been adopted during the derivation. The theory is applicable to FG shells with material properties either being constant or changing continuously in each layer. A cylindrical bending example of a homogeneous substrate with a thin SiC–Al functionally graded coating under sinusoidal pressure on the top surface has been used to validate the present model. Excellent agreement between the present model and the 3-D exact solution have been observed for the 3-D distribution of displacements and stresses through the shell thickness. Moreover, the present model is valid for large displacements and global rotations and can capture all the geometric nonlinearity of a shell when the strains are small. The present paper has built on the second author’s previous work [15,21] with the following new contributions: 1. The present work has the capability of analyzing multilayer FGM cylindrical shell with material properties as functions of transverse locations or constants while in previous work these properties are treated as constants for each layer. 2. Interface continuity conditions are explicitly derived and solved to obtain the multilayer solutions. 3. Simplifications have been made in deriving B, C, D matrices. 4. Explicit analytical solutions for the second-order approximation of warping functions have been provided.
Acknowledgments This research is supported by the Fundamental Research Funds for the Central Universities, China, under Grant No. 0218005202017 and China Postdoctoral Science Foundation, under Grant No. 2011M501386. The views and conclusions contained herein are those of the writers and should not be interpreted as
necessarily representing the official policies or endorsement, either expressed or implied, of NSFC. References [1] Muller E, Drasar C, Schilz J, Kaysser WA. Functionally graded materials for sensor and energy applications. Mater Sci Eng A – Struct Mater Prop Microstruct Process 2003;362(1–2):17–39. [2] Mishnaevsky LL. Functionally gradient metal matrix composites: numerical analysis of the microstructure-strength relationships. Compos Sci Technol 2006;66(11–12):1873–87. [3] Shaw LL. Thermal residual stresses in plates and coatings composed of multilayered and functionally graded materials. Composites Part B – Eng 1998;29(3):199–210. [4] Kashtalyan M, Menshykova M. Three-dimensional elastic deformation of a functionally graded coating/substrate system. Int J Solids Struct 2007;44(16):5272–88. [5] Das M, Barut A, Madenci E, Ambur DR. A triangular plate element for thermoelastic analysis of sandwich panels with a functionally graded core. Int J Numer Meth Eng 2006;68(9):940–66. [6] Pan E, Heyliger PR. Exact solutions for magneto–electro-elastic laminates in cylindrical bending. Int J Solids Struct 2003;40(24):6859–76. [7] Vel SS, Batra RC. Exact solution for thermoelastic deformations of functionally graded thick rectangular plates. AIAA J 2002;40(7):1421–33. [8] Qatu MS, Sullivan RW, Wang W. Recent research advances on the dynamic analysis of composite shells: 2000–2009. Compos Struct 2010;93(1):14–31. [9] Maalawi KY. A generalized formulation for radially graded composite cylinders/ rings with maximized buckling stability limits. J Eng Appl Sci 2010;23(5):129–48. [10] Khazaeinejad P, Najafizadeh MM, Jenabi J, Isvandzibaei MR. On the buckling of functionally graded cylindrical shells under combined external pressure and axial compression. J Press Ves Technol–Trans ASME 2010;132(6):1551–7. [11] Shariyat M. Dynamic buckling of suddenly loaded imperfect hybrid FGM cylindrical shells with temperature-dependent material properties under thermo–electro-mechanical loads. Int J Mech Sci 2008;50(12):1561–71. [12] Yas MH, Shakeri M, Heshmati M, Mohammadi S. Layer-wise finite element analysis of functionally graded cylindrical shell under dynamic load. J Mech Sci Technol 2011;25(3):597–604. [13] Bian ZG, Chen WQ, Lim CW, Zhang N. Analytical solutions for single- and multi-span functionally graded plates in cylindrical bending. Int J Solids Struct 2005;42(24–25):6433–56. [14] Soldatos KP, Watson P. A method for improving the stress analysis performance of one- and two-dimensional theories for laminated composites. Acta Mech 1997;123(1–4):163–86. [15] Yu WB, Hodges DH, Volovoi VV. Asymptotic generalization of Reissner– Mindlin theory: accurate three-dimensional recovery for composite shells. Comput Meth Appl Mech Eng 2002;191(44):5087–109. [16] Yu WB, Hodges DH, Volovoi VV. Asymptotically accurate 3-D recovery from Reissner-like composite plate finite elements. Comput Struct 2003;81(7):439–54. [17] Hodges DH, Atilgan AR, Danielson DA. A geometrically nonlinear-theory of elastic plates. J Appl Mech–Trans ASME 1993;60(1):109–16. [18] Danielson DA, Hodges DH. Nonlinear beam kinematics by decomposition of the rotation tensor. J Appl Mech–Trans ASME 1987;54(2):258–62. [19] Atilgan AR, Hodges DH. On the strain-energy of laminated composite plates. Int J Solids Struct 1992;29(20):2527–43. [20] Gilhooley DF, Batra RC, Xiao JR, McCarthy MA, Gillespie JW. Analysis of thick functionally graded plates by using higher-order shear and normal deformable plate theory and MLPG method with radial basis functions. Compos Struct 2007;80(4):539–52. [21] Yu WB. Mathematical construction of a Reissner–Mindlin plate theory for composite laminates. Int J Solids Struct 2005;42(26):6680–99.