An isogeometric finite element approach for thermal bending and buckling analyses of laminated composite plates

Accepted Manuscript Isogeometric finite element approach for thermal bending and buckling analyses of laminated composite plate Loc V. Tran, M. Abdel Wahab, Seung-Eock Kim PII: DOI: Reference:

S0263-8223(17)30714-6 http://dx.doi.org/10.1016/j.compstruct.2017.07.056 COST 8711

To appear in:

Composite Structures

Received Date: Revised Date: Accepted Date:

4 March 2017 9 July 2017 18 July 2017

Please cite this article as: Tran, L.V., Wahab, M.A., Kim, S-E., Isogeometric finite element approach for thermal bending and buckling analyses of laminated composite plate, Composite Structures (2017), doi: http://dx.doi.org/ 10.1016/j.compstruct.2017.07.056

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.

Isogeometric finite element approach for thermal bending and buckling analyses of laminated composite plate Loc V. Tran1, M. Abdel Wahab2,3, Seung-Eock Kim1* 1

Department of Civil and Environmental Engineering, Sejong University, 209 Neungdong-ro, Gwangjin-gu, Seoul 05006, Republic of Korea

2

Institute of Research and Development, Duy Tan University, 03 Quang Trung, Da Nang, Vietnam

3

Department of Mechanical Construction and Production, Faculty of Engineering and Architecture, Ghent University, 9000, Ghent – Belgium

Abstract Temperature rise in a plate structure produces non-zero transverse normal strain. Thus, a six-variable quasi-3D model with one additional variable in transverse displacement of higher-order shear deformation theory (HSDT) is developed in this paper to take into account the effects of transverse shears and normal strain in a laminated composite plate. The governing equation is discretized by isogeometric analysis (IGA), which naturally fulfills the C1-continuity requirement of the plate model. Due to the presence of bending-extension coupling, two kinds of thermal plate issues are considered – thermal buckling and thermal bending phenomena. Several numerical examples are provided to show the accuracy of the present method compared to reference results. Furthermore, it has been confirmed that the transverse normal strain cannot be discarded, especially for thick plates under a temperature environment. Keywords: Laminated composite plate, isogeometric analysis, thermo-elastic, higherorder shear deformation theory (HSDT) and quasi-3D model.

*

Corresponding author. Tel.: +82 2 3408 3291; fax: +82 2 3408 3332.

Email address: [email protected] (Loc V. Tran), [email protected](S.-E. Kim).

1

1

Introduction Laminated composite plates are made by stacking several lamina layers in which fiber

orientation follows a predetermined direction [1] in order to provide efficient material property, e.g., durability, high strength-to-weight and stiffness-to-weight ratios. Thus, they are commonly used in many weight-sensitive structures such as aerospace structures, aircraft structures, high-speed vehicle frames, etc. During operation, friction between the vehicle surface and the atmosphere produces thermal loading. It leads to harmful effects due to stress concentration, cracking and debonding at the interface between two distinct layers, or buckling resulting in large deformations and reduction of structural stiffness [2, 3]. For this reason, the thermal effect has become a necessary consideration in design and analysis. Numerous researchers have reported this problem. For instance, Noor and Burton [4-6] used threedimensional (3D) elasticity to investigate the thermal buckling response of laminated composite and sandwich plates. It is well known that such an exact 3D approach has the most potential to obtain exact solutions. However, it is unable to solve practical problems with complex (or even slightly complicated) geometries under arbitrary boundary conditions and takes high computational costs because of using 3D elements in modelling. In order to reduce the computational cost yet reach acceptable accuracy, the 3D theory can be simplified to 2D plate theories by making suitable assumptions based on an equivalent single layer (ESL) theory, which is classified as the following: the classical laminate plate theory (CLPT) [7, 8], the first order shear deformation theory (FSDT) [9-11], and the higher order shear deformable theory (HSDT) [12-14]. Among them, CLPT merely provides acceptable results for thin plates, while FSDT is suitable for moderate and thick plate. However, FSDT describes incorrectly shear strains/stresses by a constant through plate thickness. It is hence required to amend the unrealistic shear strain energy by the shear correction factors (SCFs), which are dispersed through each problem and not easily predicted in an accurate manner [15]. To eliminate these inefficiencies from CLPT and FSDT, HSDTs including higher-order terms in-plane displacements [16-19] are adopted. They gain better results and a more exact shear stress distribution as a nonlinear curve with zero values at the lateral surfaces. As known, these plate theories are due to the assumptions of σ z = 0 and ε z = 0 . However, temperature rise through the plate thickness, T(z) produces an associated non-zero transverse normal

2

strain, ε z = α zT ( z) regardless of whether T(z) distributes constantly, linearly or

nonlinearly. Fazzolari and Carrera [20] figured out that accounting for the transverse normal deformation is mandatory to accurately determine the thermal buckling parameter of thick plates. In addition, these 2D plate models may exhibit some disadvantages, e.g., by violating condition of inter-laminate continuous shear stress [21-23], or producing misleading results for thick plates subjected to thermal loading [24]. To more accurately predict the local behavior, e.g., inter-laminar continuous shear stresses in multi-layered plates, layer-wise (LW) theories [25, 26] and zigzag (ZZ) theories [27] have been proposed. LW requires a large number of unknowns, which increase rapidly according to increase in number of layers, whereas ZZ with two additional in in-plane displacements improves results of mechanically loaded plates and also produces inaccurate solution for thermally loaded plates [28] because of eliminating the transverse normal strain. Quasi-3D theories have been therefore developed to remove the inconsistency completely in 2D plate models, in which the transverse displacement is expanded with higher-order variations in respect to the thickness ratio. As a result, the thickness stretching effect can be captured. Numerous quasi-3D theories have been exploited for composite plates based on the Carrera Unified Formulation (CUF) [29], such as Carrera et al. [30] with 15 unknown variables, Ganapathi and Makhecha [31] with 13 variables, Kant and Shiyekar [32] and Matsunaga [33] with 12 variables, Ferreira and his coworkers [34-37] using 9 variables, and so on. In another way, Zenkour [38] proposed a refined trigonometric quasi-3D theory (RPT) with four unknowns, in which the transverse displacement is split into bending and shear parts. The RPT produces acceptable results for isotropic and functionally graded plates [39, 40] but shows quite a considerable error for multilayered plates, as concluded by both Kant [41, 42] and Tran et al. [43]. Being different from the aforementioned works, the present quasi-3D model is obtained by enhancing HSDT [17] with one additional variable in the thickness direction. The idea is identical to that of Zenkour [44] except using a cubic distributed function instead. Consequently, this model considers shear deformation as well as thickness stretching effect. The inter-laminar shear stress continuity can be managed by deriving from the equilibrium instead of the constitutive equations. The number of displacement variables is equal to six but the present model requires C1-continuity in approximation, which sets obstacles in the standard finite element formulations.

3

Fortunately, Hughes et al. recently proposed an efficient numerical method – so-called isogeometric analysis [45, 46], which uses the same basis functions as B-Spline or NURBS in describing the geometry and constructing the finite approximation for analysis. As well known, the basis functions, which come from computer aid design (CAD) representation, possess higher-order smoothness and continuity, e.g., Cp-1continuity by using pth NURBS order. Up to now, IGA has been rapidly developed and successfully applied in various fields [47-52] as well as functionally graded plate under a thermal environment [53-55]. Nevertheless, according to the author’s opinion, there are limited studies on developing isogeometric finite element approach for the thermal behavior of composite plates. Therefore, the main motivation of this paper is to propose an efficient approach based on conjunction IGA with quasi-3D model to study laminated composite plates under thermal loading in the hope of filling this research gap. Two kinds of thermal problems are herein considered. Thermal stress occurs in multilayered composite plates due to the presence of bending-extension coupling because they experience an inconstant temperature field or are stacked unsymmetrically. Otherwise, true bifurcation buckling occurs in symmetric composite plates subjected to a uniform temperature rise. Firstly, the present formulation is verified by comparison with available references from literature. The significant effects of aspect ratio, length-to-thickness ratio, fiber orientation and boundary condition on the behavior of composite plates are then analyzed. 2

Quasi-3D model for laminate composite plate formulation

Let us consider a laminated composite plate with length L, width W, and total T

thickness h depicted in Figure 1. The displacements of an arbitrary point u = {u , v, w}

