Aerothermoelastic analysis of composite laminates with variable fiber spacing

Computational Materials Science 91 (2014) 83–90

Contents lists available at ScienceDirect

Computational Materials Science journal homepage: www.elsevier.com/locate/commatsci

Aerothermoelastic analysis of composite laminates with variable fiber spacing Shih-Yao Kuo ⇑ Department of Air Transportation Management, Aletheia University, Tainan City 721, Taiwan, ROC

a r t i c l e

i n f o

Article history: Received 10 February 2014 Received in revised form 26 March 2014 Accepted 21 April 2014

Keywords: Variable fiber spacing Thermal postbuckling Vibration Flutter Composite laminate

a b s t r a c t This study presents the effects of variable fiber spacing on the thermal postbuckling, vibration, and flutter behaviors of composite laminates subjected to aerodynamic force and thermal stress. A 54 degree-offreedom high order triangular plate element is developed based on the Von Karman large deflection assumptions and quasi-steady supersonic aerodynamic theory. The numerical results reveal that the redistribution of fibers can dramatically increase the critical buckling temperature, and efficiently increase the natural frequencies and flutter boundary. The sequence of the natural mode may be altered and the frequency coalescence pair may change. The stability boundary is different and some specific phenomena are discussed. Ó 2014 Elsevier B.V. All rights reserved.

1. Introduction Many studies have examined the buckling, vibration, and flutter behaviors of composite laminates, and thus the related phenomena are now reasonably well understood. Such works generally consider the effects of fiber angle, plate shape, boundary conditions, material properties, and stacking sequence, while few consider laminates with variable fiber spacing (VFS) [1–11]. Martin and Leissa [1] investigated the plane stress problem of a rectangular composite sheet with variable fiber content. The stresses obtained from this analysis were treated as inputs for the first vibration and buckling study of composites with VFS [2]. They also presented exact solutions for the stress, strain and displacement fields for four types of problems with arbitrary fiber spacing [3]. Kuo and Shiau [10] studied the buckling and vibration of composite laminates with VFS using the finite element method. More recently, Kuo [11] presented the first flutter analysis of composite laminates with VFS. All of the above-mentioned works adopted small deflection theory to carry out the buckling, vibration, and flutter analyses of laminates with VFS. High-speed aircraft structural panels are subjected not only to aerodynamic loading, but also to aerodynamic heating, and this rise in temperature may buckle the plate and reduce its load-carrying capacity. Plates in a buckled state may be expected to survive under dynamic disturbances, and thus how such conditions affect composite laminates is of considerable interest. The postbuckling ⇑ Tel.: +886 6 5703100x7431. E-mail address: [email protected] http://dx.doi.org/10.1016/j.commatsci.2014.04.045 0927-0256/Ó 2014 Elsevier B.V. All rights reserved.

of laminates [12–16], vibration of buckled laminates [17–20], and flutter of buckled laminates [21–30] have been studied by many scholars. However, to the best of the author’s knowledge, no study has yet investigated the influence of VFS on the aerothermoelastic behavior of laminates. This study develops a 54 degree-of-freedom high order triangular plate element based on the Von Karman large deflection assumptions and quasi-steady supersonic aerodynamic theory. The formulation of the location dependent linear stiffness matrix, first-order nonlinear stiffness matrix, second-order nonlinear stiffness matrix, consistent mass matrix, thermal loading due to nonhomogeneous material properties, geometrical stiffness matrix due to thermal effects, aerodynamic matrix, and aerodynamic damping matrix are derived. This research attempts to study the effects of VFS on the thermal postbuckling, vibration, and flutter behaviors of composite laminates subjected to aerodynamic force and thermal stress in the frequency and time domains.

2. Formulations Consider a rectangular composite laminate with length a, width b, and thickness h. Supersonic airflow with air density q1, flow velocity V1, Mach number M1, and aerodynamic pressure DP is assumed passing over the top surface of the plate along the x-axis. The rise in temperature on the top surface (Tt) is assumed to be higher than that on the bottom surface (Tb) if the temperature gradient (Tg) is considered. The temperature distribution DT(z) may vary across the panel thickness, but not along the panel, i.e.

