- Email: [email protected]
Dynamic response of functionally graded material shells with a discrete double directors shell element
Accepted Manuscript Dynamic response of functionally graded material shells with a discrete double directors shell element A. Frikha, M. Wali, A. Hajlaoui, F. Dammak PII: DOI: Reference:
S0263-8223(16)30304-X http://dx.doi.org/10.1016/j.compstruct.2016.07.021 COST 7621
To appear in:
Composite Structures
Received Date: Revised Date: Accepted Date:
15 April 2016 30 June 2016 2 July 2016
Please cite this article as: Frikha, A., Wali, M., Hajlaoui, A., Dammak, F., Dynamic response of functionally graded material shells with a discrete double directors shell element, Composite Structures (2016), doi: http://dx.doi.org/ 10.1016/j.compstruct.2016.07.021
This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
Dynamic response of functionally graded material shells with a discrete double directors shell element. A. Frikhaa,∗, M. Walia , A. Hajlaouia , F. Dammaka a Mechanical
Modeling and Manufacturing Laboratory (LA2MP), National School of Engineers of Sfax, B.P 1173-3038, Sfax, University of Sfax, Tunisia
Abstract The purpose of the paper is to compute the dynamic behavior of functionally graded material (FGM) shell structures subjected to time-varying excitation using 3D-shell model based on a discrete double directors shell element. The third-order shear deformation theory is introduced in the present method to remove the shear correction factor and improve the accuracy of transverse shear stresses. Material properties of the shell are assumed to be graded in the thickness direction by varying the volume fraction of the ceramic and the metallic constituents using general four-parameter power-law distribution. The transient excitation is defined in the time domain and known at each time. The damping material is neglected and the time derivative is approximated by Newmark method. Numerical results for deflection and stresses are presented for plates and spherical caps. The effect of an imposed force on the response of the FGM shell is discussed. The numerical examples prove a good accuracy and reliability compared to the few results available in literature. Keywords: FGM shell structures, Dynamic response, Third-order theory, Double director vector, Newmark method.
1. Introduction The improved properties of composite materials have led to their use in space and aerospace applications. However, the abrupt change of the properties across the interface between different materials in conventional composite materials causes high interlaminar stresses leading to delamination. In addition, large plastic deformation at the interface may activate the initiation and the propagation of cracks in the conventional composite material. Functionally graded materials (FGMs), which are special kind of composite with a gradual transition of material properties from one material to another, are made to overcome the problems associated to the discontinuity in conventional composite. The most known FGMs are made of a transition alloys from metal at one surface to ceramic at the opposite surface. This kind of FGMs was introduced as ultra high temperature resistant materials for nuclear reactor, chemical plants, heat engine components and aerospace vehicles. The thermal resistance of FGM is due to low thermal conductivity of ceramic [1, 2]. The static behavior of FGM plates was studied in literature with analytical and numerical methods. Thin FGM plate can be studied by a classical plate theory [3, 4]. For thick FG plates, some models take into account the transverse shear effect by using the First-order Shear Deformation Theory (FSDT) [5, 6]. In the FSDT theories, transverse shear is assumed to be constant through the shell thickness and thus require the computation of shear correction coefficients, [7–9]. In fact, shear correction coefficients are problem dependent and cumbersome. This limitation of the FSDT forced the development of High-Order Shear Deformation theory (HSDT) which includes the consideration of realistic parabolic variation of transverse ∗ Corresponding
author Email address: [email protected] (A. Frikha)
Preprint submitted to Elsevier
July 21, 2016
shear stress through the shell thickness in the isotropic cases. There is no need of a shear correction factor when using a HSDT but equations of motion are more complicated to obtain than those of the FSDT. The HSDT are introduced to develop a plate or shell elements, [10–14] among authors. An other possibility to introduce a HSDT is by using the enhanced assumed strain formulation (EAS) in solid-shell elements with quadratic transverse shear enhancement, [15, 16]. Recently, a unified formulation developed by Carrera and his co-workers (CUF), which can generate any refined theory, is developed in static and free vibration for laminate composites and FGM shells [17–21]. FGM shell structures are subjected to severe dynamic loads. Therefore, an assessment of natural frequency and transient response of structures seems required. Based on the Mindlin’s first-order shear deformation theory, free vibration of functionally graded plates and cylinder was investigated by Hosseini-Hashemi et al [22]. Using a discrete double directors shell element, free vibration analysis of FGM shell structures was studied by Wali et al. [12]. Based on layer-wise finite element, the dynamic response of functionally graded cylindrical shell was studied by Yas et al. [23]. Using a refined 8-node shell element, a forced vibration of FGM to arbitrary loading is investigated in [24]. The transient response of FGM plate is analyzed by [5]. The objective of this work is to present a general formulation for forced vibration of FGM shell using discrete double directors shell model. The formulation is investigated for the transient response of FGM shells. Square plate and spherical cap are used to compare the accuracy of the formulation. Transient deflection and axial stress are analyzed along the thickness. 2. Double directors shell model In this section, the geometry and kinematic of double directors shell model are described. The reference surface of the shell is assumed to be smooth, continuous and differentiable. Variables associated with the undeformed state will be denoted by upper-case letters and by a lower-case letter when referred to the deformed configuration. 2.1. Double directors shell kinematic assumption
Figure 1: Double directors shell model.
Parameterizations, which define material points of the shell, are carried out in terms of curvilinear coordinates ξ = (ξ 1 , ξ 2 , ξ 3 = z). The position vectors of any material point (q), whose normal projection on mid-surface is the material point (p), in the initial states C0 are given by: 2
X q (ξ 1 , ξ 2 , z) = X p (ξ 1 , ξ 2 ) + zD(ξ 1 , ξ 2 ),
z ∈ [−h/2, h/2] ,
(1)
where h is the thickness. X p and D are respectively a point of the reference surface and the initial shell director. The covariant basis (G1 , G2 , G3 ) is obtained from the position vector by (G1 , G2 , G3 )= (∂X q /∂ξ 1 , ∂X q /∂ξ 2 , ∂X q /∂z), which yields in base vectors relative to the initial state : Gα = Aα + zD ,α ; G3 = D , α = 1, 2.
(2)
The surface element dA in the initial state is given by: dA =
√
AdAξ ,
√
A = kA1 ∧ A2 k , dAξ = dξ 1 dξ 2 .
(3)
The covariant reference metric tensor G at a material point ξ is defined by: G = [Gi · Gj ] , i, j = 1, 2, 3.
(4) p
For use, the geometrical variable Det(G) and dV are related by dV = (G)dξ 1 dξ 2 dz, where Det(G) = p later p (G) = (|Gij |). With the assumption of a double directors shell model, the position vector of any point q, in the deformed configuration is given by (show Fig. 1) : xq (ξ 1 , ξ 2 , z) = xp (ξ 1 , ξ 2 ) + f1 (z)d1 (ξ 1 , ξ 2 ) + f2 (z)d2 (ξ 1 , ξ 2 ).
(5)
In the deformed state, the base vectors are: 0
0
g α = aα + f1 (z)d1,α + f2 (z)d2,α ; g 3 = f1 (z)d1 + f2 (z)d2 .
(6)
The metric tensor components in the deformed configuration Ct are separated into the in-plane and out-ofplane part components. With some approximations, the metric tensor can be written as: g ≈ aαβ + f1 (z)b1αβ + f2 (z)b2αβ αβ 0 0 gα3 ≈ f1 (z)c1α + f2 (z)c2α , (7) gij = g i · g j , 2 0 0 g33 ≈ f1 + f2 d
where aαβ , bkαβ and ckα (k = 1, 2) represent the covariant metric surface, the first curvature tensors and the shear, respectively. The parameter d denotes the thickness stretching. Taking into account d1 · d1 ≈ d2 · d2 ≈ d1 · d2 , these components can be computed as : aαβ = aα · aβ , ckα = aα · dk , bkαβ = aα · dk,β + aβ · dk,α , d = d1 · d1 , k = 1, 2.
(8)
Similar expressions for the in-plane and out-of-plane components of the metric tensor can be obtained in the case of the initial configuration C0 . Using the kinematic assumption, Eq. (7), the linearized strains can be written as follows [11, 12]: αβ = eαβ + f1 (z)χ1αβ + f2 (z)χ2αβ , α, β = 1, 2 . (9) 0 0 2α3 = f1 (z)γα1 + f2 (z)γα2 In matrix notation, the vectors e11 e22 e= 2e12
of membrane, bending and shear strains are given by: k χk11 γ1 χk22 , k = 1, 2. , γk = , χk = γ2k k 2χ12
The variation of the strain can be written, in the initial configuration, as: 3
(10)
δeαβ = 1/2(Aα · δx,β + Aβ · δx,α ), δγαk = Aα · δdk + δx,α · dk , δχkαβ = 1/2(Aα · δdk,β + Aβ · δdk,α + δx,α · dk,β + δx,β · dk,α )
k = 1, 2.
(11)
Or in matrix notation:
δe = B m · δx,
δχk = B bmk δx + B bbk δdk ,
δγ k = B smk δx + B sbk δdk ,
k = 1, 2,
(12)
where the matrix differential operators, relative to the initial configuration (Linear theory), are given by:
Bm
T d0k,1 ∂ξ∂ 1 AT1 ∂ξ∂ 1 T , AT2 ∂ξ∂ 2 = B bbk = d0k,2 ∂ξ∂ 2 , B bmk = T T T ∂ T ∂ A1 ∂ξ2 + A2 ∂ξ1 d0k,1 ∂ξ∂ 2 + d0k,2 ∂ξ∂ 1 " # T T d0k ∂ξ∂ 1 A1 B smk = , d0k = D, k = 1, 2. = , B T sbk AT2 d0k ∂ξ∂ 2
(13)
2.2. Third-order double director shell model The two Functions f1 (z) and f2 (z), used in Eq. (5), are related by: f1 (z) = z − f2 (z), f2 (z) =
4z 3 . 3h2
(14)
This gives the following transverse shear strain: 0
0
2α3 = f1 (z)γα1 + f2 (z)γα2 ,
0
0
f1 (z) = 1 − f2 (z).
(15)
The zero condition of the transverse shear stress on the top and bottom face of the shell, is simply obtained by imposing: 0 0 γα2 = 0, 2α3 = 1 − f2 (z) γα1 , f2 (±h/2) = 1. (16) The kinematic constrain γα2 = 0 will be imposed in a discrete form in the finite element approximation. 3. Weak form The numerical solution with the finite element method is based on the weak form of equilibrium equations. The three dimensional form of the latter is given as : % Z ! 2 X 1 k M k · δχ dA − Gext (Φ, δΦ) = 0, N · δe + T 1 · δγ + (17) G(Φ, δΦ) = A
k=1
where δΦ = (δu, δd1 , δd2 ) is an arbitrary variation. Gext is the external virtual work. δe, δχk and δγ 1 are the variations of shell strains. N , M k and T 1 are the membrane, bending and shear stresses resultants, which can be written in matrix form as : 11 1 N Mk11 T1 N 22 Mk22 N= , Mk = , T1 = , k = 1, 2. (18) T12 12 12 N Mk
Their components are defined as follows :
N αβ =
Z
h/2
σ αβ −h/2
p
G/A dz,
Mkαβ =
Z
h/2
fk (z)σ αβ
−h/2
4
p G/A dz,
T1α =
Z
h/2
−h/2
0
f1 (z)σ α3
p G/A dz, (19)
where σ αβ are the components of the stress tensor. The generalized resultant of stress R and strain vectors Σ are defined by: e N 1 χ M1 . (20) , Σ= R= χ2 M2 1 γ T 1 11×1 11×1 The weak form of the equilibrium equation can then be rewritten as: Z δΣT · RdA − Gext (Φ, δΦ) = 0. G(Φ, δΦ) =
(21)
A
Using strain-displacement relations, Eq. (12), the strain variation can then be given by: Bm 0 0 B bm1 B bb1 0 . δΣ = B δΦ, B = B bm2 0 B bb2 B sm1 B sb1 0
(22)
4. Constitutive equations In the case of an elastic isotropic constitutive model, the stress resultant R is expressed as: R = H T Σ,
(23)
where H T is the linear elastic modulus. Its expression is: H 11 H 12 H 13 H 12 H 22 H 23 HT = H 13 H 23 H 33 0 0 0
0 0 , 0 H 44
where:
(H 11 , H 12 , H 13 , H 22 , H 23 , H 33 ) =
Z
h/2
−h/2
1, f1 , f2 , f12 , f1 f2 , f22 Hdz.
(24)
(25)
f1 and f2 are given by Eq. (14). In Eq. (25), H 44 is given by: H 44 =
Z
h/2 −h/2
0 2 f1 H τ dz,
(26)
where H and H τ , in Eqs. (24) and (25), are in plane and out-of-plane linear elastic sub-matrices, which can be expressed in a Cartesian system as: 1 ν(z) 0 E(z) E(z) 1 0 , Hτ = ν(z) 1 0 , (27) H= 1 − ν 2 (z) 2(1 + ν(z)) 0 1 0 0 (1 − ν(z))/2
where E(z) and ν(z) are the Young’s modulus and the Poisson’s ratio respectively. As, a functionally graded material is considered, these material properties are assumed to vary through the shell thickness according to the rule of mixture written as: P (z) = (Pc − Pm )Vc + Pm , 5
(28)
where P (z), Pm and Pc denote respectively the effective material property, the properties of the metal and ceramic. Vc is the volume fraction, which varies according to the two general four-parameter power-laws distribution [25]. c p 1 z 1 z + + , (29) F GMI (a, b, c, p) : Vc (z) = 1 − a +b 2 h 2 h c p 1 z 1 z , (30) − − +b F GMII (a, b, c, p) : Vc (z) = 1 − a 2 h 2 h where p is the power-law index. The parameters a, b and c determine the material variation profile through the functionally graded shell thickness. Instead of the rule of mixture, Eqs. (28-30), another analytical method for estimating the effective properties can be used as Mori-Tanaka or self-consistent methods [26]. 5. Finite element formulation In this section, the numerical implementation of the presented shell theoretical formulation based upon a four node shell element is established. Using the isoparametric concept, the mid-surface (X) and surface displacement field (u = x − X) are approximated by: X=
4 X
N I XI ,
u=
I=1
4 X
N I uI ,
(31)
I=1
where X I ∈ R3 and uI ∈ R3 are position and displacement at node I. N I are the standard bilinear isoparametric shape functions. The first director vector d1 is approximated with the same functions d1 =
4 X
N I d01I ,
(32)
I=1
where d01I denote the first directors at the nodal points in the initial configuration C0 . 5.1. Local Cartesian system Let n0 = A1 ∧ A2 / kA1 ∧ A2 k be the normal field to the mid-surface in the initial state C0 . A local Cartesian system with base vectors n01 , n02 , n0 , can be defined by means of an orthogonal transformation. The Jacobian transformation J from the basis n01 , n02 to {A1 , A2 } is defined by: 0 n1 · A1 n02 · A1 J= . (33) n01 · A2 n02 · A2 Since the formulation is developed in local Cartesian coordinate, the derivatives of the shape functions need to be transformed as follows: ( I ) I N,1 N ,1 −1 . (34) = [J ] I I N,2 N ,2 5.2. Membrane and first bending strain field The shell membrane part of the problem is considered. The strain-displacement relation is: δe = B m Φn , where B m is the discrete membrane strain-displacement operator defined as:
6
(35)
B Im =
B Imm
0 0
,
B Imm
T I n01 N ,1 T I = . n02 N ,2 T T I I 0 0 n1 N ,2 + n2 N ,1
(36)
For the first bending part, the strain-displacement relation is given by: δχ1 = B 1 · δΦn ,
(37)
where B 1 is the discrete first bending strain-displacement operator: T T I I d01,1 N ,1 n01 N ,1 T T I I I I = B I1 = B I1m B I1b 0 , B1m , B1b = d01,2 N ,2 n02 N ,2 , T T T T I I I I d01,1 N ,2 + d01,2 N ,1 n01 N ,2 + n02 N ,1 where d01,a =
P4
(38)
I
I=1
N ,a d01I , a = 1, 2.
5.3. Construction of the assumed natural transverse shear strain field A typical isoparametric finite element is considered as depicted in Fig. 2. A, B, C and D denote the mid-points of the element boundaries set. Following [27], the assumed natural transverse shear strain field is expressed by: δγ11 (1 − η)δγ11 (B) + (1 + η)δγ11 (D) 1 , (39) δγξ = = (1 − ξ)δγ21 (A) + (1 + ξ)δγ21 (C) δγ21 where γ21 (A), γ11 (B), γ21 (C) et γ11 (D) are the strains at points defined in Fig. 2.
Figure 2: Assumed strain construction of the isoparametric shell element.
δγ 1ξ is related to δΦn by: δγ 1ξ = B sξ δΦn ,
(40)
where B sξ is the discrete shear strain-displacement operator expressed as: Bsξ =
"
T
T
T
T
N1,1 d01B
N2,1 AT 1B
0
N2,1 d01B
N2,1 AT 1B
0
N3,1 d01D
N3,1 AT 1D
0
N4,1 d01D
N3,1 AT 1D
0
N1,2 d01A
N4,2 AT 2A
0
N2,2 d01C
N3,2 AT 2C
0
N3,2 d01C
N3,2 AT 2C
0
N4,2 d01A
N4,2 AT 2A
0
T
T
T
T
#
. (41)
In the local Cartesian system, the transverse shear strains are: δγ 1 = B s · δΦn ,
B s = J −1 B sξ .
7
(42)
5.4. Second bending strain field: discrete constraints The shear part relative to the second director vector d2 will be vanished in a discrete form. A quadratic interpolation is selected as proposed in the works of [11, 28]: d2 =
4 X
N I d02I +
8 X
PK αK t0K ,
δd2 =
K=5
I=1
4 X
N I δd02I +
I=1
8 X
PK δαK t0K ,
(43)
K=5
where (I) and (K) represent respectively a node of the element and the mid-point of the element boundaries. δαK are variables associated to δd2 on the element boundaries. The vector t0K is unit and its direction is defined by the position of the nodes couple (I, J) : t0K = (xJ − xI )/LK ,
LK = kxJ − xI k ,
(44)
where LK is the I − J side length. The shape functions PK are quadratic [11]. By introducing the vanishing shearing hypothesis, on top and bottom faces, over the element boundaries under integral form, the transverse strain is given for side (I, J) by: Z
which leads to: δd2 =
4 X
N I δd2I +
I=1
J I
2 δγsz ds = 0,
(45)
8 X 3 1 1 t0K . − PK 0 0 2 2 (δd + δd ) · t L (δu + δu ) · d I J K I J K K K=5
(46)
In a matrix form, δd2 becomes: δd2 =
4 X I=1
M Id δuI + M Ir δd2I ,
(47)
where matrixes M Id and M Ir are given by the following expressions: M Id = PK tdIK + PM tdIM ,
tdIK =
3 0 t ⊗ d0K , 2LK K
(48)
3 0 t ⊗ t0K . (49) 4 K The (K) and (M ) are the two mid-side of every side of the quadrilateral, which are bound to the node (I). Finally, the second bending strain is expressed as: M Ir = N I I + PK ttIK + PM ttIM ,
ttIK =
δχ2 = B 2 · δΦn ,
(50)
where B 2 is the discrete second bending strain-displacement operator: B I2 = B I2m and B I2b are given by:
B I2m
B I2m
0 B I2b
T T d02,1 N ,I1 +n01 M Id,1 T T = , d02,2 N ,I2 +n02 M Id,2 T T T T 0 0 I I I 0 I 0 d2,1 N ,2 +n1 M d,2 + d2,2 N ,1 +n2 M d,1
.
(51)
B I2b =
T n01
T n01 · M Ir,1 T n02 · M Ir,2 . T I I 0 · M r,2 + n2 · M r,1
Finally, the continuum generalized strain δΣ, given by Eq. (17) will be in an approximate form as: 8
(52)
δe δχ1 = B · δΦn , δΣ = 2 δχ1 δγ
Bm B1 B= B2 . Bs
Using Eq.(23), the internal virtual work becomes: Z Gint = δΣT · RdA = δΦTn KΦn , A
where K is the stiffness matrix given by K = 5.5. Nodal transformation
R
A
(53)
(54)
B T H T BdA.
In all equations, δdk and δdk,α are the variation of the directors and there derivatives. These variations can be written either in spatial description: δdk = δθ k ∧ dk = Λk δθ k ,
ek , Λk = −d
k = 1, 2,
(55)
ek is the skew-symmetric tensor such that d ek dk = 0 , or in material description: where d
e 3, δdk = Qk δθ k E 3 = Λk δθ k , Λk = Qk E (56) t where we assumed that dk = Qk E 3 and E 3 = 0 0 1 . A spatial description leads to a shell problem with 9 DOF/node and the material description leads to a shell problem with 7 DOF/node, where the transformation Λ take the following form: Λk =
−t2k
t1k
3×2
.
(57)
All numerical examples given in section 6, are obtained with the material description and 7DOF/node. 5.6. Mass matrix Several different schemes for deriving element mass matrices exist in the literature. One finds the consistent mass matrix and different lumped mass matrix constructions. The variation of the kinetic energy is written as: Z (58) Gext = δxT q ρx¨q dv, v
¨ q denote virtual displacement where ρ being the material density and vary along the thickness, δxq and x and acceleration of any point in the shell. By using the displacement field, Eq. (5), the variation of the kinetic energy becomes: Z Z h/2 1 f1 f2 ¨ Gext = δΦT ρΦdA, ρ f1 f12 f1 f2 dz, (59) ρ= A −h/2 f2 f1 f2 f22
the finite element discretization of Gext can then be written in the following form: ¨ Gext = δΦT n M Φn .
(60)
M is the mass matrix. This matrix can be either the diagonal mass matrix or the consistent mass matrix. To construct a diagonal mass matrix, several alternatives were investigated by [29], in a free vibration analysis of Mindlins plate. In the present numerical implementation, the efficient scheme considered in [12] is used where the diagonal mass matrix is given, with 7DOF/node, as follows: 9
M ii
with :
M1i
=
M1i
0
M1i M2i
0
M2i M3i M3i
,
(61)
R ρ N Ni dA A ρkk dA , k = 1, 2, 3. ( ρkk Ni Ni ) dA
(62)
The governing equation for direct transient analysis can then be written as: n o ¨ + [K] {Φ} = {F} . [M] Φ
(63)
M ki =
R
A R kkP i
The above equation becomes transient when the mechanical force, {F} varies with time. It is solved by using Newmark’s method. 5.7. Newmark algorithm The Newmark’s direct integration scheme is used to approximate the time derivatives and to solve the governing dynamic equation. This n o n oimplicit method allows to compute the solution at time t + ∆t from ¨t . The Newmark’s formula to process the solution for damped system known vectors {Ut }, U˙ t and U is given by: K {Ut+∆t } = {Rt+∆t } . (64) n o n o ¨0 is computed by For a given initial displacement {U0 } and velocity U˙ 0 , the initial acceleration U solving Eq. (63). K is update by using: K = [K] +
1 γ [M] + [C] . β∆t2 β∆t
(65)
In our study, β and γ are equal respectively to 0.25 and 0.5. For each time step t + ∆t, it is required to proceed as follows: • to compute the residu {Rt+∆t } by using : n o 1 1 n˙ o 1 ¨t {Rt+∆t } = {F t+∆t } + [M] {U } + − 1 U U + + t t β∆t2 β∆t 2β n o n o γ γ γ ¨t U˙ t + ∆t U , {U t } + −1 −1 [C] β∆t β 2β
• to update K ,
10
(66)
• to solve Eq. (64) for {Ut+∆t }, n o n o ¨ t+∆t and the velocity U˙ t+∆t at time t + ∆t by using : • to compute the acceleration U n o n o 1 n o 1 ¨t , ˙ ¨ t+∆t = 1 U − β U {U } − {U } − ∆t − U t t+∆t t β ∆t2 2
(67)
n o n o n o n o ¨ t + γ∆t U ¨ t+∆t . U˙ t+∆t = U˙ t + ∆t (1 − γ) U
(68)
and
6. Numerical results To show the dynamic response of FGM shells with a discrete double directors shell element, the authors offer numerical results for plate and spherical caps. The model was validated in each structure for FGM or isotropic material with literature. The effect of functionally graded material is then investigated for each structure. In this section, the dynamic response of FGM shell, using the developed SHO4 element, is presented. The abbreviation SHO4 is used to identify the present linear discrete double directors shell element. 6.1. Square plate The forced vibration was performed on a square plate of side length l = 0.2 m and thickness h = 0.01 m, simply supported on all its edges. The plate is subjected to suddenly applied uniform pressure loading equal to q0 = −106 N/m2 [30]. The non-dimensional parameters are center deflection, w = wEm h/(q0 l2 ), stress, p 2 2 2 σ = σh /(|q0 | l ), and time, t = t Em /(l ρm ), where Em is the Young modulus of metal. The analysis was conducted using Aluminum and Zirconia (ceramic). The material properties of Aluminum and Zirconia are given in Tab.1. Table 1: Material properties of FGM plate [30].
Parameter Aluminum Zirconia
Young’s modulus (GPa) 70 151
Poisson’s ratio 0.3 0.3
Density (kg/m3 ) 2707 3000
Only one-half of the square plate was modeled. The mesh consists of 128 nodes and 192 elements. One layer of SHO4 elements was used in the thickness direction. Fig.3 shows the dynamic response of simply supported FGM square plate under suddenly applied uniform pressure loading. The non-dimensional center deflection of the plate is considered here, and the present results for different values of the power law exponent n are compared with the results from literature [5, 30]. The discrepancy between results obtained from two methods is closely small. It should be noted that a non linear model is used in [5, 30]. The non-linearity effect is more significant for lower stiffness plate corresponding to metal. As known in literature, non-linear analysis gives a reduced deflection than linear one. The gap between the present and the reference results for metal can be explained by the taking into account the non-linearity effect in [5, 30]. The magnitude of vibration is the maximum for metallic plates and a minimum for the ceramic plate. The magnitude of vibration increases smoothly with the grading index p. The bending rigidity of ceramic, corresponding to the lower magnitude, is higher than the bending rigidity of metal, corresponding to the highest magnitude. It is seen that the frequency of vibration of the ceramic plates is much higher than that of metallic plates since the maximum magnitude of ceramic plate is located before the maximum magnitude of the metallic plate along the time axis. 11
Figure 3: Temporal variation of center deflection of simply supported FGM square plate under suddenly applied uniform pressure of intensity q0 = −106 N/m2 , compared with [30].
The computational domain for the dynamic analysis with clamped boundary conditions is the same as that of Fig. 3. The analysis was performed for the same two material combinations. Fig. 4 shows the temporal response of the plates under suddenly applied mechanical loading of −106 N/m2 . The deflection of clamped FGM plates is much lower than the simply supported case. The effective stiffness of clamped plate is higher than the simply supported plates.
Figure 4: Non-dimensional center deflection of the of clamped FGM square plate under suddenly applied uniform pressure of intensity q0 = −106 N/m2 .
Fig. 5 shows the temporal variation of axial stress in the center of the FGM plates under suddenly applied mechanical loading of −106 N/m2 by varying the grading index p. Left column of Fig. 5 (a, c, e) depicts the axial stress respectively in z = −h/2, z = 0 and z = h/2 of clamped FGM plates. Right column of Fig. 5 (b, d, f) depicts the axial stress respectively in z = −h/2, z = 0 and z = h/2 of simply supported FGM plates. Usually, the axial stress modulus of clamped plates are lower than simply supported. Since the applied pressure is negative, the axial stress on the bottom, h = −h/2, (resp. top, h = h/2) is mostly positive (resp. negative). At z = −h/2, the axial stress modulus decreases between the grading index p = 0, corresponding to ceramic material, and near to p = 0.2 and then increases to metallic material. The maximum axial stress modulus of metallic and ceramic plates are similar for both clamped and simply supported plates (see Fig. 5 (a-b)). At z = 0, it is shown that the axial stress is zero for isotropic material (metallic and ceramic) for clamped (Fig 5 (c)) and simply supported (Fig 5 (d)) plates, according to the plate theory. The axial stress 12
modulus increases between p = 0, corresponding to ceramic material, and the grading index near p = 2 then decreases to metallic material. At z = h/2, the axial modulus increases between the grading index p = 0 and near to p = 2 and then decreases to metallic material.
Figure 5: Non-dimensional axial stress of the center of FGM square plate under suddenly applied uniform pressure of intensity q0 = −106 N/m2 for (a, c, e): simply supported FGM plate, (b, d, f): clamped FGM plate, (a, b): z = −h/2, (c, d): z = 0, (e, f): z = h/2.
Fig. 6 shows the plots of the axial stress through the thickness at non-dimensional time t = 5 of simply supported and clamped plate under uniform pressure applied on the top surface. The axial stress is compressive for z = h/2 and tensile for z = −h/2 for both simply supported and clamped plates. The axial stress varies linearly versus the thickness for isotropic material. For the different grading index chosen, the plate corresponding to p = 2 yielding to the maximum compressive stress at the top surface. At the bottom surface, the plate corresponding to p = 0.2 yielding to the minimum tensile stress. These phenomena, observed at p = 0.2 and p = 2, are in good agreement with [5]. 6.2. Spherical cap Dynamic response of clamped spherical cap with isotropic material subjected to concentrated apex load is presented. The characteristics of the shell and the concentrated load are shown in Fig. 7. Its parameters are given in Tab. 2: 13
Figure 6: Non-dimensional axial stress in the center of FGM square plate under suddenly applied uniform pressure of q0 = −106 N/m2 for (a) simply supported FGM plate, (b) clamped FGM plate.
Table 2: Parameters of spherical cap [31].
Parameter Value
R 4.76
θ 10.9◦
E 107
ν 0.3
h 0.01576
H 0.0859
ρ 0.000245
F 100
Only one-quarter of the workpiece was modeled. The mesh consists of 217 nodes and 192 elements. One layer of SHO4 elements was used in the thickness direction with Simpson integration (see Fig. 8).
Figure 7: Characteristics of the spherical cap.
A concentrated non-dimensional load of F = 100 is applied between 0 and 250 µs. The response of the cap, represented by the normalized vertical displacement at shell apex, is compared to linear dynamic analysis obtained by [31] for homogeneous material. The present results are found to be in good agreement with those reported by [31]. Since the angle of spherical cap is small with a great radius, the geometry is close to plate. In order to show the effect of the radius on dynamic behavior of spherical shell, an hemisphere (θ = 90◦ ) with the same parameters than Tab. 2 is considered. The spherical cap is subjected to suddenly applied uniform pressure loading equal to q0 = −106 N/m2 . The non-dimensional parameters are center deflection, w = p wEm h/(q0 R2 ), stress, σ = σh2 /(|q0 | R2 ), and time, t = t Em /(R2 ρm ), where Em is the Young’s modulus of metal. The analysis was conducted using Aluminum and Zirconia (ceramic). The Young’s modulus, Poisson’s ratio and the density for Aluminum and Zirconia are given in Tab. 1. Fig. 9 shows the transient non-dimensional center deflection of clamped FGM hemisphere under suddenly applied uniform pressure loading. The magnitude of vibration is the maximum for metallic hemisphere and a minimum for the ceramic case. The transient magnitude of FGM spherical cap is intermediate to that of metal and ceramic. It increases smoothly with the grading index p. The bending rigidity of ceramic, corresponding to the lower magnitude, is higher than the bending rigidity of metal, corresponding to the 14
Figure 8: Transient center deflection of isotropic spherical cap.
higher magnitude. It is seen that the frequency of vibration of the ceramic hemisphere is higher than that of metallic case since the maximum magnitude of ceramic plate is located before the maximum magnitude of the metallic plate along the time axis.
Figure 9: Non-dimensional center deflection of clamped FGM hemisphere under suddenly applied uniform pressure q0 = −106 N/m2 .
Fig. 10 shows the temporal variation of axial stress in the center of the FGM plates under suddenly applied mechanical loading of 106 N/m2 by varying the grading index p. Fig. 10 (a-c) depicts the axial stress respectively in z = −h/2, z = 0 and z = h/2 of clamped FGM hemisphere. Since the applied pressure on the hemisphere is positive, the membrane effect yields positive the axial stress at z = −h/2, z = 0 and z = h/2 (see Fig. 10). At z = −h/2, the axial stress modulus decreases between the grading index p = 0, corresponding to ceramic material, and near to p = 0.2 and then increases to metallic material. At z = h/2, the axial modulus increases between the grading index p = 0 and near to p = 2 and then decreases to metallic material. The same phenomenon is observed for plates and verified with [5]. Fig. 11 shows the non-dimensional axial stress along the thickness at non-dimensional time t = 1 for a clamped FGM hemisphere subjected to a uniform pressure 106 N/m2 . It is shown that the ceramic spherical cap yields to a higher membrane stiffness than metallic one. The variation of axial stress along the thickness is small for isotropic material. This variation is great for FGM spherical cap.
15
Figure 10: Non-dimensional axial stress of the center of clamped FGM hemisphere under suddenly applied uniform pressure of intensity q0 = 106 N/m2 for (a) z = −h/2, (b) z = 0 and (c) z = h/2.
16
Figure 11: Non-dimensional axial stress in the center of clamped FGM hemisphere under suddenly applied uniform pressure q0 = −106 N/m2 .
Conclusions Finite element model based on a discrete double-directors shell element is presented to study the dynamic response of 3D-FGM-shell structures subjected to mechanical loads. The grading function is assumed along the thickness direction. Dimensionless deflection and stress of functionally graded plates and spherical caps are computed by the present model. To validate the results of the present study and demonstrate its accuracy, the results are compared with the existing ones in the literature for plates and spherical caps. A very good agreement among the results confirms the high accuracy of the current approach. Transient response of deflection and axial stress are analyzed by varying the grading index. [1] D. P. H. Hasselman and G. E. Youngblood. Enhanced thermal stress resistance of structural ceramics with thermal conductivity gradient. Journal of The American Ceramic Society, 61 (1,2):49–53, 1978. [2] M. Niino and S. Maeda. Recent development status of functionally gradient materials. International Iron and Steel Institute, 30:699–703, 1990. [3] S.-H. Chi and Y.-L. Chung. Mechanical behavior of functionally graded material plates under transverse load. part 1: Analysis. International Journal of Solids and Structures, 43:3657–3674, 2006. [4] S.-H. Chi and Y.-L. Chung. Mechanical behavior of functionally graded material plates under transverse load. part 2: Numerical results. International Journal of Solids and Structures, 43:3675–3691, 2006. [5] G. N. Praveen and J. N. Reddy. Nonlinear transient thermoelastic analysis of functionally graded ceramic-metal plates. International Journal of Solids and Structures, 35:4457–4476, 1997. [6] J.N. Reddy, C.M. Wang, and S. Kitipornchai. Axisymmetric bending of functionally graded circular and annular plates. European Journal of Mechanics - A/Solids, 18:185–199, 1999. [7] T. K. Nguyen and G. Bonnet. First-order shear deformation plate models for functionally graded materials. Composite Structures, 83:25–36, 2008. [8] M. Shariyat and M. M. Alipour. A novel shear correction factor for stress and modal analyses of annular fgm plates with non-uniform inclined tractions and nonuniform elastic foundations. International Journal of Mechanical Sciences, 87:60–71, 2014. [9] A. Hajlaoui, A. Jarraya, K. El Bikri, and F. Dammak. Buckling analysis of functionally graded materials structures with enhanced solid-shell elements and transverse shear correction. Composite Structures, 132:87–97, 2015. [10] F. Tornabene. Free vibration analysis of functionally graded conical, cylindrical shell and annular plate structures with a four-parameter power-law distribution. Computer Methods in Applied Mechanics and Engineering, 198:2911–2935, 2009. [11] M. Wali, A. Hajlaoui, and F. Dammak. Discrete double directors shell element for the functionally graded material shell structures analysis. Computer Methods in Applied Mechanics and Engineering, 278:388–403, 2014. [12] M. Wali, T. Hentati, A. Jarraya, and F. Dammak. Free vibration analysis of fgm shell structures with a discrete double directors shell element. Composite Structures, 125:295–303, 2015. [13] Mallikarjuna and T. Kant. Effect of cross-sectional warping of anisotropic sandwich laminates due to dynamic loads using a refined theory and c0 finite elements. International Journal for Numerical Methods in Engineering, 35:2031–2047, 1992. [14] T. Kant and J. R. Kommineni. Geometrically non-linear transient analysis of laminated composite and sandwich shells with a refined theory and c0 finite elements. Computers & Structures, 52:1243–1259, 1994. [15] N. D. Quy and A. Matzenmiller. A solid-shell element with enhanced assumed strains for higher order shear deformations in laminates. Technishe Mechanik, 28:334–355, 2008. [16] A. Hajlaoui, M. Wali, M. Ben Jdidia, and F. Dammak. An improved enhanced solid shell element for static and buckling analysis of shell structures. Mechanics & Industry, 17(5):510, 2016.
17
[17] A.M.A. Neves, A.J.M. Ferreira, E. Carrera, M. Cinefra, C.M.C. Roque, R.M.N. Jorge, and C.M.M. Soares. Free vibration analysis of functionally graded shells by a higher-order shear deformation theory and radial basis functions collocation, accounting for through-the-thickness deformation. European Journal of Mechanics - A/Solids, 37:24–34, 2013. [18] M. Cinefra, S. Belouettar, M. Soave, and E. Carrera. Variable kinematic models applied to free vibration analysis of functionally graded materials shells. European Journal of Mechanics - A/Solids, 29:1078–1087, 2010. [19] M. Cinefra, E. Carrera, L. Della Croce, and C. Chinosi. Refined shell elements for the analysis of functionally graded structures. Composite Structures, 94:415–422, 2012. [20] M. Cinefra and S. Valvano. A variable kinematic doubly-curved mitc9 shell element for the analysis of laminated composites. Mechanics of advanced materials and structures, 23(11):1312–1325, 2016. [21] M. Cinefra and E. Carrera. Shell finite elements with different through-the-thickness kinematics for the linear analysis of cylindrical multilayered structures. International Journal for Numerical Methods in Engineering, 93:160–182, 2013. [22] Sh. Hosseini-Hashemi, J. N. Fadaee, and M. Es’haghi. A novel approach for in-plane/out plane frequency analysis of functionally graded circular/annular plates. International Journal of Mechanical Sciences, 52:1025–1035, 2010. [23] M. H. Yas, M. Shakeri, M. Heshmati, and S. Mohammadi. Layer-wise finite element analysis of functionally graded cylindrical shell under dynamic load. Journal of Mechanical Science and technology, 25(3):597–604, 2011. [24] W. Y. Jung and S. C. Han. Transient analysis of fgm and laminated composite structures using a refined 8-node ans shell element. Composites: Part B, 56:372–383, 2014. [25] Z. Su, G. Jin, S. Shi, T. Ye, and X. Jia. A unified solution for vibration analysis of functionally graded cylindrical, conical shells and annular plates with general boundary conditions. International Journal of Mechanical Sciences, 80:62–80, 2014. [26] S. S. Vel and R. C. Batra. Three-dimensional exact solution for the vibration of functionally graded rectangular plates. Journal of Sound and Vibration, 272:703–730, 2004. [27] K. J. Bathe and E. Dvorkin. A four-node plate bending element based on mindlin/reissner plate theory and a mixed interpolation. International Journal for Numerical Methods in Engineering, 21:367–383, 1985. [28] F. Dammak, S. Abid, A. Gakwaya, and G. Dhatt. A formulation of the non linear discrete kirchhoff quadrilateral shell element with finite rotations and enhanced strains. Revue Europeenne des Elements Finis, 14:7–31, 2005. [29] T. A. Rock and E. A. Hinton. A finite element method for the free vibration of plates allowing for transverse shear deformation. Composite Structures, 6:37–44, 1976. [30] J. N. Reddy. Analysis of functionally graded plates. International Journal for Numerical Methods in Engineering, 47:663–684, 2000. [31] L. A. Duarte Filho and M. A. Armando. Geometrically nonlinear static and dynamic analysis of shells and plates using the eight-node hexahedral element with one-point quadrature. Finite Elements in Analysis and Design, 40:1297–1315, 2004.
18
S0263-8223(16)30304-X http://dx.doi.org/10.1016/j.compstruct.2016.07.021 COST 7621
To appear in:
Composite Structures
Received Date: Revised Date: Accepted Date:
15 April 2016 30 June 2016 2 July 2016
Please cite this article as: Frikha, A., Wali, M., Hajlaoui, A., Dammak, F., Dynamic response of functionally graded material shells with a discrete double directors shell element, Composite Structures (2016), doi: http://dx.doi.org/ 10.1016/j.compstruct.2016.07.021
This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
Dynamic response of functionally graded material shells with a discrete double directors shell element. A. Frikhaa,∗, M. Walia , A. Hajlaouia , F. Dammaka a Mechanical
Modeling and Manufacturing Laboratory (LA2MP), National School of Engineers of Sfax, B.P 1173-3038, Sfax, University of Sfax, Tunisia
Abstract The purpose of the paper is to compute the dynamic behavior of functionally graded material (FGM) shell structures subjected to time-varying excitation using 3D-shell model based on a discrete double directors shell element. The third-order shear deformation theory is introduced in the present method to remove the shear correction factor and improve the accuracy of transverse shear stresses. Material properties of the shell are assumed to be graded in the thickness direction by varying the volume fraction of the ceramic and the metallic constituents using general four-parameter power-law distribution. The transient excitation is defined in the time domain and known at each time. The damping material is neglected and the time derivative is approximated by Newmark method. Numerical results for deflection and stresses are presented for plates and spherical caps. The effect of an imposed force on the response of the FGM shell is discussed. The numerical examples prove a good accuracy and reliability compared to the few results available in literature. Keywords: FGM shell structures, Dynamic response, Third-order theory, Double director vector, Newmark method.
1. Introduction The improved properties of composite materials have led to their use in space and aerospace applications. However, the abrupt change of the properties across the interface between different materials in conventional composite materials causes high interlaminar stresses leading to delamination. In addition, large plastic deformation at the interface may activate the initiation and the propagation of cracks in the conventional composite material. Functionally graded materials (FGMs), which are special kind of composite with a gradual transition of material properties from one material to another, are made to overcome the problems associated to the discontinuity in conventional composite. The most known FGMs are made of a transition alloys from metal at one surface to ceramic at the opposite surface. This kind of FGMs was introduced as ultra high temperature resistant materials for nuclear reactor, chemical plants, heat engine components and aerospace vehicles. The thermal resistance of FGM is due to low thermal conductivity of ceramic [1, 2]. The static behavior of FGM plates was studied in literature with analytical and numerical methods. Thin FGM plate can be studied by a classical plate theory [3, 4]. For thick FG plates, some models take into account the transverse shear effect by using the First-order Shear Deformation Theory (FSDT) [5, 6]. In the FSDT theories, transverse shear is assumed to be constant through the shell thickness and thus require the computation of shear correction coefficients, [7–9]. In fact, shear correction coefficients are problem dependent and cumbersome. This limitation of the FSDT forced the development of High-Order Shear Deformation theory (HSDT) which includes the consideration of realistic parabolic variation of transverse ∗ Corresponding
author Email address: [email protected] (A. Frikha)
Preprint submitted to Elsevier
July 21, 2016
shear stress through the shell thickness in the isotropic cases. There is no need of a shear correction factor when using a HSDT but equations of motion are more complicated to obtain than those of the FSDT. The HSDT are introduced to develop a plate or shell elements, [10–14] among authors. An other possibility to introduce a HSDT is by using the enhanced assumed strain formulation (EAS) in solid-shell elements with quadratic transverse shear enhancement, [15, 16]. Recently, a unified formulation developed by Carrera and his co-workers (CUF), which can generate any refined theory, is developed in static and free vibration for laminate composites and FGM shells [17–21]. FGM shell structures are subjected to severe dynamic loads. Therefore, an assessment of natural frequency and transient response of structures seems required. Based on the Mindlin’s first-order shear deformation theory, free vibration of functionally graded plates and cylinder was investigated by Hosseini-Hashemi et al [22]. Using a discrete double directors shell element, free vibration analysis of FGM shell structures was studied by Wali et al. [12]. Based on layer-wise finite element, the dynamic response of functionally graded cylindrical shell was studied by Yas et al. [23]. Using a refined 8-node shell element, a forced vibration of FGM to arbitrary loading is investigated in [24]. The transient response of FGM plate is analyzed by [5]. The objective of this work is to present a general formulation for forced vibration of FGM shell using discrete double directors shell model. The formulation is investigated for the transient response of FGM shells. Square plate and spherical cap are used to compare the accuracy of the formulation. Transient deflection and axial stress are analyzed along the thickness. 2. Double directors shell model In this section, the geometry and kinematic of double directors shell model are described. The reference surface of the shell is assumed to be smooth, continuous and differentiable. Variables associated with the undeformed state will be denoted by upper-case letters and by a lower-case letter when referred to the deformed configuration. 2.1. Double directors shell kinematic assumption
Figure 1: Double directors shell model.
Parameterizations, which define material points of the shell, are carried out in terms of curvilinear coordinates ξ = (ξ 1 , ξ 2 , ξ 3 = z). The position vectors of any material point (q), whose normal projection on mid-surface is the material point (p), in the initial states C0 are given by: 2
X q (ξ 1 , ξ 2 , z) = X p (ξ 1 , ξ 2 ) + zD(ξ 1 , ξ 2 ),
z ∈ [−h/2, h/2] ,
(1)
where h is the thickness. X p and D are respectively a point of the reference surface and the initial shell director. The covariant basis (G1 , G2 , G3 ) is obtained from the position vector by (G1 , G2 , G3 )= (∂X q /∂ξ 1 , ∂X q /∂ξ 2 , ∂X q /∂z), which yields in base vectors relative to the initial state : Gα = Aα + zD ,α ; G3 = D , α = 1, 2.
(2)
The surface element dA in the initial state is given by: dA =
√
AdAξ ,
√
A = kA1 ∧ A2 k , dAξ = dξ 1 dξ 2 .
(3)
The covariant reference metric tensor G at a material point ξ is defined by: G = [Gi · Gj ] , i, j = 1, 2, 3.
(4) p
For use, the geometrical variable Det(G) and dV are related by dV = (G)dξ 1 dξ 2 dz, where Det(G) = p later p (G) = (|Gij |). With the assumption of a double directors shell model, the position vector of any point q, in the deformed configuration is given by (show Fig. 1) : xq (ξ 1 , ξ 2 , z) = xp (ξ 1 , ξ 2 ) + f1 (z)d1 (ξ 1 , ξ 2 ) + f2 (z)d2 (ξ 1 , ξ 2 ).
(5)
In the deformed state, the base vectors are: 0
0
g α = aα + f1 (z)d1,α + f2 (z)d2,α ; g 3 = f1 (z)d1 + f2 (z)d2 .
(6)
The metric tensor components in the deformed configuration Ct are separated into the in-plane and out-ofplane part components. With some approximations, the metric tensor can be written as: g ≈ aαβ + f1 (z)b1αβ + f2 (z)b2αβ αβ 0 0 gα3 ≈ f1 (z)c1α + f2 (z)c2α , (7) gij = g i · g j , 2 0 0 g33 ≈ f1 + f2 d
where aαβ , bkαβ and ckα (k = 1, 2) represent the covariant metric surface, the first curvature tensors and the shear, respectively. The parameter d denotes the thickness stretching. Taking into account d1 · d1 ≈ d2 · d2 ≈ d1 · d2 , these components can be computed as : aαβ = aα · aβ , ckα = aα · dk , bkαβ = aα · dk,β + aβ · dk,α , d = d1 · d1 , k = 1, 2.
(8)
Similar expressions for the in-plane and out-of-plane components of the metric tensor can be obtained in the case of the initial configuration C0 . Using the kinematic assumption, Eq. (7), the linearized strains can be written as follows [11, 12]: αβ = eαβ + f1 (z)χ1αβ + f2 (z)χ2αβ , α, β = 1, 2 . (9) 0 0 2α3 = f1 (z)γα1 + f2 (z)γα2 In matrix notation, the vectors e11 e22 e= 2e12
of membrane, bending and shear strains are given by: k χk11 γ1 χk22 , k = 1, 2. , γk = , χk = γ2k k 2χ12
The variation of the strain can be written, in the initial configuration, as: 3
(10)
δeαβ = 1/2(Aα · δx,β + Aβ · δx,α ), δγαk = Aα · δdk + δx,α · dk , δχkαβ = 1/2(Aα · δdk,β + Aβ · δdk,α + δx,α · dk,β + δx,β · dk,α )
k = 1, 2.
(11)
Or in matrix notation:
δe = B m · δx,
δχk = B bmk δx + B bbk δdk ,
δγ k = B smk δx + B sbk δdk ,
k = 1, 2,
(12)
where the matrix differential operators, relative to the initial configuration (Linear theory), are given by:
Bm
T d0k,1 ∂ξ∂ 1 AT1 ∂ξ∂ 1 T , AT2 ∂ξ∂ 2 = B bbk = d0k,2 ∂ξ∂ 2 , B bmk = T T T ∂ T ∂ A1 ∂ξ2 + A2 ∂ξ1 d0k,1 ∂ξ∂ 2 + d0k,2 ∂ξ∂ 1 " # T T d0k ∂ξ∂ 1 A1 B smk = , d0k = D, k = 1, 2. = , B T sbk AT2 d0k ∂ξ∂ 2
(13)
2.2. Third-order double director shell model The two Functions f1 (z) and f2 (z), used in Eq. (5), are related by: f1 (z) = z − f2 (z), f2 (z) =
4z 3 . 3h2
(14)
This gives the following transverse shear strain: 0
0
2α3 = f1 (z)γα1 + f2 (z)γα2 ,
0
0
f1 (z) = 1 − f2 (z).
(15)
The zero condition of the transverse shear stress on the top and bottom face of the shell, is simply obtained by imposing: 0 0 γα2 = 0, 2α3 = 1 − f2 (z) γα1 , f2 (±h/2) = 1. (16) The kinematic constrain γα2 = 0 will be imposed in a discrete form in the finite element approximation. 3. Weak form The numerical solution with the finite element method is based on the weak form of equilibrium equations. The three dimensional form of the latter is given as : % Z ! 2 X 1 k M k · δχ dA − Gext (Φ, δΦ) = 0, N · δe + T 1 · δγ + (17) G(Φ, δΦ) = A
k=1
where δΦ = (δu, δd1 , δd2 ) is an arbitrary variation. Gext is the external virtual work. δe, δχk and δγ 1 are the variations of shell strains. N , M k and T 1 are the membrane, bending and shear stresses resultants, which can be written in matrix form as : 11 1 N Mk11 T1 N 22 Mk22 N= , Mk = , T1 = , k = 1, 2. (18) T12 12 12 N Mk
Their components are defined as follows :
N αβ =
Z
h/2
σ αβ −h/2
p
G/A dz,
Mkαβ =
Z
h/2
fk (z)σ αβ
−h/2
4
p G/A dz,
T1α =
Z
h/2
−h/2
0
f1 (z)σ α3
p G/A dz, (19)
where σ αβ are the components of the stress tensor. The generalized resultant of stress R and strain vectors Σ are defined by: e N 1 χ M1 . (20) , Σ= R= χ2 M2 1 γ T 1 11×1 11×1 The weak form of the equilibrium equation can then be rewritten as: Z δΣT · RdA − Gext (Φ, δΦ) = 0. G(Φ, δΦ) =
(21)
A
Using strain-displacement relations, Eq. (12), the strain variation can then be given by: Bm 0 0 B bm1 B bb1 0 . δΣ = B δΦ, B = B bm2 0 B bb2 B sm1 B sb1 0
(22)
4. Constitutive equations In the case of an elastic isotropic constitutive model, the stress resultant R is expressed as: R = H T Σ,
(23)
where H T is the linear elastic modulus. Its expression is: H 11 H 12 H 13 H 12 H 22 H 23 HT = H 13 H 23 H 33 0 0 0
0 0 , 0 H 44
where:
(H 11 , H 12 , H 13 , H 22 , H 23 , H 33 ) =
Z
h/2
−h/2
1, f1 , f2 , f12 , f1 f2 , f22 Hdz.
(24)
(25)
f1 and f2 are given by Eq. (14). In Eq. (25), H 44 is given by: H 44 =
Z
h/2 −h/2
0 2 f1 H τ dz,
(26)
where H and H τ , in Eqs. (24) and (25), are in plane and out-of-plane linear elastic sub-matrices, which can be expressed in a Cartesian system as: 1 ν(z) 0 E(z) E(z) 1 0 , Hτ = ν(z) 1 0 , (27) H= 1 − ν 2 (z) 2(1 + ν(z)) 0 1 0 0 (1 − ν(z))/2
where E(z) and ν(z) are the Young’s modulus and the Poisson’s ratio respectively. As, a functionally graded material is considered, these material properties are assumed to vary through the shell thickness according to the rule of mixture written as: P (z) = (Pc − Pm )Vc + Pm , 5
(28)
where P (z), Pm and Pc denote respectively the effective material property, the properties of the metal and ceramic. Vc is the volume fraction, which varies according to the two general four-parameter power-laws distribution [25]. c p 1 z 1 z + + , (29) F GMI (a, b, c, p) : Vc (z) = 1 − a +b 2 h 2 h c p 1 z 1 z , (30) − − +b F GMII (a, b, c, p) : Vc (z) = 1 − a 2 h 2 h where p is the power-law index. The parameters a, b and c determine the material variation profile through the functionally graded shell thickness. Instead of the rule of mixture, Eqs. (28-30), another analytical method for estimating the effective properties can be used as Mori-Tanaka or self-consistent methods [26]. 5. Finite element formulation In this section, the numerical implementation of the presented shell theoretical formulation based upon a four node shell element is established. Using the isoparametric concept, the mid-surface (X) and surface displacement field (u = x − X) are approximated by: X=
4 X
N I XI ,
u=
I=1
4 X
N I uI ,
(31)
I=1
where X I ∈ R3 and uI ∈ R3 are position and displacement at node I. N I are the standard bilinear isoparametric shape functions. The first director vector d1 is approximated with the same functions d1 =
4 X
N I d01I ,
(32)
I=1
where d01I denote the first directors at the nodal points in the initial configuration C0 . 5.1. Local Cartesian system Let n0 = A1 ∧ A2 / kA1 ∧ A2 k be the normal field to the mid-surface in the initial state C0 . A local Cartesian system with base vectors n01 , n02 , n0 , can be defined by means of an orthogonal transformation. The Jacobian transformation J from the basis n01 , n02 to {A1 , A2 } is defined by: 0 n1 · A1 n02 · A1 J= . (33) n01 · A2 n02 · A2 Since the formulation is developed in local Cartesian coordinate, the derivatives of the shape functions need to be transformed as follows: ( I ) I N,1 N ,1 −1 . (34) = [J ] I I N,2 N ,2 5.2. Membrane and first bending strain field The shell membrane part of the problem is considered. The strain-displacement relation is: δe = B m Φn , where B m is the discrete membrane strain-displacement operator defined as:
6
(35)
B Im =
B Imm
0 0
,
B Imm
T I n01 N ,1 T I = . n02 N ,2 T T I I 0 0 n1 N ,2 + n2 N ,1
(36)
For the first bending part, the strain-displacement relation is given by: δχ1 = B 1 · δΦn ,
(37)
where B 1 is the discrete first bending strain-displacement operator: T T I I d01,1 N ,1 n01 N ,1 T T I I I I = B I1 = B I1m B I1b 0 , B1m , B1b = d01,2 N ,2 n02 N ,2 , T T T T I I I I d01,1 N ,2 + d01,2 N ,1 n01 N ,2 + n02 N ,1 where d01,a =
P4
(38)
I
I=1
N ,a d01I , a = 1, 2.
5.3. Construction of the assumed natural transverse shear strain field A typical isoparametric finite element is considered as depicted in Fig. 2. A, B, C and D denote the mid-points of the element boundaries set. Following [27], the assumed natural transverse shear strain field is expressed by: δγ11 (1 − η)δγ11 (B) + (1 + η)δγ11 (D) 1 , (39) δγξ = = (1 − ξ)δγ21 (A) + (1 + ξ)δγ21 (C) δγ21 where γ21 (A), γ11 (B), γ21 (C) et γ11 (D) are the strains at points defined in Fig. 2.
Figure 2: Assumed strain construction of the isoparametric shell element.
δγ 1ξ is related to δΦn by: δγ 1ξ = B sξ δΦn ,
(40)
where B sξ is the discrete shear strain-displacement operator expressed as: Bsξ =
"
T
T
T
T
N1,1 d01B
N2,1 AT 1B
0
N2,1 d01B
N2,1 AT 1B
0
N3,1 d01D
N3,1 AT 1D
0
N4,1 d01D
N3,1 AT 1D
0
N1,2 d01A
N4,2 AT 2A
0
N2,2 d01C
N3,2 AT 2C
0
N3,2 d01C
N3,2 AT 2C
0
N4,2 d01A
N4,2 AT 2A
0
T
T
T
T
#
. (41)
In the local Cartesian system, the transverse shear strains are: δγ 1 = B s · δΦn ,
B s = J −1 B sξ .
7
(42)
5.4. Second bending strain field: discrete constraints The shear part relative to the second director vector d2 will be vanished in a discrete form. A quadratic interpolation is selected as proposed in the works of [11, 28]: d2 =
4 X
N I d02I +
8 X
PK αK t0K ,
δd2 =
K=5
I=1
4 X
N I δd02I +
I=1
8 X
PK δαK t0K ,
(43)
K=5
where (I) and (K) represent respectively a node of the element and the mid-point of the element boundaries. δαK are variables associated to δd2 on the element boundaries. The vector t0K is unit and its direction is defined by the position of the nodes couple (I, J) : t0K = (xJ − xI )/LK ,
LK = kxJ − xI k ,
(44)
where LK is the I − J side length. The shape functions PK are quadratic [11]. By introducing the vanishing shearing hypothesis, on top and bottom faces, over the element boundaries under integral form, the transverse strain is given for side (I, J) by: Z
which leads to: δd2 =
4 X
N I δd2I +
I=1
J I
2 δγsz ds = 0,
(45)
8 X 3 1 1 t0K . − PK 0 0 2 2 (δd + δd ) · t L (δu + δu ) · d I J K I J K K K=5
(46)
In a matrix form, δd2 becomes: δd2 =
4 X I=1
M Id δuI + M Ir δd2I ,
(47)
where matrixes M Id and M Ir are given by the following expressions: M Id = PK tdIK + PM tdIM ,
tdIK =
3 0 t ⊗ d0K , 2LK K
(48)
3 0 t ⊗ t0K . (49) 4 K The (K) and (M ) are the two mid-side of every side of the quadrilateral, which are bound to the node (I). Finally, the second bending strain is expressed as: M Ir = N I I + PK ttIK + PM ttIM ,
ttIK =
δχ2 = B 2 · δΦn ,
(50)
where B 2 is the discrete second bending strain-displacement operator: B I2 = B I2m and B I2b are given by:
B I2m
B I2m
0 B I2b
T T d02,1 N ,I1 +n01 M Id,1 T T = , d02,2 N ,I2 +n02 M Id,2 T T T T 0 0 I I I 0 I 0 d2,1 N ,2 +n1 M d,2 + d2,2 N ,1 +n2 M d,1
.
(51)
B I2b =
T n01
T n01 · M Ir,1 T n02 · M Ir,2 . T I I 0 · M r,2 + n2 · M r,1
Finally, the continuum generalized strain δΣ, given by Eq. (17) will be in an approximate form as: 8
(52)
δe δχ1 = B · δΦn , δΣ = 2 δχ1 δγ
Bm B1 B= B2 . Bs
Using Eq.(23), the internal virtual work becomes: Z Gint = δΣT · RdA = δΦTn KΦn , A
where K is the stiffness matrix given by K = 5.5. Nodal transformation
R
A
(53)
(54)
B T H T BdA.
In all equations, δdk and δdk,α are the variation of the directors and there derivatives. These variations can be written either in spatial description: δdk = δθ k ∧ dk = Λk δθ k ,
ek , Λk = −d
k = 1, 2,
(55)
ek is the skew-symmetric tensor such that d ek dk = 0 , or in material description: where d
e 3, δdk = Qk δθ k E 3 = Λk δθ k , Λk = Qk E (56) t where we assumed that dk = Qk E 3 and E 3 = 0 0 1 . A spatial description leads to a shell problem with 9 DOF/node and the material description leads to a shell problem with 7 DOF/node, where the transformation Λ take the following form: Λk =
−t2k
t1k
3×2
.
(57)
All numerical examples given in section 6, are obtained with the material description and 7DOF/node. 5.6. Mass matrix Several different schemes for deriving element mass matrices exist in the literature. One finds the consistent mass matrix and different lumped mass matrix constructions. The variation of the kinetic energy is written as: Z (58) Gext = δxT q ρx¨q dv, v
¨ q denote virtual displacement where ρ being the material density and vary along the thickness, δxq and x and acceleration of any point in the shell. By using the displacement field, Eq. (5), the variation of the kinetic energy becomes: Z Z h/2 1 f1 f2 ¨ Gext = δΦT ρΦdA, ρ f1 f12 f1 f2 dz, (59) ρ= A −h/2 f2 f1 f2 f22
the finite element discretization of Gext can then be written in the following form: ¨ Gext = δΦT n M Φn .
(60)
M is the mass matrix. This matrix can be either the diagonal mass matrix or the consistent mass matrix. To construct a diagonal mass matrix, several alternatives were investigated by [29], in a free vibration analysis of Mindlins plate. In the present numerical implementation, the efficient scheme considered in [12] is used where the diagonal mass matrix is given, with 7DOF/node, as follows: 9
M ii
with :
M1i
=
M1i
0
M1i M2i
0
M2i M3i M3i
,
(61)
R ρ N Ni dA A ρkk dA , k = 1, 2, 3. ( ρkk Ni Ni ) dA
(62)
The governing equation for direct transient analysis can then be written as: n o ¨ + [K] {Φ} = {F} . [M] Φ
(63)
M ki =
R
A R kkP i
The above equation becomes transient when the mechanical force, {F} varies with time. It is solved by using Newmark’s method. 5.7. Newmark algorithm The Newmark’s direct integration scheme is used to approximate the time derivatives and to solve the governing dynamic equation. This n o n oimplicit method allows to compute the solution at time t + ∆t from ¨t . The Newmark’s formula to process the solution for damped system known vectors {Ut }, U˙ t and U is given by: K {Ut+∆t } = {Rt+∆t } . (64) n o n o ¨0 is computed by For a given initial displacement {U0 } and velocity U˙ 0 , the initial acceleration U solving Eq. (63). K is update by using: K = [K] +
1 γ [M] + [C] . β∆t2 β∆t
(65)
In our study, β and γ are equal respectively to 0.25 and 0.5. For each time step t + ∆t, it is required to proceed as follows: • to compute the residu {Rt+∆t } by using : n o 1 1 n˙ o 1 ¨t {Rt+∆t } = {F t+∆t } + [M] {U } + − 1 U U + + t t β∆t2 β∆t 2β n o n o γ γ γ ¨t U˙ t + ∆t U , {U t } + −1 −1 [C] β∆t β 2β
• to update K ,
10
(66)
• to solve Eq. (64) for {Ut+∆t }, n o n o ¨ t+∆t and the velocity U˙ t+∆t at time t + ∆t by using : • to compute the acceleration U n o n o 1 n o 1 ¨t , ˙ ¨ t+∆t = 1 U − β U {U } − {U } − ∆t − U t t+∆t t β ∆t2 2
(67)
n o n o n o n o ¨ t + γ∆t U ¨ t+∆t . U˙ t+∆t = U˙ t + ∆t (1 − γ) U
(68)
and
6. Numerical results To show the dynamic response of FGM shells with a discrete double directors shell element, the authors offer numerical results for plate and spherical caps. The model was validated in each structure for FGM or isotropic material with literature. The effect of functionally graded material is then investigated for each structure. In this section, the dynamic response of FGM shell, using the developed SHO4 element, is presented. The abbreviation SHO4 is used to identify the present linear discrete double directors shell element. 6.1. Square plate The forced vibration was performed on a square plate of side length l = 0.2 m and thickness h = 0.01 m, simply supported on all its edges. The plate is subjected to suddenly applied uniform pressure loading equal to q0 = −106 N/m2 [30]. The non-dimensional parameters are center deflection, w = wEm h/(q0 l2 ), stress, p 2 2 2 σ = σh /(|q0 | l ), and time, t = t Em /(l ρm ), where Em is the Young modulus of metal. The analysis was conducted using Aluminum and Zirconia (ceramic). The material properties of Aluminum and Zirconia are given in Tab.1. Table 1: Material properties of FGM plate [30].
Parameter Aluminum Zirconia
Young’s modulus (GPa) 70 151
Poisson’s ratio 0.3 0.3
Density (kg/m3 ) 2707 3000
Only one-half of the square plate was modeled. The mesh consists of 128 nodes and 192 elements. One layer of SHO4 elements was used in the thickness direction. Fig.3 shows the dynamic response of simply supported FGM square plate under suddenly applied uniform pressure loading. The non-dimensional center deflection of the plate is considered here, and the present results for different values of the power law exponent n are compared with the results from literature [5, 30]. The discrepancy between results obtained from two methods is closely small. It should be noted that a non linear model is used in [5, 30]. The non-linearity effect is more significant for lower stiffness plate corresponding to metal. As known in literature, non-linear analysis gives a reduced deflection than linear one. The gap between the present and the reference results for metal can be explained by the taking into account the non-linearity effect in [5, 30]. The magnitude of vibration is the maximum for metallic plates and a minimum for the ceramic plate. The magnitude of vibration increases smoothly with the grading index p. The bending rigidity of ceramic, corresponding to the lower magnitude, is higher than the bending rigidity of metal, corresponding to the highest magnitude. It is seen that the frequency of vibration of the ceramic plates is much higher than that of metallic plates since the maximum magnitude of ceramic plate is located before the maximum magnitude of the metallic plate along the time axis. 11
Figure 3: Temporal variation of center deflection of simply supported FGM square plate under suddenly applied uniform pressure of intensity q0 = −106 N/m2 , compared with [30].
The computational domain for the dynamic analysis with clamped boundary conditions is the same as that of Fig. 3. The analysis was performed for the same two material combinations. Fig. 4 shows the temporal response of the plates under suddenly applied mechanical loading of −106 N/m2 . The deflection of clamped FGM plates is much lower than the simply supported case. The effective stiffness of clamped plate is higher than the simply supported plates.
Figure 4: Non-dimensional center deflection of the of clamped FGM square plate under suddenly applied uniform pressure of intensity q0 = −106 N/m2 .
Fig. 5 shows the temporal variation of axial stress in the center of the FGM plates under suddenly applied mechanical loading of −106 N/m2 by varying the grading index p. Left column of Fig. 5 (a, c, e) depicts the axial stress respectively in z = −h/2, z = 0 and z = h/2 of clamped FGM plates. Right column of Fig. 5 (b, d, f) depicts the axial stress respectively in z = −h/2, z = 0 and z = h/2 of simply supported FGM plates. Usually, the axial stress modulus of clamped plates are lower than simply supported. Since the applied pressure is negative, the axial stress on the bottom, h = −h/2, (resp. top, h = h/2) is mostly positive (resp. negative). At z = −h/2, the axial stress modulus decreases between the grading index p = 0, corresponding to ceramic material, and near to p = 0.2 and then increases to metallic material. The maximum axial stress modulus of metallic and ceramic plates are similar for both clamped and simply supported plates (see Fig. 5 (a-b)). At z = 0, it is shown that the axial stress is zero for isotropic material (metallic and ceramic) for clamped (Fig 5 (c)) and simply supported (Fig 5 (d)) plates, according to the plate theory. The axial stress 12
modulus increases between p = 0, corresponding to ceramic material, and the grading index near p = 2 then decreases to metallic material. At z = h/2, the axial modulus increases between the grading index p = 0 and near to p = 2 and then decreases to metallic material.
Figure 5: Non-dimensional axial stress of the center of FGM square plate under suddenly applied uniform pressure of intensity q0 = −106 N/m2 for (a, c, e): simply supported FGM plate, (b, d, f): clamped FGM plate, (a, b): z = −h/2, (c, d): z = 0, (e, f): z = h/2.
Fig. 6 shows the plots of the axial stress through the thickness at non-dimensional time t = 5 of simply supported and clamped plate under uniform pressure applied on the top surface. The axial stress is compressive for z = h/2 and tensile for z = −h/2 for both simply supported and clamped plates. The axial stress varies linearly versus the thickness for isotropic material. For the different grading index chosen, the plate corresponding to p = 2 yielding to the maximum compressive stress at the top surface. At the bottom surface, the plate corresponding to p = 0.2 yielding to the minimum tensile stress. These phenomena, observed at p = 0.2 and p = 2, are in good agreement with [5]. 6.2. Spherical cap Dynamic response of clamped spherical cap with isotropic material subjected to concentrated apex load is presented. The characteristics of the shell and the concentrated load are shown in Fig. 7. Its parameters are given in Tab. 2: 13
Figure 6: Non-dimensional axial stress in the center of FGM square plate under suddenly applied uniform pressure of q0 = −106 N/m2 for (a) simply supported FGM plate, (b) clamped FGM plate.
Table 2: Parameters of spherical cap [31].
Parameter Value
R 4.76
θ 10.9◦
E 107
ν 0.3
h 0.01576
H 0.0859
ρ 0.000245
F 100
Only one-quarter of the workpiece was modeled. The mesh consists of 217 nodes and 192 elements. One layer of SHO4 elements was used in the thickness direction with Simpson integration (see Fig. 8).
Figure 7: Characteristics of the spherical cap.
A concentrated non-dimensional load of F = 100 is applied between 0 and 250 µs. The response of the cap, represented by the normalized vertical displacement at shell apex, is compared to linear dynamic analysis obtained by [31] for homogeneous material. The present results are found to be in good agreement with those reported by [31]. Since the angle of spherical cap is small with a great radius, the geometry is close to plate. In order to show the effect of the radius on dynamic behavior of spherical shell, an hemisphere (θ = 90◦ ) with the same parameters than Tab. 2 is considered. The spherical cap is subjected to suddenly applied uniform pressure loading equal to q0 = −106 N/m2 . The non-dimensional parameters are center deflection, w = p wEm h/(q0 R2 ), stress, σ = σh2 /(|q0 | R2 ), and time, t = t Em /(R2 ρm ), where Em is the Young’s modulus of metal. The analysis was conducted using Aluminum and Zirconia (ceramic). The Young’s modulus, Poisson’s ratio and the density for Aluminum and Zirconia are given in Tab. 1. Fig. 9 shows the transient non-dimensional center deflection of clamped FGM hemisphere under suddenly applied uniform pressure loading. The magnitude of vibration is the maximum for metallic hemisphere and a minimum for the ceramic case. The transient magnitude of FGM spherical cap is intermediate to that of metal and ceramic. It increases smoothly with the grading index p. The bending rigidity of ceramic, corresponding to the lower magnitude, is higher than the bending rigidity of metal, corresponding to the 14
Figure 8: Transient center deflection of isotropic spherical cap.
higher magnitude. It is seen that the frequency of vibration of the ceramic hemisphere is higher than that of metallic case since the maximum magnitude of ceramic plate is located before the maximum magnitude of the metallic plate along the time axis.
Figure 9: Non-dimensional center deflection of clamped FGM hemisphere under suddenly applied uniform pressure q0 = −106 N/m2 .
Fig. 10 shows the temporal variation of axial stress in the center of the FGM plates under suddenly applied mechanical loading of 106 N/m2 by varying the grading index p. Fig. 10 (a-c) depicts the axial stress respectively in z = −h/2, z = 0 and z = h/2 of clamped FGM hemisphere. Since the applied pressure on the hemisphere is positive, the membrane effect yields positive the axial stress at z = −h/2, z = 0 and z = h/2 (see Fig. 10). At z = −h/2, the axial stress modulus decreases between the grading index p = 0, corresponding to ceramic material, and near to p = 0.2 and then increases to metallic material. At z = h/2, the axial modulus increases between the grading index p = 0 and near to p = 2 and then decreases to metallic material. The same phenomenon is observed for plates and verified with [5]. Fig. 11 shows the non-dimensional axial stress along the thickness at non-dimensional time t = 1 for a clamped FGM hemisphere subjected to a uniform pressure 106 N/m2 . It is shown that the ceramic spherical cap yields to a higher membrane stiffness than metallic one. The variation of axial stress along the thickness is small for isotropic material. This variation is great for FGM spherical cap.
15
Figure 10: Non-dimensional axial stress of the center of clamped FGM hemisphere under suddenly applied uniform pressure of intensity q0 = 106 N/m2 for (a) z = −h/2, (b) z = 0 and (c) z = h/2.
16
Figure 11: Non-dimensional axial stress in the center of clamped FGM hemisphere under suddenly applied uniform pressure q0 = −106 N/m2 .
Conclusions Finite element model based on a discrete double-directors shell element is presented to study the dynamic response of 3D-FGM-shell structures subjected to mechanical loads. The grading function is assumed along the thickness direction. Dimensionless deflection and stress of functionally graded plates and spherical caps are computed by the present model. To validate the results of the present study and demonstrate its accuracy, the results are compared with the existing ones in the literature for plates and spherical caps. A very good agreement among the results confirms the high accuracy of the current approach. Transient response of deflection and axial stress are analyzed by varying the grading index. [1] D. P. H. Hasselman and G. E. Youngblood. Enhanced thermal stress resistance of structural ceramics with thermal conductivity gradient. Journal of The American Ceramic Society, 61 (1,2):49–53, 1978. [2] M. Niino and S. Maeda. Recent development status of functionally gradient materials. International Iron and Steel Institute, 30:699–703, 1990. [3] S.-H. Chi and Y.-L. Chung. Mechanical behavior of functionally graded material plates under transverse load. part 1: Analysis. International Journal of Solids and Structures, 43:3657–3674, 2006. [4] S.-H. Chi and Y.-L. Chung. Mechanical behavior of functionally graded material plates under transverse load. part 2: Numerical results. International Journal of Solids and Structures, 43:3675–3691, 2006. [5] G. N. Praveen and J. N. Reddy. Nonlinear transient thermoelastic analysis of functionally graded ceramic-metal plates. International Journal of Solids and Structures, 35:4457–4476, 1997. [6] J.N. Reddy, C.M. Wang, and S. Kitipornchai. Axisymmetric bending of functionally graded circular and annular plates. European Journal of Mechanics - A/Solids, 18:185–199, 1999. [7] T. K. Nguyen and G. Bonnet. First-order shear deformation plate models for functionally graded materials. Composite Structures, 83:25–36, 2008. [8] M. Shariyat and M. M. Alipour. A novel shear correction factor for stress and modal analyses of annular fgm plates with non-uniform inclined tractions and nonuniform elastic foundations. International Journal of Mechanical Sciences, 87:60–71, 2014. [9] A. Hajlaoui, A. Jarraya, K. El Bikri, and F. Dammak. Buckling analysis of functionally graded materials structures with enhanced solid-shell elements and transverse shear correction. Composite Structures, 132:87–97, 2015. [10] F. Tornabene. Free vibration analysis of functionally graded conical, cylindrical shell and annular plate structures with a four-parameter power-law distribution. Computer Methods in Applied Mechanics and Engineering, 198:2911–2935, 2009. [11] M. Wali, A. Hajlaoui, and F. Dammak. Discrete double directors shell element for the functionally graded material shell structures analysis. Computer Methods in Applied Mechanics and Engineering, 278:388–403, 2014. [12] M. Wali, T. Hentati, A. Jarraya, and F. Dammak. Free vibration analysis of fgm shell structures with a discrete double directors shell element. Composite Structures, 125:295–303, 2015. [13] Mallikarjuna and T. Kant. Effect of cross-sectional warping of anisotropic sandwich laminates due to dynamic loads using a refined theory and c0 finite elements. International Journal for Numerical Methods in Engineering, 35:2031–2047, 1992. [14] T. Kant and J. R. Kommineni. Geometrically non-linear transient analysis of laminated composite and sandwich shells with a refined theory and c0 finite elements. Computers & Structures, 52:1243–1259, 1994. [15] N. D. Quy and A. Matzenmiller. A solid-shell element with enhanced assumed strains for higher order shear deformations in laminates. Technishe Mechanik, 28:334–355, 2008. [16] A. Hajlaoui, M. Wali, M. Ben Jdidia, and F. Dammak. An improved enhanced solid shell element for static and buckling analysis of shell structures. Mechanics & Industry, 17(5):510, 2016.
17
[17] A.M.A. Neves, A.J.M. Ferreira, E. Carrera, M. Cinefra, C.M.C. Roque, R.M.N. Jorge, and C.M.M. Soares. Free vibration analysis of functionally graded shells by a higher-order shear deformation theory and radial basis functions collocation, accounting for through-the-thickness deformation. European Journal of Mechanics - A/Solids, 37:24–34, 2013. [18] M. Cinefra, S. Belouettar, M. Soave, and E. Carrera. Variable kinematic models applied to free vibration analysis of functionally graded materials shells. European Journal of Mechanics - A/Solids, 29:1078–1087, 2010. [19] M. Cinefra, E. Carrera, L. Della Croce, and C. Chinosi. Refined shell elements for the analysis of functionally graded structures. Composite Structures, 94:415–422, 2012. [20] M. Cinefra and S. Valvano. A variable kinematic doubly-curved mitc9 shell element for the analysis of laminated composites. Mechanics of advanced materials and structures, 23(11):1312–1325, 2016. [21] M. Cinefra and E. Carrera. Shell finite elements with different through-the-thickness kinematics for the linear analysis of cylindrical multilayered structures. International Journal for Numerical Methods in Engineering, 93:160–182, 2013. [22] Sh. Hosseini-Hashemi, J. N. Fadaee, and M. Es’haghi. A novel approach for in-plane/out plane frequency analysis of functionally graded circular/annular plates. International Journal of Mechanical Sciences, 52:1025–1035, 2010. [23] M. H. Yas, M. Shakeri, M. Heshmati, and S. Mohammadi. Layer-wise finite element analysis of functionally graded cylindrical shell under dynamic load. Journal of Mechanical Science and technology, 25(3):597–604, 2011. [24] W. Y. Jung and S. C. Han. Transient analysis of fgm and laminated composite structures using a refined 8-node ans shell element. Composites: Part B, 56:372–383, 2014. [25] Z. Su, G. Jin, S. Shi, T. Ye, and X. Jia. A unified solution for vibration analysis of functionally graded cylindrical, conical shells and annular plates with general boundary conditions. International Journal of Mechanical Sciences, 80:62–80, 2014. [26] S. S. Vel and R. C. Batra. Three-dimensional exact solution for the vibration of functionally graded rectangular plates. Journal of Sound and Vibration, 272:703–730, 2004. [27] K. J. Bathe and E. Dvorkin. A four-node plate bending element based on mindlin/reissner plate theory and a mixed interpolation. International Journal for Numerical Methods in Engineering, 21:367–383, 1985. [28] F. Dammak, S. Abid, A. Gakwaya, and G. Dhatt. A formulation of the non linear discrete kirchhoff quadrilateral shell element with finite rotations and enhanced strains. Revue Europeenne des Elements Finis, 14:7–31, 2005. [29] T. A. Rock and E. A. Hinton. A finite element method for the free vibration of plates allowing for transverse shear deformation. Composite Structures, 6:37–44, 1976. [30] J. N. Reddy. Analysis of functionally graded plates. International Journal for Numerical Methods in Engineering, 47:663–684, 2000. [31] L. A. Duarte Filho and M. A. Armando. Geometrically nonlinear static and dynamic analysis of shells and plates using the eight-node hexahedral element with one-point quadrature. Finite Elements in Analysis and Design, 40:1297–1315, 2004.
18