can be expressed through unknown variables, which are functions of x and y only, as follows u ( x, y, z ) = u0 − zw0, x + f ( z )φx

h  −h ≤z≤  v( x, y, z ) = v0 − zw0, x + f ( z )φ y ,  2 2  w( x, y, z ) = w0 + g ( z )φ z

(1)

where {u0 , v0 , w0 } denote the axial displacements at the middle-plane Ω and {φ x , φ y } are the rotations of the normal to middle-plane in the y-z, x-z planes, respectively. Herein, the

4

normal deformation effect is considered by introducing an additional term φz and function

g ( z ) = f ′( z ) , which describes the distribution of the thickness stretching effect through plate thickness. By eliminating these terms, the generalized shear deformation plate theory is obtained as in Ref. [17]. According to Carrera’s recommendation [56], “Plate theories with at least a quadratic transverse displacement field in the z direction are required to capture transverse normal strain caused by the linear distribution of temperature across the thickness”, the distributed function is chosen as a cubic function,

f ( z ) = z − 4 3h2 z 3 [18].

Figure 1. Laminated geometry with displacement components. According to the displacement field in Eq. (1), the compatibility relationship between strain and displacement is given by ε = [ε x ε y γ xy ]T = ε m + z κ1 + f ( z )κ 2

ε z = g ′( z )φ z

(2) T

γ = γ xz γ yz  = f ′( z )ε s + g ( z)κ s where  u0, x   w, xx   φ x, x        ε m =  v0, y  , κ1 = −  w, yy  , κ 2 =  φ y , y  , u0, y + v0, x   2 w, xy  φ x , y + φ y , x       

φx  φz , x  εs =   , κs =   φy  φz , y  It is noted that subscript comma denotes the partial derivative.

5

(3)

For a kth orthotropic elastic lamina, the constitutive equation in local coordinate is derived from Hooke’s law as (k )

σ 1  σ   2 τ 12    σ 3  τ 31    τ 23 

 Q11 Q12 Q  21 Q22 0 0 =  Q31 Q32  0 0  0  0

0 0

0 0

0 0

Q66 0

0

0 Q33

0

0 0

0  0  0   0  0   Q44 

0 Q55 0 0

(k)

(k)

ε1 − α1∆T  ε − α ∆T  2  2  γ 12    ε 3 − α 3 ∆T  γ 31    γ 23 

(4)

where Qij are the elastic coefficients in the material coordinate [1], αi is the thermal coefficient of expansion in the principal i-direction, and ∆ T is the temperature change. Using the transformation formula, the constitutive equations for the kth layer with fiber orientation θ in the structural reference system (x,y,z) can be described as (k )

σ x    σ y  τ  xy   σ  z τ   xz  τ yz 

 Q11 Q12  Q21 Q22 Q Q26 =  16  Q31 Q32  0  0  0 0

Q16 Q13

0

Q26 Q23

0

Q66 Q36

0

Q36 Q33

0

0 0

0 0

0   0  0   0   Q45  Q44 

Q55 Q45

(k )

(k )

ε x − α x ∆T    ε y − α y ∆T  γ − α ∆T  xy xy   − ε α z ∆T   z γ   xz  γ yz 

(5)

in which the transformed material constants are given as below [1] Q11   c 4    Q12   c 2 s 2 Q   4 s  22   = 3 c Q16   s Q   cs3  26   2 2 Q66   c s Q44  Q45  Q55 

  Q   11  2 2 4 2 2  Q12  2c s c 4c s   , cs( s 2 − c 2 ) −cs 3 2cs( s 2 − c 2 )  Q22  cs(c 2 − s 2 ) − c3 s 2cs(c 2 − s 2 )  Q66   −2c 2 s 2 c2 s2 (c 2 − s 2 ) 2  2c 2 s 2 c4 + s 4

s4 c2 s2

Q13   c 2 s 2     Q44 −Q36  =  − cs cs   Q Q23   s 2 c 2   55

4c 2 s 2 −4c 2 s 2

Q13  , Q23 

6

Q33 = Q33

(6)

α x   c 2 s2     2  α1  c2    , α z = α3 α y  =  s    2cs −2cs  α 2   α xy  

(7)

where c = cosθ and s = sin θ . In case of HSDT, in which the transverse normal stress and strain are neglected, the stress-strain relation in Eq. (5) is reduced to (k )

σ x    σ y    τ xy    τ xz  τ yz   

 Q '11 Q '12   Q '21 Q '22 =  Q '61 Q '62  0  0  0 0 

Q '16 Q '26

0 0

Q '66 0 0

0 Q '55 Q '45

0   0  0   Q '54  Q '44 

(k )

ε x − α x ∆T    ε y − α y ∆T    γ xy − α xy ∆T    γ xz  γ yz   

(k )

(8)

where Q 'ij are the plane stress-reduced stiffness coefficients mentioned in the aforementioned works [22]. The stress resultants are defined as follows

 Nij   1  h/ 2      M ij  = ∫− h /2 σ ij  z  dz i, j := x, y    f (z)    Pij  h/ 2 Qα z   f ′( z )    = ∫− h /2 τ α z   dz g (z)   Rα z 

α := x, y

and

(9)

Nz = ∫

h/ 2

− h /2

g ′( z )σ z dz

By choosing f ( z ) = z − 4 3h2 z 3 , both functions f ′( z) and g ( z) get a zero value at

z = ± h / 2 . It means that the traction-free boundary condition is automatically satisfied at the top and bottom plate surfaces. Substituting Eq. (5) into Eq. (9), the stress resultants can be written as N   A M   B   σˆ =   =  P   E  N z   AT3

B D

E F

F

H

BT3

E3T

A3  ε m   Nth    B 3   κ1  M th  ˆ ˆ −     = Dεˆ − N th E3   κ 2   Pth   H 3  φ z  Oth 

7

(10)

Q   A s τˆ =   =  s R   B

B s  ε s  ˆ    = Ds εˆ s Ds   κ s 

(11)

in which h/ 2

( A , B , D , E , F , H ) = ∫ (1, z, z , f ( z ), zf ( z ), f ij

ij

ij

ij

ij

ij

2

− h /2

h /2

( A3i , B3i , E3i ) = ∫− h/ 2 (1, z, f ( z ) ) g ′( z )Q3i dz and ( Aijs , Bijs , Dijs ) = ∫

h /2

( z) ) Qij dz i, j := 1,2,6

i := 1,2,6 , H 3 = ∫

h/ 2

−h/ 2

([ f ′( z)] , f ′( z) g ( z),[ g ( z)] ) Q dz 2

− h/ 2

2

2

ij

( g ′( z ) )

2

Q33 dz

(12)

i, j := 4,5

And thermal stress resultants are expressed as

{N th

α x  α  h /2  y Pth } = ∫ Q  {1 z − h/ 2 α xy  α   z 

M th

Oth = ∫

h/ 2

− h /2

g ′( z )Q 3i α x

f ( z )} ∆Tdz

(13) T

α y α xy α z  ∆Tdz

The virtual strain energy is defined as

δ U ε = ∫ ( N T δε m + M T δκ1 + P T δ κ 2 + N z δφ z + QT δ ε s + R T δ κ s ) dΩ Ω

(14)

By substituting Eq.(9) into Eq.(14), the variation of strain energy can be written as

(

)

T ˆ δ εˆ + εˆ T D ˆ δεˆ dΩ − N δU ε = ∫ εˆT D s s s ∫ ˆ thδ εˆdΩ Ω

Ω

(15)

Under only thermal loading, the total potential energy must be equal to the strain energy.

δΠ = δU ε = 0 3 3.1

(16)

Isogeometric finite element for plate formulation The basic function A knot vector Ξ = {ξ1 , ξ 2 ,..., ξ m + p +1 } is a non-decreasing sequence of parameter values

ξi , i = 1,..., m + p , in which ξi ∈ R + is the ith knot, p is the polynomial order and m is

8

number of the basis functions. In an open knot, the first and the last knots are repeated p+1 times and often get normalized values of ξ1 = 0 and ξ m + p +1 = 1 . Furthermore, the knot vector is divided into m-p of inner knot spans, some of which may possibly have zero length if a knot value appears more than once. The B-splines basis function is C∞ continuous inside a knot span and Cp-k continuous at a knot repeated k times. As a result, it achieves Cp-1-continuity at a single knot and C-1-continuity at the first and last knots. Using the Cox-de Boor algorithm, the univariate B-spline basis functions Nip (ξ ) are defined recursively [57] on the corresponding knot vector with the unit step function at the beginning. N i p (ξ ) =