84

S.-Y. Kuo / Computational Materials Science 91 (2014) 83–90

DTðzÞ ¼ T u þ

Tg Tt þ Tb Tt  Tb þ z¼ z h 2 h

ð1Þ

The displacement field may be expressed as:

8 9 > =

8 9 > < u0 > =

8 9 > < w0;x > = v > ¼ > v 0 >  z> w0;y > > : ; : ; : ; w w0 0

a1 ðfÞ ¼ ð2Þ

where u0 and v0 are the in-plane displacement components in the mid-plane of the plate, and w0 is the out-of-plane displacement component in the mid-plane of the plate. Adopting the Von Karman model in order to consider moderately large deflections, the kinematic relation can be determined as

feg ¼ fe0 g þ zfjg þ fdg

ð3Þ

where {e0}, {j}, and {d} are the strain, plate curvature, and large deflection strain in the mid-plane, respectively. The plate is assumed to consist of N layers of orthotropic sheets bonded together. The fibers in each layer are aligned parallel to the longitudinal direction, but distributed unevenly in the transverse direction. The fiber volume fraction, Vf, is thus a function of the nondimensional coordinate f = 2y/b having its origin at the plate center. For the fiber distribution function [10]

  n V f ðfÞ ¼ ðV f Þin  ðV f Þout ð1  f2 Þ þ ðV f Þout ;

n ¼ 1; 2; 3

ð4Þ

where (Vf)in is the fiber volume fraction at the plate center (f = 0) and (Vf)out is the fiber volume fraction at the edges (f = ±1), the volume fraction index n controls the variation of the volume fraction. It is seen that the more fibers are centralized in the central portion of the lamina for the higher volume fraction (Vf)in and higher order n, as shown in Fig. 1. With this VFS, the elastic moduli E1, E2, m12, G12, mass density q, and the coefficients of thermal expansion a1 and a2 for the composite material are also functions of f. The formulas used for the calculation of these effective engineering constants are based on the following rules of mixture.



E1 ðfÞ ¼ Ef V f ðfÞ þ Em 1  V f ðfÞ E2 ðfÞ ¼



Ef Em Ef ð1  V f ðfÞÞ þ Em V f ðfÞ

m12 ðfÞ ¼ mf V f ðfÞ þ mm ð1  V f ðfÞÞ G12 ðfÞ ¼

Gf Gm Gf ð1  V f ðfÞÞ þ Gm V f ðfÞ

qðfÞ ¼ qf V f ðfÞ þ qm ð1  V f ðfÞÞ

ð5dÞ

ð5fÞ ð5gÞ

It is known that thermal stresses are not caused by external loads, but are the consequences of restrained thermal distortion. The stress–strain relation of the composite laminate subjected to temperature rise DT is given by:

 ðfeg  DTfag Þ frgk ¼ ½Q k k

ð6Þ

  is the transformed reduced stiffness matrix, and k where ½Q denotes the kth layer of the laminate. The force and moment resultants of the laminate are defined as:

fNg ¼ ½Afe0 g þ ½Afdg þ ½Bfjg  fNDT g

ð7aÞ

fMg ¼ ½Bfe0 g þ ½Dfjg  fMDT g

ð7bÞ

where [A], [B], and [D] are the extensional, coupling, and bending stiffness matrices, while {NDT} and {MDT} are the thermal forces and thermal moments induced by the rise in temperature. The total strain energy U, kinetic energy due to the transverse oscillation of the plate T, and work done by the nonconservative aerodynamic loading Wnc can be expressed as:

Z fe0 gT ½Afe0 gdA þ 1=2 fjgT ½DfjgdA Z þ 1=2 ðfdgT ½Afe0 g þ fe0 gT ½AfdgÞdA Z Z þ 1=2 fdgT ½AfdgdA  1=2 fdgT fNDT gdA Z    1=2 fe0 gT fNDT g þ fjgT fM DT g dA

U ¼1=2

T¼

ð5cÞ

