- Email: [email protected]
A refined three-node triangular element based on the HW variational theorem for multilayered composite plates
Accepted Manuscript A refined three-node triangular element based on the HW variational theorem for multilayered composite plates Wu Zhen, Ma Rui, Chen Wanji PII: DOI: Reference:
S0263-8223(16)32001-3 http://dx.doi.org/10.1016/j.compstruct.2016.11.040 COST 8000
To appear in:
Composite Structures
Received Date: Revised Date: Accepted Date:
29 September 2016 5 November 2016 11 November 2016
Please cite this article as: Zhen, W., Rui, M., Wanji, C., A refined three-node triangular element based on the HW variational theorem for multilayered composite plates, Composite Structures (2016), doi: http://dx.doi.org/10.1016/ j.compstruct.2016.11.040
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.
A refined three-node triangular element based on the HW variational theorem for multilayered composite plates Wu Zhen*, Ma Rui, Chen Wanji Key Laboratory of Liaoning Province for Composite Structural Analysis of Aerocraft and Simulation, Shenyang Aerospace University, Shenyang 110136, China
Abstract: For the three-node triangular elements, the displacement parameters are linear variation within one element, so that the second and the third derivatives of displacement parameters are close to zero. Therefore, the three-node triangular element will encounter difficulty to accurately predict transverse shear stresses by integrating the three-dimensional equilibrium equation through the thickness direction. Thus, the few three-node triangular plate elements can accurately yield the transverse shear stresses of multilayered composite plates. A higher-order zig-zag model accounting for the zero transverse shear strain conditions at the top and the bottom surfaces is firstly presented in this work. Based on the proposed model, a refined three-node triangular plate element satisfying the requirement of C1 weak-continuity conditions in the inter-element is developed to accurately produce the distributions of displacements and stresses through the thickness of multilayered composite plates. To obtain the improved transverse shear stresses, the simplified three-dimensional equilibrium equation has been a priori used. It is significant that the higher-order derivatives of displacement parameters have been taken out from the expression of transverse shear stresses by employing the three-field Hu-Washizu (HW) variational principle, which is convenient for the finite element implementation. Performance of the present element is tested by comparing with three-dimensional elasticity solutions as well as the reference results evaluated using other models. A detailed study is conducted to highlight the influences of boundary condition, number of layers and aspect ratio on the static response of multilayered composite plates, and numerical results can show the accuracy and range of applicability for the proposed element. Keywords: Three-node triangular element; Higher-order zig-zag model; Three-field Hu-Washizu variational principle; Multilayered composite plates; Transverse shear stress.
†
E-mail address: [email protected] (Wu Z.) 1
1. Introduction Due to their high strength-to-weight ratio, stiffness-to-weight ratio and superior fatigue characteristics, the fiber-reinforced composite materials are widely used as the fundamental members of aerospace, such as parts of wing, fuselage as well as satellite body etc. In order to apply the advantage of the high specific strength of such materials, it needs to accurately predict structural behaviors of these laminates. In the early stages of the development of models, the three-dimensional models [1,2] are usually employed to capture the significant stresses arising at the interfaces of the laminated
composites.
Bogdanovich and
Yushanov [3] also developed
a
three-dimensional displacement-assumed variational-analysis approach for bending analysis of laminated composite plates. Three-dimensional analysis can serve as benchmark solutions to assess the performance of approximate laminated plate models, whereas three-dimensional analysis is limited because it needs the computationally expensive cost as some layers much thinner compared to others [4]. Since three-dimensional analytical solutions for laminated composite structures are limited, the two-dimensional plate models have been developed to study the practical problems by using finite element methods and some analytical procedure in very simple cases. Based on Kirchhoff's hypothesis, the classical laminated plate theory is employed for the analysis of laminated composites. Owning to neglect of transverse shear and normal strains in the laminate, the classical laminated plate theory is not adaptable to analysis of multilayered composite plates with relatively soft transverse shear modulus [5]. In order to consider the influence of transverse shear strain, a number of first-order shear deformation theories [6-9] have been proposed. Nevertheless, experience of many investigators [9-11] has shown that the first-order theory is inadequate for prediction of the interlaminar stresses which are used in predicting the delamination onset in laminated composite structures. Thus, various higher-order shear deformation models including the effects of transverse shear deformation have been developed for the determination of transverse shear stresses in laminated composites. Based on the higher-order model considering transverse normal deformation, a closed form result for laminated composite plates has been presented by Lo et al. [12,13]. Marur and Kant [14] employed three higher-order models considering the warping of the cross-section for the free vibration analysis of laminated composite 2
and sandwich beams. Subsequently, the higher-order models are extended to the static and dynamic analysis of laminated composite and sandwich plates [15,16]. In terms of the third-order model, Babu and Kant [17] presented a quadrilateral element for buckling analysis of skew laminated composite and sandwich panels. The proposed element was also used to study the thermal buckling behaviors of skew laminated composite and sandwich plates [18]. On the basis of Reddy’s theory [19], the higher-order quadrilateral elements have been proposed for free vibration analysis of isotropic, laminated composite and sandwich plates by assuming interpolations for the shear strain [20]. Moreover, the proposed elements are employed for the buckling and vibration analysis of laminated composite and sandwich plates subjected to initial stresses [21]. By means of Reddy’s theory, Sheikh and Chakrabarti [22] proposed a six-node triangular plate element for bending analysis of laminated composite plates and skew plates. Based on the third-order zig-zag theory [23], a three-node triangular plate element has been proposed to study the mechanical, thermal and electric behaviors [24]. To obtain accurately transverse shear stresses, the integration of three-dimensional equilibrium equation through the thickness of laminates has to be used. Once the equilibrium equation approach is applied, the second and the third derivatives of displacement parameters will be involved into the three-dimensional equilibrium equation. For a three-node triangular element, the displacement parameters are linear variation within one element, so that the second and the third derivatives of displacement parameters are close to zero. Thus, the equilibrium equation approach will encounter difficulty to accurately yield the transverse shear stresses for the three-node triangular element. Therefore, the stress smoothing technique within whole domain has to be used in regular meshes [24]. In order to accurately produce transverse stresses by integrating local equilibrium equation, an eight-node quadrilateral element with quadratic approximation has been proposed for bending analysis of laminated composite plates [25]. However, such element contains the second-order derivatives of transverse displacement as nodal unknowns, which are not encouraged in a practical analysis [26]. Versino et al. [27] proposed a triangular plate element for the bending analysis of homogeneous, multilayered composite and sandwich plates, whereas the transverse shear stresses obtained from the proposed element can not satisfy the continuity conditions at interfaces. Employing the principle of virtual displacements (PVD) and the Reissner mixed variational theorem (RMVT), Carrera and Demasi [28] derived the finite element formulations for the 3
bending analysis of laminated composite plates [29]. For the Reissner mixed variational theorem (RMVT), transverse stress variables are considered as the nodal parameters of finite element, which will increase the computational cost. In addition, it has been indicated that the third-order displacement models are only adequate in predicting the global response such as natural frequency and buckling stress, but such models are inadequate to predict the local response such as displacement and stress distributions in each ply of laminates [30]. Moreover, Matsunaga [31] also pointed out that it seems to be a hasty conclusion to study the global and the local behaviors of laminated composite plates using the third-order models, so that a nine-order model has been proposed for the prediction of displacement and stress distribution in each ply of laminated composite and sandwich plates. In review article [32], Carrera showed that Lekhnitskii [33] was the first to present a zig-zag model for the analysis of laminated composite beam. In addition, Ambartsumian [34] also proposed an independent zig-zag model. Differing from previous models, this paper proposes a third-order zig-zag model accounting for the zero transverse shear strain conditions at the top and the bottom free surfaces. Based on the proposed model, a refined three-node triangular plate element satisfying the requirement of C1 weak-continuity conditions in the inter-element is developed to accurately predict the displacement and stress distributions in each ply of laminates. To obtain the improved transverse shear stresses, the simplified three-dimensional equilibrium equation has been a priori used. It is significant that the higher-order derivatives of displacement parameters have been eliminated from the expression of transverse shear stresses by employing the three-field HW variational principle [35,36], which is convenient for the present model's finite element implementation. Performance of the proposed element is tested by comparing with three-dimensional elasticity solutions as well as the reference results computed using other models. A detailed study is conducted to highlight the influences of boundary condition, number of layers and aspect ratio on the static response of multilayered composite plates, and numerical results show the accuracy and range of applicability for the proposed element. Moreover, a number of examples also show that the present third-order displacement model can accurately yield the displacement and stress distribution in each ply of laminates.
2. Higher-order zig-zag theory (HZT) satisfying the bounding surface 4
free traction condition Owning to the different mechanical properties at each layer, the in-plane displacements of laminated composite structures show a zig-zag shape distribution through the thickness direction [37]. In order to model the zig-zag pattern of the in-plane displacements along the thickness, the zig-zag function proposed by Murakami [38] is introduced to displacement field of the third-order model, which can be written as
u = u0 + zu1 + z 2u2 + z 3u3 + M k ( z ) uz v = v0 + zv1 + z 2v2 + z 3v3 + M k ( z ) v z
(1)
w = w0 where k is the number of layers in laminates; u0, v0 and w0 are the displacements of a generic point in the middle surface; u1 and v1 is the rotations of normal to the middle surface about y axis and x axis, respectively; u2, u3, v2 and v3 are the higher-order terms in Taylor's series expansions; uz and vz are the generalized displacement variables associated with the zig-zag function Mk proposed by Murakami [38]. The zig-zag function Mk can be written as k
M k ( z ) = ( −1) ζ k
in which ζ k = a k z − bk , a k =
(2) z + zk 2 , bk = k +1 , and zi has been defined in z k +1 − z k z k +1 − z k
Figure 1. Employing the zero shear strain boundary conditions at the top and bottom surfaces of the plate, the higher-order zig-zag theory can be obtained, which can be given by ∂w0 ∂x ∂w v = v0 + Φ1v1 + Φ 2 vz + Φ3 0 ∂y w = w0 u = u0 + Φ1u1 + Φ 2u z + Φ3
(3)
where
5
Φ1 = z −
4z3 3h 2
z 2 2 z 3 − 3z1 z 2 M1 ( z1 ) 2 z 3 − 3z1 z 2 M n ( zn +1 ) Φ2 = M ( z) − + − 2 12 zn2+1 ∂z ∂z 2 z1 12 zn +1 3 4z Φ3 = − 2 3h k
(4)
For linear elasticity, strain of the present model HZT can be expressed as ε ε= I ε T
(5)
in which T
ε I = {ε x , ε y , γ xy } = B I EI
(6)
T
ε T = {γ xz , γ yz } = BT ET where
BI X BI = 0 0
0 BI Y
B BT = TX 0
0 BTY
EI = {E IX
0 0 B I XY
0
T
EIXY }
EIY
ET = {ETX
(7)
(8)
T
ETY }
in which B IX = B IY = [1 Φ1 Φ 2
Φ3]
B IXY = [1 Φ1 Φ 2 1 Φ1 Φ 2 ∂Φ BTX = B TY = 1 ∂z ∂u EIX = 0 ∂x ∂v EIY = 0 ∂y
∂u EIXY = 0 ∂y
∂Φ 2 ∂z
1+
∂u z ∂x
∂ 2 w0 ∂x 2
∂v1 ∂y
∂vz ∂y
∂ 2 w0 ∂y 2
∂u z ∂y
∂v0 ∂x
(9)
∂Φ 3 ∂z
∂u1 ∂x
∂u1 ∂y
2Φ 3 ]
(10) ∂v1 ∂x
∂vz ∂x
6
∂ 2 w0 ∂x∂y
ETX = u1 uz
∂w0 ∂x
ETY = v1 v z
∂w0 ∂y
Stress vector at the kth ply of the cross-ply composite plate with reference to a common structural axis system can be written as σ k = Qk ε k
where σ k = {σ xk
(11) T
T
σ yk τ xyk τ xzk τ yzk } , ε k = {ε xk ε yk γ xyk γ xzk γ yzk } .
Q11k Q12k k k Q21 Q22 Qk = 0 0 0 0 0 0
0 0
0 0
Q33k 0
0 k Q44
0
0
0 0 0 0 Q55k
k
(12)
where Qijk is the transformed elastic stiffness coefficients at the kth ply referred to the (x, y, z) coordinate system.
3. Finite element formulation In equations (6) and (10), it is found that strain components of the present model possess the first and the second derivatives of transverse displacement w, so the C0 and the C1 continuity transverse displacement functions are required. 3.1 Transverse displacement function satisfying C0 continuity in the inter-element The nine-parameter nonconforming element BCIZ [39] satisfying the C0 continuity condition can be given by
w = Fq
(13)
where
F = {F1
F2
F3 }
q = {q 1
q2
q 3}
(14)
T
in which Fi = Fi , Fxi , Fyi
(15)
T
q i = {wi , wxi , w yi }
7
where Fi = Li + L2i L j + L2i Lk − Li L2j − Li L2k ,
(
)
(
)
Fxi = c k L2i L j − c j L2i Lk + c k − c j Li L j Lk / 2 , Fyi = b j L2i Lk − bk L2i L j + b j − bk Li L j Lk / 2 , Li =
a i + bi x + ci y , ai = x j y k − xk y j , 2∆
bi = y j − y k , c i = x k − x j , (i = 1 ~ 3) . The second order derivatives of transverse displacement w can be written as ∂2w 2 ∂x
T
∂2w = Bq ∂x∂y
∂2w ∂y 2
where B = [ B1 B 2
(16)
B3 ] .
∂ 2 Fi 2 ∂x ∂ 2 Fi Bi = 2 ∂y ∂2F i ∂ x ∂ y
∂ 2 Fxi ∂x 2 ∂ 2 Fxi ∂y 2 ∂ 2 Fxi ∂x∂ y
∂ 2 F yi ∂x 2 ∂ 2 F yi , (i = 1 ~ 3) . ∂y 2 ∂ 2 F yi ∂ x ∂ y
3.2 Transverse displacement function satisfying C1 weak -continuity Following the refined element method [40], the displacement function wˆ satisfies the C1 weak-continuity condition can be given by ˆ wˆ = w + P ( Bc − B0 ) q = Fq
(17)
where w is the transverse displacement function satisfying the C0 continuity condition which B0 =
is
presented
in equation
(13);
P = x 2 / 2
y 2 / 2 xy ;
1 B dxdy , and ∆ is the area of element. ∆∫
In equation (17), it is seen that the matrix B c needs to be given to obtain the shape function Fˆ of element transverse displacement. In terms of the refined element method [40], the matrix B c can be expressed as
B c = [B c1
Bc2
B c3 ]
(18) 8
where l1m1 − l3m3 Bc1 = l3m3 −l1m1 m12 − m32
(l12 y21 + l32 y13 ) 1 (m12 y21 + m32 y13 ) 2 1 2
1 2
(l12 x12 + l32 x31)
(l12 x12 + l32 x31 ) 1 (m12 x12 + m32 x31 ) 2 2 2 1 2 (m1 y21 + m3 y13 ) 1 2
(19)
in which li and mi are cosines of the vector normal to the ith boundary. xij=xi-xj, yij=yi-yj, where xi and yi are the coordinates of node i. Matrices B c 2 and Bc 3 can be respectively presented by the permutation of the subscript. Following the equation (17), the shape function Fˆ can be given by Fˆ = F + P ( B c − B 0 )
(20)
3.3 Displacement functions of three-node triangular element
The primary displacement parameters in triangular element can be expressed in terms of nodal variables and shape functions as follows 3
3
3
i =1
i =1
i =1
3
3
3
i =1
i =1
i =1
u0 = ∑ Li u0i , u1 = ∑ Li u1i , u z = ∑ Li uzi
v0 = ∑ Li v0i , v1 = ∑ Li v1i , vz = ∑ Li vzi
(21)
3
w=
∑ (F w i
0 ,i
+ F xi w 0 , xi + F yi w 0 , yi )
0 ,i
+ Fˆ xi w 0 , xi + Fˆ yi w 0 , yi )
i =1 3
wˆ =
∑ ( Fˆ w i
i =1
where Li is area coordinate; the expressions of Fi , Fxi and Fyi can be found in equation (15); the expressions of Fˆi , Fˆxi and Fˆyi may refer to equation (20). 3.4 The strain matrix
In terms of the linear strain-displacement relationships, the strain for the kth layer can be written as follows ε = ∂u = Bδ e = [ B1 B 2
B3 ] δ e
(22)
T
where δe = δ1e δ2e δ3e , δei are displacement variables at the ith node; uk = {u v
T
w} . The nine displacement variables at each node can be written as
9
δie = u0i
v0 i
∂ ∂x ∂= 0 0
w0i
u1i
u zi
w0, xi
∂ ∂z
∂ ∂y
∂ ∂y ∂ ∂x
0
0
0 ∂ ∂z ∂ ∂y
0
∂Li ∂x 0 ∂L Bi = i ∂y 0 0
0 ∂Li ∂y ∂Li ∂x
0 ∂ ∂x
Φ1
B1,3 B2,3
B4,3
0
B5,3
vzi
w0, yi .
T
(23)
∂Li ∂x
Φ2
0 ∂Li ∂y ∂Φ1 Li ∂z
Φ1
B3,3
0
v1i
0
∂Li ∂x
B1,6
0
B2,6
∂Li ∂y ∂Φ 2 Li ∂z
Φ2
0
B3,6
0 ∂Li ∂y ∂L Φ1 i ∂x
∂Li ∂y ∂L Φ2 i ∂x
Φ1
B4,6 B5,6
0 Φ2
0 Li
0
∂Φ1 ∂z
Li
∂Φ 2 ∂z
B1,9 B2,9 B3,9 B4,9 B5,9
(24)
where
B1,3 = Φ 3
∂ 2 Fˆi ; ∂x 2
B1,6 = Φ 3
∂ 2 Fˆxi ; ∂x 2
B1,9 = Φ 3
B2,3 = Φ 3
∂ 2 Fˆi ; ∂y 2
B2,6 = Φ 3
∂ 2 Fˆxi ; ∂y 2
B2,9 = Φ 3
∂ 2 Fˆi B3,3 = 2Φ 3 ; ∂x∂y
B3,6
∂ 2 Fˆxi = 2Φ 3 ; ∂x∂y
∂ 2 Fˆyi
∂x 2
∂ 2 Fˆyi
∂y 2
B3,9 = 2Φ3
;
;
∂ 2 Fˆyi
∂x∂y
;
B4,3 = (1 +
∂Φ 3 ∂Fi ) ; ∂z ∂x
B4,6 = (1 +
∂Φ 3 ∂Fxi ) ; ∂z ∂x
B4,9 = (1 +
∂Φ 3 ∂Fyi ) ; ∂z ∂x
B5,3 = (1 +
∂Φ 3 ∂Fi ) ; ∂z ∂y
B5,6 = (1 +
∂Φ 3 ∂Fxi ) ; ∂z ∂y
B5,9 = (1 +
∂Φ 3 ∂Fyi ) . ∂z ∂y
where i=1, 2 and 3. Based on the global higher-order models, transverse shear stresses obtained from constitutive equations are not continuous at interfaces. Thus, three-dimensional equilibrium equation is required to improve the accuracy of transverse shear stresses, which is usually called as the postprocessing method based on three-dimensional equilibrium equations [41,42]. 10
Neglecting the body forces, the 3-D equilibrium equation at kth ply can be expressed as follows k ∂σ xk ∂τ xy ∂τ xzk + + =0 ∂x ∂y ∂z
∂τ xyk ∂x
+
∂σ yk ∂y
+
∂τ yzk ∂z
(25)
=0
By integrating equation (25) along the z-coordinate direction and enforcing the transverse shear free condition at the bottom surface, the transverse shear stress components satisfying continuity conditions at interfaces can be expressed as z
τˆxzk ( x, y, z ) = − ∫ ∂σ xk / ∂x + ∂τ xyk / ∂y dz z1
(26)
z
τˆyzk ( x, y , z ) = − ∫ ∂τ xyk / ∂x + ∂σ yk / ∂y dz z1
As the in-plane stresses expressed using displacement parameters of the model HZT are substituted in equation (26), a number of higher-order partial derivatives of displacement parameters are to be involved in the expression of transverse shear stresses. In order to obtain simple expression of transverse shear stresses, Iurlaro et al. [43] suggested that transverse shear stresses are expressed using two separate states of cylindrical bending in equation (26). The same approach is to be adopted in the present work, so the expression of transverse shear stresses can be written as
τˆxzk ( x, y , z ) = −Q zx U x
(27)
τˆyzk ( x, y , z ) = −Q zy U y where
{∫ Q dz, ∫ Φ Q dz, ∫ Φ Q dz, ∫ Φ Q dz} = {∫ Q dz, ∫ Φ Q dz, ∫ Φ Q dz, ∫ Φ Q dz}
Q zx = Q zy
z
z1 z
z1
k 11
k 22
z
z1
1
k 11
1
k 22
z
z1
z
2
z1
z
2
z1
z
k 11
k 11
3
z1
z
k 22
z1
T
3
k 22
(28)
∂ 2u ∂ 2u ∂ 2u ∂ 3w0 U x = 20 , 21 , 2z , ∂x 3 ∂x ∂x ∂x
T
∂ 2v ∂ 2v ∂ 2v ∂ 3w0 U y = 20 , 21 , 2z , ∂y 3 ∂y ∂y ∂y
Using the transverse shear free conditions on the top surface of plate (z=zn+1), the following equation can be given by
11
Q zx U x = 0
(29)
Q zy U y = 0
in which
{∫ = {∫
Q zx = Q zy
z n +1
z1 z n +1
z1
Q11k dz , ∫
z n +1
Q dz , ∫
zn +1
z1
k 22
z1
Φ1Q11k dz, ∫
z n +1
z1
Φ1Q dz, ∫ k 22
Φ 2Q11k dz , ∫
z n +1
z1
zn +1
z1
Φ 2Q dz , ∫ k 22
zn +1
z1
} Φ Q dz}
Φ 3Q11k dz 3
(30)
k 22
By means of the equation (30), the second derivative of displacement parameters u0 and v0 can be expressed as follows ∂ 2u0 ˆ ˆ = Q zx U x ∂x 2 ∂ 2v0 ˆ ˆ = Q zy U y ∂y 2
(31)
in which ˆ = {−Q ( 2 ) / Q (1) , − Q ( 3) / Q (1) , − Q ( 4 ) / Q (1)} Q zx zx zx zx zx zx zx ˆ = {−Q ( 2 ) / Q (1) , − Q ( 3) / Q (1) , − Q ( 4 ) / Q (1)} Q zy zy zy zy zy zy zy T
2 2 3 ˆ = ∂ u1 , ∂ u z , ∂ w0 U x 2 2 ∂x 3 ∂x ∂x
(32)
T
2 2 3 ˆ = ∂ v1 , ∂ vz , ∂ w0 U y 2 2 ∂y 3 ∂y ∂y
Inserting equation (31) in equation (27), transverse shear stress components can be rewritten as
τˆxzk ( x, y, z ) = Q zx Uˆ x
(33)
ˆy τˆyzk ( x, y , z ) = Q zy U
in which Q zα (1) = − Q zα ( 2 ) = − Q zα ( 3) = −
(Qˆ (Qˆ (Qˆ
zα
(1) Q zα (1) +Q zα ( 2 ) )
zα
( 2 ) Q zα (1) +Q zα ( 3) )
zα
( 3) Q zα (1) +Q zα ( 4 ) )
(34)
where, α is x or y. In terms of the nodal variables and shape functions of three-node triangular plate element, the equation (33) can be rewritten as
12
k τˆxz ˆ e k = Bδ τˆyz
(35)
where ˆ Bˆ = Bˆ 1 B 2
ˆ B 3
0 0 Bˆ 13 Bˆ 14 Bˆ i = ˆ 0 0 B 23 0
(36) Bˆ 15
0
ˆ B 16 ˆ B
26
0 ˆ B
27
0 ˆ B
28
ˆ B 19 ˆB 29
(37)
where 3 2 2 3 ˆ (3) ∂ Fi , Bˆ = Q ˆ (1) ∂ Li , Bˆ = Q ˆ (2) ∂ Li , Bˆ = Q ˆ (3) ∂ Fxi , Bˆ 13 = Q zx 14 zx 15 zx 16 zx ∂x 3 ∂x 2 ∂x 2 ∂x 3 3
ˆ (3) ∂ Fyi ; Bˆ 17 = Q zx ∂x 3 3 3 2 2 ˆ (3) ∂ Fi , Bˆ = Q ˆ (3) ∂ Fxi , Bˆ = Q ˆ (1) ∂ Li , Bˆ = Q ˆ (2) ∂ Li , Bˆ 23 = Q zy 26 zy 27 zy 28 zy ∂y 3 ∂y 3 ∂y 2 ∂y 2 3
ˆ (3) ∂ Fyi . Bˆ 29 = Q zy ∂y 3 The area coordinates Li are linear functions within one element, so that the second derivatives of area coordinates Li are zero. Equation (37) can be rewritten as ˆ 0 0 Bˆ 13 0 0 B 16 Bˆ i = ˆ ˆ 0 0 B 23 0 0 B 26
ˆ 0 0 B 19 ˆ 0 0 B 29
(38)
From equations (35) to (38), it is observed that it is difficult to produce accurately transverse shear stresses using the equilibrium equation based on the three-node triangular element, as only displacement parameters w0, w0,x and w0,y take part in computation of transverse shear stresses. In order to obtain accurately transverse shear stresses based on the equilibrium equation approach, the higher-order derivatives of displacement parameters in the equilibrium equation have to be eliminated by applying the three-field Hu-Washizu (HW) variational principle [35,36].
3.5 Transverse shear stresses obtained from the HW variational principle Transverse shear stress components τˆxzk and τˆyzk in equation (33) can satisfy the interlaminar continuity conditions, so that the transverse shear stress components 13
τ xzk and τ yzk obtained from constitutive equations will be respectively replaced by the τˆxzk and τˆyzk . Thus, the corresponding strain vector can be given by T
ε = {ε x
ε y γ xy γˆxz γˆ yz }
where γˆxz =
(39)
1 1 τˆxz and γˆ yz = k τˆyz . k Q44 Q55
By means of the HW variational principle discarding body forces, the functional suitable for laminated composite plates [44] can be expressed as follows 1 Π ( u, ε, σ ) = ∫ εT Qε + σ ( ε c ( u ) − ε ) d Ω − Π ext ( u ) 2 Ω
(40)
in which u, ε and σ are independent variables subjected to variation, and Ω is region occupied by a whole plate. T
∂u ∂v ∂u ∂v ∂u ∂w ∂v ∂w εc = , , + , + , + ∂x ∂y ∂y ∂x ∂z ∂x ∂z ∂y
(41)
Π ext = ∫ qu dS
(42)
S
in which q is traction vector. The first variation of the functional presented in equation (40) is used, so that the following equation can be given by
δΠ ( u, ε, σ ) = ∫ ( εT Q − σ ) δ ε dΩ + ∫ δ σ ( ε c − ε ) dΩ + ∫ σδ ε c ( u ) dΩ − δΠ ext ( u ) Ω
Ω
Ω
(43)
=0 By employing integration by parts, the following equation can be given by
∫ σδ ε ( u ) dΩ − δΠ ( u ) = ∫ Dσδ udΩ + ∫ D σδ udS ext
c
(44)
S
Ω
Ω
S
where D and DS are differential operators. By substituting the equation (44) into the equation (43), the equation (43) can be rewritten as
δΠ ( u, ε, σ ) = ∫ ( εT Q − σ ) δ ε dΩ + ∫ δ σ ( ε c − ε ) dΩ + ∫ Dσδ udΩ + ∫ DS σδ udS Ω
Ω
Ω
S
(45)
=0 On the basis of the first integral in equation (45), the following equation can be
14
written as σ = εT Q = {σ I
σT }
(46)
in which σ I = {σ x , σ y , τ xy } , and σ T = { τˆxz , τˆyz } . Inserting the equations (46) into the second integral of equation (45), the following equation can be given by
∫ δ σ (B E I
I
I
− ε I ) dΩ = 0
Ω
(47)
∫ δ σ T ( BT ET − εˆ T ) d Ω = 0 Ω
T
T
in which ε I = {ε x , ε y , γ xy } , and εˆ T = {γˆxz , γˆ yz } . Inserting the equation (33) in the second equation of (47), the higher-order ˆ and U ˆ can be written as follows derivatives of displacement parameters U x y
ˆ =Q E U IX TX x ˆ =Q E U IY TY y
(48)
where
zn +1 1 QIX = ∫ k QTzxQ zx dz z 1 Q44
−1
zn +1 1 QIY = ∫ k QTzy Q zy dz z 1 Q55
−1
∫
zn +1
∫
zn +1
z1
QTzx B TX dz
(49) z1
QTzy B TY dz
Substituting equation (48) in equation (33), the final expressions of transverse shear stress components at kth ply can be expressed as follows
τˆxzk ( x, y , z ) = ΦTX ETX
(50)
τˆyzk ( x, y , z ) = ΦTY ETY where ΦTX = Q zxQ IX and ΦTY = Q zy QIY .
Equation (50) shows that the higher-order derivatives of displacement parameters have been taken out from transverse shear stress components. 3.6 The stiffness matrix of triangular plate element Following the equations (44) and (46), the stiffness matrix of three-node triangular plate element can be presented. Equation (44) within ith element can be 15
written as follows
∫ σδ ε ( u ) d Ω − δΠ ( u ) = c
ext i
i
e h /2
∫ ∫ {σ Ai
− h /2
k x
k y
k xy
k xz
, σ ,τ ,τˆ ,τˆ
k yz
}δ {ε
k x
k y
k xy
k xz
, ε , δγ , δγ , δγ
k T yz
}
(51) dzdxdy + δ Wi = 0
where Ai is area the ith element. Inserting equations (21), (47) and (50) in equation (51), the element stiffness matrix K e within ith element can be readily written as follows Ke = ∫
Ai
∫
h /2
− h /2
BT QBdzdxdy
(52)
where strain matrix B has been defined in equation (22), and strain matrix B can be written as follows
B = B1 B 2
B 3
(53)
in which ∂Li ∂x 0 Bi = ∂Li ∂y 0 0
0
∂Li ∂y ∂Li ∂x 0 0
B1,3
Φ1
∂Li ∂x
Φ2
∂Li ∂x
B1,6 B2,6
B2,3
0
0
B3,3 B4,3 B5,3
∂L Φ1 i ∂y B4,4 0
∂L Φ2 i ∂y B4,5 0
B3,6 B4,6 B5,6
0
∂Li ∂y ∂L Φ1 i ∂x 0 B5,7
Φ1
0
∂Li ∂y ∂L Φ2 i ∂x 0 B5,8
Φ2
B1,9 B2,9 B3,9 B4,9 B5,9
where ∂F B4,3 = Skx Φ TX ( 3) i , B4,4 = Skx Φ TX (1) Li , B4,5 = Skx Φ TX ( 2 ) Li , ∂x
∂F ∂F B4,6 = SkxΦ TX ( 3) xi , B4,9 = Skx Φ TX ( 3) yi ; ∂x ∂x
∂F ∂F B5,3 = Sky Φ TY ( 3) i , B5,6 = Sky Φ TY ( 3) xi , B5,7 = Sky Φ TY (1) Li , ∂y ∂y ∂F B5,8 = Sky Φ TY ( 2 ) Li , B5,9 = Sky Φ TY ( 3) yi . ∂y
16
(54)
where Skx =
1 1 , Sky = k , and the ΦTX ( i ) , ΦTY ( i ) have been defined in k Q44 Q55
equation (50), and i=1, 2, 3. According to the following equation, the vector of nodal displacement δe can be given by
∑K δ
e e
=P
(55)
in which P is loading vector which can be presented by δ Wi in equation (51). The results are obtained from the refined three-node triangular element based on the HW variational theorem, which be called as RTTHW in the present work. For the sake of comparison, the finite element results are also computed based on the principle of minimum potential energy (PMP), namely RTTPMP.
4. Numerical examples and discussion To test the performance of the triangular element RTTHW and its range of applicability, the bending of multilayered composite plates is studied by solving numerically a large number of problems covering a wide range of lamination, geometric parameters and boundary conditions of the plates. The three-dimensional elasticity solutions of laminated composite plates presented by Pagano [1,2] are employed as the benchmark solutions. In addition, the quasi-3D solution obtained from the nine-order model HSDT-98 proposed by Matsunaga [31] are employed to assess the performance of the present element. The in-plane displacement field of the model HSDT-98 consists of 9th-order polynomial in the global thickness coordinate z and the transverse deflection is denoted by an 8th-order polynomial in z. On the other hand, a comparison with the RZTM model [43] is presented for all examples to evaluate the improvement of the proposed triangular plate element. Owning to the symmetry of geometry, material properties and loading conditions, only one-quarter of the laminated composite plate has been used. The mesh configuration can be found in Figure 2. The acronyms used in all Figures and Tables will be simply explained as follows Exact:
Exact three-dimensional elasticity solution presented by Pagano [1].
Quasi-3D: The quasi-3D solution is obtained from the nine-order model HSDT-98 proposed by Matsunaga [31]. 17
RTTHW: The present result is obtained from the refined three-node triangular
element based on the HW variational theorem. RZTM:
Based on the Reissner's mixed variational principle, analytical result is obtained from the refined zigzag theory (RZT) [43].
RHSDT: Analytical result is obtained from Reddy's higher-order shear deformation
theory [19]. FSDT:
Analytical result is obtained from the first-order shear deformation theory (FSDT).
Example 1 ：Cylindrical bending of the laminated composite plate strip subjected to a
sinusoidal transverse loading q=q0sin(π x/a) is analyzed. The ply material properties of all laminated composite plates [1,2] are given by EL=172.5GPa, ET=6.9Gpa, GLT= 3.45GPa, GTT=1.38GPa, vLT=0.25, vTT=0.25 where L and T denote the direction parallel to the fiber and the transverse direction for orthotropic materials. The dimensionless parameters of the displacements and stresses are expressed as follows u=
ET u ( 0, z ) 100ET h3 w ( a / 2, z ) σ ( a / 2, z ) τ ( 0, z ) , w= , σx = x , τ xz = xz . 4 q0 h q0 a q0 q0
For all examples, a and b are respectively the length and width of laminated plates, h is the thickness of lamination. Figure 3 shows the convergence rate of the transverse displacements with mesh refinement for 11-ply plate (a/h=4). It is seen that the present results can rapidly converge to three-dimensional elasticity solutions which are calculated by the present authors following Pagano’s methods [1]. Numerical results show that the present results with 16×16 mesh density can converge to the exact solutions. In Figure 4, it is found that the present element is free from the shear locking problem even in the analysis of the 11-ply plate with a length-to-thickness ratio (a/h) of 1000. In Table 1, displacements and stresses at the strategic point of composite plates obtained from different models are compared with three-dimensional elasticity solution (Exact) [1]. It ought to be shown that the RTTPMP and RZT in Table 1 are the results obtained directly from the principle of minimum potential energy without using the transverse shear stresses τˆxzk and τˆyzk . 18
Nevertheless, RTTHW and RZTM are the results obtained respectively from the Hu-Washizu (HW) variational principle and the mixed Reissner's variational principle considering the effects of the transverse shear stresses τˆxzk and τˆyzk . Numerical results show that the maximum percentage error of transverse displacement of eight-ply plate rapidly decreases from 15.50 to 1.41 for the model RZT as effects of the transverse shear stresses τˆxzk and τˆyzk are considered. The maximum percentage error of in-plane stress of eight-ply plate decreases from 15.29 to 9.82 for the model RZT, so the transverse shear stresses τˆxzk and τˆyzk has significant impact on displacements and stresses of the model RZT. Nevertheless, effects of the transverse shear stresses τˆxzk and τˆyzk on displacements and stresses of the model RTTHW are slight. As the number of layers is increased to 50, the maximum percentage errors of in-plane displacement and in-plane stress respectively decrease from 24.36 and 24.26 to 23.71 and 23.62 for the model RZT, which show that the transverse shear stresses
τˆxzk and τˆyzk has less effects on displacements and in-plane stresses of multilayered composite plates for the model RZT. On the other hand, the transverse shear stresses
τˆxzk and τˆyzk still impact actively the accuracy of displacements and stresses obtained from the RTTHW. In Figures 5-7, distributions of displacements and stresses through the thickness of symmetric composite plates obtained from different models are compared with three-dimensional elasticity solutions. It is found that the present results RTTHW for the displacements and the stresses are in excellent agreement with those of three-dimensional elasticity theory. The model RZTM can describe a zig-zag shape distribution of in-plane displacement through the thickness direction of composite plate whereas in-plane displacement obtained from the model RZTM is obviously diverge from the three-dimensional elasticity solutions. In addition, the transverse shear stresses obtained from the model RZTM satisfy the continuity conditions at the interfaces, but the model RZTM overestimates transverse shear stresses in comparison with the exact solutions. It is also observed that Reddy's higher-order shear deformation theory (RHSDT) is unable to describe a zig-zag shape distribution of in-plane displacements along the thickness, and interlaminar continuity of transverse shear stresses is also unable to be satisfied. For a six-ply antisymmetric plate, the displacement and the stresses obtained from the RTTHW are compared with 19
three-dimensional elasticity solutions and results obtained from other models in Figure 8. Again it is observed that the present RTTHW can accurately predict displacement and stresses. Nevertheless, the in-plane displacement on the top surface obtained from the model RZTM departs noticeably from the exact solution. Example 2: Simply-supported multilayered composite plates subjected to a doubly
sinusoidal transverse loading q=q0sin( π x/a)sin( π y/b) are analyzed. The displacements and stresses are normalized as follows (σ x , σ y )(0,0, z )h 2 τ xy (a / 2, b / 2, z )h 2 E2 h 2u (a / 2,0, z ) , , , u= (σ x , σ y ) = τ xy = a3 q0 a 2 q0 a 2
τ xz =
τ (0, b / 2, z )h τ xz (a / 2,0, z )h , τ yz = yz q0 a q0 a
Numerical results of in-plane stresses and transverse stresses are given in Table 2. It is found that the maximum percentage error of the RTTHW is 5.38 whereas the maximum percentage error of the RZTM is 34.72. Therefore, the results obtained using the RTTHW are in close agreement with three-dimensional elasticity solution [2]. Figure 9 shows the through the thickness variations of displacements and stresses for a simply-supported six-ply square plate. It shows that the results obtained from RTTHW are in good agreement with the three-dimensional solutions whereas the values predicted by the RZTM differ from three-dimensional solutions. Figure 10 shows the through the thickness variation of in-plane displacement and stresses for an eight-ply square plate. It shows that the displacement and stress values obtained from the RTTHW are in excellent agreement with exact solutions [45]. To test the range of applicability of the RTTHW, the results of parametric studies for multilayered rectangular plates are presented. These researches were conducted to provide some insight into the influences of variation in the lamination and geometric parameters of composite plates on the accuracy of the RTTHW and RZTM. Figure 11 shows the distribution of displacements and stresses of moderately thick six-ply rectangular plate (a/h=4, a/b=3), and the results are also compared with quasi-3D solutions [31]. For the displacement and stress distribution, the present results RTTHW are in excellent agreement with those of the quasi-3D theory. In Figure 12, it is observed that the results obtained using the RTTHW are quite acceptable even in the case of a very thick laminated composite plates (a/h=2, a/b=3). Compared with the quasi-3D solutions, the distribution of transverse shear stress satisfies the stress 20
boundary conditions on the upper and the lower surfaces of plate and continuity condition at interfaces. However, results obtained using the model RZTM are less accurate in comparison with the quasi-3D solution for thick composite plates. Figure 13 shows the distribution of displacements and stresses of 15-ply rectangular plate (a/h=4, a/b=3). A close agreement between the RTTHW and the quasi-3D solution is noticed, which shows that accuracy of the RTTHW is not affected by the number of layers in laminate.
5. Conclusions Based on the higher-order zig-zag model (HZT) satisfying the free surface conditions, a refined three-node triangular plate element (RTTHW) has been presented. In order to fulfill the continuity condition of transverse shear stresses at interfaces, transverse shear stresses are usually computed by integrating three-dimensional equilibrium equation through the thickness of laminates. Once integration of three-dimensional equilibrium equation is employed, the second and the third derivatives of displacement parameters will be involved into three-dimensional equilibrium equation. For the three-node triangular element, displacement parameters are linear variation within one element, so that the second and the third derivatives of displacement parameters are close to zero. Thus, three-dimensional equilibrium equation will encounter difficulty to accurately predict transverse shear stresses for the three-node triangular element. To obtain improved transverse shear stresses, the simplified three-dimensional equilibrium equation has been a priori used. It is significant that the higher-order derivatives of displacement parameters have been taken out from the expression of transverse shear stresses by employing the three-field HW variational principle. Applying the present element RTTHW in conjunction with other published models to investigate the typical bending problems of multilayered composite plates, the following conclusions can be drawn: (1) The proposed model HZT can model a zig-zag shape distribution of in-plane displacement along the thickness. Moreover, the present element RTTHW based on the model HZT can predict the more accurate in-plane displacement in comparison with the model RZTM. (2) Based on the element RTTHW, transverse shear stresses can be accurately computed without using any postprocessing approaches. Moreover, accurate transverse shear stresses have active impact on accuracy of displacements and in-plane stresses. (3) Even in the case of a very thick laminated composite plates 21
(a/h=2), the results obtained using the element RTTHW are quite acceptable. Moreover, the accuracy of the RTTHW is not affected by the number of layers and aspect ratio (a/b) in laminate.
Acknowledgement The work described in this paper was supported by the National Natural Sciences Foundation of China [No. 11272217, 11402152, 11572204].
References [1] Pagano NJ, Exact solutions for composite laminates in cylindrical bending. J Compos Mater, 1969; 3: 398-411. [2] Pagano NJ, Exact solutions for rectangular bidirectional composites. J Compos Mater, 1970; 4:20-34. [3] Bogdanovich AE, Yushanov SP, Three-dimensional variational analysis of Pagano's problems for laminated composite plates. Compos Sci Tech, 2000; 60:2407-2425. [4] Mendonca PTR, Barcellos CS, Torres DAF, Robust Ck/C0 generalized FEM approximations for higher-order conformity requirement: Application to Reddy's HSDT model for anisotropic laminated plates. Compos Struct, 2013; 96:332-345. [5] Matsunaga H, A comparison between 2-D single-layer and 3-D layerwise theories for computing interlaminar stresses of laminated composite and sandwich plates subjected to thermal loadings. Compos Struct, 2004; 64:161-177. [6] Whitney JM, The effect of transverse shear deformation on the ending of laminated plates. J Compos Mat, 1969; 3:534-847. [7] Rolfes R, Rohwer K, Improved transverse shear stresses in composite finite elements based on first-order shear deformation theory. Int J Num Meth Eng, 1997; 40:51-60. [8] Liew KM, He XQ, Tan MJ, Lim HK, Dynamic analysis of laminated composite plates with piezoelectric sensor/actuator patches using the FSDT mesh-free method. Int J Mech Sci, 2004; 46:411-431. [9] Noor AK, Malik M, An assessment of five modeling approaches for thermo-mechanical stress analysis of laminated composite panels. Comput Mech, 22
2000; 25:43-58. [10] Kant T, Swaminathan K, Analytical solutions for free vibration of laminated composite and sandwich plates based on a higher-order refined theory. Compos Struct, 2001; 53:73-85. [11] Matsunaga H, Interlaminar stress analysis of laminated composite and sandwich circular arches subjected to thermal/mechanical loading. Compos Struct, 2003; 60:345-358. [12] Lo KH, Christensen PM, Wu EM, A higher order theory of plate deformation, Part 1: Homogeneous plates. J Appl Mech, 1977; 44:663-668. [13] Lo KH, Christensen PM, Wu EM, A higher order theory of plate deformation, Part 2: laminated plates. J Appl Mech, 1977; 44:669-676. [14] Marur SR, Kant T, Free vibration analysis of fiber reinforced composite beams using higher order theories and finite element modelling. J Sound Vib, 1996; 194(3):337-351. [15] Kant T, Swaminathan K, Free vibration of isotropic, orthotropic, and multilayer plates based on higher order refined theories. J Sound Vib, 2001; 241(2):319-327. [16] Kant T, Swaminathan K, Analytical solutions for the static analysis of laminated composite and sandwich plates based on a higher-order refined theory. Compos Struct, 2002; 56:329-344. [17] Babu CS, Kant T, Two shear deformable finite element models for buckling analysis of skew fibre-reinforced composite and sandwich panels. Compos Struct, 1999; 46:115-124. [18] Kant T, Babu CS, Thermal buckling analysis of skew fibre-reinforced composite and sandwich plates using shear deformable finite element models. Compos Struct, 2000; 49:77-85. [19] Reddy JN, A simple higher-order theory for laminated composite plates. J Appl Mech, 1984; 51(12):745-752. [20] Nayak AK, Moy SSJ, Shenoi RA, Free vibration analysis of composite sandwich plates based on Reddy's higher-order theory. Composites: Part B, 2002; 33:505-519.
23
[21] Nayak AK, Moy SSJ, Shenoi RA, A higher order finite element theory for buckling and vibration analysis of initially stressed composite sandwich plates. J Sound Vib, 2005; 286:763-780. [22] Sheikh AH, Chakrabarti A, A new plate bending element based on higher-order shear deformation theory for the analysis of composite plates. Finite Element in Analysis and Design, 2003; 39:883-903. [23] Cho M, Oh J, Higher order zig-zag theory for fully coupled thermoelectric-mechanical smart composite plates. Int J Solids Struct, 2004; 41:1331-1356. [24] Oh J, Cho M, A finite element based on cubic zig-zag plate theory for the prediction of thermo-electric-mechanical behaviors. Int J Solids Struct, 2004; 41:1357-1375. [25] Icardi U, Eight-noded zig-zag element for deflection and stress analysis of plates with general lay-up. Composite Part B, 1998; 29:425-441. [26] Chakrabarti A, Sheikh AH, Analysis of laminated sandwich plates based on interlaminar shear stress continuous plate theory. J Eng Mech, 2004; 131(4):377-384. [27] Versino D, Gherlone M, Mattone M, Di Sciuva M, Tessler A, C0 triangular elements based on the refined zigzag theory for multilayer composite and sandwich plates. Composites: Part B, 2013; 44:218-230. [28] Carrera E, Demasi L, Classical and advanced multilayered plate elements based upon PVD and RMVT. Part 1: Derivation of finite element matrices. Int J Numer Meth Engng, 2002; 55:191-231. [29] Carrera E, Demasi L, Classical and advanced multilayered plate elements based upon PVD and RMVT. Part 2: Numerical implementations. Int J Numer Meth Engng, 2002; 55:253-291. [30] Wu CP, Chen WY, Vibration and stability of laminated plates based on a local higher order plate theory. J Sound Vib, 1994; 177(4):503-520. [31] Matsunaga H, Assessment of a global higher-order deformation theory for laminated composite and sandwich plates. Compos Struct, 2002; 56:279-291. [32] Carrera E. Historical review of Zig-Zag theories for multilayered plates and 24
shells, Appl Mech Rev, 2003; 56:287-308. [33] Lekhnitskii SG, Strength calculation of composite beams. Vestnik inzhen i tekhnikov, 1935, No. 9. [34] Ambartsumian SA, On a theory of bending of anisotropic plates. Investiia Akad Nauk SSSR, Ot Tekh Nauk, 1958, No. 4. [35] Hu HC, On some variational principle in the theory of elasticity and plasticity. Scientia Sinica Beijing, 1955; 4:33-54. [36] Washizu K, On some variational principle in the theory of elasticity and plasticity. Technical Report 25-18, Aeroelastic and Structures Research Laboratory, Massachusetts Institute of Technology, March 1955. [37] Carrera E, On the use of the Murakami’s zig-zag function in the modeling of layered plates and shells. Comput Struct, 2004; 82:541-554. [38] Murakami H, Laminated composite plate theory with improved in-plane response. J Appl Mech, 1986; 53:661-666. [39] Bazeley GP, Cheung YK, Irons BM, Zienkiewiz OC, Triangular elements in bending conforming and non-conforming solution, Proc. Conf. Matrix Methods in Structural Mechanics, Air Force Ins. Tech., Wright-Patterson A.F. Base, Ohio,
1965; 547-576. [40] Cheung YK, Chen WJ, Refined Nine-parameter Triangular Thin Plate Bending Element by Using Refined Direct Stiffness Method, Int J Numer Methods Eng, 1995; 38:283-298. [41] Cho M, Kim JS, Four-noded finite element post-process method using a displacement field of higher-order laminated composite plate theory. Compos Struct, 1996; 61(2):283-290. [42] Oh J, Cho M, Higher-order zig-zag theory for smart composite shells under mechanical-thermo-electric loading. Int J Solids Struct, 2007; 44:100-127. [43] Iurlaro L, Gherlone M, Di Sciuva M, Tessler A, Refined Zigzag Theory for laminated composite and sandwich plates derived from Reissner's mixed variational Theorem. Compos Struct, 2015; 133:809-817. [44] Rah K, Van Paepegem W, Habraken AM, Degrieck, A mixed solid-shell element for the analysis of laminated composites. Int J Numer Meth Engng, 2012;
25
89:805-828. [45] Bogdanovich AE, Deepak BP, Three-dimensional analysis of thick composite plates with multiple layers. Composite Part B, 1997; 28:345-357.
26
q t ( x, y )
z y
z N +1 zN
h/ 2 h/ 2
x
z1
Fig. 1 A laminated composite plate subjected to transverse loading
y
b 2
x a/2
Fig. 2 A quarter of the plate divided into a 16 × 16 mesh
27
3.23 3.22 Exact RTTHW
3.21 3.2
w
3.19 3.18 3.17 3.16 3.15 3.14
12 × 2
4 2× 4
8 ×3 8
12 4× 12
16 ×5 16
24 ×6 24
Mesh density
Fig. 3 Convergence rate of transverse displacement for 11-ply plate (a/h=4)
3.5 Exact RTTHW (16x16)
3
2.5
w
2
1.5
1
0.5
14
2 10
3 20
4 50
5 100
6 500
1000 7
a/h
Fig. 4 Transverse displacements of 11-ply square plates with different length-to-thickness ratios (a/h)
28
0.5
0.5
0.4
Exac t RTTHW (16x16) RZTM RHSDT FSDT
0.3 0.2
z h
0.4 0.3 0.2
0.5 Exact RTTHW (16x16) RZTM RHSDT FSDT
0.4 0.3 0.2
0.1
0.1
0.1
0
0
0
-0.1
-0.1
-0.1
-0.2
-0.2
-0.2
-0.3
-0.3
-0.4 -0.5 -1.5
-1
-0.5
0
-0.3
-0.4
u ( a / 2,0, z )
0.5
1
1.5
-0.5 -30
Exac t RTTHW (16x16) RZTM RHSDT FSDT
-0.4 τ xz ( a / 2, 0, z )
σ x ( 0, 0, z )
-20
-10
0
10
20
30
-0.5 -0.5
0
0.5
1
1. 5
2
2.5
3
Fig. 5 Comparison of displacements and stresses through thickness of 11-ply plate under cylindrical bending (a/h=4) 0.5
0.5
0.4
0.2
z h
0.4
Exact RTTHW (16x16) RZTM FSDT
0.3
0.3 0.2
0.5 Exact RTTHW (16x16) RZTM FSDT
0.4 0.3 0.2
0.1
0.1
0.1
0
0
0
-0.1
-0.1
-0.1
-0.2
-0.2
-0.2
-0.3
-0.3
-0.3 -0.4 -0.5 -1.5
-0.4
u ( a / 2,0, z )
-1
-0.5
0
0.5
1
1.5
-0.5 -30
-0.4
σ x ( 0, 0, z )
-20
-10
0
10
20
30
-0.5 -0.5
Exact RTTHW (16x16) RZTM FSDT
τ xz ( a / 2,0, z )
0
0.5
1
Fig. 6 Comparison of displacements and stresses through thickness of 15-ply plate under cylindrical bending (a/h=4)
29
1.5
2
0.5
0.5
0.4 0.3 0.2
z h
0.4
Exact RTTHW (16x16) RZTM FSDT
0.3 0.2
0.5 Exact RTTHW (16x16) RZTM FSDT
0.4 0.3 0.2
0.1
0.1
0.1
0
0
0
-0.1
-0.1
-0.1
-0.2
-0.2
-0.2
-0.3
-0.3
-0.3 -0.4
-0.4
u ( a / 2,0, z )
-0.5 -1.5
-1
-0.5
0
0.5
1
1.5
-0.5 -30
-0.4
σ x ( 0, 0, z )
-20
-10
0
10
20
30
Exact RTTHW (16x16) RZTM FSDT
τ xz ( a / 2,0, z )
-0.5 -0.5
0
0.5
1
1.5
2
Fig. 7 Comparison of displacements and stresses through thickness of 21-ply plate under cylindrical bending (a/h=4) 0.5
0.5
0.4
Exact RTTHW (16x16) RZTM RHSDT FSDT
0.3 0.2
z h
0.4 0.3 0.2
0.5 0.4
Exact RTTHW (16x16) RZTM RHSDT FSDT
0.2
0.1
0.1
0.1
0
0
0
-0.1
-0.1
-0.1
-0.2
-0.2
-0.2
-0.3
-0.3
-0.3
-0.4
-0.4
-0.5 -2.5
u ( a / 2, 0, z )
-2
-1.5
-1
-0.5
0
0.5
1
1.5
-0.5 -30
-20
-10
Exact RTTHW (16x16) RZTM RHSDT FSDT
0.3
σ x ( 0,0, z )
-0.4
0
-0.5 0
10
20
30
τ xz ( a / 2,0, z )
0.5
1
1.5
2
2.5
Fig. 8 Comparison of displacements and stresses through thickness of six-ply plate under cylindrical bending (a/h=4)
30
3
3.5
0.5
0.5 Quasi-3D RTTHW (16x16) RZTM RHSDT FSDT
0.4 0.3 0.2
z h
0.4 0.3 0.2
0.4
0.2
0.1
0.1
0
0
0
-0.1
-0.1
-0.1
-0.2
-0.2
-0.2
-0.3
-0.3
-0.4
-0.4
u ( a / 2, 0, z )
-0.015
-0.01
-0.005
0
0.005
0.01
0.5 0.4 0.3 0.2
-0.5 -0.8
-0.3
-0.6
-0.4
-0.2
0
0.2
0.4
0.6
Quasi-3D RTTHW (16x16) RZTM RHSDT FSDT
0.4 0.3 0.2
0.1
0.1
0.8
-0.5 -0.8
0
-0.2
-0.2
-0.2
-0.3
-0.3
-0.3
-0.02
-0.01
0
0.01
-0.4 0.02
0.03
0.04
0.2
0.4
0.6
Quasi-3D RTTHW (16x16) RZTM RHSDT FSDT
0.1
-0.1
-0.03
0
0.2
0
-0.5 0
-0.2
0.3
-0.1
τ xy ( a / 2, b / 2, z )
-0.4
0.4
0
-0.5 -0.04
-0.6
0.5
-0.1
-0.4
σ y ( 0, 0, z )
-0.4
σ x ( 0, 0, z )
0.5 Quasi-3D RTTHW (16x16) RZTM RHSDT FSDT
Quasi-3D RTTHW (16x16) RZTM RHSDT FSDT
0.3
0.1
-0.5 -0.02
z h
0.5 Quasi-3D RTTHW (16x16) RZTM RHSDT FSDT
0.05
0.1
0.15
0. 2
0.25
0.3
τ xz ( a / 2,0, z )
-0.4
0.35
-0.5 0
0. 4
τ yz ( 0, b / 2, z ) 0.05
0.1
0.15
0.2
0.25
0.3
0.35
Fig. 9 Comparison of displacements and stresses through thickness of six-ply plate subjected to double sinusoidal load (a/h=4)
31
0.4
0.8
0.5
0.5
0.4 0.3 0.2
z h
0.4
Quasi-3D RTTHW (16x16) RZTM RHSDT FSDT
0.2
0.1
0.1
0
0
-0.1
-0.1
-0.2
-0.2
-0.3
-0.3
-0.4 -0.5 -0.015
-0.4
u ( a / 2,0, z ) -0.01
-0.005
0
0.005
0.01
0.5 0.4 0.3 0.2
z h
σ x (0, 0, z )
-0.5 -0.8
-0.6
-0.4
-0.2
0
0.2
0.4
0.6
0.8
0.5 0.4
Quasi-3D RTTHW (16x16) RZTM RHSDT FSDT
0.3
Exact RTTHW (16x16) RZTM RHSDT FSDT
0.2
0.1
0.1
0
0
-0.1
-0.1
-0.2
-0.2
-0.3
-0.3
τ xy ( a / 2, b / 2, z )
-0.4 -0.5 -0.04
Exact RTTHW (16x16) RZTM RHSDT FSDT
0.3
-0.03
-0.02
-0.01
0
τ xz ( a / 2,0, z )
-0.4
0.01
0.02
0.03
0.04
-0.5
-0.05
0
0.05
0.1
0.15
0.2
0.25
0.3
0.35
Fig. 10 Comparison of displacements and stresses through thickness of eight-ply plate subjected to double sinusoidal load (a/h=4)
0.5
0.5
0.4
0.4
0.3 0.2
z h
0.1
0.3
Quasi-3D RTTHW (16x16) RZTM RHSDT FSDT
0.2 0.1
0
0
-0.1
-0.1
-0.2
-0.2
-0.3
-0.3
-0.4 -0.5 -0.04
-0.4
u ( a / 2, 0, z ) -0.03
-0.02
Quasi-3D RTTHW (16x16) RZTM RHSDT FSDT
-0.01
0
0.01
0.02
32
-0.5 -1.5
σ x (0, 0, z )
-1
-0.5
0
0.5
1
1.5
0.5
0.5
0.4 0.3 0.2
z h
Quasi-3D RTTHW (16x16) RZTM RHSDT FSDT
0.4
Quasi-3D RTTHW (16x16) RZTM RHSDT FSDT
0.3 0.2
0.1
0.1
0
0
-0.1
-0.1
-0.2
-0.2
-0.3
-0.3
τ xy ( a / 2, b / 2, z )
-0.4 -0.5 -0.04
-0.03
-0.02
-0.01
0
τ xz ( a / 2,0, z )
-0.4
0.01
0.02
0.03
0.04
-0.5 0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
Fig. 11 Comparison of displacements and stresses through thickness of six-ply rectangular plate subjected to double sinusoidal load (a/h=4, b/a=3) 0.5
0.5
0.4
0.4
0.3
Quasi-3D RTTHW (16x16) RZTM RHSDT FSDT
0.2
z h
0.1
0.3 0.2 0.1
0
0
-0.1
-0.1
-0.2
-0.2
-0.3
-0.3
-0.4
-0.4
u ( a / 2, 0, z )
-0.5 -0.08
-0.06
-0.04
-0.02
0
0.02
0.04
0.5
-0.5 -0.6
σ y ( 0,0, z ) -0.4
-0.2
0
0.2
0.4
0.6
0.8
0.5
0.4
0.4
Quasi-3D RTTHW (16x16) RZTM RHSDT FSDT
0.3 0.2
z h
Quasi-3D RTTHW (16x16) RZTM RHSDT FSDT
0.3 0.2
0.1
0.1
0
0
-0.1
-0.1
-0.2
-0.2
-0.3
-0.3
τ xy ( a / 2, b / 2, z )
-0.4 -0.5 -0.08
-0.06
-0.04
-0.02
0
0.02
-0.4 0.04
0.06
0.08
-0.5 -0.2
Quasi-3D RTTHW (16x16) RZTM RHSDT FSDT
τ xz ( a / 2,0, z ) -0.1
0
0.1
0.2
0.3
0.4
0.5
Fig. 12 Comparison of displacements and stresses through thickness of six-ply rectangular thick plate subjected to double sinusoidal load (a/h=2, b/a=3)
33
0.6
0.5
0.5
0.4
0.2
z h
0.4
Quasi-3D RTTHW (16x16) RZTM RHSDT FSDT
0.3
0.2
0.1
0.1
0
0
-0.1
-0.1
-0.2
-0.2
-0.3
-0.3
-0.4
0
0.005
0.01
0.015
0.02
0.5 0.4 0.3 0.2
z h
-0.5 -1.5
-1
-0.5
0
0.5
1
1.5
0.5 0.4
Quasi-3D RTTHW (16x16) RZTM RHSDT FSDT
0.3
Quasi-3D RTTHW (16x16) RZTM RHSDT FSDT
0.2
0.1
0.1
0
0
-0.1
-0.1
-0.2
-0.2
-0.3
-0.3
τ xy ( a / 2, b / 2, z )
-0.4 -0.5 -0.04
σ x (0, 0, z )
-0.4
u ( a / 2, 0, z )
-0.5 -0.02 -0.015 -0.01 -0.005
Quasi-3D RTTHW (16x16) RZTM RHSDT FSDT
0.3
-0.03
-0.02
-0.01
0
0.01
τ xz ( a / 2,0, z )
-0.4 0.02
0.03
0.04
-0.5 -0.1
0
0.1
0.2
0.3
0.4
0.5
0.6
Fig. 13 Comparison of displacements and stresses through thickness of 15-ply rectangular plate subjected to double sinusoidal load (a/h=4, b/a=3)
34
0.7
Table 1 Effect of transverse shear stresses on displacements and in-plane stresses (a/h=4) w ( a / 2, 0 )
u ( 0, − h / 2 )
σ x ( a / 2, − h / 2 )
Exact
3.7391
1.1524
-22.6552
RTTHW(16×16)
3.7456 (0.17)
1.1930 (3.52)
-23.7562 (4.86)
RTTPMP(16×16) 3.6373 (2.72)
1.1501 (0.20)
-22.9051 (1.10)
RZTM
3.7920 (1.41)
1.0379 (9.94)
-20.4303 (9.82)
RZT
3.1597 (15.50)
0.9749 (15.40)
-19.1911 (15.29)
Exact
3.5812
1.1797
-23.1923
RTTHW(16×16)
3.5752 (0.17)
1.2261 (3.93)
-24.4078 (5.24)
RTTPMP(16×16) 3.4410 (3.91)
1.1204 (5.03)
-22.3111 (3.80)
RZTM
3.6519 (1.97)
0.9865 (16.38)
-19.4185 (16.27)
RZT
3.1606 (11.74)
0.9599 (18.63)
-18.8949 (18.53)
Exact
3.5379
1.2325
-24.2313
RTTHW(16×16)
3.5286 (0.26)
1.2844 (4.21)
-25.5560 (5.47)
RTTPMP(16×16) 3.3591 (5.05)
1.1256 (8.67)
-22.4092 (7.52)
RZTM
3.6166 (2.22)
0.9668 (21.56)
-19.0315 (21.46)
RZT
3.1613 (10.64)
0.9540 (22.60)
-18.7797 (22.50)
Exact
3.5310
1.2591
-24.7538
RTTHW(16×16)
3.5192 (0.33)
1.3129 (4.27)
-26.1186 (5.51)
RTTPMP(16×16) 3.3333 (5.60)
1.1309 (10.18)
-22.5118 (9.06)
RZTM
3.6097 (2.23)
0.9605 (23.71)
-18.9079 (23.62)
RZT
3.1615 (10.46)
0.9524 (24.36)
-18.7473 (24.26)
Number of Models layers 8
16
32
50
The number in bracket is the % errors of various models relative to three-dimensional elasticity solution.
35
Table 2 Comparison of results obtained from various models with three-dimensional elasticity solutions (a/h=4)
σx
τ xy
τ xz
τ yz
(a/2,b/2,h/2)
(0,0,h/2)
(0,b/2,0)
(a/2,0,0)
Exact
0.8010
-0.0511
0.256
0.217
RTTHW(16×16)
0.8175 (2.06)
-0.0525 (2.74)
0.254 (0.78)
0.212 (2.30)
RZTM
0.6873 (14.19)
-0.0488 (4.50)
0.258 (0.78)
0.219 (0.92)
FSDT
0.4369 (45.46)
-0.0369 (27.79)
0.120 (53.13)
0.130 (40.09)
Exact
0.684
-0.0337
0.223
0.223
RTTHW(16×16)
0.702 (2.63)
-0.0342 (1.48)
0.235 (5.38)
0.224 (0.45)
RZTM
0.577 (15.64)
-0.0220 (34.72)
0.257 (15.25)
0.230 (3.14)
FSDT
0.466 (31.87)
-0.0224 (33.53)
0.244 (9.42)
0.083 (62.78)
Number of layers
Models
3
9
The number in bracket is the % errors of various models relative to three-dimensional elasticity solution.
36
S0263-8223(16)32001-3 http://dx.doi.org/10.1016/j.compstruct.2016.11.040 COST 8000
To appear in:
Composite Structures
Received Date: Revised Date: Accepted Date:
29 September 2016 5 November 2016 11 November 2016
Please cite this article as: Zhen, W., Rui, M., Wanji, C., A refined three-node triangular element based on the HW variational theorem for multilayered composite plates, Composite Structures (2016), doi: http://dx.doi.org/10.1016/ j.compstruct.2016.11.040
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.
A refined three-node triangular element based on the HW variational theorem for multilayered composite plates Wu Zhen*, Ma Rui, Chen Wanji Key Laboratory of Liaoning Province for Composite Structural Analysis of Aerocraft and Simulation, Shenyang Aerospace University, Shenyang 110136, China
Abstract: For the three-node triangular elements, the displacement parameters are linear variation within one element, so that the second and the third derivatives of displacement parameters are close to zero. Therefore, the three-node triangular element will encounter difficulty to accurately predict transverse shear stresses by integrating the three-dimensional equilibrium equation through the thickness direction. Thus, the few three-node triangular plate elements can accurately yield the transverse shear stresses of multilayered composite plates. A higher-order zig-zag model accounting for the zero transverse shear strain conditions at the top and the bottom surfaces is firstly presented in this work. Based on the proposed model, a refined three-node triangular plate element satisfying the requirement of C1 weak-continuity conditions in the inter-element is developed to accurately produce the distributions of displacements and stresses through the thickness of multilayered composite plates. To obtain the improved transverse shear stresses, the simplified three-dimensional equilibrium equation has been a priori used. It is significant that the higher-order derivatives of displacement parameters have been taken out from the expression of transverse shear stresses by employing the three-field Hu-Washizu (HW) variational principle, which is convenient for the finite element implementation. Performance of the present element is tested by comparing with three-dimensional elasticity solutions as well as the reference results evaluated using other models. A detailed study is conducted to highlight the influences of boundary condition, number of layers and aspect ratio on the static response of multilayered composite plates, and numerical results can show the accuracy and range of applicability for the proposed element. Keywords: Three-node triangular element; Higher-order zig-zag model; Three-field Hu-Washizu variational principle; Multilayered composite plates; Transverse shear stress.
†
E-mail address: [email protected] (Wu Z.) 1
1. Introduction Due to their high strength-to-weight ratio, stiffness-to-weight ratio and superior fatigue characteristics, the fiber-reinforced composite materials are widely used as the fundamental members of aerospace, such as parts of wing, fuselage as well as satellite body etc. In order to apply the advantage of the high specific strength of such materials, it needs to accurately predict structural behaviors of these laminates. In the early stages of the development of models, the three-dimensional models [1,2] are usually employed to capture the significant stresses arising at the interfaces of the laminated
composites.
Bogdanovich and
Yushanov [3] also developed
a
three-dimensional displacement-assumed variational-analysis approach for bending analysis of laminated composite plates. Three-dimensional analysis can serve as benchmark solutions to assess the performance of approximate laminated plate models, whereas three-dimensional analysis is limited because it needs the computationally expensive cost as some layers much thinner compared to others [4]. Since three-dimensional analytical solutions for laminated composite structures are limited, the two-dimensional plate models have been developed to study the practical problems by using finite element methods and some analytical procedure in very simple cases. Based on Kirchhoff's hypothesis, the classical laminated plate theory is employed for the analysis of laminated composites. Owning to neglect of transverse shear and normal strains in the laminate, the classical laminated plate theory is not adaptable to analysis of multilayered composite plates with relatively soft transverse shear modulus [5]. In order to consider the influence of transverse shear strain, a number of first-order shear deformation theories [6-9] have been proposed. Nevertheless, experience of many investigators [9-11] has shown that the first-order theory is inadequate for prediction of the interlaminar stresses which are used in predicting the delamination onset in laminated composite structures. Thus, various higher-order shear deformation models including the effects of transverse shear deformation have been developed for the determination of transverse shear stresses in laminated composites. Based on the higher-order model considering transverse normal deformation, a closed form result for laminated composite plates has been presented by Lo et al. [12,13]. Marur and Kant [14] employed three higher-order models considering the warping of the cross-section for the free vibration analysis of laminated composite 2
and sandwich beams. Subsequently, the higher-order models are extended to the static and dynamic analysis of laminated composite and sandwich plates [15,16]. In terms of the third-order model, Babu and Kant [17] presented a quadrilateral element for buckling analysis of skew laminated composite and sandwich panels. The proposed element was also used to study the thermal buckling behaviors of skew laminated composite and sandwich plates [18]. On the basis of Reddy’s theory [19], the higher-order quadrilateral elements have been proposed for free vibration analysis of isotropic, laminated composite and sandwich plates by assuming interpolations for the shear strain [20]. Moreover, the proposed elements are employed for the buckling and vibration analysis of laminated composite and sandwich plates subjected to initial stresses [21]. By means of Reddy’s theory, Sheikh and Chakrabarti [22] proposed a six-node triangular plate element for bending analysis of laminated composite plates and skew plates. Based on the third-order zig-zag theory [23], a three-node triangular plate element has been proposed to study the mechanical, thermal and electric behaviors [24]. To obtain accurately transverse shear stresses, the integration of three-dimensional equilibrium equation through the thickness of laminates has to be used. Once the equilibrium equation approach is applied, the second and the third derivatives of displacement parameters will be involved into the three-dimensional equilibrium equation. For a three-node triangular element, the displacement parameters are linear variation within one element, so that the second and the third derivatives of displacement parameters are close to zero. Thus, the equilibrium equation approach will encounter difficulty to accurately yield the transverse shear stresses for the three-node triangular element. Therefore, the stress smoothing technique within whole domain has to be used in regular meshes [24]. In order to accurately produce transverse stresses by integrating local equilibrium equation, an eight-node quadrilateral element with quadratic approximation has been proposed for bending analysis of laminated composite plates [25]. However, such element contains the second-order derivatives of transverse displacement as nodal unknowns, which are not encouraged in a practical analysis [26]. Versino et al. [27] proposed a triangular plate element for the bending analysis of homogeneous, multilayered composite and sandwich plates, whereas the transverse shear stresses obtained from the proposed element can not satisfy the continuity conditions at interfaces. Employing the principle of virtual displacements (PVD) and the Reissner mixed variational theorem (RMVT), Carrera and Demasi [28] derived the finite element formulations for the 3
bending analysis of laminated composite plates [29]. For the Reissner mixed variational theorem (RMVT), transverse stress variables are considered as the nodal parameters of finite element, which will increase the computational cost. In addition, it has been indicated that the third-order displacement models are only adequate in predicting the global response such as natural frequency and buckling stress, but such models are inadequate to predict the local response such as displacement and stress distributions in each ply of laminates [30]. Moreover, Matsunaga [31] also pointed out that it seems to be a hasty conclusion to study the global and the local behaviors of laminated composite plates using the third-order models, so that a nine-order model has been proposed for the prediction of displacement and stress distribution in each ply of laminated composite and sandwich plates. In review article [32], Carrera showed that Lekhnitskii [33] was the first to present a zig-zag model for the analysis of laminated composite beam. In addition, Ambartsumian [34] also proposed an independent zig-zag model. Differing from previous models, this paper proposes a third-order zig-zag model accounting for the zero transverse shear strain conditions at the top and the bottom free surfaces. Based on the proposed model, a refined three-node triangular plate element satisfying the requirement of C1 weak-continuity conditions in the inter-element is developed to accurately predict the displacement and stress distributions in each ply of laminates. To obtain the improved transverse shear stresses, the simplified three-dimensional equilibrium equation has been a priori used. It is significant that the higher-order derivatives of displacement parameters have been eliminated from the expression of transverse shear stresses by employing the three-field HW variational principle [35,36], which is convenient for the present model's finite element implementation. Performance of the proposed element is tested by comparing with three-dimensional elasticity solutions as well as the reference results computed using other models. A detailed study is conducted to highlight the influences of boundary condition, number of layers and aspect ratio on the static response of multilayered composite plates, and numerical results show the accuracy and range of applicability for the proposed element. Moreover, a number of examples also show that the present third-order displacement model can accurately yield the displacement and stress distribution in each ply of laminates.
2. Higher-order zig-zag theory (HZT) satisfying the bounding surface 4
free traction condition Owning to the different mechanical properties at each layer, the in-plane displacements of laminated composite structures show a zig-zag shape distribution through the thickness direction [37]. In order to model the zig-zag pattern of the in-plane displacements along the thickness, the zig-zag function proposed by Murakami [38] is introduced to displacement field of the third-order model, which can be written as
u = u0 + zu1 + z 2u2 + z 3u3 + M k ( z ) uz v = v0 + zv1 + z 2v2 + z 3v3 + M k ( z ) v z
(1)
w = w0 where k is the number of layers in laminates; u0, v0 and w0 are the displacements of a generic point in the middle surface; u1 and v1 is the rotations of normal to the middle surface about y axis and x axis, respectively; u2, u3, v2 and v3 are the higher-order terms in Taylor's series expansions; uz and vz are the generalized displacement variables associated with the zig-zag function Mk proposed by Murakami [38]. The zig-zag function Mk can be written as k
M k ( z ) = ( −1) ζ k
in which ζ k = a k z − bk , a k =
(2) z + zk 2 , bk = k +1 , and zi has been defined in z k +1 − z k z k +1 − z k
Figure 1. Employing the zero shear strain boundary conditions at the top and bottom surfaces of the plate, the higher-order zig-zag theory can be obtained, which can be given by ∂w0 ∂x ∂w v = v0 + Φ1v1 + Φ 2 vz + Φ3 0 ∂y w = w0 u = u0 + Φ1u1 + Φ 2u z + Φ3
(3)
where
5
Φ1 = z −
4z3 3h 2
z 2 2 z 3 − 3z1 z 2 M1 ( z1 ) 2 z 3 − 3z1 z 2 M n ( zn +1 ) Φ2 = M ( z) − + − 2 12 zn2+1 ∂z ∂z 2 z1 12 zn +1 3 4z Φ3 = − 2 3h k
(4)
For linear elasticity, strain of the present model HZT can be expressed as ε ε= I ε T
(5)
in which T
ε I = {ε x , ε y , γ xy } = B I EI
(6)
T
ε T = {γ xz , γ yz } = BT ET where
BI X BI = 0 0
0 BI Y
B BT = TX 0
0 BTY
EI = {E IX
0 0 B I XY
0
T
EIXY }
EIY
ET = {ETX
(7)
(8)
T
ETY }
in which B IX = B IY = [1 Φ1 Φ 2
Φ3]
B IXY = [1 Φ1 Φ 2 1 Φ1 Φ 2 ∂Φ BTX = B TY = 1 ∂z ∂u EIX = 0 ∂x ∂v EIY = 0 ∂y
∂u EIXY = 0 ∂y
∂Φ 2 ∂z
1+
∂u z ∂x
∂ 2 w0 ∂x 2
∂v1 ∂y
∂vz ∂y
∂ 2 w0 ∂y 2
∂u z ∂y
∂v0 ∂x
(9)
∂Φ 3 ∂z
∂u1 ∂x
∂u1 ∂y
2Φ 3 ]
(10) ∂v1 ∂x
∂vz ∂x
6
∂ 2 w0 ∂x∂y
ETX = u1 uz
∂w0 ∂x
ETY = v1 v z
∂w0 ∂y
Stress vector at the kth ply of the cross-ply composite plate with reference to a common structural axis system can be written as σ k = Qk ε k
where σ k = {σ xk
(11) T
T
σ yk τ xyk τ xzk τ yzk } , ε k = {ε xk ε yk γ xyk γ xzk γ yzk } .
Q11k Q12k k k Q21 Q22 Qk = 0 0 0 0 0 0
0 0
0 0
Q33k 0
0 k Q44
0
0
0 0 0 0 Q55k
k
(12)
where Qijk is the transformed elastic stiffness coefficients at the kth ply referred to the (x, y, z) coordinate system.
3. Finite element formulation In equations (6) and (10), it is found that strain components of the present model possess the first and the second derivatives of transverse displacement w, so the C0 and the C1 continuity transverse displacement functions are required. 3.1 Transverse displacement function satisfying C0 continuity in the inter-element The nine-parameter nonconforming element BCIZ [39] satisfying the C0 continuity condition can be given by
w = Fq
(13)
where
F = {F1
F2
F3 }
q = {q 1
q2
q 3}
(14)
T
in which Fi = Fi , Fxi , Fyi
(15)
T
q i = {wi , wxi , w yi }
7
where Fi = Li + L2i L j + L2i Lk − Li L2j − Li L2k ,
(
)
(
)
Fxi = c k L2i L j − c j L2i Lk + c k − c j Li L j Lk / 2 , Fyi = b j L2i Lk − bk L2i L j + b j − bk Li L j Lk / 2 , Li =
a i + bi x + ci y , ai = x j y k − xk y j , 2∆
bi = y j − y k , c i = x k − x j , (i = 1 ~ 3) . The second order derivatives of transverse displacement w can be written as ∂2w 2 ∂x
T
∂2w = Bq ∂x∂y
∂2w ∂y 2
where B = [ B1 B 2
(16)
B3 ] .
∂ 2 Fi 2 ∂x ∂ 2 Fi Bi = 2 ∂y ∂2F i ∂ x ∂ y
∂ 2 Fxi ∂x 2 ∂ 2 Fxi ∂y 2 ∂ 2 Fxi ∂x∂ y
∂ 2 F yi ∂x 2 ∂ 2 F yi , (i = 1 ~ 3) . ∂y 2 ∂ 2 F yi ∂ x ∂ y
3.2 Transverse displacement function satisfying C1 weak -continuity Following the refined element method [40], the displacement function wˆ satisfies the C1 weak-continuity condition can be given by ˆ wˆ = w + P ( Bc − B0 ) q = Fq
(17)
where w is the transverse displacement function satisfying the C0 continuity condition which B0 =
is
presented
in equation
(13);
P = x 2 / 2
y 2 / 2 xy ;
1 B dxdy , and ∆ is the area of element. ∆∫
In equation (17), it is seen that the matrix B c needs to be given to obtain the shape function Fˆ of element transverse displacement. In terms of the refined element method [40], the matrix B c can be expressed as
B c = [B c1
Bc2
B c3 ]
(18) 8
where l1m1 − l3m3 Bc1 = l3m3 −l1m1 m12 − m32
(l12 y21 + l32 y13 ) 1 (m12 y21 + m32 y13 ) 2 1 2
1 2
(l12 x12 + l32 x31)
(l12 x12 + l32 x31 ) 1 (m12 x12 + m32 x31 ) 2 2 2 1 2 (m1 y21 + m3 y13 ) 1 2
(19)
in which li and mi are cosines of the vector normal to the ith boundary. xij=xi-xj, yij=yi-yj, where xi and yi are the coordinates of node i. Matrices B c 2 and Bc 3 can be respectively presented by the permutation of the subscript. Following the equation (17), the shape function Fˆ can be given by Fˆ = F + P ( B c − B 0 )
(20)
3.3 Displacement functions of three-node triangular element
The primary displacement parameters in triangular element can be expressed in terms of nodal variables and shape functions as follows 3
3
3
i =1
i =1
i =1
3
3
3
i =1
i =1
i =1
u0 = ∑ Li u0i , u1 = ∑ Li u1i , u z = ∑ Li uzi
v0 = ∑ Li v0i , v1 = ∑ Li v1i , vz = ∑ Li vzi
(21)
3
w=
∑ (F w i
0 ,i
+ F xi w 0 , xi + F yi w 0 , yi )
0 ,i
+ Fˆ xi w 0 , xi + Fˆ yi w 0 , yi )
i =1 3
wˆ =
∑ ( Fˆ w i
i =1
where Li is area coordinate; the expressions of Fi , Fxi and Fyi can be found in equation (15); the expressions of Fˆi , Fˆxi and Fˆyi may refer to equation (20). 3.4 The strain matrix
In terms of the linear strain-displacement relationships, the strain for the kth layer can be written as follows ε = ∂u = Bδ e = [ B1 B 2
B3 ] δ e
(22)
T
where δe = δ1e δ2e δ3e , δei are displacement variables at the ith node; uk = {u v
T
w} . The nine displacement variables at each node can be written as
9
δie = u0i
v0 i
∂ ∂x ∂= 0 0
w0i
u1i
u zi
w0, xi
∂ ∂z
∂ ∂y
∂ ∂y ∂ ∂x
0
0
0 ∂ ∂z ∂ ∂y
0
∂Li ∂x 0 ∂L Bi = i ∂y 0 0
0 ∂Li ∂y ∂Li ∂x
0 ∂ ∂x
Φ1
B1,3 B2,3
B4,3
0
B5,3
vzi
w0, yi .
T
(23)
∂Li ∂x
Φ2
0 ∂Li ∂y ∂Φ1 Li ∂z
Φ1
B3,3
0
v1i
0
∂Li ∂x
B1,6
0
B2,6
∂Li ∂y ∂Φ 2 Li ∂z
Φ2
0
B3,6
0 ∂Li ∂y ∂L Φ1 i ∂x
∂Li ∂y ∂L Φ2 i ∂x
Φ1
B4,6 B5,6
0 Φ2
0 Li
0
∂Φ1 ∂z
Li
∂Φ 2 ∂z
B1,9 B2,9 B3,9 B4,9 B5,9
(24)
where
B1,3 = Φ 3
∂ 2 Fˆi ; ∂x 2
B1,6 = Φ 3
∂ 2 Fˆxi ; ∂x 2
B1,9 = Φ 3
B2,3 = Φ 3
∂ 2 Fˆi ; ∂y 2
B2,6 = Φ 3
∂ 2 Fˆxi ; ∂y 2
B2,9 = Φ 3
∂ 2 Fˆi B3,3 = 2Φ 3 ; ∂x∂y
B3,6
∂ 2 Fˆxi = 2Φ 3 ; ∂x∂y
∂ 2 Fˆyi
∂x 2
∂ 2 Fˆyi
∂y 2
B3,9 = 2Φ3
;
;
∂ 2 Fˆyi
∂x∂y
;
B4,3 = (1 +
∂Φ 3 ∂Fi ) ; ∂z ∂x
B4,6 = (1 +
∂Φ 3 ∂Fxi ) ; ∂z ∂x
B4,9 = (1 +
∂Φ 3 ∂Fyi ) ; ∂z ∂x
B5,3 = (1 +
∂Φ 3 ∂Fi ) ; ∂z ∂y
B5,6 = (1 +
∂Φ 3 ∂Fxi ) ; ∂z ∂y
B5,9 = (1 +
∂Φ 3 ∂Fyi ) . ∂z ∂y
where i=1, 2 and 3. Based on the global higher-order models, transverse shear stresses obtained from constitutive equations are not continuous at interfaces. Thus, three-dimensional equilibrium equation is required to improve the accuracy of transverse shear stresses, which is usually called as the postprocessing method based on three-dimensional equilibrium equations [41,42]. 10
Neglecting the body forces, the 3-D equilibrium equation at kth ply can be expressed as follows k ∂σ xk ∂τ xy ∂τ xzk + + =0 ∂x ∂y ∂z
∂τ xyk ∂x
+
∂σ yk ∂y
+
∂τ yzk ∂z
(25)
=0
By integrating equation (25) along the z-coordinate direction and enforcing the transverse shear free condition at the bottom surface, the transverse shear stress components satisfying continuity conditions at interfaces can be expressed as z
τˆxzk ( x, y, z ) = − ∫ ∂σ xk / ∂x + ∂τ xyk / ∂y dz z1
(26)
z
τˆyzk ( x, y , z ) = − ∫ ∂τ xyk / ∂x + ∂σ yk / ∂y dz z1
As the in-plane stresses expressed using displacement parameters of the model HZT are substituted in equation (26), a number of higher-order partial derivatives of displacement parameters are to be involved in the expression of transverse shear stresses. In order to obtain simple expression of transverse shear stresses, Iurlaro et al. [43] suggested that transverse shear stresses are expressed using two separate states of cylindrical bending in equation (26). The same approach is to be adopted in the present work, so the expression of transverse shear stresses can be written as
τˆxzk ( x, y , z ) = −Q zx U x
(27)
τˆyzk ( x, y , z ) = −Q zy U y where
{∫ Q dz, ∫ Φ Q dz, ∫ Φ Q dz, ∫ Φ Q dz} = {∫ Q dz, ∫ Φ Q dz, ∫ Φ Q dz, ∫ Φ Q dz}
Q zx = Q zy
z
z1 z
z1
k 11
k 22
z
z1
1
k 11
1
k 22
z
z1
z
2
z1
z
2
z1
z
k 11
k 11
3
z1
z
k 22
z1
T
3
k 22
(28)
∂ 2u ∂ 2u ∂ 2u ∂ 3w0 U x = 20 , 21 , 2z , ∂x 3 ∂x ∂x ∂x
T
∂ 2v ∂ 2v ∂ 2v ∂ 3w0 U y = 20 , 21 , 2z , ∂y 3 ∂y ∂y ∂y
Using the transverse shear free conditions on the top surface of plate (z=zn+1), the following equation can be given by
11
Q zx U x = 0
(29)
Q zy U y = 0
in which
{∫ = {∫
Q zx = Q zy
z n +1
z1 z n +1
z1
Q11k dz , ∫
z n +1
Q dz , ∫
zn +1
z1
k 22
z1
Φ1Q11k dz, ∫
z n +1
z1
Φ1Q dz, ∫ k 22
Φ 2Q11k dz , ∫
z n +1
z1
zn +1
z1
Φ 2Q dz , ∫ k 22
zn +1
z1
} Φ Q dz}
Φ 3Q11k dz 3
(30)
k 22
By means of the equation (30), the second derivative of displacement parameters u0 and v0 can be expressed as follows ∂ 2u0 ˆ ˆ = Q zx U x ∂x 2 ∂ 2v0 ˆ ˆ = Q zy U y ∂y 2
(31)
in which ˆ = {−Q ( 2 ) / Q (1) , − Q ( 3) / Q (1) , − Q ( 4 ) / Q (1)} Q zx zx zx zx zx zx zx ˆ = {−Q ( 2 ) / Q (1) , − Q ( 3) / Q (1) , − Q ( 4 ) / Q (1)} Q zy zy zy zy zy zy zy T
2 2 3 ˆ = ∂ u1 , ∂ u z , ∂ w0 U x 2 2 ∂x 3 ∂x ∂x
(32)
T
2 2 3 ˆ = ∂ v1 , ∂ vz , ∂ w0 U y 2 2 ∂y 3 ∂y ∂y
Inserting equation (31) in equation (27), transverse shear stress components can be rewritten as
τˆxzk ( x, y, z ) = Q zx Uˆ x
(33)
ˆy τˆyzk ( x, y , z ) = Q zy U
in which Q zα (1) = − Q zα ( 2 ) = − Q zα ( 3) = −
(Qˆ (Qˆ (Qˆ
zα
(1) Q zα (1) +Q zα ( 2 ) )
zα
( 2 ) Q zα (1) +Q zα ( 3) )
zα
( 3) Q zα (1) +Q zα ( 4 ) )
(34)
where, α is x or y. In terms of the nodal variables and shape functions of three-node triangular plate element, the equation (33) can be rewritten as
12
k τˆxz ˆ e k = Bδ τˆyz
(35)
where ˆ Bˆ = Bˆ 1 B 2
ˆ B 3
0 0 Bˆ 13 Bˆ 14 Bˆ i = ˆ 0 0 B 23 0
(36) Bˆ 15
0
ˆ B 16 ˆ B
26
0 ˆ B
27
0 ˆ B
28
ˆ B 19 ˆB 29
(37)
where 3 2 2 3 ˆ (3) ∂ Fi , Bˆ = Q ˆ (1) ∂ Li , Bˆ = Q ˆ (2) ∂ Li , Bˆ = Q ˆ (3) ∂ Fxi , Bˆ 13 = Q zx 14 zx 15 zx 16 zx ∂x 3 ∂x 2 ∂x 2 ∂x 3 3
ˆ (3) ∂ Fyi ; Bˆ 17 = Q zx ∂x 3 3 3 2 2 ˆ (3) ∂ Fi , Bˆ = Q ˆ (3) ∂ Fxi , Bˆ = Q ˆ (1) ∂ Li , Bˆ = Q ˆ (2) ∂ Li , Bˆ 23 = Q zy 26 zy 27 zy 28 zy ∂y 3 ∂y 3 ∂y 2 ∂y 2 3
ˆ (3) ∂ Fyi . Bˆ 29 = Q zy ∂y 3 The area coordinates Li are linear functions within one element, so that the second derivatives of area coordinates Li are zero. Equation (37) can be rewritten as ˆ 0 0 Bˆ 13 0 0 B 16 Bˆ i = ˆ ˆ 0 0 B 23 0 0 B 26
ˆ 0 0 B 19 ˆ 0 0 B 29
(38)
From equations (35) to (38), it is observed that it is difficult to produce accurately transverse shear stresses using the equilibrium equation based on the three-node triangular element, as only displacement parameters w0, w0,x and w0,y take part in computation of transverse shear stresses. In order to obtain accurately transverse shear stresses based on the equilibrium equation approach, the higher-order derivatives of displacement parameters in the equilibrium equation have to be eliminated by applying the three-field Hu-Washizu (HW) variational principle [35,36].
3.5 Transverse shear stresses obtained from the HW variational principle Transverse shear stress components τˆxzk and τˆyzk in equation (33) can satisfy the interlaminar continuity conditions, so that the transverse shear stress components 13
τ xzk and τ yzk obtained from constitutive equations will be respectively replaced by the τˆxzk and τˆyzk . Thus, the corresponding strain vector can be given by T
ε = {ε x
ε y γ xy γˆxz γˆ yz }
where γˆxz =
(39)
1 1 τˆxz and γˆ yz = k τˆyz . k Q44 Q55
By means of the HW variational principle discarding body forces, the functional suitable for laminated composite plates [44] can be expressed as follows 1 Π ( u, ε, σ ) = ∫ εT Qε + σ ( ε c ( u ) − ε ) d Ω − Π ext ( u ) 2 Ω
(40)
in which u, ε and σ are independent variables subjected to variation, and Ω is region occupied by a whole plate. T
∂u ∂v ∂u ∂v ∂u ∂w ∂v ∂w εc = , , + , + , + ∂x ∂y ∂y ∂x ∂z ∂x ∂z ∂y
(41)
Π ext = ∫ qu dS
(42)
S
in which q is traction vector. The first variation of the functional presented in equation (40) is used, so that the following equation can be given by
δΠ ( u, ε, σ ) = ∫ ( εT Q − σ ) δ ε dΩ + ∫ δ σ ( ε c − ε ) dΩ + ∫ σδ ε c ( u ) dΩ − δΠ ext ( u ) Ω
Ω
Ω
(43)
=0 By employing integration by parts, the following equation can be given by
∫ σδ ε ( u ) dΩ − δΠ ( u ) = ∫ Dσδ udΩ + ∫ D σδ udS ext
c
(44)
S
Ω
Ω
S
where D and DS are differential operators. By substituting the equation (44) into the equation (43), the equation (43) can be rewritten as
δΠ ( u, ε, σ ) = ∫ ( εT Q − σ ) δ ε dΩ + ∫ δ σ ( ε c − ε ) dΩ + ∫ Dσδ udΩ + ∫ DS σδ udS Ω
Ω
Ω
S
(45)
=0 On the basis of the first integral in equation (45), the following equation can be
14
written as σ = εT Q = {σ I
σT }
(46)
in which σ I = {σ x , σ y , τ xy } , and σ T = { τˆxz , τˆyz } . Inserting the equations (46) into the second integral of equation (45), the following equation can be given by
∫ δ σ (B E I
I
I
− ε I ) dΩ = 0
Ω
(47)
∫ δ σ T ( BT ET − εˆ T ) d Ω = 0 Ω
T
T
in which ε I = {ε x , ε y , γ xy } , and εˆ T = {γˆxz , γˆ yz } . Inserting the equation (33) in the second equation of (47), the higher-order ˆ and U ˆ can be written as follows derivatives of displacement parameters U x y
ˆ =Q E U IX TX x ˆ =Q E U IY TY y
(48)
where
zn +1 1 QIX = ∫ k QTzxQ zx dz z 1 Q44
−1
zn +1 1 QIY = ∫ k QTzy Q zy dz z 1 Q55
−1
∫
zn +1
∫
zn +1
z1
QTzx B TX dz
(49) z1
QTzy B TY dz
Substituting equation (48) in equation (33), the final expressions of transverse shear stress components at kth ply can be expressed as follows
τˆxzk ( x, y , z ) = ΦTX ETX
(50)
τˆyzk ( x, y , z ) = ΦTY ETY where ΦTX = Q zxQ IX and ΦTY = Q zy QIY .
Equation (50) shows that the higher-order derivatives of displacement parameters have been taken out from transverse shear stress components. 3.6 The stiffness matrix of triangular plate element Following the equations (44) and (46), the stiffness matrix of three-node triangular plate element can be presented. Equation (44) within ith element can be 15
written as follows
∫ σδ ε ( u ) d Ω − δΠ ( u ) = c
ext i
i
e h /2
∫ ∫ {σ Ai
− h /2
k x
k y
k xy
k xz
, σ ,τ ,τˆ ,τˆ
k yz
}δ {ε
k x
k y
k xy
k xz
, ε , δγ , δγ , δγ
k T yz
}
(51) dzdxdy + δ Wi = 0
where Ai is area the ith element. Inserting equations (21), (47) and (50) in equation (51), the element stiffness matrix K e within ith element can be readily written as follows Ke = ∫
Ai
∫
h /2
− h /2
BT QBdzdxdy
(52)
where strain matrix B has been defined in equation (22), and strain matrix B can be written as follows
B = B1 B 2
B 3
(53)
in which ∂Li ∂x 0 Bi = ∂Li ∂y 0 0
0
∂Li ∂y ∂Li ∂x 0 0
B1,3
Φ1
∂Li ∂x
Φ2
∂Li ∂x
B1,6 B2,6
B2,3
0
0
B3,3 B4,3 B5,3
∂L Φ1 i ∂y B4,4 0
∂L Φ2 i ∂y B4,5 0
B3,6 B4,6 B5,6
0
∂Li ∂y ∂L Φ1 i ∂x 0 B5,7
Φ1
0
∂Li ∂y ∂L Φ2 i ∂x 0 B5,8
Φ2
B1,9 B2,9 B3,9 B4,9 B5,9
where ∂F B4,3 = Skx Φ TX ( 3) i , B4,4 = Skx Φ TX (1) Li , B4,5 = Skx Φ TX ( 2 ) Li , ∂x
∂F ∂F B4,6 = SkxΦ TX ( 3) xi , B4,9 = Skx Φ TX ( 3) yi ; ∂x ∂x
∂F ∂F B5,3 = Sky Φ TY ( 3) i , B5,6 = Sky Φ TY ( 3) xi , B5,7 = Sky Φ TY (1) Li , ∂y ∂y ∂F B5,8 = Sky Φ TY ( 2 ) Li , B5,9 = Sky Φ TY ( 3) yi . ∂y
16
(54)
where Skx =
1 1 , Sky = k , and the ΦTX ( i ) , ΦTY ( i ) have been defined in k Q44 Q55
equation (50), and i=1, 2, 3. According to the following equation, the vector of nodal displacement δe can be given by
∑K δ
e e
=P
(55)
in which P is loading vector which can be presented by δ Wi in equation (51). The results are obtained from the refined three-node triangular element based on the HW variational theorem, which be called as RTTHW in the present work. For the sake of comparison, the finite element results are also computed based on the principle of minimum potential energy (PMP), namely RTTPMP.
4. Numerical examples and discussion To test the performance of the triangular element RTTHW and its range of applicability, the bending of multilayered composite plates is studied by solving numerically a large number of problems covering a wide range of lamination, geometric parameters and boundary conditions of the plates. The three-dimensional elasticity solutions of laminated composite plates presented by Pagano [1,2] are employed as the benchmark solutions. In addition, the quasi-3D solution obtained from the nine-order model HSDT-98 proposed by Matsunaga [31] are employed to assess the performance of the present element. The in-plane displacement field of the model HSDT-98 consists of 9th-order polynomial in the global thickness coordinate z and the transverse deflection is denoted by an 8th-order polynomial in z. On the other hand, a comparison with the RZTM model [43] is presented for all examples to evaluate the improvement of the proposed triangular plate element. Owning to the symmetry of geometry, material properties and loading conditions, only one-quarter of the laminated composite plate has been used. The mesh configuration can be found in Figure 2. The acronyms used in all Figures and Tables will be simply explained as follows Exact:
Exact three-dimensional elasticity solution presented by Pagano [1].
Quasi-3D: The quasi-3D solution is obtained from the nine-order model HSDT-98 proposed by Matsunaga [31]. 17
RTTHW: The present result is obtained from the refined three-node triangular
element based on the HW variational theorem. RZTM:
Based on the Reissner's mixed variational principle, analytical result is obtained from the refined zigzag theory (RZT) [43].
RHSDT: Analytical result is obtained from Reddy's higher-order shear deformation
theory [19]. FSDT:
Analytical result is obtained from the first-order shear deformation theory (FSDT).
Example 1 ：Cylindrical bending of the laminated composite plate strip subjected to a
sinusoidal transverse loading q=q0sin(π x/a) is analyzed. The ply material properties of all laminated composite plates [1,2] are given by EL=172.5GPa, ET=6.9Gpa, GLT= 3.45GPa, GTT=1.38GPa, vLT=0.25, vTT=0.25 where L and T denote the direction parallel to the fiber and the transverse direction for orthotropic materials. The dimensionless parameters of the displacements and stresses are expressed as follows u=
ET u ( 0, z ) 100ET h3 w ( a / 2, z ) σ ( a / 2, z ) τ ( 0, z ) , w= , σx = x , τ xz = xz . 4 q0 h q0 a q0 q0
For all examples, a and b are respectively the length and width of laminated plates, h is the thickness of lamination. Figure 3 shows the convergence rate of the transverse displacements with mesh refinement for 11-ply plate (a/h=4). It is seen that the present results can rapidly converge to three-dimensional elasticity solutions which are calculated by the present authors following Pagano’s methods [1]. Numerical results show that the present results with 16×16 mesh density can converge to the exact solutions. In Figure 4, it is found that the present element is free from the shear locking problem even in the analysis of the 11-ply plate with a length-to-thickness ratio (a/h) of 1000. In Table 1, displacements and stresses at the strategic point of composite plates obtained from different models are compared with three-dimensional elasticity solution (Exact) [1]. It ought to be shown that the RTTPMP and RZT in Table 1 are the results obtained directly from the principle of minimum potential energy without using the transverse shear stresses τˆxzk and τˆyzk . 18
Nevertheless, RTTHW and RZTM are the results obtained respectively from the Hu-Washizu (HW) variational principle and the mixed Reissner's variational principle considering the effects of the transverse shear stresses τˆxzk and τˆyzk . Numerical results show that the maximum percentage error of transverse displacement of eight-ply plate rapidly decreases from 15.50 to 1.41 for the model RZT as effects of the transverse shear stresses τˆxzk and τˆyzk are considered. The maximum percentage error of in-plane stress of eight-ply plate decreases from 15.29 to 9.82 for the model RZT, so the transverse shear stresses τˆxzk and τˆyzk has significant impact on displacements and stresses of the model RZT. Nevertheless, effects of the transverse shear stresses τˆxzk and τˆyzk on displacements and stresses of the model RTTHW are slight. As the number of layers is increased to 50, the maximum percentage errors of in-plane displacement and in-plane stress respectively decrease from 24.36 and 24.26 to 23.71 and 23.62 for the model RZT, which show that the transverse shear stresses
τˆxzk and τˆyzk has less effects on displacements and in-plane stresses of multilayered composite plates for the model RZT. On the other hand, the transverse shear stresses
τˆxzk and τˆyzk still impact actively the accuracy of displacements and stresses obtained from the RTTHW. In Figures 5-7, distributions of displacements and stresses through the thickness of symmetric composite plates obtained from different models are compared with three-dimensional elasticity solutions. It is found that the present results RTTHW for the displacements and the stresses are in excellent agreement with those of three-dimensional elasticity theory. The model RZTM can describe a zig-zag shape distribution of in-plane displacement through the thickness direction of composite plate whereas in-plane displacement obtained from the model RZTM is obviously diverge from the three-dimensional elasticity solutions. In addition, the transverse shear stresses obtained from the model RZTM satisfy the continuity conditions at the interfaces, but the model RZTM overestimates transverse shear stresses in comparison with the exact solutions. It is also observed that Reddy's higher-order shear deformation theory (RHSDT) is unable to describe a zig-zag shape distribution of in-plane displacements along the thickness, and interlaminar continuity of transverse shear stresses is also unable to be satisfied. For a six-ply antisymmetric plate, the displacement and the stresses obtained from the RTTHW are compared with 19
three-dimensional elasticity solutions and results obtained from other models in Figure 8. Again it is observed that the present RTTHW can accurately predict displacement and stresses. Nevertheless, the in-plane displacement on the top surface obtained from the model RZTM departs noticeably from the exact solution. Example 2: Simply-supported multilayered composite plates subjected to a doubly
sinusoidal transverse loading q=q0sin( π x/a)sin( π y/b) are analyzed. The displacements and stresses are normalized as follows (σ x , σ y )(0,0, z )h 2 τ xy (a / 2, b / 2, z )h 2 E2 h 2u (a / 2,0, z ) , , , u= (σ x , σ y ) = τ xy = a3 q0 a 2 q0 a 2
τ xz =
τ (0, b / 2, z )h τ xz (a / 2,0, z )h , τ yz = yz q0 a q0 a
Numerical results of in-plane stresses and transverse stresses are given in Table 2. It is found that the maximum percentage error of the RTTHW is 5.38 whereas the maximum percentage error of the RZTM is 34.72. Therefore, the results obtained using the RTTHW are in close agreement with three-dimensional elasticity solution [2]. Figure 9 shows the through the thickness variations of displacements and stresses for a simply-supported six-ply square plate. It shows that the results obtained from RTTHW are in good agreement with the three-dimensional solutions whereas the values predicted by the RZTM differ from three-dimensional solutions. Figure 10 shows the through the thickness variation of in-plane displacement and stresses for an eight-ply square plate. It shows that the displacement and stress values obtained from the RTTHW are in excellent agreement with exact solutions [45]. To test the range of applicability of the RTTHW, the results of parametric studies for multilayered rectangular plates are presented. These researches were conducted to provide some insight into the influences of variation in the lamination and geometric parameters of composite plates on the accuracy of the RTTHW and RZTM. Figure 11 shows the distribution of displacements and stresses of moderately thick six-ply rectangular plate (a/h=4, a/b=3), and the results are also compared with quasi-3D solutions [31]. For the displacement and stress distribution, the present results RTTHW are in excellent agreement with those of the quasi-3D theory. In Figure 12, it is observed that the results obtained using the RTTHW are quite acceptable even in the case of a very thick laminated composite plates (a/h=2, a/b=3). Compared with the quasi-3D solutions, the distribution of transverse shear stress satisfies the stress 20
boundary conditions on the upper and the lower surfaces of plate and continuity condition at interfaces. However, results obtained using the model RZTM are less accurate in comparison with the quasi-3D solution for thick composite plates. Figure 13 shows the distribution of displacements and stresses of 15-ply rectangular plate (a/h=4, a/b=3). A close agreement between the RTTHW and the quasi-3D solution is noticed, which shows that accuracy of the RTTHW is not affected by the number of layers in laminate.
5. Conclusions Based on the higher-order zig-zag model (HZT) satisfying the free surface conditions, a refined three-node triangular plate element (RTTHW) has been presented. In order to fulfill the continuity condition of transverse shear stresses at interfaces, transverse shear stresses are usually computed by integrating three-dimensional equilibrium equation through the thickness of laminates. Once integration of three-dimensional equilibrium equation is employed, the second and the third derivatives of displacement parameters will be involved into three-dimensional equilibrium equation. For the three-node triangular element, displacement parameters are linear variation within one element, so that the second and the third derivatives of displacement parameters are close to zero. Thus, three-dimensional equilibrium equation will encounter difficulty to accurately predict transverse shear stresses for the three-node triangular element. To obtain improved transverse shear stresses, the simplified three-dimensional equilibrium equation has been a priori used. It is significant that the higher-order derivatives of displacement parameters have been taken out from the expression of transverse shear stresses by employing the three-field HW variational principle. Applying the present element RTTHW in conjunction with other published models to investigate the typical bending problems of multilayered composite plates, the following conclusions can be drawn: (1) The proposed model HZT can model a zig-zag shape distribution of in-plane displacement along the thickness. Moreover, the present element RTTHW based on the model HZT can predict the more accurate in-plane displacement in comparison with the model RZTM. (2) Based on the element RTTHW, transverse shear stresses can be accurately computed without using any postprocessing approaches. Moreover, accurate transverse shear stresses have active impact on accuracy of displacements and in-plane stresses. (3) Even in the case of a very thick laminated composite plates 21
(a/h=2), the results obtained using the element RTTHW are quite acceptable. Moreover, the accuracy of the RTTHW is not affected by the number of layers and aspect ratio (a/b) in laminate.
Acknowledgement The work described in this paper was supported by the National Natural Sciences Foundation of China [No. 11272217, 11402152, 11572204].
References [1] Pagano NJ, Exact solutions for composite laminates in cylindrical bending. J Compos Mater, 1969; 3: 398-411. [2] Pagano NJ, Exact solutions for rectangular bidirectional composites. J Compos Mater, 1970; 4:20-34. [3] Bogdanovich AE, Yushanov SP, Three-dimensional variational analysis of Pagano's problems for laminated composite plates. Compos Sci Tech, 2000; 60:2407-2425. [4] Mendonca PTR, Barcellos CS, Torres DAF, Robust Ck/C0 generalized FEM approximations for higher-order conformity requirement: Application to Reddy's HSDT model for anisotropic laminated plates. Compos Struct, 2013; 96:332-345. [5] Matsunaga H, A comparison between 2-D single-layer and 3-D layerwise theories for computing interlaminar stresses of laminated composite and sandwich plates subjected to thermal loadings. Compos Struct, 2004; 64:161-177. [6] Whitney JM, The effect of transverse shear deformation on the ending of laminated plates. J Compos Mat, 1969; 3:534-847. [7] Rolfes R, Rohwer K, Improved transverse shear stresses in composite finite elements based on first-order shear deformation theory. Int J Num Meth Eng, 1997; 40:51-60. [8] Liew KM, He XQ, Tan MJ, Lim HK, Dynamic analysis of laminated composite plates with piezoelectric sensor/actuator patches using the FSDT mesh-free method. Int J Mech Sci, 2004; 46:411-431. [9] Noor AK, Malik M, An assessment of five modeling approaches for thermo-mechanical stress analysis of laminated composite panels. Comput Mech, 22
2000; 25:43-58. [10] Kant T, Swaminathan K, Analytical solutions for free vibration of laminated composite and sandwich plates based on a higher-order refined theory. Compos Struct, 2001; 53:73-85. [11] Matsunaga H, Interlaminar stress analysis of laminated composite and sandwich circular arches subjected to thermal/mechanical loading. Compos Struct, 2003; 60:345-358. [12] Lo KH, Christensen PM, Wu EM, A higher order theory of plate deformation, Part 1: Homogeneous plates. J Appl Mech, 1977; 44:663-668. [13] Lo KH, Christensen PM, Wu EM, A higher order theory of plate deformation, Part 2: laminated plates. J Appl Mech, 1977; 44:669-676. [14] Marur SR, Kant T, Free vibration analysis of fiber reinforced composite beams using higher order theories and finite element modelling. J Sound Vib, 1996; 194(3):337-351. [15] Kant T, Swaminathan K, Free vibration of isotropic, orthotropic, and multilayer plates based on higher order refined theories. J Sound Vib, 2001; 241(2):319-327. [16] Kant T, Swaminathan K, Analytical solutions for the static analysis of laminated composite and sandwich plates based on a higher-order refined theory. Compos Struct, 2002; 56:329-344. [17] Babu CS, Kant T, Two shear deformable finite element models for buckling analysis of skew fibre-reinforced composite and sandwich panels. Compos Struct, 1999; 46:115-124. [18] Kant T, Babu CS, Thermal buckling analysis of skew fibre-reinforced composite and sandwich plates using shear deformable finite element models. Compos Struct, 2000; 49:77-85. [19] Reddy JN, A simple higher-order theory for laminated composite plates. J Appl Mech, 1984; 51(12):745-752. [20] Nayak AK, Moy SSJ, Shenoi RA, Free vibration analysis of composite sandwich plates based on Reddy's higher-order theory. Composites: Part B, 2002; 33:505-519.
23
[21] Nayak AK, Moy SSJ, Shenoi RA, A higher order finite element theory for buckling and vibration analysis of initially stressed composite sandwich plates. J Sound Vib, 2005; 286:763-780. [22] Sheikh AH, Chakrabarti A, A new plate bending element based on higher-order shear deformation theory for the analysis of composite plates. Finite Element in Analysis and Design, 2003; 39:883-903. [23] Cho M, Oh J, Higher order zig-zag theory for fully coupled thermoelectric-mechanical smart composite plates. Int J Solids Struct, 2004; 41:1331-1356. [24] Oh J, Cho M, A finite element based on cubic zig-zag plate theory for the prediction of thermo-electric-mechanical behaviors. Int J Solids Struct, 2004; 41:1357-1375. [25] Icardi U, Eight-noded zig-zag element for deflection and stress analysis of plates with general lay-up. Composite Part B, 1998; 29:425-441. [26] Chakrabarti A, Sheikh AH, Analysis of laminated sandwich plates based on interlaminar shear stress continuous plate theory. J Eng Mech, 2004; 131(4):377-384. [27] Versino D, Gherlone M, Mattone M, Di Sciuva M, Tessler A, C0 triangular elements based on the refined zigzag theory for multilayer composite and sandwich plates. Composites: Part B, 2013; 44:218-230. [28] Carrera E, Demasi L, Classical and advanced multilayered plate elements based upon PVD and RMVT. Part 1: Derivation of finite element matrices. Int J Numer Meth Engng, 2002; 55:191-231. [29] Carrera E, Demasi L, Classical and advanced multilayered plate elements based upon PVD and RMVT. Part 2: Numerical implementations. Int J Numer Meth Engng, 2002; 55:253-291. [30] Wu CP, Chen WY, Vibration and stability of laminated plates based on a local higher order plate theory. J Sound Vib, 1994; 177(4):503-520. [31] Matsunaga H, Assessment of a global higher-order deformation theory for laminated composite and sandwich plates. Compos Struct, 2002; 56:279-291. [32] Carrera E. Historical review of Zig-Zag theories for multilayered plates and 24
shells, Appl Mech Rev, 2003; 56:287-308. [33] Lekhnitskii SG, Strength calculation of composite beams. Vestnik inzhen i tekhnikov, 1935, No. 9. [34] Ambartsumian SA, On a theory of bending of anisotropic plates. Investiia Akad Nauk SSSR, Ot Tekh Nauk, 1958, No. 4. [35] Hu HC, On some variational principle in the theory of elasticity and plasticity. Scientia Sinica Beijing, 1955; 4:33-54. [36] Washizu K, On some variational principle in the theory of elasticity and plasticity. Technical Report 25-18, Aeroelastic and Structures Research Laboratory, Massachusetts Institute of Technology, March 1955. [37] Carrera E, On the use of the Murakami’s zig-zag function in the modeling of layered plates and shells. Comput Struct, 2004; 82:541-554. [38] Murakami H, Laminated composite plate theory with improved in-plane response. J Appl Mech, 1986; 53:661-666. [39] Bazeley GP, Cheung YK, Irons BM, Zienkiewiz OC, Triangular elements in bending conforming and non-conforming solution, Proc. Conf. Matrix Methods in Structural Mechanics, Air Force Ins. Tech., Wright-Patterson A.F. Base, Ohio,
1965; 547-576. [40] Cheung YK, Chen WJ, Refined Nine-parameter Triangular Thin Plate Bending Element by Using Refined Direct Stiffness Method, Int J Numer Methods Eng, 1995; 38:283-298. [41] Cho M, Kim JS, Four-noded finite element post-process method using a displacement field of higher-order laminated composite plate theory. Compos Struct, 1996; 61(2):283-290. [42] Oh J, Cho M, Higher-order zig-zag theory for smart composite shells under mechanical-thermo-electric loading. Int J Solids Struct, 2007; 44:100-127. [43] Iurlaro L, Gherlone M, Di Sciuva M, Tessler A, Refined Zigzag Theory for laminated composite and sandwich plates derived from Reissner's mixed variational Theorem. Compos Struct, 2015; 133:809-817. [44] Rah K, Van Paepegem W, Habraken AM, Degrieck, A mixed solid-shell element for the analysis of laminated composites. Int J Numer Meth Engng, 2012;
25
89:805-828. [45] Bogdanovich AE, Deepak BP, Three-dimensional analysis of thick composite plates with multiple layers. Composite Part B, 1997; 28:345-357.
26
q t ( x, y )
z y
z N +1 zN
h/ 2 h/ 2
x
z1
Fig. 1 A laminated composite plate subjected to transverse loading
y
b 2
x a/2
Fig. 2 A quarter of the plate divided into a 16 × 16 mesh
27
3.23 3.22 Exact RTTHW
3.21 3.2
w
3.19 3.18 3.17 3.16 3.15 3.14
12 × 2
4 2× 4
8 ×3 8
12 4× 12
16 ×5 16
24 ×6 24
Mesh density
Fig. 3 Convergence rate of transverse displacement for 11-ply plate (a/h=4)
3.5 Exact RTTHW (16x16)
3
2.5
w
2
1.5
1
0.5
14
2 10
3 20
4 50
5 100
6 500
1000 7
a/h
Fig. 4 Transverse displacements of 11-ply square plates with different length-to-thickness ratios (a/h)
28
0.5
0.5
0.4
Exac t RTTHW (16x16) RZTM RHSDT FSDT
0.3 0.2
z h
0.4 0.3 0.2
0.5 Exact RTTHW (16x16) RZTM RHSDT FSDT
0.4 0.3 0.2
0.1
0.1
0.1
0
0
0
-0.1
-0.1
-0.1
-0.2
-0.2
-0.2
-0.3
-0.3
-0.4 -0.5 -1.5
-1
-0.5
0
-0.3
-0.4
u ( a / 2,0, z )
0.5
1
1.5
-0.5 -30
Exac t RTTHW (16x16) RZTM RHSDT FSDT
-0.4 τ xz ( a / 2, 0, z )
σ x ( 0, 0, z )
-20
-10
0
10
20
30
-0.5 -0.5
0
0.5
1
1. 5
2
2.5
3
Fig. 5 Comparison of displacements and stresses through thickness of 11-ply plate under cylindrical bending (a/h=4) 0.5
0.5
0.4
0.2
z h
0.4
Exact RTTHW (16x16) RZTM FSDT
0.3
0.3 0.2
0.5 Exact RTTHW (16x16) RZTM FSDT
0.4 0.3 0.2
0.1
0.1
0.1
0
0
0
-0.1
-0.1
-0.1
-0.2
-0.2
-0.2
-0.3
-0.3
-0.3 -0.4 -0.5 -1.5
-0.4
u ( a / 2,0, z )
-1
-0.5
0
0.5
1
1.5
-0.5 -30
-0.4
σ x ( 0, 0, z )
-20
-10
0
10
20
30
-0.5 -0.5
Exact RTTHW (16x16) RZTM FSDT
τ xz ( a / 2,0, z )
0
0.5
1
Fig. 6 Comparison of displacements and stresses through thickness of 15-ply plate under cylindrical bending (a/h=4)
29
1.5
2
0.5
0.5
0.4 0.3 0.2
z h
0.4
Exact RTTHW (16x16) RZTM FSDT
0.3 0.2
0.5 Exact RTTHW (16x16) RZTM FSDT
0.4 0.3 0.2
0.1
0.1
0.1
0
0
0
-0.1
-0.1
-0.1
-0.2
-0.2
-0.2
-0.3
-0.3
-0.3 -0.4
-0.4
u ( a / 2,0, z )
-0.5 -1.5
-1
-0.5
0
0.5
1
1.5
-0.5 -30
-0.4
σ x ( 0, 0, z )
-20
-10
0
10
20
30
Exact RTTHW (16x16) RZTM FSDT
τ xz ( a / 2,0, z )
-0.5 -0.5
0
0.5
1
1.5
2
Fig. 7 Comparison of displacements and stresses through thickness of 21-ply plate under cylindrical bending (a/h=4) 0.5
0.5
0.4
Exact RTTHW (16x16) RZTM RHSDT FSDT
0.3 0.2
z h
0.4 0.3 0.2
0.5 0.4
Exact RTTHW (16x16) RZTM RHSDT FSDT
0.2
0.1
0.1
0.1
0
0
0
-0.1
-0.1
-0.1
-0.2
-0.2
-0.2
-0.3
-0.3
-0.3
-0.4
-0.4
-0.5 -2.5
u ( a / 2, 0, z )
-2
-1.5
-1
-0.5
0
0.5
1
1.5
-0.5 -30
-20
-10
Exact RTTHW (16x16) RZTM RHSDT FSDT
0.3
σ x ( 0,0, z )
-0.4
0
-0.5 0
10
20
30
τ xz ( a / 2,0, z )
0.5
1
1.5
2
2.5
Fig. 8 Comparison of displacements and stresses through thickness of six-ply plate under cylindrical bending (a/h=4)
30
3
3.5
0.5
0.5 Quasi-3D RTTHW (16x16) RZTM RHSDT FSDT
0.4 0.3 0.2
z h
0.4 0.3 0.2
0.4
0.2
0.1
0.1
0
0
0
-0.1
-0.1
-0.1
-0.2
-0.2
-0.2
-0.3
-0.3
-0.4
-0.4
u ( a / 2, 0, z )
-0.015
-0.01
-0.005
0
0.005
0.01
0.5 0.4 0.3 0.2
-0.5 -0.8
-0.3
-0.6
-0.4
-0.2
0
0.2
0.4
0.6
Quasi-3D RTTHW (16x16) RZTM RHSDT FSDT
0.4 0.3 0.2
0.1
0.1
0.8
-0.5 -0.8
0
-0.2
-0.2
-0.2
-0.3
-0.3
-0.3
-0.02
-0.01
0
0.01
-0.4 0.02
0.03
0.04
0.2
0.4
0.6
Quasi-3D RTTHW (16x16) RZTM RHSDT FSDT
0.1
-0.1
-0.03
0
0.2
0
-0.5 0
-0.2
0.3
-0.1
τ xy ( a / 2, b / 2, z )
-0.4
0.4
0
-0.5 -0.04
-0.6
0.5
-0.1
-0.4
σ y ( 0, 0, z )
-0.4
σ x ( 0, 0, z )
0.5 Quasi-3D RTTHW (16x16) RZTM RHSDT FSDT
Quasi-3D RTTHW (16x16) RZTM RHSDT FSDT
0.3
0.1
-0.5 -0.02
z h
0.5 Quasi-3D RTTHW (16x16) RZTM RHSDT FSDT
0.05
0.1
0.15
0. 2
0.25
0.3
τ xz ( a / 2,0, z )
-0.4
0.35
-0.5 0
0. 4
τ yz ( 0, b / 2, z ) 0.05
0.1
0.15
0.2
0.25
0.3
0.35
Fig. 9 Comparison of displacements and stresses through thickness of six-ply plate subjected to double sinusoidal load (a/h=4)
31
0.4
0.8
0.5
0.5
0.4 0.3 0.2
z h
0.4
Quasi-3D RTTHW (16x16) RZTM RHSDT FSDT
0.2
0.1
0.1
0
0
-0.1
-0.1
-0.2
-0.2
-0.3
-0.3
-0.4 -0.5 -0.015
-0.4
u ( a / 2,0, z ) -0.01
-0.005
0
0.005
0.01
0.5 0.4 0.3 0.2
z h
σ x (0, 0, z )
-0.5 -0.8
-0.6
-0.4
-0.2
0
0.2
0.4
0.6
0.8
0.5 0.4
Quasi-3D RTTHW (16x16) RZTM RHSDT FSDT
0.3
Exact RTTHW (16x16) RZTM RHSDT FSDT
0.2
0.1
0.1
0
0
-0.1
-0.1
-0.2
-0.2
-0.3
-0.3
τ xy ( a / 2, b / 2, z )
-0.4 -0.5 -0.04
Exact RTTHW (16x16) RZTM RHSDT FSDT
0.3
-0.03
-0.02
-0.01
0
τ xz ( a / 2,0, z )
-0.4
0.01
0.02
0.03
0.04
-0.5
-0.05
0
0.05
0.1
0.15
0.2
0.25
0.3
0.35
Fig. 10 Comparison of displacements and stresses through thickness of eight-ply plate subjected to double sinusoidal load (a/h=4)
0.5
0.5
0.4
0.4
0.3 0.2
z h
0.1
0.3
Quasi-3D RTTHW (16x16) RZTM RHSDT FSDT
0.2 0.1
0
0
-0.1
-0.1
-0.2
-0.2
-0.3
-0.3
-0.4 -0.5 -0.04
-0.4
u ( a / 2, 0, z ) -0.03
-0.02
Quasi-3D RTTHW (16x16) RZTM RHSDT FSDT
-0.01
0
0.01
0.02
32
-0.5 -1.5
σ x (0, 0, z )
-1
-0.5
0
0.5
1
1.5
0.5
0.5
0.4 0.3 0.2
z h
Quasi-3D RTTHW (16x16) RZTM RHSDT FSDT
0.4
Quasi-3D RTTHW (16x16) RZTM RHSDT FSDT
0.3 0.2
0.1
0.1
0
0
-0.1
-0.1
-0.2
-0.2
-0.3
-0.3
τ xy ( a / 2, b / 2, z )
-0.4 -0.5 -0.04
-0.03
-0.02
-0.01
0
τ xz ( a / 2,0, z )
-0.4
0.01
0.02
0.03
0.04
-0.5 0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
Fig. 11 Comparison of displacements and stresses through thickness of six-ply rectangular plate subjected to double sinusoidal load (a/h=4, b/a=3) 0.5
0.5
0.4
0.4
0.3
Quasi-3D RTTHW (16x16) RZTM RHSDT FSDT
0.2
z h
0.1
0.3 0.2 0.1
0
0
-0.1
-0.1
-0.2
-0.2
-0.3
-0.3
-0.4
-0.4
u ( a / 2, 0, z )
-0.5 -0.08
-0.06
-0.04
-0.02
0
0.02
0.04
0.5
-0.5 -0.6
σ y ( 0,0, z ) -0.4
-0.2
0
0.2
0.4
0.6
0.8
0.5
0.4
0.4
Quasi-3D RTTHW (16x16) RZTM RHSDT FSDT
0.3 0.2
z h
Quasi-3D RTTHW (16x16) RZTM RHSDT FSDT
0.3 0.2
0.1
0.1
0
0
-0.1
-0.1
-0.2
-0.2
-0.3
-0.3
τ xy ( a / 2, b / 2, z )
-0.4 -0.5 -0.08
-0.06
-0.04
-0.02
0
0.02
-0.4 0.04
0.06
0.08
-0.5 -0.2
Quasi-3D RTTHW (16x16) RZTM RHSDT FSDT
τ xz ( a / 2,0, z ) -0.1
0
0.1
0.2
0.3
0.4
0.5
Fig. 12 Comparison of displacements and stresses through thickness of six-ply rectangular thick plate subjected to double sinusoidal load (a/h=2, b/a=3)
33
0.6
0.5
0.5
0.4
0.2
z h
0.4
Quasi-3D RTTHW (16x16) RZTM RHSDT FSDT
0.3
0.2
0.1
0.1
0
0
-0.1
-0.1
-0.2
-0.2
-0.3
-0.3
-0.4
0
0.005
0.01
0.015
0.02
0.5 0.4 0.3 0.2
z h
-0.5 -1.5
-1
-0.5
0
0.5
1
1.5
0.5 0.4
Quasi-3D RTTHW (16x16) RZTM RHSDT FSDT
0.3
Quasi-3D RTTHW (16x16) RZTM RHSDT FSDT
0.2
0.1
0.1
0
0
-0.1
-0.1
-0.2
-0.2
-0.3
-0.3
τ xy ( a / 2, b / 2, z )
-0.4 -0.5 -0.04
σ x (0, 0, z )
-0.4
u ( a / 2, 0, z )
-0.5 -0.02 -0.015 -0.01 -0.005
Quasi-3D RTTHW (16x16) RZTM RHSDT FSDT
0.3
-0.03
-0.02
-0.01
0
0.01
τ xz ( a / 2,0, z )
-0.4 0.02
0.03
0.04
-0.5 -0.1
0
0.1
0.2
0.3
0.4
0.5
0.6
Fig. 13 Comparison of displacements and stresses through thickness of 15-ply rectangular plate subjected to double sinusoidal load (a/h=4, b/a=3)
34
0.7
Table 1 Effect of transverse shear stresses on displacements and in-plane stresses (a/h=4) w ( a / 2, 0 )
u ( 0, − h / 2 )
σ x ( a / 2, − h / 2 )
Exact
3.7391
1.1524
-22.6552
RTTHW(16×16)
3.7456 (0.17)
1.1930 (3.52)
-23.7562 (4.86)
RTTPMP(16×16) 3.6373 (2.72)
1.1501 (0.20)
-22.9051 (1.10)
RZTM
3.7920 (1.41)
1.0379 (9.94)
-20.4303 (9.82)
RZT
3.1597 (15.50)
0.9749 (15.40)
-19.1911 (15.29)
Exact
3.5812
1.1797
-23.1923
RTTHW(16×16)
3.5752 (0.17)
1.2261 (3.93)
-24.4078 (5.24)
RTTPMP(16×16) 3.4410 (3.91)
1.1204 (5.03)
-22.3111 (3.80)
RZTM
3.6519 (1.97)
0.9865 (16.38)
-19.4185 (16.27)
RZT
3.1606 (11.74)
0.9599 (18.63)
-18.8949 (18.53)
Exact
3.5379
1.2325
-24.2313
RTTHW(16×16)
3.5286 (0.26)
1.2844 (4.21)
-25.5560 (5.47)
RTTPMP(16×16) 3.3591 (5.05)
1.1256 (8.67)
-22.4092 (7.52)
RZTM
3.6166 (2.22)
0.9668 (21.56)
-19.0315 (21.46)
RZT
3.1613 (10.64)
0.9540 (22.60)
-18.7797 (22.50)
Exact
3.5310
1.2591
-24.7538
RTTHW(16×16)
3.5192 (0.33)
1.3129 (4.27)
-26.1186 (5.51)
RTTPMP(16×16) 3.3333 (5.60)
1.1309 (10.18)
-22.5118 (9.06)
RZTM
3.6097 (2.23)
0.9605 (23.71)
-18.9079 (23.62)
RZT
3.1615 (10.46)
0.9524 (24.36)
-18.7473 (24.26)
Number of Models layers 8
16
32
50
The number in bracket is the % errors of various models relative to three-dimensional elasticity solution.
35
Table 2 Comparison of results obtained from various models with three-dimensional elasticity solutions (a/h=4)
σx
τ xy
τ xz
τ yz
(a/2,b/2,h/2)
(0,0,h/2)
(0,b/2,0)
(a/2,0,0)
Exact
0.8010
-0.0511
0.256
0.217
RTTHW(16×16)
0.8175 (2.06)
-0.0525 (2.74)
0.254 (0.78)
0.212 (2.30)
RZTM
0.6873 (14.19)
-0.0488 (4.50)
0.258 (0.78)
0.219 (0.92)
FSDT
0.4369 (45.46)
-0.0369 (27.79)
0.120 (53.13)
0.130 (40.09)
Exact
0.684
-0.0337
0.223
0.223
RTTHW(16×16)
0.702 (2.63)
-0.0342 (1.48)
0.235 (5.38)
0.224 (0.45)
RZTM
0.577 (15.64)
-0.0220 (34.72)
0.257 (15.25)
0.230 (3.14)
FSDT
0.466 (31.87)
-0.0224 (33.53)
0.244 (9.42)
0.083 (62.78)
Number of layers
Models
3
9
The number in bracket is the % errors of various models relative to three-dimensional elasticity solution.
36