ξi + p +1 − ξ ξ − ξi N p −1 (ξ ) + N p −1 (ξ ) ξi + p − ξi i ξ i + p +1 − ξi +1 i +1

p ≥1

(17) at p = 0,

1 if ξ i < ξ < ξi +1 Ni0 (ξ ) =  otherwise 0

Note that when evaluating these functions, the ratio in the form 0/0 is assumed to be zero. Their kth derivative thereafter is defined as  dkk−−11 Nip −1(ξ ) d kk−−11 N ip+−11 ( ξ )  dk p  N i (ξ ) = p  dξ − dξ  ξi + p − ξ i dξ k ξi + p +1 − ξ i +1   

(18)

An example of univariate B-spline basis functions based on an open, uniform knot vector Ξ = {0, 0,0, 0, 13 , 23 ,1,1,1,1} is illustrated in Figure 2. It can be seen that the cubic functions are C2 continuous and their third derivatives are constant values in each knot interval. Thus, they naturally satisfy the C1 -continuity requirement of the plate model and have the ability to produce the third order derivatives required for calculating shear stress in subsection 3.5.

9

10

1 N 31

N 36 N 33

N 32

0.6

5

N 35

N 34

dN36

dN32

0.8

dN34 0

0.4

dN33 -5

0.2

dN35

dN31

0

-10 0

0.2

0.4

0.6

0.8

1

0

0.2

a) Cubic basis functions

0.4

0.6

0.8

1

b) First derivative

60 d2N 36

d2N 31

40

300 3

d3N2 d3N3

200

6

20

3

d3N3

100

0

3

d3N34

d3N5

d2N 33

-40

d2N 34

-100

d3N33

-200

d3N31

3 4

3

d3N3

d3N32

d3N33

d3N6

0

-20

d3N

2

d3N3 4

-60

d2N 32

d2N 35

3

d3N5

-80 -300

0

0.2

0.4

0.6

0.8

1

c) Second derivative

0

0.2

0.4

0.6

0.8

1

d) Third derivative

Figure 2. Basis functions and their derivatives corresponding to Ξ = {0, 0, 0, 0, 13 , 23 ,1,1,1,1} . By utilizing a tensor product of univariate B-splines, multivariate B-spline basis functions are generated d

( )

N ip ( ξ ) = ∏ N iαpα ξα

(19)

α =1

where parametric d = 1, 2, 3 represents 1D, 2D and 3D spaces, respectively. Figure 3 gives an illustration of bivariate B-splines basic based on the tensor product of two knot vectors Ξ = {0,0, 0, 15 , 52 , 53 , 53 , 54 ,1,1,1} and Ψ = {0, 0, 0, 0, 13 , 13 , 32 , 32 ,1,1,1,1} according to the parametric dimensions ξ and η.

10

Figure 3. Bivariate B-spline basis functions. B-spline is convenient for modelling geometries with curved boundaries but it is unable to exactly represent some conic shapes, e.g., circles, ellipses, and spheres. Thus, NURBS is introduced as a generalization of B-spline in the form of rational functions, i.e.

N ip ( ξ ) ζ i R (ξ ) = ∑ N jp ( ξ ) ζ j p i

(20)

j

where ζ i > 0 is the individual weight corresponding to B-splines basis functions N ip (ξ ) . It is observed that with a constant weight value, NURBS will transform to a B-spline. 3.2

Geometric description By linear combination of the basis function and the control points Pi , geometry is

constructed as S ( ξ ) = ∑ N ip ( ξ ) Pi i

(21)

Figure 4 reveals a comparison of constructing a 180 degree arc with two elements by quadratic NURBS and the Lagrange basis function. It is observed that when using four control points at the red dots, NURBS describes geometry exactly while Lagrange interpolation functions using five nodes just produce an approximated geometry. Furthermore, at the inter-element boundary (the arc center), different from the approximated geometry by Lagrange function, NURBS geometry is smooth due to higher continuity of the basis function.

11

1.05 1 0.95 −0.2

−0.1

0

0.1

0.2

x

Figure 4. Model of a half circle with two elements by quadratic NURBS basis (red curve) and Lagrange quadratic function (blue curve), in which red dots are the control points and blue squares are the nodal points. For 2D geometry, as shown in Figure 5, a circle is represented by only one quadratic element based on a bivariate NURBS constructed from a tensor product of two knot vectors Ξ = {0,0,0,1,1,1} and Ψ = {0,0,0,1,1,1} and nine control points [16].

Figure 5.Circular geometry is described in the coarsest mesh by one quadratic NUBRS with 9 control points denoted by red dots.

12

3.3

Discrete system equation Similar to the traditional finite element method, isogeometric analysis also invokes the

isoparametric concept, in which NURBS from the geometric description is utilized to construct a finite approximation of the displacements u h ( ξ ) = ∑ RA ( ξ )q A

(22)

A

where q A denotes the vector of nodal degrees of freedom associated with the control T

point PA , which has 5 degrees of freedom (DoFs) q A = u0 A v0 A w0 A φxA φ yA  for HSDT T

and 6 DoFs q A = u0 A v0 A w0 A φxA φ yA φzA  for the quasi-3D model. By substituting Eq. (22) into Eq.(2), the in-plane and shear strains can be rewritten as: εˆ T = ∑ ( B mA )  A=1 

T

b1 T A

b2 T A

b3 T A

(B ) (B ) (B )

εˆ Ts = ∑  ( BsA0 )  A=1 

T

T

 q = Bˆ q b  A (23)

T

T ( BsA1 )  q A = Bˆ sq

in which  RA, x  B = 0  RA, y  m A

B

b2 A

0

RA, y RA, x

0 0 0 0  0 0 0 0 , 0 0 0 0 

 0 0 0 RA, x  0 = 0 0 0  0 0 0 RA, y 

0 0 0 B sA0 =  0 0 0

RA 0

0

RA, y RA, x 0 RA

 0 0 RA, xx  B = −  0 0 RA, yy  0 0 2 RA, xy  b1 A

0  0 , 0 

0 , 0 

0 0 0  0 0 0 0 0 0 

BbA3 = [ 0 0 0 0 0 RA ]

(24)

 0 0 0 0 0 RA, x  B sA1 =    0 0 0 0 0 RA, y 

In case of HSDT, strain matrices have five columns and are detailed in [23, 54]. By substituting Eq. (23) into Eq.(16), we obtain the system of linear static equation as follows

Kq = Fth

(25)

where the global stiffness matrix K is given by

13

ˆ ˆ +B ˆTD ˆ B ˆ K = ∫ Bˆ Tb DB b s s s dΩ

(26)

Ω

and Fth is thermal load vector. ˆ dΩ Fth = ∫ Bˆ Tb N th

(27)

Ω

3.4

Thermal buckling analysis Normally, temperature rise produces thermal moment resultants, which make the

plates deflect. However, these thermal moments can vanish in some special cases. For example, due to the symmetrically mid-plane configuration of the plate, the initial moment under uniformly temperature rise in Eq. (13) equals to zero. In another case, that is clamped edges, the supports are capable of handling the extra moments [58]. Thus, the initially perfect plate is still flat with no transverse deflection. In these cases, thermal inplane compression increases according to temperature rise. As thermal force reaches a predetermined magnitude, the structure suddenly experiences a large deflection. This is the so-called a thermal buckling phenomenon. At the adjacent equilibrium state, the

N virtual work done by the produced thermal force N 0 =  x  N xy

δV =

N xy   is N y 

T 1  w, x w, y  N 0 δ w, x δ w, y  dΩ ∫ Ω 2

(28)

Hence, the total potential energy is given by

δΠ = δ U ε + δ V = 0

(29)

Following the above assembled procedure, the stability equation is expressed as

(K + λK g )q = 0

(30)

where K g is the initial stress stiffness matrix T

K g = ∫ ( B g ) N 0 B g dΩ Ω

in which

14

(31)

 0 0 RA , x B gA =   0 0 RA , y

g (0) RA, x  0 0 g (0) RA, y  0 0

(32)

Herein, a value of g(z) at z=0 means that one Gaussian point is used for computing through the plate thickness. The thermal load factor λ is determined as the smallest eigenvalue. The critical buckling temperature is then calculated by

Tcr = λ∆T 3.5

(33)