Ef af V f ðfÞ þ Em am ð1  V f ðfÞÞ Ef V f ðfÞ þ Em ð1  V f ðfÞÞ

a2 ðfÞ ¼ af V f ðfÞ þ am ð1  V f ðfÞÞð1 þ mm Þ

ð5aÞ ð5bÞ

ð5eÞ

1 2

Z

Z

W nc ¼

qh

Z

ð8aÞ

 2 dw dA dt

ð8bÞ

D p wdA

ð8cÞ

For high supersonic Mach numbers (M1 > 1.6), the aerodynamic pressure loading is assumed to follow quasi-steady supersonic aerodynamic theory, as follows:

! dw M 21  2 1 dw DP ¼  qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi þ dx M 21  1 V 1 dt M21  1

q1 V 21

ð9Þ

For accuracy and simplicity, the three displacement functions for the in-plane and out-of-plane deformations are assumed to have the same form. The strain energy, kinetic energy, and work done by the aerodynamic loading can be expressed as: DT

U ¼1=2fqgT ð½k þ 1=3½n1  þ 1=6½n2 Þfqg  1=2DTfqgT ½kN fqg  DTfqgT ff DT g T ¼ 1=2

T

dq dq ½m dt dt



q1 V 21 q V 1 ðM21  2Þ T dq ffi fqgT ½afqg þ 1 W nc ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi fqg ½c 3=2 2 dt 2  1Þ ðM M1  1 1 Fig. 1. Fiber distribution function.

ð10aÞ

ð10bÞ

ð10cÞ

where [k], [n1], and [n2] are the element linear, first-order, and DT second-order nonlinear stiffness matrices, respectively; ½kN  is the

85

S.-Y. Kuo / Computational Materials Science 91 (2014) 83–90

geometric stiffness matrix due to thermal effects; [m] is the mass matrix; [a] and [c] are the element aerodynamic and aerodynamic damping matrices, respectively; {fDT} is the thermal loading vector; and {q} contains the nodal degrees of freedom. The following nondimensional parameters and constants are introduced when assembling all the element matrices and load vectors

q1 V 21

3

3

E21 h k ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 0 ; D011 ¼ ; 12ðE1  m212 E2 Þ M21  1 D11 sffiffiffiffiffiffiffiffiffiffiffi q a D011 l ¼ 1 ; s ¼ x0 t; x0 ¼ qh qha4 b

þ

a4 D011

sffiffiffiffiffiffiffiffiffi a3 g 3 ½CfQ_ g b

ð½K þ 1=2½N1  þ 1=3½N2   DT½K NDT Þ þ k

a4 b

!

½A fQg 3

DTa n DT o ¼ 0 F D11 4

ð11Þ

in which the capital letters represent the overall structural equations, while the dots indicate differentiation with respect to nondimensional time s. Because aerodynamic damping g is usually small, it is advantageous to ignore g when analyzing the behavior in the frequency domain [29]. Let the total response of the plate be the sum of displacement due to static deformation and that due to dynamic deformation, Eq. (11) can then be written into two sets of equations as:

½K þ 1=2½N1  þ 1=3½N2   DT½K DNT  þ k

D011 3

b

Distribution

Mesh

DTcr (°C)

xn

kcr

Vf = 2/3

Vf = 1  f2

44 66 88 66 88 10  10 12  12

24.53 24.53 24.53 52.17 53.86 54.60 54.93 (+123.9%)

10.64 10.64 10.64 12.90 13.04 13.09 13.11 (+23.2%)

250.58 250.64 250.64 337.62 340.88 342.61 343.55 (+37.1%)

Vf = 8/15 Vf = (1  f2)2

66 12  12

16.15 34.61 (+114.3%)

9.57 12.04 (+25.8%)

200.40 295.71 (+47.6%)

Vf = 16/35 Vf = (1  f2)3

66 12  12

12.72 24.83 (+95.2%)

8.94 11.37 (+27.2%)

172.58 265.48 (+53.8%)

kl g¼ ; M1

The governing equation in matrix form for the aerothermoelastic analysis of the composite laminate is obtained as follows:

1 €g þ ½MfQ qh

Table 1 Critical buckling temperature, first natural frequency, and flutter boundary.

! ½A fQ s g ¼ DTfF DT g

mesh will give satisfying results for a square composite laminate with uniform fiber spacing. In the following, a 6  6 mesh is used for the cross-ply laminates with uniform fiber distribution (Vf = 2/3, Vf = 8/15, and Vf = 16/35) and a 12  12 mesh is used for the laminates with VFS (V f ¼ 1  f2 ; Vf = (1  f2)2, Vf = (1  f2)3). If more fibers are distributed in the central portion of the plate, then this can efficiently increase the critical buckling temperature, first mode natural frequency, and flutter boundary. In particular, the critical buckling temperature is increased significantly due to the redistribution of fibers. Fig. 2 provides a clear picture of the effects of the fiber volume fraction (Vf) on the material properties. The stiffness of the plate is mainly provided by the graphite fibers in the longitudinal direction (E11), and the thermal stress is mainly induced by the epoxy matrix in the transverse direction (a2). If more fibers are distributed in the central portion of the plate, then this can strengthen the plate in an efficient manner. On the other hand, higher thermal stress exists

ð12aÞ ja ½M þ ½K þ ½N1  þ ½N2   DT½K DNT  þ k

D011 3

b

! ½A fQ d g ¼ f0g ð12bÞ

where ja is the nondimensional eigenvalue that is equal to the square of the nondimensional natural frequencies. Eq. (12a) is a set of nonlinear algebraic equations that yields the static deflection of the thermally buckled laminate, and Eq. (12b) is a set of linear equations for the determination of the dynamic behavior of the thermally buckled laminate. The nonlinear equation in Eq. (12a) is solved by the Newton–Raphson method. After the equilibrium state is obtained, Eq. (12b) becomes a standard eigenvalue problem. In addition, the Newmark method is also employed to obtain the dynamic response of the system in the time domain. 3. Results and discussion The correctness of the aerothermoelastic analysis of composite laminates has been confirmed in previous studies [11,29], and thus the element used in this work is an acceptable one. Table 1 shows the critical buckling temperature (DTcr), first natural frequency (xn) and flutter boundary (kcr ) of the simply supported square [0/90/0/90]s graphite-epoxy laminate with different fiber distribution functions (V f ¼ 1  f2 ; Vf = (1  f2)2, Vf = (1  f2)3) and the corresponding uniform fiber distribution with same fiber volume fraction (V f ¼ 2=3; Vf = 8/15, Vf = 16/35). Since the 54 degree-offreedom element is efficient and accurate, even a coarse 4  4

Fig. 2. Fiber volume fraction and material properties.

86

S.-Y. Kuo / Computational Materials Science 91 (2014) 83–90

along the edges of the plate where the in-plane displacement is fixed. These two factors make the effect of the fiber distribution on the critical buckling temperature more significant. Therefore, the critical bucking temperature of the laminate with Vf = 1  f2 (54.93 °C) is much higher than the one with Vf = 2/3(24.53 °C). As shown in Fig. 2, the density distribution at the central portion of the plate is also higher than at the outer portion due to VFS. This explains why the increase in the natural frequencies of the laminate with Vf = 1  f2 is only 23.2% in Table 1. If more fibers are distributed in the central portion of the plate, then this can efficiently increase the buckling load and natural frequencies. However, if more fibers are distributed in the outer portion of the plate, then this may also increase the critical buckling load [10]. This is why the increase in the critical buckling temperature is relatively limited for the laminate with the higher volume fraction index n. On the other hand, the increase in the natural frequencies is more significant for the laminate with a higher volume fraction index. Fig. 3 shows the static and dynamic behaviors of the cross-ply laminate with temperature vs. nondimensional static deflections (w/h) and natural frequencies (x/x0). Before the critical buckling temperature, the plate remains flat. As the critical buckling temperature is reached, the plate will buckle into mode (1, 1). The frequency of the natural mode (1, 1) with a similar mode shape will drop to zero. When the plates are subjected to an uneven temperature gradient (Tg = 1 for Tt = 1.5 DT and Tb = 0.5 DT; Tg = 2 for Tt = 2 DT and Tb = 0), the plate will be deformed into the (1, 1) shape as the temperature is increased from the reference temperature, and none of the natural frequencies will drop to zero as the temperature is increased to the critical buckling temperature due to the presence of the thermal moment in the plate. In addition, the natural frequencies do not change abruptly at the critical buckling temperature, because the plate is already deformed before buckling. The 3D thermal buckling mode shapes and contour plots for the laminates with VFS (Vf = 1  f2) and corresponding uniform fiber