Shear stresses Concerning the calculation of the transverse shear stresses, it is pointed out that the

traction free condition at the lateral surfaces and the continuous shear stress condition at the layer interfaces are usually violated in plate models, including FSDT [59], HSDT [22, 43] and quasi-3D. Therefore, to satisfy these conditions, the transverse shear stresses are derived from the equilibrium equations and integrated the equations through the thickness as follows

τ xz = − ∫ (σ x , x + τ xy , y ) dz h

(34)

τ yz = − ∫ (σ y , y + τ xy , x ) dz h

4

Numerical results This section focuses on the study of thermal behaviors of laminated composite plates,

in which we assume that the plate is constrained on all edges by: - Simply supported (SS) condition is divided into two cases: movable and immovable in the in-plane directions. Movable edge (SS1):

v0 = w0 = φ y = φ z = 0  u0 = w0 = φ x = φz = 0

Immovable edge (SS2):

u0 = v0 = w0 = φ y = φ z = 0  u0 = v0 = w0 = φ x = φz = 0

- Clamped supported condition

15

on x = 0, L on y = 0, W

on x = 0, L on y = 0, W

(35)

(36)

u =v =w=0

(37)

on all edges

Some sets of the material property are used for numerical investigations: - Material I [28, 60]

EL / ET = 25, GLT / ET = 0.5, GTT / ET = 0.2,

ν LT = ν TT = 0.25, αT / α L = 1125, αT / α 0 = 1 - Material II [4, 14]

EL / ET = 15, GLT / ET = 0.5, GTT / ET = 0.3356,

ν LT = 0.3, ν TT = 0.49, α L / α0 = 0.015, αT / α 0 = 1 - Material III E2 = 200GPa , E1 = 2.45 E2 , G12 = G13 = 0.48 E2 , G23 = 0.2 E2 ,

ν 12 = 0.23, ν 13 = ν 23 = 0.25, α1 = α 2 = α 3 = 10 −6 / C where L and T refer to the directions parallel and perpendicular to the fibers and α0 is the normalized factor for the thermal expansion coefficient. Depending on the value of the load vector, plate problems are categorized into two main groups: thermal bending and thermal buckling phenomena. Numerical examples are solved by using cubic NURBS with full ( p + 1) × ( q + 1) integrated Gaussian points. 4.1

Thermal bending It is assumed that the laminated composite plate is subjected to a steady-state

temperature rise. Due to the same thermal conductivity of all layers in the transverse direction, the temperature varies linearly though the plate thickness. For convenience, the obtained results of displacement and stresses are presented following the normalized form (u , v ) =

4.1.1

1

Lα 0T1

(u , v), w =

h 2

L α 0T1

w, σ ij =

1

ELα 0T1

σ ij

Square composite plate with lower thermal anisotropy subjected to gradient temperature rise

Let us first consider a simply supported (SS1) square plate undergone a uniformly distributed temperature (UDT) ∆T ( x, y , z ) = Th1 z

16

or a sinusoidally one (SDT)

∆T ( x, y, z ) = Th1 z sin( πL x)sin( Wπ y ) . Material I is set with a lower level of thermal

anisotropy αT / α L = 3 . Table 1 presents the effect of the length-to-thickness ratio on nondimensional thermal deflection of orthotropic and four-layered composite plates. Herein, the obtained results based on IGA in combination with four-variable RPT [39], fivevariable HSDT [17] and six-variable quasi-3D model are compared to available ones by Cetkovic [61] using layer wise (LW), FSDT and CLPT. It can be seen that HSDT (εz=0) and present quasi-3D (εz ≠ 0) achieve good agreement of deflection for all cases. Meanwhile, RPT, which also considers the transverse normal deformation (εz ≠0), significantly underestimates the results with differences of over 15% in the case of L/h = 10. The inaccurate prediction of RPT when studying the multilayered plates was pointed out in the literature [42, 43]. Thus, in the next examples we prefer to use HSDT and quasi-3D models. It is also observed that an increase in the length-to-thickness ratio L/h seems to have an insignificant effect on the transverse displacement. Table 1: Thermal central deflection w∗ = 10 w  L2 ,W2 ,± h2  of orthotropic and four-layered plates with small thermal anisotropy αT / α L = 3 under gradient temperature distribution.

L/h 10 Orthotropic SDL RPT 0.8904 HSDT 1.0438 IGA Quasi-3D 1.0785 LW [61] 1.0420 FSDT 1.0440 CLPT 1.0312 Laminate [0/90/90/0] RPT 0.8904 IGA HSDT 1.0455 Quasi-3D 1.0806 LW [61] 1.0778 FSDT 1.0421 CLPT 1.0312

20

100

20

100

1.0293 1.0313 1.0316 1.0313 1.0313 1.0312

10 UDL 1.2504 1.4605 1.4963 1.4568 1.4603 1.4331

0.9881 1.0346 1.0433 1.0340 1.0346 1.0312

1.3830 1.4412 1.4503 1.4401 1.4409 1.4331

1.4315 1.4338 1.4342 1.4338 1.4334 1.4331

0.9881 1.0353 1.0440 1.0708 1.0343 1.0312

1.0293 1.0313 1.0317 1.068 1.0313 1.0312

1.3604 1.5511 1.5849 1.6013 1.5452 1.5316

1.4863 1.5380 1.5463 1.5912 1.5357 1.5316

1.5300 1.5321 1.5325 1.586 1.5318 1.5316

17

4.1.2

Three layer [0°/90°/0°] composite plate under a sinusoidally distributed temperature

Next, thermal elastic bending of a three layer [0°/90 °/0°] square plate subjected to gradient temperature ∆T ( x, y, z ) =

2T1 h

z sin( πL x)sin( Wπ y ) is investigated. The material

property is the same as that of the previous example except that higher level of thermal anisotropy αT / α L = 1125 is set. The normalized displacements and stresses are tabulated in Table 2, in which the present results are compared with 3D exact solution [62], that using layer-wise (LW) model [61], higher order shear and normal deformation theory (HOSNT) with 12 DoFs [32] and zigzag (ZZ) theory [28] for thick (L/h = 4) to thin (L/h = 100) plates. Being different from the above example, for highly anisotropic plates, the transverse displacement is sensitive to changes in the length-to-thickness ratio and there is considerable difference between the obtained results. Ignoring the transverse normal strain effect, these plate theories: HSDT, ZZ, LW produce inaccurate results compared to the 3D exact solution, especially for thick plates. For instance, for L/h = 4 the error of udisplacement is up to 20% for HSDT, 21% for ZZ, and 10% for LW, while that of transverse deflection is 39.6%, 38.6%, and 28.1% using HSDT, ZZ, and LW, respectively. This issue motivated us to develop six-variable quasi-3D model to study composite plates subjected to thermal loading. It is worth noting that using six unknowns in total, a half variables compared to HOSNT, the present quasi-3D model gains accurate results, especially for the displacements as shown in Figure 6. As seen, non-linear variation is observed for the in-plane displacement especially for thick laminated plates. Furthermore, the transverse displacement varies through the plate thickness as a curve and is identical with 3D exact solution. Calculating following to the equilibrium equations [63] instead of using the constitutive relation, the obtained shear stresses satisfy the traction free condition at the lateral surfaces and the inter-laminate continuity condition at the layer interfaces as shown in Figure 7. The effect of the L/h ratio on the axial and shear stress of three-layer laminate plate is revealed in Figure 8. Table 2: Thermal displacements and stresses of three-layer [0 °/90°/0°] square plate for different L/h ratios and thermal loading. Method εz= 0 IGA εz ≠ 0

L/h

u (± h2 ) *

v (± h2 )

w( ± h2 )

σ x (± h2 )

σ y ( ± h2 )

τ xz (± h6 )

τ yz ( ± h6 )

4

14.413 18.671

62.348 66.976

25.781 42.237

869.017 1204.322

920.332 900.641

90.908 80.438

-137.678 -138.244

18

LW(εz= 0) [61] HOSNT [32] ZZ(εz= 0) [28] 3D exact [62] εz= 0 IGA εz ≠ 0 LW(εz= 0) [61] HOSNT[32] ZZ(εz= 0) [28] 3D exact [62] εz= 0 IGA εz ≠ 0 LW(εz= 0) [61] HOSNT [32] ZZ(εz= 0) [28] 3D exact [62] εz= 0 IGA εz ≠ 0 LW(εz= 0) [61] HOSNT[32] ZZ(εz= 0) [28] 3D exact [62] εz= 0 IGA εz ≠ 0 HOSNT[32] 3D exact [62] CLPT [62]

10

20

50

100

16.272 18.178 14.29 18.11 15.539 16.252 15.732 15.539 15.90 16.610 15.834 16.014 15.919 16.128 15.99 16.170 15.928 15.957 15.964 16.00 15.99 16.020 15.942 15.949 15.983 16.000 15.99

88.962 75.9776 77.43 81.83 27.796 28.525 32.837 27.796 31.23 31.950 19.185 19.366 20.292 19.8194 20.16 20.340 16.478 16.507 16.405 16.6067 16.68 16.710 16.080 16.087 16.135 16.170 15.99

30.683 42.02 26.22 42.69 14.111 16.848 16.855 14.111 14.65 17.390 11.252 11.941 11.447 11.953 11.43 12.120 10.356 10.466 10.389 10.468 10.39 10.500 10.224 10.252 10.244 10.260 10.18

792.572 1189.39 879 1183 935.724 991.859 950.314 935.724 969.70 1026.000 953.497 967.621 968.381 979.754 967.50 982.000 959.155 961.388 972.764 967.113 965.20 965.256 959.996 960.524 965.256 965.400 964.60

904.129 871.19 873 856.1 1028.313 1025.144 1019.610 1028.313 1017 1014 1055.215 1054.448 1048.400 1052.76 1052 1051 1063.671 1063.575 1048.400 1063.14 1063 1063 1064.916 1064.917 1065.000 1065 1065

97.007 88.618 87.97 84.81 61.435 60.776 64.923 61.435 60.88 60.54 33.820 33.742 36.008 34.174 34.03 33.97 13.919 13.913 14.858 14.08 14.09 14.08 6.991 6.990 7.075 7.073 -

*z value is given in parentheses while (x, y) values are (0, W/2) for u , τ xz ; (L/2, 0) for v , τ yz and (L/2, W/2) for w , σ x , σ y , respectively.

19

-136.27 -129 -128.7 -66.700 -66.755 -66.700 -65.99 -66.010 -34.802 -34.811 -34.862 -34.76 -34.760 -14.103 -14.104 -14.13 -14.13 14.130 -7.065 -7.065 -7.080 -7.080 -

0.5

0.5 L/h

0.3

3D exact

0.2

Quasi 3D

0.1

HSDT

20 4 4 20 4 20

Thickness, z/h

Thickness, z/h

0.4

0 -0.1

0

-0.2 -0.3 -0.4 -0.5 -20

-15

-10

-5

0

5

10

15

-0.5 10

20

15

20

Inplane displacement,¯ u(0, W/2)

25

30

35

40

45

Displacement,w(L/2, ¯ W/2)

Figure 6. Comparison of in-plane displacement u (left) and transverse deflection w (right) in three-layer [0°/90°/0°] plates subjected to gradient temperature with 3D elastic solution [62].

0.5

0.5

3D, L/h=4 3D, L/h=20

Thickness, z/h

Thickness, z/h

L/h=4

L/h=20

3D, L/h=4 3D, L/h=20

0

L/h=20 0

-0.5 -140

L/h=4

-120

-100

-80

-60

-40

-20

0

Shear stress,¯ τyz (0, W/2)

-0.5 0

20

40

60

80

100

120

140

Shear stress,¯ τxz (L/2, 0) W Figure 7. Shear stress distribution of τ xz ( L2 ,0) (left) and τ yz (0, 2 ) (right) obtained from 3D elastic

solution [62], equilibrium equation (solid line), and constitutive relation (dash line).

20

0.5

L/h=4 L/h=10 L/h=20 L/h=50 L/h=4 L/h=10 L/h=20 L/h=50

Quasi 3D

0

HSDT

-0.5 -1500

-1000

-500

0

Thickness, z/h

Thickness, z/h

0.5

500

1000

0

-0.5 -2000

1500

-1000

Axial stress,¯ σx (L/2, W/2)

0

1000

2000

Axial stress,¯ σy (L/2, W/2)

0.5

0.5 0.4 0.3

Thickness, z/h

Quasi 3D

0.1 0

HSDT

-0.1 -0.2

Thickness, z/h

L/h

0.2

4 10 20 50 4 10 50 20

0

-0.3 -0.4 -0.5 0

20

40

60

80

100

120

-0.5 -140

-120

-100

-80

-60

-40

-20

0

Shear stress,¯τyz (0, W/2)

Shear stress,¯ τxz (L/2, 0)

Figure 8. Stress distributions of three-layer [0°/90°/0°] plate subjected to gradient temperature via length-to-thickness ratio L/h.

4.2

Thermal buckling For the symmetric laminated plate subjected to uniform temperature rise, only

membrane forces are generated. Thus, bifurcation buckling physically occurs. Calculation of the buckling temperature requires two steps. Firstly, the linear static equation (25) is solved to determine the thermal in-plane resultants following Eq. (10). These stress resultants are then used to compute the initial stress stiffness matrix K g , which is subsequently used in the eigenvalue equation (30) to find the thermal load factor λ. In all

21

the examples in this subsection, material II is set for each lamina, unless indicated otherwise. The critical buckling temperature is normalized as Tcr = α 0Tcr . 4.2.1

Single layer simply supported plate

Let us consider a one layer orthotropic (material II) or isotropic ( E = 106 ,ν = 0.3 ,

α 0 = 10−6 ) square plate. The plate with a length-to-thickness ratio L/h ranging from 4 to 100 is under simply supported (SS2) condition and uniform temperature distribution. Table 3 summarizes the obtained results in comparison with 3D elasticity solution proposed by Noor and Burton [4], the analytical solutions by Matsunaga [33] and Cetkovic [64], and those of FEM [14] and mesh-free method (MM) [65]. Again, quasi3D model provides accurate results, while HSDT overestimates the values of critical temperature as the plate becomes thicker. The error percentages compared with 3D exact increase from 1.25% to 4.53% for isotropic plate, and 2.27% to 4.79% for orthotropic plate, respectively, in the range of 10 to 4 for L/h ratio. Figure 9 reveals the mode shape at buckling state for each type of plate. As also concluded by Matsunaga [33], the isotropic plate is buckled in the form of a half sine wave while mode shape of the orthotropic plate is a full sine wave shape in direction of the reinforced fiber. Table 3: Comparison of critical temperature for single layer plates 3D exact [4, 33] Isotropic L/h

Quasi-3D [14] 9 DoFs 11 DoFs

4

5.60e-2

-

-

5

3.99e-2

-

-

10

1.18e-2

-

-

20

3.11e-3

-

-

100

1.26e-4

-

-

-

-

0.1417 -1.32 5.78e-2

0.1441 0.35 5.82e-2

LW [64]

MM[65] HSDT

εz= 0

IGA

εz ≠ 0

5.85e-2 4.52 4.14e-2 3.71 1.20e-2 1.27 3.12e-3 0.35 1.27e-4 0.08

5.86e-2 4.64 4.13e-2 3.61 1.20e-2 1.27 3.12e-3 0.35 1.26e-4 -0.63

5.854e-2 4.53 4.132e-2 3.55 1.198e-2 1.25 3.119e-3 0.33 1.265e-4 0.05

5.737e-2 2.44 4.074e-2 2.10 1.193e-2 0.85 3.116e-3 0.23 1.265e-4 0.05

0.1868 5.12 0.1504 4.74 5.92e-2

0.1878 5.68 0.1506 4.87 5.92e-2

0.1862 4.79 0.1499 4.40 5.913e-2

0.1785 0.44 0.1452 1.14 5.848e-2

Orthotropic 4

0.1777

5

0.1436

10

5.78e-2

22

-0.07

0.61

20

1.739e-2

-

-

100

7.46e-4

7.47e-4 0.09

7.47e-4 0.11

2.44 1.75e-2 0.81 7.47e-4 0.04

2.35 1.75e-2 0.75 7.46e-4 0.00

2.27 1.752e-2 0.75 7.466e-4 0.04

1.15 1.747e-2 0.44 7.465e-4 0.03

Upper line: The critical temperature Tcr , lower line: the error as compared to 3D elastic solution [4]. Noted that ●e-i means ●×10-i .

Figure 9. Thermal buckling mode shapes of isotropic (left) and orthotropic (right) square plates with L/h=10.