distribution (Vf = 2/3) are shown in Fig. 4. It is seen that the ‘‘plateau’’ buckling mode shape for the laminate with VFS is no longer the half sine wave form. Fig. 5 shows the frequency coalescence plot of the laminates with Vf = 1  f2 and Vf = 2/3. The coalescence occurs between the 1st mode (1, 1) and 3rd mode (2, 1). If more fibers are distributed in the central portion of the plate, then this can efficiently increase the flutter boundary. Fig. 6 shows the first and third vibration mode shapes, and the flutter mode shape of the plate. The redistribution of fibers makes the natural mode and flutter mode slightly different. Fig. 7 shows the map in the k, DT/DTcr space which identifies the stability boundaries for the laminates, where k is nondimensional aerodynamic pressure and DT/DTcr is nondimensional temperature rise. The extra axis (DT=DT 0cr ) is also provided for convenience, where the critical buckling temperature DT 0cr ¼ 54.93 °C for the laminate with Vf = 1  f2. Four types of panel behavior, namely flat, buckled, limit-cycle, and chaotic motions, are shown. In addition, the effects of the temperature gradient are also presented. For a plate subjected to a temperature gradient,

Fig. 4. Buckling mode shape and contour plot.

Fig. 3. Thermal postbuckling deflection and natural frequency.

Fig. 5. Frequency coalescence plot.

S.-Y. Kuo / Computational Materials Science 91 (2014) 83–90

87

Fig. 6. Natural mode shape and flutter mode shape.

Fig. 8. Frequency coalescence plot (Vf = 1  f2).

Fig. 7. Stability boundary (Vf = 1  f2 and Vf = 2/3).

the thermal moment makes the stability boundary between the flat and buckled regions disappear. Moreover, the temperature gradient increases the overall stiffness of the plate, which in turn stabilizes the plate and increases the flutter boundary. It is interesting to note that the flutter boundary of the laminate with Vf = 1  f2 is no longer a straight line, and the slope of the flutter boundary is not continuous if the plate is subjected to the temperature gradient Tg = 2. Fig. 8 shows the frequency coalescence plots of the laminates with Vf = 1  f2 at four different temperatures. The coalescence

occurs between the 1st mode (1, 1) and 3rd mode (2, 1) when DT=DT 0cr < 1.26, while the coalescence pair becomes the 2nd mode (1, 2) and 4th mode (2, 2) when DT=DT 0cr > 1.26. Fig. 9 shows the first four vibration mode shapes of the laminate at various airflow speeds and DT=DT 0cr = 1.26. As the airflow speed increases, these mode shapes gradually change. Finally, the mode shapes of the coalescence pairs become the same for k P kcr . The 1st and 3rd modes thus become the same. In addition, the 2nd and 4th modes also become the same. In other words, these two pairs coalesce simultaneously as DT/DTcr = 1.26. The change in the coalescence pair causes a change in the slope of the flutter boundary. A similar phenomenon is also observed for the plates subjected to

88

S.-Y. Kuo / Computational Materials Science 91 (2014) 83–90

Fig. 9. Contour plot of natural mode and flutter mode (DT=DT 0cr = 1.26).