4.2.2

Symmetric cross-ply laminated plate

In order to verify the accuracy of the present formulation, let us consider symmetric cross-ply three-layer [0°/90°/0°] and four-layer [0°/90°/90°/0°] square composite plates under immovable simply supported constraint. Table 4 presents the comparison of the obtained critical temperature of three-layer [0°/90°/0°] composite plates using present models and that from 3D [4], layer-wise [64], and HSDT [33, 65] models from the literature. As observed, quasi-3D, considering the thickness-stretch effect, gains closer results than HSDT model compared with 3D exact solution by Noor [4]. This conclusion is clearly observed for thick plates, e.g., L/h <10. Table 4: The critical temperature Tcr of a three-layer [0°/90°/0°] laminated plate via different L/h ratios. L/h LW[64] 3D exact[4] HSDT[33]

4 0.2272 0.214 0.2133

5 0.1848 0.1763 0.1752

10 0.7628e-1 0.7467e-1 0.7442e-1

23

20 0.2316e-1 0.2308e-1 0.2297e-1

100 0.9964e-3 0.9961e-3 0.9960e-3

HSDT[65] εz= 0 IGA εz ≠ 0

0.2253 0.2201 0.2160

0.1828 0.1800 0.1775

0.7433e-1 0.7540e-1 0.7505e-1

0.2308e-1 0.2307e-1 0.2304e-1

0.9917e-3 0.9915e-3 0.9915e-3

Next, the effect of boundary conditions (BCs) on the thermal critical parameter Tcr of the four-layer [0 °/90°/90°/0°] laminated plate with the two ratios of L/h = 10 and 100 is investigated. Table 5 illustrates the comparison of critical temperature obtained by cubic NURBS-based IGA and FEM-Q16 [13] with the same polynomial order. In this example, each plate edge is imposed by simply supported (S) or clamped (C) constraints. It is observed that the solution dramatically depends on the boundary conditions. Changing the BC from CCCC to SCSC, SSSS, or SFSF means that more constraints are released. As a consequence, the magnitude of temperature is decreased accordingly. Figure 10 plots the thermal buckling mode shape via the change of BC. Table 5: Effect of the BC on the critical temperature Tcr of a symmetric four-layer [0°/90°/90°/0°] laminated plate. L/h

100*

10

Boundary condition FEM-Q16 HSDT [13] FSDT εz= 0 IGA εz ≠ 0 FEM-Q16 HSDT FSDT [13] εz= 0 IGA εz ≠ 0

SSSS 0.0996 0.0997 0.0996 0.0996 0.0757 0.0770 0.0766 0.0762

CSCS 0.144 0.1441 0.1441 0.1444 0.1044 0.1069 0.1069 0.1066

SCSC 0.253 0.2532 0.2531 0.2536 0.1305 0.1344 0.1354 0.1345

* Result of Tcr for thinner plate (L/h=100) is scaled up by 100.

24

CCCC 0.3348 0.3352 0.3353 0.3364 0.1601 0.1655 0.1657 0.1646

SSSS

CSCS

SCSC

CCCC

Figure 10. Buckling mode shapes of the symmetric cross-ply [0°/90°/90°/0°] laminate plate under various boundary conditions with L/h = 10.

4.2.3

Symmetric angle-ply [θ/-θ/θ/(-θ/θ)n] laminated plate

After the validations, in this subsection, we extend the present formulation for a simply supported angle-ply [θ/-θ/θ/(-θ/θ)n ] composite plate (L/h=10) with a range of n from 0 to 3. Figure 11 plots the variations of critical buckling temperature via fiber orientation for the angle-ply plate under SS1, SS2 and fully clamped BCs. It is observed that the critical buckling temperature strongly depends on the value of the orientation angle θ and is symmetric with respect to θ=45°. With immovable BCs, an increase in θ leads to a monotonic increase in Tcr because for the angle-ply plates, the stiffness increases accordingly to increase in the orientation angle of fibers. However, there is a different observation in the case of movable BC. The critical temperature decreases as θ increases from 0° to around 20°-25° depending on number of layers and subsequently

25

increases according to increase in the orientation angle θ. The difference in the nature of

Tcr with respect to θ is caused by obtained in-plane stress resultants. In the case of immovable BCs, they remain uniform throughout the plate and equal to thermal stress resultants, whereas these stress resultants obtained in the SS1 plate exhibit a non-uniform distribution and increase to maximum value of around 20°-25° [10, 14]. Furthermore, increase in the number of lamina makes the plate become stiffer and it is therefore able to resist higher temperature before buckling. However, its effect seems to be insignificant when the number of lamina is over five as shown in Figure 12. Figure 13 further investigates the influence of the aspect ratio (L/W) and the orientation angle θ on the critical buckling temperature of a simply supported (SS2) plate with L/h = 10. It is seen that for a square plate, θ = 45° leads to highest thermal buckling value. However, for a rectangular plate, the highest value is given at θ = 60°. By increasing in the aspect ratio, a higher buckling temperature is obtained and more half-sine waves are revealed in the buckling mode shape as plotted in Figure 14.

26

0.18

n=0 n=1 0.16

n=2

Critical temperature

cr

(T )

n=3 0.14

0.12

0.10

0.08

0.06

0

15

30

45

60

75

90

Fiber orientation

0.250

n = 0 n = 1

0.225

n = 2

Critical temperature

cr

(T )

n=3 0.200

0.175

0.150

0.125

0.100 0

15

30

45

60

75

90

Fiber orientation

Figure 11. Variations of critical buckling temperature via fiber orientation for an angleply θ/-θ/θ/[- θ/θ]n plate under: (upper) immovable: simply supported SS2 (dash line) and fully clamped (solid line) and (lower) movable simply supported SS1 contraints.

27

25

Critical temperature

cr

(100T )

20

n = 0 n = 1

15

n = 2 n = 3 10

5

0 0

20

40

60

80

100

L/h ratio

Figure 12. Variation of critical temperature via length-to-thickness ratio, L/h for a fivelayer angle-ply plate with fiber orientation θ = 60°.

30

Critical temperature

cr

(100T )

60

o

45

25

o

75

20

90

o

o

30

o

15

15

10

0

o

o

5 1.0

1.5

2.0

2.5

3.0

Aspect ratio L/W

3.5

4.0

Figure 13. Effect of aspect ratioon critical temperature of a five-layer [θ/-θ/θ/-θ/θ] plate.

28

L/W = 1

L/W = 1.5

L/W = 2

L/W = 3

Figure 14. Increase in wavelength of buckling mode shape of a five-layer [60/-60/60/-60/60] plate via increase in aspect ratio, L/W.

4.2.4

Complicated shaped plate

In order to illustrate the ability of IGA in modelling complicated shapes, a square plate with a heart cutout is studied. As seen in Figure 15a, the complicated shape plate can be described exactly in a coarse mesh with four NURBS patches ( 4 × 1×1 elements). By applying the refinement technique, we can obtain a finer mesh with 4 × 5 × 5 cubic elements (Figure 15b). Firstly, to validate the present model, free vibration of an isotropic plate with material and geometrical parameters following to Liu’s work [66] is studied. This problem is solved by eigenvalue analysis in the same way as the buckling problem, Eq. (30). Table 6 shows the normalized fundamental frequencies ωˆ =(ω 2 ρ hL4 / D )1/4 , where D is the flexural rigidity. Herein, the bending strip method [67] is adopted to maintain C1-continuity at the patch interfaces. With this treatment, the present model enhances the solutions based on IGA in combination with four-unknown shear and

29

normal deformation theory [39] and achieves identical results to the literature [66, 68, 69]. Furthermore, the first six mode shapes are represented smoothly in Figure 16. Finally, let us consider a three-layer composite plate with the aforementioned shape subjected to uniform temperature rise. The thickness is h=0.06 m, and material III is chosen. The effect of fiber orientation on the critical temperature Tcr is revealed in Table 7 for various boundary conditions. The first four buckling mode shapes of SSSS, SCSC, and CCCC, respectively, with fiber orientation θ = 45° are also plotted in Figure 17.

10

10

9

8

8 7

6 y

y

6

4

5 4 3

2

2 1

0 0

2

4

6

8

0

10

0

x

2

4

6

8

10

x

a) Coarsest mesh

b) Finer mesh

Figure 15. Four-patch plate with heart cutout, in which the element boundary is in blue while red dots denote the control points.

Table 6: Normalized frequencies of a full clamped isotropic plate with heart cutout. Mode number 1 2 3 4 5 6

MM CLPT [66] 7.548 10.764 11.113 11.328 12.862 13.300

CLPT [68] 7.621 9.81 9.948 11.135 11.216 12.482

IGA FSDT [69] 7.453 9.825 9.845 10.964 11.165 12.381

30

RPT [39] 6.9177 9.5876 9.5913 10.5567 11.2211 11.8972

present 7.4263 9.8944 9.9536 11.0327 11.3621 12.5902

Figure 16. First six mode shapes of a fully clamped plate with a heart cutout.

Table 7: Critical temperature of symmetric three-layered plate with heart cutout via different boundary conditions (BCs). BC SSSS CSCS CCCC

[0°/0°/0°] 32.4844 72.2126 102.1204

[15°/-15°/15°] 33.8294 71.6637 103.554

[30°/-30°/30°] 36.8027 68.1465 104.9144

31

[45°/-45°/45°] 37.899 60.8116 107.1027

[0°/90°/0°] 31.9618 67.5953 101.5320

SSSS

CSCS

CCCC

Mode 1

Mode 2

Mode 3

Mode 4

Figure 17. First four buckling mode shapes of a three-layered [45°/-45°/45°] plate with a heart cutout via different boundary conditions.

5

Conclusions In this paper, two types of plate theories (HSDT and quasi-3D models) in conjunction

with isogeometric approach are proposed for thermal bending and buckling analyses of laminated composite plates, in which the latter model is an expansion of the former with one additional variable in the transverse displacement. Thus, with six unknown variables, the quasi-3D model accounts not only for the shear deformation but also for the thickness stretching effect. With salient features utilizing the higher-order smoothness and continuity basis function, IGA naturally fulfills the C1 continuous requirement of these plate models without any additional variables. It helps the model to reduce more unknown variables compared to LW theory or other C0 continuous quasi-3D models but still achieves accurate results. Through various numerical results, some concluding remarks can be drawn:

32

-

HSDT obtains acceptable results only for moderate and thin plates and produces more errors as plates become thicker, especially for highly anisotropic material. In contract, the quasi-3D model always provides good results compared to 3D elastic solutions. It is concluded that transverse normal strain ε z and stress σ z are very important and cannot be eliminated for thick composite plates under a thermal environment.

-

Due to thermal field and layer arrangement, temperature rise produces the extension-bending effect, which causes thermal bending. Otherwise, apart from one special case, where a symmetric composite plate subjected to uniform temperature rise, the bifurcation buckling does exist.

-

Thermal buckling parameter increases according to increase in the aspect ratio L/W, the thickness-to-length ratio h/L, number of lamina, and stiffer constraints at plate boundaries.

-

For angle-ply laminated plates, the critical buckling temperature strongly depends on value of the orientation angle. The maximum value is observed at

θ=45° for square plate and θ=60° for rectangular one. The present work is restricted to purely thermal loading. In practical applications, plates can be enforced by combined electrical-thermal-mechanical loads. For those cases in which the inter-laminar stresses play a predominant role, including ZZ in the quasi-3D is a very promising way to accurately predict the local behavior.

Acknowledgement This work was supported by the National Research Foundation of Korea (NRF) grant funded by the Korean government (MSIP) (No.2011-0030040) References [1] Reddy JN. Mechanics of laminated composite plates-theory and analysis. New York: CRC Press 2nd Edit, 2004. [2] Wang X, Dai H. Thermal buckling for local delamination near the surface of laminated cylindrical shells and delaminated growth. Journal of thermal stresses. 2003;26:423-42. [3] Wang X, Lu G, Xiao D. Non-linear thermal buckling for local delamination near the surface of laminated cylindrical shell. International journal of mechanical sciences. 2002;44:947-65.

33

[4] Noor AK, Burton WS. Three-dimensional solutions for thermal buckling of multilayered anisotropic plates. Journal of engineering mechanics. 1992;118:683-701. [5] Noor AK. Stability of multilayered composite plates. Fibre Science and Technology. 1975;8:81-9. [6] Burton WS. Three-dimensional solutions for the free vibrations and buckling of thermally stressed multilayered angle-ply composite plates. Journal of Applied Mechanics. 1992;59:869. [7] Chen L-W, Chen L-Y. Thermal buckling of laminated cylindrical plates. Composite structures. 1987;8:189-205. [8] Wu C, Tauchert T. Thermoelastic analysis of laminated plates. I: Symmetric specially orthotropic laminates. Journal of thermal stresses. 1980;3:247-59. [9] Chen W, Lin P, Chen L. Thermal buckling behavior of thick composite laminated plates under nonuniform temperature distribution. Computers & structures. 1991;41:63745. [10] Prabhu M, Dhanaraj R. Thermal buckling of laminated composite plates. Computers & structures. 1994;53:1193-204. [11] Rolfes R, Noor A, Sparr H. Evaluation of transverse thermal stresses in composite plates based on first-order shear deformation theory. Computer Methods in Applied Mechanics and Engineering. 1998;167:355-68. [12] Kant T, Khare R. Finite element thermal stress analysis of composite laminates using a higher-order theory. Journal of thermal stresses. 1994;17:229-55. [13] Kant T, Babu C. Thermal buckling analysis of skew fibre-reinforced composite and sandwich plates using shear deformable finite element models. Composite structures. 2000;49:77-85. [14] Babu CS, Kant T. Refined higher order finite element models for thermal buckling of laminated composite and sandwich plates. Journal of thermal stresses. 2000;23:111-30. [15] Ferreira A, Castro LM, Bertoluzza S. A high order collocation method for the static and vibration analysis of composite plates using a first-order theory. Composite Structures. 2009;89:424-32. [16] Tran LV, Ferreira AJM, Nguyen-Xuan H. Isogeometric analysis of functionally graded plates using higher-order shear deformation theory. Composites Part B: Engineering. 2013;51:368-83. [17] Thai CH, Kulasegaram S, Tran LV, Nguyen-Xuan H. Generalized shear deformation theory for functionally graded isotropic and sandwich plates based on isogeometric approach. Computers & Structures. 2014;141:94-112. [18] Reddy JN. A Simple Higher-Order Theory for Laminated Composite Plates. Journal of Applied Mechanics. 1984;51:745-52. [19] Nguyen-Xuan H, Tran LV, Thai CH, Kulasegaram S, Bordas SPA. Isogeometric analysis of functionally graded plates using a refined plate theory. Composites Part B: Engineering. 2014;64:222-34.

34

[20] Fazzolari FA, Carrera E. Thermo-Mechanical Buckling Analysis of Anisotropic Multilayered Composite and Sandwich Plates by Using Refined Variable-Kinematics Theories. Journal of thermal stresses. 2013;36:321-50. [21] Thai CH, Ferreira AJM, Bordas SPA, Rabczuk T, Nguyen-Xuan H. Isogeometric analysis of laminated composite and sandwich plates using a new inverse trigonometric shear deformation theory. European Journal of Mechanics - A/Solids. 2014;43:89-108. [22] Tran LV, Nguyen-Thoi T, Thai CH, Nguyen-Xuan H. An Edge-Based Smoothed Discrete Shear Gap Method Using the C0-Type Higher-Order Shear Deformation Theory for Analysis of Laminated Composite Plates. Mechanics of Advanced Materials and Structures. 2015;22:248-68. [23] Tran LV, Lee J, Nguyen-Van H, Nguyen-Xuan H, Wahab MA. Geometrically nonlinear isogeometric analysis of laminated composite plates based on higher-order shear deformation theory. International Journal of Non-Linear Mechanics. 2015;72:4252. [24] Sayyad AS, Ghugal YM, Mhaske BA. A four-variable plate theory for thermoelastic bending analysis of laminated composite plates. Journal of thermal stresses. 2015;38:90425. [25] Carrera E. Evaluation of layerwise mixed theories for laminated plates analysis. AIAA journal. 1998;36:830-9. [26] Ferreira A. Analysis of composite plates using a layerwise theory and multiquadrics discretization. Mechanics of Advanced Materials and Structures. 2005;12:99-112. [27] Murakami H. Laminated composite plate theory with improved in-plane responses. J Appl Mech(Trans ASME). 1986;53:661-6. [28] Ali J, Bhaskar K, Varadan T. A new theory for accurate thermal/mechanical flexural analysis of symmetric laminated plates. Composite structures. 1999;45:227-32. [29] Carrera E. Developments, ideas, and evaluations based upon Reissner’s Mixed Variational Theorem in the modeling of multilayered plates and shells. Applied Mechanics Reviews. 2001;54:301-29. [30] Carrera E, Miglioretti F, Petrolo M. Accuracy of refined finite elements for laminated plate analysis. Composite structures. 2011;93:1311-27. [31] Ganapathi M, Makhecha D. Free vibration analysis of multi-layered composite laminates based on an accurate higher-order theory. Composites Part B: Engineering. 2001;32:535-43. [32] Kant T, Shiyekar S. An assessment of a higher order theory for composite laminates subjected to thermal gradient. Composite structures. 2013;96:698-707. [33] Matsunaga H. Thermal buckling of cross-ply laminated composite and sandwich plates according to a global higher-order deformation theory. Composite structures. 2005;68:439-54. [34] Ferreira AJM, Roque CMdC, Carrera E, Cinefra M, Polit O. Analysis of sandwich plates by radial basis functions collocation, according to Murakami's Zig-Zag theory. Journal of Sandwich Structures & Materials. 2012;14:505-24.