temperature gradient Tg = 2 and a higher rise in temperature. The coalescence occurs between the 5th (1, 3) and 6th (2, 3) modes. It can be seen that the coalescence pair should be between the (1, n) and (2, n) modes, where n = 1, 2, 3. Figs. 10–12 show the time history responses and phase plane plots of the laminates with Vf = 1  f2. The complete time history response is displayed during the nondimensional time interval [0, 10], and the steady state phase plane plot is shown by deleting the transient response during the time interval [0, 4]. These time history responses and phase plane plots verify the behaviors of the laminate in the time domain. If the dynamic pressure k is fixed at 300, the time history responses and phase plane plots for the plate motion at different DT=DT 0cr are shown in Fig. 10. The locus in the phase plane plot for the plate subjected to a uniform temperature change is symmetric to the origin. Once the flutter begins, the displacement and velocity of limit cycle oscillation will increase with the temperature. Fig. 11 shows the effects of the temperature gradient on the responses of the laminates. The time history response moves upward, and the plate will buckle due to the

thermal moments. The buckled plate corresponds to the nonlinear static equilibrium point in the phase plane. The higher the temperature gradient is, the larger the deflection. When the temperature rise DT=DT 0cr is fixed at 3, the time history responses and phase plane plots of chaotic motion at different dynamic pressure k are shown in Fig. 12. Finally, the effects of the fiber distribution function on the aerothermoelastic behavior of the laminate are shown in Figs. 13–15. Fig. 13 shows the frequency coalescence plot of the laminates with three different fiber distribution functions. The coalescence occurs between the 2nd (1, 2) and 4th modes (2, 2) for the higher volume fraction index n = 3 (Vf = (1  f2)3). In summary, the higher modes may lead to coalescence for the laminates with a more uneven fiber distribution, or plates subjected to a higher temperature gradient. Fig. 14 shows the stability boundaries for the laminates with different fiber distribution functions. The axis (DT=DT 00cr ) is provided for convenience where the critical buckling temperature DT 00cr ¼ 34.61 °C for the laminate with Vf = (1  f2)2. The slopes of

S.-Y. Kuo / Computational Materials Science 91 (2014) 83–90

Fig. 10. Time history response and phase plane plot (k = 300).

89

Fig. 12. Time history response and phase plane plot (DT=DT 0cr = 3).

Fig. 13. Frequency coalescence plot (DT/DTcr = 0).

Fig. 11. Time history response and phase plane plot (DT=DT 0cr = 3, k = 300).

the flutter boundaries of the laminate with Vf = (1  f2)2 subjected to Tg = 2 and Tg = 1 are not continuous. The phenomenon of a closed section for chaotic motion is found for the laminate with Vf = (1  f2)2. Fig. 15 shows the frequency coalescence plots of the laminates with Vf = (1  f2)2 at four different temperatures. When DT=DT 00cr = 1.4, the plate is in a buckled condition for the lower dynamic pressure. The eigenvalue of mode (1, 1) will drop to zero because the buckle pattern changes and the plate will have a flat condition when the dynamic pressure reaches the critical value (k=85).

As the dynamic pressure reaches the flutter boundary (kcr = 117), the two real eigenvalues become a complex conjugate pair, and the plate will be in a limit cycle condition. When DT=DT 00cr = 1.7, the two real eigenvalues become a complex conjugate pair when the dynamic pressure reaches the flutter boundary (kcr =96). The eigenvalue of mode (1, 1) will drop to zero because of the change in the buckle pattern when the dynamic pressure reaches the critical value (k=101). It is interesting to note that the eigenvalues do not merge for the plate in a chaotic condition when DT/DTcr = 2.2. If more fibers are distributed in the central portion of the plate, then this can efficiently increase the flutter boundary and critical buckling temperature. However, the aerothermoelastic behavior is dominated by mode shapes and becomes quite complicated. The stability boundary of the laminate with VFS does not simply

90

S.-Y. Kuo / Computational Materials Science 91 (2014) 83–90

(1) If more fibers are distributed in the central portion of the plate, then this can dramatically increase the buckling temperature. Accordingly, the magnitude of thermal postbuckling deflection can be reduced significantly. (2) If more fibers are distributed in the central portion of the plate, then this can efficiently increase the natural frequencies and flutter boundary. (3) The fiber distribution can change the stiffness, mass and thermal stress of the plate. Therefore, the sequence of the natural mode may be altered and the frequency coalescence pair may change. (4) The stability boundary of the laminate with VFS is not as the same as the one with uniform fiber spacing. Some particular phenomena are found, including two pairs coalesce simultaneously, coalescence between higher modes, and a closed section for chaotic motion.