35

[35] Maturi D, Ferreira AJ, Zenkour AM, Mashat DS. Analysis of Laminated Shells by Murakami’s Zig-Zag Theory and Radial Basis Functions Collocation. Journal of Applied Mathematics. 2013;2013. [36] Natarajan S, Ferreira AJM, Bordas SPA, Carrera E, Cinefra M. Analysis of composite plates by a unified formulation-cell based smoothed finite element method and field consistent elements. Composite Structures. 2013;105:75-81. [37] Natarajan S, Ferreira AJM, Nguyen-Xuan H. Analysis of cross-ply laminated plates using isogeometric analysis and unified formulation. Curved and Layered Structures2014. [38] Zenkour AM. A simple four-unknown refined theory for bending analysis of functionally graded plates. Applied Mathematical Modelling. 2013;37:9041-51. [39] Thai CH, Zenkour A, Wahab MA, Nguyen-Xuan H. A simple four-unknown shear and normal deformations theory for functionally graded isotropic and sandwich plates based on isogeometric analysis. Composite Structures. 2016;139:77-95. [40] Nguyen N-T, Hui D, Lee J, Nguyen-Xuan H. An efficient computational approach for size-dependent analysis of functionally graded nanoplates. Computer Methods in Applied Mechanics and Engineering. 2015;297:191-218. [41] Kant T, Swaminathan K. Analytical solutions for free vibration of laminated composite and sandwich plates based on a higher-order refined theory. Composite Structures. 2001;53:73-85. [42] Kant T, Swaminathan K. Analytical solutions for the static analysis of laminated composite and sandwich plates based on a higher order refined theory. Composite Structures. 2002;56:329-44. [43] Tran LV, Thai CH, Le HT, Gan BS, Lee J, Nguyen-Xuan H. Isogeometric analysis of laminated composite plates based on a four-variable refined plate theory. Engineering Analysis with Boundary Elements. 2014;47:68-81. [44] Zenkour A. Benchmark trigonometric and 3-D elasticity solutions for an exponentially graded thick rectangular plate. Archive of Applied Mechanics. 2007;77:197-214. [45] Hughes TJR, Cottrell JA, Bazilevs Y. Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement. Computer Methods in Applied Mechanics and Engineering. 2005;194:4135-95. [46] Cottrell JA, Hughes TJ, Bazilevs Y. Isogeometric analysis: toward integration of CAD and FEA: John Wiley & Sons, 2009. [47] Cottrell J, Reali A, Bazilevs Y, Hughes T. Isogeometric analysis of structural vibrations. Computer Methods in Applied Mechanics and Engineering. 2006;195:525796. [48] Manh ND, Evgrafov A, Gersborg AR, Gravesen J. Isogeometric shape optimization of vibrating membranes. Computer Methods in Applied Mechanics and Engineering. 2011;200:1343-53.

36

[49] Bazilevs Y, Hsu M-C, Scott M. Isogeometric fluid–structure interaction analysis with emphasis on non-matching discretizations, and with application to wind turbines. Computer Methods in Applied Mechanics and Engineering. 2012;249:28-41. [50] Nguyen-Thanh N, Valizadeh N, Nguyen M, Nguyen-Xuan H, Zhuang X, Areias P, et al. An extended isogeometric thin shell analysis based on Kirchhoff–Love theory. Computer Methods in Applied Mechanics and Engineering. 2015;284:265-91. [51] Tran LV, Ly HA, Lee J, Wahab MA, Nguyen-Xuan H. Vibration analysis of cracked FGM plates using higher-order shear deformation theory and extended isogeometric approach. International journal of mechanical sciences. 2015;96:65-78. [52] Herath MT, Natarajan S, Prusty BG, John NS. Isogeometric analysis and Genetic Algorithm for shape-adaptive composite marine propellers. Computer Methods in Applied Mechanics and Engineering. 2015;284:835-60. [53] Tran LV, Thai CH, Nguyen-Xuan H. An isogeometric finite element formulation for thermal buckling analysis of functionally graded plates. Finite Elements in Analysis and Design. 2013;73:65-76. [54] Tran LV, Phung-Van P, Lee J, Wahab MA, Nguyen-Xuan H. Isogeometric analysis for nonlinear thermomechanical stability of functionally graded plates. Composite structures. 2016;140:655-67. [55] Phung-Van P, Tran LV, Ferreira AJM, Nguyen-Xuan H, Abdel-Wahab M. Nonlinear transient isogeometric analysis of smart piezoelectric functionally graded material plates based on generalized shear deformation theory under thermo-electromechanical loads. Nonlinear Dynamics. 2017;87:879-94. [56] Carrera E. Transverse normal strain effect on thermal stress analysis of homogeneous and layered plates. AIAA journal. 2005;43:2232-42. [57] Les Piegl, Wayne Tiller. The NURBS book. Germany: Springer 2nd edition, 1997. [58] Kiani Y, Eslami MR. An exact solution for thermal buckling of annular FGM plates on an elastic medium. Composites Part B: Engineering. 2013;45:101-10. [59] Thai CH, Nguyen-Xuan H, Nguyen-Thanh N, Le TH, Nguyen-Thoi T, Rabczuk T. Static, free vibration, and buckling analysis of laminated composite Reissner–Mindlin plates using NURBS-based isogeometric approach. International Journal for Numerical Methods in Engineering. 2012;91:571-603. [60] Pagano N. Exact solutions for rectangular bidirectional composites and sandwich plates. Journal of composite materials. 1970;4:20-34. [61] Cetkovic M. Thermo-mechanical bending of laminated composite and sandwich plates using layerwise displacement model. Composite structures. 2015;125:388-99. [62] Bhaskar K, Varadan T, Ali J. Thermoelastic solutions for orthotropic and anisotropic composite laminates. Composites Part B: Engineering. 1996;27:415-20. [63] Nguyen TN, Thai CH, Nguyen-Xuan H. On the general framework of high order shear deformation theories for laminated composite plate structures: A novel unified approach. International journal of mechanical sciences. 2016;110:242-55. [64] Cetkovic M. Thermal buckling of laminated composite plates using layerwise displacement model. Composite structures. 2016;142:238-53.

37

[65] Singh S, Singh J, Shukla KK. Buckling of laminated composite plates subjected to mechanical and thermal loads using meshless collocations. Journal of Mechanical Science and Technology. 2013;27:327-36. [66] Liu G, Chen X. A mesh-free method for static and free vibration analyses of thin plates of complicated shape. Journal of sound and vibration. 2001;241:839-55. [67] Kiendl J, Bazilevs Y, Hsu MC, Wüchner R, Bletzinger KU. The bending strip method for isogeometric analysis of Kirchhoff–Love shell structures comprised of multiple patches. Computer Methods in Applied Mechanics and Engineering. 2010;199:2403-16. [68] Shojaee S, Izadpanah E, Valizadeh N, Kiendl J. Free vibration analysis of thin plates by using a NURBS-based isogeometric approach. Finite Elements in Analysis and Design. 2012;61:23-34. [69] Yin S, Hale JS, Yu T, Bui TQ, Bordas SP. Isogeometric locking-free plate element: a simple first order shear deformation theory for functionally graded plates. Composite Structures. 2014;118:121-38.

38