Acknowledgement Fig. 14. Stability boundary (Vf = 1  f2 and Vf = (1  f2)2).

The financial support this study received from the National Science Council through Grant NSC 102-2221-E-156-001 is appreciated. References [1] [2] [3] [4] [5] [6] [7] [8] [9]

Fig. 15. Frequency coalescence plot (Vf = (1  f2) 2).

expand linearly due to the increases in critical buckling temperature and flutter boundary.

4. Conclusions A high precision triangular plate element is developed in this work for the aerothermoelastic analysis of cross-ply laminate with VFS. Based on the results, the following conclusions can be drawn:

[10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30]

A.F. Martin, A.W. Leissa, Int. J. Numer. Methods Eng. 28 (1989) 1813–1825. A.W. Leissa, A.F. Martin, Compos. Struct. 14 (1990) 339–357. A.W. Leissa, A.F. Martin, J. Elasticity 23 (1990) 97–112. L.C. Shiau, Y.H. Chue, Compos. Struct. 19 (1991) 107–205. L.C. Shiau, G.C. Lee, Compos. Struct. 24 (1993) 107–115. S.A. Meftah, R. Yeghnem, A. Tounsi, E.A. Adda, Constr. Build. Mater. 21 (2007) 1661–1671. S.A. Meftah, R. Yeghnem, A. Tounsi, E.A. Adda, Mater. Des. 29 (2008) 1955–1964. M. Bouremana, A. Tounsi, A. Kaci, I. Mechab, Mater. Des. 30 (2009) 2532–2537. M.A. Benatta, I. Mechab, A. Tounsi, E.A. Adda, Comput. Mater. Sci. 44 (2008) 765–773. S.Y. Kuo, L.C. Shiau, Compos. Struct. 90 (2009) 196–200. S.Y. Kuo, Compos. Struct. 93 (2011) 2533–2540. K.K. Raju, G.V. Rao, Compos. Struct. 12 (1989) 149–154. L.W. Chen, L.Y. Chen, Compos. Struct. 19 (1991) 267–283. C.A. Meyers, C.A. Hyer, J. Therm. Stresses 14 (1991) 519–540. G. Singh, G.V. Rao, N.G.R. Iyengar, AIAA J. 32 (1993) 1336–1338. L.C. Shiau, S.Y. Kuo, S.Y. Chang, J. Mech. 27 (2011) 559–566. T.Y. Yang, A.D. Han, AIAA J. 21 (1983) 758–766. J.E. Locke, J. Aircraft 31 (1994) 696–702. L. Librescu, W. Lin, M.P. Nemeth, J.H. Starnes Jr., J. Spacecraft Rockets 33 (1996) 285–291. L.C. Shiau, S.Y. Kuo, S.Y. Chang, Compos. Struct. 93 (2011) 2678–2684. A.D. Han, T.Y. Yang, AIAA J. 21 (1983) 1453–1461. D.Y. Xue, C. Mei, AIAA J. 31 (1993) 154–162. T. Prakash, M. Ganapathi, Compos. Struct. 72 (2006) 10–18. H.M. Navazi, H. Haddadpour, Compos. Struct. 80 (2007) 580–587. H.H. Ibrahim, M. Tawfik, M. Al-Ajmi, J. Aircraft 44 (2007) 1610–1618. K.J. Sohn, J.H. Kim, Compos. Struct. 88 (2009) 380–387. S.L. Lee, J.H. Kim, Compos. Struct. 92 (2010) 422–429. M.A. Kouchakzadeh, M. Rasekh, H. Haddadpour, Compos. Struct. 92 (2010) 2906–2915. S.Y. Kuo, L.C. Shiau, C.H. Lai, Smart Mater. Struct. 21 (2012) 35020–35030. Z.G. Song, F.M. Li, Compos. Struct. 106 (2013) 653–660.