Hybrid-Trefftz formulation for analysis of anisotropic and symmetric laminated plates

Accepted Manuscript Hybrid-Trefftz formulation for analysis of anisotropic and symmetric laminated plates Mohammad Karkon PII: DOI: Reference:

S0263-8223(15)00794-1 http://dx.doi.org/10.1016/j.compstruct.2015.08.098 COST 6803

To appear in:

Composite Structures

Please cite this article as: Karkon, M., Hybrid-Trefftz formulation for analysis of anisotropic and symmetric laminated plates, Composite Structures (2015), doi: http://dx.doi.org/10.1016/j.compstruct.2015.08.098

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.

Hybrid-Trefftz formulation for analysis of anisotropic and symmetric laminated plates Mohammad Karkon (Ph.D) Civil Engineering Department, Ferdowsi University of Mashhad [email protected]

Abstract: In this paper, two novel and efficient elements are proposed for analysis of anisotropic and symmetric laminated thick plates based on hybrid-Trefftz formulation. In this approach, two independent displacement fields are defined for internal and boundary of the elements. The internal field is selected in such a manner that the governing equation of thick plates could be satisfied. At the first step, the governing differential equations of anisotropic and symmetric laminated thick plates are derived. Then, the analytical solution of these equations is considered. By using these analytical functions, the internal displacement field of proposed elements is achieved. Also, the boundary field is related to the nodal degrees of freedom by the boundary interpolation functions. For this aim, the interpolation functions of tow-node laminated Timoshenko beam element are calculated for boundary field. By utilizing these displacement fields and hybrid-Trefftz functional, a triangular element THT-7 and quadrilateral element QHT-11 are proposed. These elements have 9 and 12 degrees of freedom, respectively. The results of the numerical testes are proven the absence of the shear locking and high accuracy and efficiency of the proposed elements for analysis of anisotropic and symmetric laminated plates. Keywords: Finite element, Hybrid-Trefftz method, thick laminated plate, quadrilateral element, triangular element.

1. Introduction Composite laminated plates due to superior characteristics are widely used in aerospace, civil, mechanical and in other fields of modern technology. Several theories have been proposed to analyze such structures. In the classical laminate plate theory (CLPT) which is based on the assumptions of Kirchhoff’s plate theory, the shear deformation is neglected [1]. Although, this method is appropriate for thin plate's analysis, the error of responses obtained from this theory rises by increasing the thickness of plates. Two types of theory, first-order shear deformation theory (FSDT) and the higher order shear deformation theory (HSDT) have been developed based on shear deformations affect. In the FSDT is assumed that transverse shear strain, to be constant through the laminate thickness [2, 3]. Therefore, this theory provides accurately appreciate responses for thin and moderately thick plates. The HSDT theory that has been proposed by Reddy [4], of course provides a slightly improved accuracy, but at the cost of a significant increase in computation time [5].

1

Taking advantage of the analytical solution of differential equation in the finite element formulation, causes the element accuracy is increased. Unlike the finite element method (FEM) based on the Ritz method, the boundary element method (BEM) provides using the analytical solution of governing differential equation (fundamental solution) for analyzing the laminated plates. For the sake of inserting the analytical solution of differential equation in the finite element formulation, the hybrid-Trefftz functional is used. This functional includes two independent displacement fields in internal domain and boundary of element. This technique, like the Ritz’s process, is suitable to solve solid mechanic problems [6]. However, dissimilar to the Ritz’s method, the used polynomial function for the internal field establishes the governing equation of the structure. In other words, homogeneous and particular solution of the governing equation is used to set up the finite element formulation. The Hybrid-Trefftz formulation was first employed to assess the distorted meshes of the thin plates by Jirousek and Leon in 1970s [7]. In this method, weak form of compatibility between elements is achieved by applying independent boundary displacement field and using variational principles. Also, The Main advantage of this formulation is that it allows for the integration on the boundary, which is not only easier than integration on the surface, but also makes it possible to create polygon elements [8]. Because of high accuracy and efficiency of hybridTrefftz method, it has been widely used to analyze the isotropic thin and thick plates [9-14]. But so far, due to the complexity, this method is not extended for analysis of laminated thick plates. Up to now, few researches have been done about analyzing the anisotropic and laminated plates based on the analytical solution. In 1990 Jirousek and Diaye developed the hybridTrefftz finite element formulation for analysis of thin orthotropic plates [15]. They are proposed a simple method for generating the set of trial functions that satisfies in the governing differential equation of orthotropic thin plate. Petrolito by using the analytical solution of isotropic plate proposed a triangular element based on the hybrid-Trefftz formulation for static, free vibration and stability analysis of thick orthotropic plates [16, 17]. Shi and Bezine presented a boundary element analysis of plate bending problems based on Kirchhoff assumptions [18]. They used the real variable fundamental solution proposed by Wu and Altiero [19]. Rajamohan and Raamachandran to avoid domain discretization, presented a boundary element method, called the charge simulation method, for analysis of anisotropic thin-plate bending problems [20]. Dong et al. developed a Trefftz boundary collocation method to analyze the bending of anisotropic thin plates. In this method, using the complicated real variable fundamental solution is avoided [21]. In 2006, the radial integration method was used by Albuquerque et al., to obtain a boundary element formulation without any domain integral for analysis of general anisotropic plate bending problems [22]. Reis et al. developed a boundary element method for stress analysis of composite laminate thin plates [23]. Wünsche et al. presented a hypersingular boundary element method (BEM) for analysis of thin anisotropic plates [24]. They adopted discontinuous quadratic elements for boundary of plate. The exact series solutions of orthotropic and symmetric laminated rectangular thin plates with any combination of clamped and simply supported edges are presented by Bhaskar and Kaushik [25]. These solutions are derived using the principle of virtual work. However, fewer studies have been done on the anisotropic and laminated thick plates. Wang and Huang presented a boundary element method for orthotropic thick plates [26]. This 2

formulation was extended to analyze laminated composite plates by Wang and Schweizerhof [27]. In 2006, Sladek et al. presented a meshless local Petrov-Galerkin (MLPG) method to solve static and dynamic problems of orthotropic plates described by the Reissner-Mindlin theory [28]. Also, Reis and coworkers introduced a boundary element method for analysis of orthotropic shear deformable plates based on Mindlin assumption [29]. Thai and Kim presented analytical solution of orthotropic Levy-type plates based on a two variable refined plate theory [30]. In this paper, by using Hybrid-Trefftz functional two efficient elements, triangular element (THT-7) and quadrilateral element (QHT-11), are proposed for analysis of anisotropic and symmetric laminated thick plates. These elements have 9 and 12 degrees of freedom, respectively. At first, the governing differential equations of such plates are derived based on first-order shear deformation theory (FSDT). Then the homogenous and particular solutions of the governing equations are obtained to establish the inner displacement field of the elements. Besides, the boundary displacement field is depended on the nodal degrees of freedom by using interpolation functions. For this aim, each edge of proposed elements is modeled as two-node laminated Timoshenko beam and interpolation functions of this beam are applied for boundary displacement field. For introducing the interpolation functions, the deflection, rotation and torsion fields of Timoshenko beam, are selected from cubic, quadratic and linear functions, respectively. By depicting the aforesaid fields in general coordinates, xy; the interpolation functions of the presented elements’ edges are obtained. Finally, several numerical tests are performed to investigate the robustness of these elements. The findings prove that the suggested elements are able to analyze the thick and thin laminated plates with high level of accuracy.

2. Hybrid-Trefftz functional Hybrid-Trefftz functional with independent inner displacement field ( u ) and boundary displacement field ( u
) has the following appearance [10, 14]: Π (u , u
) =

∫

Ω

In which, W

W u d Ω − ∫ bi u i d Ω − Ω

u

,

∫

ΓT

T i u
i d Γ t −

∫ T (u Γ

i

i

− u
i ) d Γ

(1)

bi , Ti and T i denote strain energy density function, body force, tractions,

and prescribed traction along boundary

ΓT

, respectively. The inner displacement field ( u )

which is utilized in this equation should satisfy the governing equation. Moreover, u
provides the compatibility and continuity between elements in the weak form. The inner displacement and tractions can be stated as follows: m

{u } = {u p } + ∑{φ j }c j j =1

= {u p } + [ Φ]{c }

(2)

{T } = {T p } + [Θ ]{c }

(3)

It should be added this equality is the sum of a homogenous and non-homogenous part. In this relation, {u p } and [Φ]{c} are particular and homogenous solutions of the governing equation, respectively. The traction resulting from {u p } is demonstrated by{T p } , and [ Θ]{c } denotes the tractions obtained from homogenous solution of the differential equation. 3

Furthermore, boundary displacement field is related to the displacement interpolation functions as follows:

{u
} =  N
 {d
}

(4)

In this equality, the nodal displacements and boundary interpolation fields are shown by {d
} and  N
 , respectively. To obtain stiffness matrix and nodal displacements of the element, the functional (1) should be stationary with respect to independent fields. The process of stationary will lead to the following formulas: δ Π (u , u
) δ u = − ∫ δT i (u i − u
i ) d Γ = 0 (5) Γ

δ Π (u , u
) δ u
=

∫T Γ

i

δ u
i d Γ − ∫ T i δ u
i d Γ = 0

(6)

ΓT

Inserting equations (2), (3) and (4) into equations (5) and (6), leads to the succeeding relations: T

δ {c } {h } + [ H

]{c } − [G ]{d } = 0

(7)

T T δ {d } {g } + [G ] {c } = 0



(8)



By deploying relationship (7), {c } can be computed as follows: −1

−1

{c } = [ H ] [G ]{d } − [H ] {h }

(9)

The parameters employed in the equations (7) to (9), are defined in below form:

{h } = ∫Γ [ Θ ]

T

{u }d Γ

(10)

p

T

{g } = ∫Γ  N


{T } d Γ − ∫ p

ΓT

T

 N
 {T }d Γ

(11)

[ H ] = ∫Γ [ Θ ] [φ ]d Γ

(12)

[G ] = ∫Γ [Θ ]

(13)

T

T

 N
 d Γ

By substituting equation (9) into equation (8), results in the coming relation: −1

−1

{ g } − [G ] [ H ] {h } + [G ] [ H ] [G ]{d
} = 0 T

T

(14)

Finally, the stiffness matrix and the nodal forces vector can be computed as the next form: T

−1

[ K ] = [G ] [ H ] [G ] T −1 {f } = [G ] [ H ] {h } − {g }

(15) (16)

Note that, the equations (15) and (16) can be applied in all solid mechanic problems. In this paper, these equations are employed for analyzing laminated thick plates.

3. Governing differential equations Governing equations of anisotropic and symmetric laminated plates are derived based on the first-order shear deformation theory (FSDT). For this aim, the equations of displacement, stress, strain, and the relationship between them, and also tractions in thick plate could be expressed as: T

{u} = {w θx θy }

(17) 4

{ε } = {κ {σ } = {M {T} = {Q

T

κ y κxy γ x γ y } = [ L]{u}

*

x

*

n

−Mnx −Mny

T

} = [ D]{ε } } = [ A]{σ }

M y Mxy Qx Qy

x

(18)

T

*

(19)

*

(20)

Figure (1).Tractions along with an arbitrary direction

Tractions of {T } are displayed in Figure (1). In addition, the strain operator matrix [ L ] and boundary tractions transformation matrix [ A ] are obtained in the following manner:

0  ∂ ∂x  0  0 ∂ ∂y  0  [ L ] =  0 ∂ ∂y ∂ ∂x    0  −1 ∂ ∂x ∂ ∂y 0 −1   0 0 0 nx  [ A ] =  −n x 0 −n y 0  0 0 −n y −n x 

(21)

ny   0  0 

(22)

n x = cos α  n y = sin α

(23)

Also, the elasticity matrix for anisotropic and symmetric laminated plates could be expressed as below:

[ D ] [ D ] =  0M 

0  [DS ]

(24)

In this relation, matrices [ DM ] and [ DS ] are the bending stiffness and transverse shear stiffness, respectively. These quantities are calculated in the following form: 5

 D11 D12 [ D M ] = − D12 D 22 D16 D 26

D16  D 26  D 66 

(25)

0 C 55 C 45   k s 1C 55 D = =  [ S ] C C  0  k s 1k s 2C 45  45 44 

(26)

ks1and k s 2 are shear correction factors. The entries of these matrices are calculated as:

Where

Dij =

k s 1k s 2C 450   0 k s 2C 44 

1 n (k ) 3 Qij z k − z k3 −1 ∑ 3 k =1

(

n

C ij0 = ∑Qij(

k)

)

( i , j = 1, 2,6)

( z k − z k −1 )

(27)

(i , j = 4,5)

(28)

k =1

Figure (2).Layer coordinates of a symmetric laminate

As shown in figure (2), the z k −1 and z k in these relations, are the coordinates of the lower and upper surfaces of the kth layer in the Z direction. The origin of the coordinate is the symmetric plane of the laminate. Also, Qij( ) are the elastic constants of the kth layer in k

laminate coordinates that are given in Appendix (A). In the present theory, the in-plane displacement components u′ and v′ , in a laminate element can be expressed as [31]: u′ ( x, y, z ) = zθ x (29) v′ ( x, y, z ) = zθ y

(30)

By substituting these relations into the strain-displacement equations of the classical theory of elasticity the following in-plane and out-of-plane strain components are obtained [31]: T

{ε } = {ε x ε y γ xy γ xz γ yz } ε x = zκ x

,

ε y = zκ y

,

(31)

γ xy = zκ xy

(32)

∂w ∂w − θx γ yz = −θy , (33) ∂x ∂y The stress-strain relationship for the kth layer (lamina) of the composite laminate has the following form:

γ xz =

6

Q11( k ) Q12( k ) Q16( k ) 0 0 ε  σ x    x  (k ) (k) (k ) σ  0 0  εy Q12 Q22 Q26  y   (k )    k) k) ( ( τ = (34) Q Q Q 0 0  xy   16  γ xy  26 66 τ     ( k )  γ xz 0 0 Q55( k ) Q45  xz   0   γ  (k ) (k )   τ yz ( k )  0 0 0 Q45 Q44  yz ( k ) The governing differential equations of laminated thick plates are obtained if the differential

equilibrium conditions are written in terms of {u } as: T T [ L ] {σ } = [ L ] [D ][ L ]{u } = − {P } The distributed load vector {P } in the present formula is defined as follows:

{P } = {q

0

T

0}

(35)

(36)

By inserting the relations (17), (21) and (24) into the equation (35), the governing equations of the symmetric laminated thick plates can be obtained as:  ∂ 2w ∂θ   ∂ 2w ∂θ  (C 55 + C 45 )  2 − x  + (C 45 + C 44 )  2 − y  = q (37)  ∂x

D11

∂x 

∂y 

∂ 2θ y ∂ 2θ y ∂ 2θ y ∂ 2θ x ∂ 2θ x ∂ 2θ x + 2 D + D + D + D + D + D ( 12 66 ) 16 66 16 26 ∂x 2 ∂x ∂y ∂y 2 ∂x 2 ∂x ∂y ∂y 2  ∂w − C 55  −θx  ∂x

D16

 ∂y

 ∂w  −θy  − C 45    ∂y

 =0 

(38)

∂ 2θ y ∂ 2θ y ∂ 2θ y ∂ 2θ x ∂ 2θ x ∂ 2θ x + D + D + D + D + 2 D + D ( ) 12 66 26 66 26 22 ∂x 2 ∂x ∂y ∂y 2 ∂x 2 ∂x ∂y ∂y 2

(39)  ∂w   ∂w  − C 45  − θ x  − C 44  −θy  = 0  ∂x   ∂y  The coupling of the governing differential equations (37) to (39) makes it difficult to generate a T-complete set of homogeneous solutions for rotations θx =

θx

w , θx and θy . To bypass this difficulty, the

and θy are supposed in the following manner:

∂w + f 1 (x , y ∂x

(40)

)

∂w + f 2 (x , y ) (41) ∂y By inserting the present relationships in the equations (38) and (39), the unknown functions f 1 and f 2 are obtained as below:

θy =

f 1 (x , y ) = f 2 (x , y ) =

1 C 44C 55 − C 452 1 C 44C 55 − C 452

( MC 44 − NC 45 )

(42)

( NC 55 − MC 45 )

(43)

Where

7

M = D11

∂ 3w ∂ 3w ∂ 3w ∂ 3w + 3 D + D + 2 D + D ( ) 16 12 66 26 ∂x 3 ∂x 2 ∂y ∂x ∂y 2 ∂y 3

(44)

N = D16

∂ 3w ∂ 3w ∂ 3w ∂ 3w + D + 2 D + 3 D + D ( ) 12 66 26 22 ∂x 3 ∂x 2 ∂y ∂x ∂y 2 ∂y 3

(45)

Therefore, the governing differential equations of anisotropic and symmetric laminated thick plates are obtained in the following appearance: ∂ 4w ∂ 4w ∂ 3w ∂ 4w ∂ 4w D11 4 D 2 D 2 D 4 D D + + + + + = q (x , y ) (46) ( ) 16 12 66 26 22 ∂x 4 ∂x 3 ∂y ∂x 2 ∂y 2 ∂x ∂y 3 ∂y 4 ∂w 1 θx = + ( MC 44 − NC 45 ) (47) ∂x C 44C 55 − C 452

θy =

∂w 1 + ( NC 55 − MC 45 ) ∂y C 44C 55 − C 452

(48)

4. Analytical solution of governing differential equation To solve the governing differential equation (46), this equation can be rewritten as[32]: L 4 L 3 L 2 L1w = 0

(49)

Where the linear differential operator L j , is defined as below: Lj =

∂ ∂ − µj ∂y ∂x

(j

= 1, 2, 3, 4 )

(50)

The parameter µ j , are the roots of the characteristic equation given by: D 22 µ 4 + 4D 26 µ 3 + 2 ( D12 + 2D 66 ) µ 2 + 4D16 µ + D11 = 0

(51)

It has been proved that the equation (51) has no real roots for any elastic homogeneous material [32]. Therefore, for any elastic bodies µ j are complex or pure imaginary numbers, and are complex parameters for anisotropic materials. Since the coefficients of µ j , there are two pairs of complex conjugates for equation (51) as follows: µ1 = α1 + i β1 , µ2 = α 2 + i β 2

µ3 = µ1

, µ4 = µ2

(52)

Where, α1 , β1 , α 2 and β 2 are real numbers and β1 both β2 are positive. Also, the complex variables z j are defined as below shape: z j = x − µj y

(j

= 1, 2, 3, 4 )

(53)

If the roots of characteristic equation (51) are not identical ( µ1 ≠ µ 2 ) , the general and complete solution of differential equation (46) can be obtained: w ( x , y ) = 2 Re U 1 ( z 1 ) + U 2 ( z 2 ) 

(54)

Also, if the roots are pairwise equal ( µ1 = µ 2 ) , the general solution has the below shape: w ( x , y ) = 2 Re  z 1U 1 ( z 1 ) + U 2 ( z 1 ) 

(55)

8

In relations (54) and (55), U 1 and U 2 are arbitrary analytic functions of the complex variable

z 1 and z 2 , respectively. The function U j ( z j , z j ) from nth order is assumed as[33]: U (j n ) ( z j , z

j

) =c z 0

n j

+ c1z nj −1z j +  + c n −1z j z jn −1 + c n z jn

(56)

Where c i ( i = 1, 2, n ) , are complex constants. Since the U j ( z j , z

j

)

should be analytical

function, the necessary and sufficient condition for the function U j ( z j , z j ) being analytical is that the Cauchy-Riemann conditions are satisfied as below[33]: ∂U j( n ) ( z j , z j ) ∂z j

(57)

=0

These conditions, give the following results: c 2 = c3 =  = c n = 0

(58)

The equation (56) can then be simplified as below: U (j n ) ( z

j

)=c z 0

n j

(59)

Therefore, function U j ( z

j

) is obtained as following shape:

∞

U j ( z j ) = ∑ c 0 z nj

(60)

n =2

By substituting equation (60) in relation (54), the analytical solution of the governing differential equation (46), for case ( µ1 ≠ µ 2 ) can be expressed in the following form: ∞

w ( x , y ) = ∑ 2 Re c 0 z 1n + d 0 z 2n  n =2

(61)

∞

n n = ∑ 2 Re ( A + i B )( x + µ1 y ) + (C + i D )( x + µ2 y )    n =2

Where, A, B, C and D are the real constants. Based on this solution, the firs eleven polynomial functions for case ( µ1 ≠ µ 2 ) , that satisfy the governing equation (46) are given in Appendix (B). Also, based on relations (55) and (60), the homogeneous solution of equation (46) for the case that roots are pairwise equal ( µ1 = µ 2 ) is given by: ∞

w ( x , y ) = ∑ 2 Re c 0′z 1z 1n −1 + d 0′z 1n  n =2 ∞

= ∑ 2 Re ( A ′ + i B ′ )( x + µ1 y )( x + µ1 y )  n =2

(62) n −1

+ (C ′ + i D ′ )( x + µ1 y )   n

In this relation, parameters A', B', C' and D' are the real numbers. The corresponding polynomial functions from 2, 3 and 4 order are reported in Appendix (B). Also, the particular solution of equation (46) under uniform distributed transverse load q has the following form[21]:

(

w p = Q D11x 4 + 4D16 x 3 y + 2 ( D12 + 2D 66 ) x 2 y 2 + 4D 26 xy 3 + D 22 y 4 Q=

(63)

q

(

2 11

2 16

)

2

2 26

8 3D + 12D + 2 ( D12 + 2D 66 ) + 12D + 3D 9

2 22

)

In the particular and homogeneous solution, there is no advantage over x and y, and therefore, the element will be rotational invariant. The minimum of the terms selected from the homogenous solutions (61) and (62), is dependent on the element’s number of degrees of freedom. Number of necessary terms (and not adequate) to prevent numerical instability and rank deficiency of the stiffness matrix is acquired by employing the coming relation [6]: m ≥ NDOF − r (64) In this equality m, NDOF and r denote the minimum number of homogenous solutions of equation (46), number of the element’s degrees of freedom and rigid body motions, respectively. It should be mentioned that using modes of the rigid body motion, product spurious strain energy in the element. Consequently, these modes should be omitted. The number of modes of rigid body motions in the plate is three, namely a transitive mode and two rotational modes in the direction of x-axis and y-axis. As a result, equation (64) can be rewritten for the plates as below: m ≥ NDOF − 3 (65) In general, the minimum number of homogenous solutions, which satisfies the condition of equation (65), provides the high level of accuracy for the element [11-13]. Therefore, the proposed elements THT-7 with 7 functions and QHT-11 with 11 Trefftz functions are the most accurate elements. After calculating the nodal displacements, for computing the moments and shear forces in the individual element, the nodal displacements vector of the element is obtained. Then by deploying the equation (9), the vector of coefficients of the internal displacement field {c} (the analytical solution that presented in Appendix (B)) in the element, is calculated. Therefore, by using the equations (18) and (19), and homogeneous and particular solutions, the moments and shear forces in the individual element can be obtained. The in-plane and out-of-plane stresses in the kth layer (lamina) of the element are also obtained by utilizing the equation (34).

5. Boundary displacement field The independent boundary displacement is related to the nodal values by interpolation functions. For this aim, the interpolation functions of laminated Timoshenko beam element as shown in figure (3), is used for boundary of proposed elements. In order to compute shape function of aforementioned beam element, a cubic displacement polynomial and a quadratic rotational field are selected. Moreover, it is assumed that shear strain has the constant value of γ 0 . Based on these assumptions, the following equations can be written: w
=

w wi (1 − s ) + j (1 + s ) + β0 l 1 − s 2 + β1ls 1 − s 2 2 2

(

)

(

θ θi (1 − s ) + j (1 + s ) + α 0 (1 − s 2 ) 2 2 γ = γ0

)

, s = 2x l − 1

(66) (67)

θ
=

(68)

10

t

s

x

j

i

θi

l wi

θj

wj Figure (3).Tow-node Timoshenko beam element

By minimizing internal strain energy of beam element with respect to γ 0 , the unknown parameters α 0 , β0 and β1 , are found as below [12]: 2 (w j −w i ) − (θi + θ j ) l 3 1  α0 = − γ 0 − Γ  2 2 

(69)

Γ=

(70)

1 1 θ i − θ j , β1 = α 0 8 6 Furthermore, the shear strain γ 0 is obtained in the following way:

(

β0 =

)

(72)

γ0 = δ Γ δ=

(71)

Qb ij 6λ , λij = Q s ij l ij2 1 + 12λ

(73)

Where, Qb ij and Q s ij are the bending stiffness constant and shear stiffness constant of the beam element, respectively. Because, the each edge of the proposed elements is modeled as laminated Timoshenko beam element, according to figure (4), for side i-j of the element, they can be determined from the following equations: 1 n Q b ij = D11ij = ∑ Q11( k ) ( z k3 − z k3 −1 ) (74) ij 3 k =1

(

)

(

Q s ij = C ij = k s 1 cos 2 α ij + k s 2 sin 2 α ij

n

) ∑ (Q ( ) ) ( z k 55

ij

k =1

k

− z k −1 )

(75)

Where

(Q ( ) ) k 11

(

k

ij

)

= Q11( ) cos 4 (θ k − αij ) + 2 Q12( ) + 2Q 66( ) cos 2 (θ k − α ij ) sin 2 (θk − α ij ) k

k

+ Q 22( ) sin 4 (θ k − α ij ) k

(Q ( ) ) k 55

ij

(k ) = Q 55( k ) cos 2 (θ k − α ij ) + Q 44 sin 2 (θ k − α ij )

11

(76) (77)

Figure (4).The degrees of freedom placed on the element’s sides

By substitution β1 , β 0 , α 0 and γ 0 into (66), the succeeding interpolation functions for laminated Timoshenko beam can be found [12]: w i    w
  N 1 N 2 N 3 N 4  θ ni   =   θ
n   N 5 N 6 N 7 N 8  w j  θ nj    1  2 + s 3 (1 − 2δ ) + s ( −3 + 2δ )   4 l N 2 =  0.5 1 − s 2 + s 3 − s ( 0.5 − δ )  4 1 N 3 =  2 − s 3 (1 − 2δ ) − s ( −3 + 2δ )  4 l N 4 =  −0.5 1 − s 2 + s 3 − s ( 0.5 − δ )   4 1 N 5 =  6 1 − s 2 ( −1 + 2δ )   4l  1 N 6 =  −1 + s ( −2 + 3s ) + 6 1 − s 2 δ   4 1 N 7 =  −6 1 − s 2 ( −1 + 2δ )   4l  1 N 8 =  −1 + s ( 2 + 3s ) + 6 1 − s 2 δ   4

(78)

N1 =

(

) (

(

(

)

) (

)

(79)

)

(

(

)

)

(

)

Moreover, the linear rotation field ( θ s ) has the succeeding form:

θ
s = N i θ si + N j θ sj N i = (1 − s ) 2 , N

(80) j

= (1 + s ) 2

(81)

Based on figure (4), the subsequent results can be obtained:

12

cos α = y ji l ij  sin α = − x ji l ij

(82)

By using the subsequent equation, x and y parameters can be expressed in terms of central coordinate s as below:  x ji = x j − x i   y ji = y j − y i (83)   x = x i + ( s + 1) x ji 2 y = y + s +1 y 2 ( ) ji i  According to figure (4), by depicting the nodal rotations in x-y space the following result is achieved:

θsm = θxm cosαij +θym sinαij  θnm = θym cosαij −θxm sinαij

m =i , j

(84)

Furthermore, the boundary fields θ
x and θ
y can be written as follows:

θ
x = θ
s cosαij −θ
n sinαij 

θy = θs sinαij +θ
n cosαij

(85)

By employing the displacement fields (78) and (80) and utilizing equations (84) and (85), the boundary displacement field for the proposed elements can be obtained as the coming shape:

{u
} = N
 {d
} {u
} = {w


{d
} = {w


θ
x i

θ
xi

(86)

θ
y

T

}

θ
yi w
j θ
xj

(87) θ
yj

T

}

(88)

N 11 N 12 N 13 N 14 N 15 N 16  N
 = N 21 N 22 N 23 N 24 N 25 N 26  (89)     N 31 N 32 N 33 N 34 N 35 N 36  All interpolation functions N ij in the matrix (89) are given explicitly in the Appendix (C).

6. Numerical tests In order to assess the accuracy of the proposed element, some numerical problems have been analyzed and its results compared with results available in the literature. At first, the patch test is performed to ensure the convergence of the proposed elements. The laminated square plates with simply supported and clamped edges under uniform distributed load with different stacking sequence are analyzed in subsections 6.2 and 6.3. To better reveal capability of proposed elements, the analysis of skew plate with different skew angles is performed in subsection 6.4. At end stage, square and rectangular laminated plates under sinusoidal are analyzed. For convenience, the following nondimensional deflection, moment and stresses will be given in terms of the following nondimensional forms:

13

w = 100 w

E2t 3 q0 a 4

,

M = 100 M

2 a a h h

σ x(3) = σ x( 3)  , ,  2  2 2 2  q0 a (3 )

h  h2 0, 0,   2  q0 a 2 

( 3) 

τ xy = τ xy

, ( 2)

1 q0a 2

(90) 2 a a h h

σ y( 2) = σ y( 2)  , ,  2  2 2 6  q0 a ( 2) 

, τ xz = τ xz

a  h a  h , τ yz( 2) = τ (yz2)  , 0, 0   0, , 0  q a 2 2   0   q0 a

(91)

Where the superscript ( k ) denotes the layer number.

6.1. Patch test The behavior of the proposed elements is assessed through the patch test. This test will check the convergence of the elements. The mesh geometry of the test is shown in Figure (5). The thickness of the plate is equal to h and the material properties are: E x = E y , G xy = G xz =G yz = E x 2.6 , ν xy = 0.3 . Also, the modulus of elasticity E x and shear correction

factor k s , are taken as 106 and 5/6, respectively. For this purpose, two patch tests with constant moment and constant shear stress are performed. For the constant moment patch test, the test is performed by enforcing the boundary displacements as below: w = −10 −3 ( x 2 + xy + y 2 )

, θ x = −10−3 ( 2x + y )

, θ y = −10−3 ( x + 2 y )

(92)

According to these quantities, the corresponding exact stresses everywhere on the plate should be equal to the following amounts: M x = M y = 238.095 h 3 , M xy = 64.103 h 3 , Q x = Q y = 0 (93) Also, the constant shear stress patch test is performed by enforcing the following boundary conditions: w = −10 −3 ( x + y ) , θ x = 10−3 , θ y = 10−3 (94) The corresponding stress solution must be everywhere as: M x = M y = 2M xy = 0 , Q x = Q y = 641.026h

(95)

In these tests, the proposed elements THT-7 and QHT-11 give the exact solutions. Based on the numerical experiences, it can be concluded that the presented elements simply pass the constant moment and shear patch test.

Figure (5).Geometry for the patch test

14

6.2. Cross-ply rectangular laminated plate under uniform load In this test, a cross-ply square laminated plate that subjected to uniformly distributed load is analyzed by using the proposed elements. At the first example, a three-layered symmetrically laminated square plate with stacking sequence  0 ,90 , 0  and simply supported at all the edges is analyzed for three different length to thickness ratios (a/h ). The shear correction factors k s 1 and k s 2 of the plate are assumed 5/6. Also, the material properties of each layer are taken as: E x = 25E y , G xy = 0.5E y , G xz = 0.5E y , G yz = 0.2E y , ν xy = 0.25

(96) Due to symmetry, one quarter of plate is analyzed. Figure (6) shows this structure with typical mesh 2×2. The nondimensional central deflection and stresses for the proposed elements along with the published results based on the FSDT theory, are presented in table (1). It is clearly seen that for all the aspect ratios studied, the proposed elements have rapid rate of convergence and variations of the aspect ratios do not have a significant influence on the accuracy of the results. This table also show that the present results, have an excellent agreement with the existing compared results.

b/2 C.L C.L

b/2

a/2

a/2

Figure (6).Square plate (a=b) with typical mesh 2×2 for quarter section

15

Table (1).Nondimensional deflection and stresses of a simply supported 3-layer symmetric cross-ply 0 , 90 , 0  laminated square plate under uniform transverse load.  

a/h

Method

10

THT-7

QHT-11

w

σx

σy

τ xy

τ xz

τ yz

3×3 6×6 12×12

0.9875 1.0032 1.0076

0.7764 0.7740 0.7735

0.3133 0.3037 0.3013

0.0412 0.0452 0.0460

0.5866 0.6068 0.6132

0.1763 0.2287 0.2604

3×3 6×6 12×12

1.0117 1.0097 1.0093

0.7923 0.7781 0.7745

0.2673 0.2912 0.2981

0.0456 0.0465 0.0465

0.5841 0.5957 0.6034

0.2295 0.2566 0.2566

0.9712 1.0190 1.0219 1.0203 1.0219 1.0211 1.0291

0.7719 0.7715 0.7719 0.7851 0.7708

0.3072 0.3079 0.3074 0.3844 0.3042

0.0514 0.0506 0.0504 0.0480 0.0456

0.7548 0.6455 0.7703 0.5367 0.6761

0.3107 0.2775 0.3311 0.1850 0.2802

Reis et al. [29] NIP-Q4 (12×12) [34] Reddy [35] NREM-NC [36] Shahbazi et al. [37] Goswami [38] MSM 7NDI [39] 20

100

THT-7

3×3 6×6 12×12

0.7404 0.7492 0.7526

0.8040 0.7996 0.7988

0.2325 0.2232 0.2217

0.0382 0.0420 0.0432

0.6003 0.6214 0.6296

0.2027 0.2298 0.2495

QHT-11

3×3 6×6 12×12

0.7544 0.7540 0.7539

0.8154 0.8028 0.7997

0.1918 0.2122 0.2188

0.0430 0.0434 0.0436

0.6031 0.6100 0.6160

0.2047 0.2382 0.2570

Reddy [35] NREM-NC [36] Shahbazi et al. [37] MSM 7NDI [39]

0.7572 0.7555 0.6697 0.7578

0.7983 0.7985 0.8072 0.7985

0.2227 0.2231 0.1927 0.2208

0.0453 0.0443 0.0428 0.0389

0.7697 0.6621 0.7903 0.6887

0.2902 0.2626 0.3044 0.2648

THT-7

3×3 6×6 12×12

0.6707 0.6690 0.6689

0.8021 0.8087 0.8077

0.2115 0.1972 0.1934

0.0406 0.0413 0.0420

0.5758 0.6115 0.6323

0.1805 0.2608 0.3090

QHT-11

3×3 6×6 12×12

0.6699 0.6696 0.6696

0.8236 0.8108 0.8081

0.1812 0.1878 0.1909

0.0435 0.0428 0.0426

0.6411 0.6294 0.6294

0.1660 0.2196 0.2496

0.6723 0.6700 0.6697 0.6655 0.6697 0.6701 0.6697

0.8072 0.8084 0.8072 0.8190 0.8072

0.1925 0.1922 0.1927 0.2923 0.1910

0.0426 0.0407 0.0428 0.0397 0.0360

0.7744 0.6715 0.7903 0.5648 0.6738

0.2842 0.2447 0.3044 0.1480 0.2397

Reis et al. [11] NIP-Q4 (12×12) [34] Reddy [35] NREM-NC [36] Shahbazi et al. [37] Goswami [38] MSM 7NDI [39]

16

In order to investigating the effect of boundary conditions and length to width ratio (a/b ) of the plate on the elements’ accuracy, the analysis of three-layer cross-ply  0 ,90 , 0  rectangular laminated plate is carried out for different length to thickness ratios (a/h ). The shear correction factors are assumed to be 5/6 and the following material properties are considered for each layer: Ex = 20.83 × 106 psi

,

E y = 10.94 ×106 psi

,

Gxy = 6.10 × 106 psi

Gxz = 3.71× 106 psi

,

G yz = 6.19 ×106 psi

,

vxy = 0.44

(97)

( )

The nondimensional central deflections and axial center stresses σ x(3 )

for the proposed

elements along with the results of Liew et al. [34], are presented in tables (2) and (3), respectively. In these tables, the simply supported, clamped and free boundary conditions are denoted by S, C, and F, respectively. Therefore, a 4-word notation such as CFCS is employed to show the boundary conditions on the four edges of the plate. The 1-4th word indicates the boundary conditions on edges x=0, y=0, x=a , and y=b respectively. These tables reveal that the proposed elements have high accuracy and the rate of convergence is quite independent of the boundary conditions. It has been also observed that neither the length to width ratio nor the thickness to span ratio has significant influences on the convergent rate. Table (2).Nondimensional central deflection of a 3-layer symmetric cross-ply 0 , 90 , 0  laminated plate under uniform transverse load.

a/b a/h Method

CCCC

CCCS

SCSC

SSSC

CFCF

CFCS

CFSF

SCSF

1

THT-7 QHT-11

1.8325 1.8449

2.1285 2.1395

2.5437 2.5561

3.1268 3.1369

2.9581 2.9641

2.7167 2.7242

4.6949 4.7046

5.0966 5.1122

Liew et al. [40]

1.8274

2.1109

-

-

2.9050

2.6729

4.6394

-

THT-7 QHT-11

1.3930 1.4021

1.6171 1.6255

2.1051 2.1153

2.6865 2.6943

2.1977 2.2039

2.0330 2.0399

3.8090 3.8187

4.4963 4.5103

Liew et al. [40]

1.3918

1.6096

-

-

2.1748

2.0131

3.7875

-

THT-7 QHT-11

0.1650 0.1656

0.2128 0.2137

0.1667 0.1671

0.2146 0.2155

3.0520 3.0651

1.5020 1.5118

4.9399 4.9594

0.8087 0.8130

Liew et al. [40]

0.1693

0.2173

0.1712

0.2194

3.0070

1.4900

4.8941

0.8209

THT-7 QHT-11

0.0981 0.0986

0.1383 0.1391

0.0990 0.0995

0.1402 0.1409

2.2801 2.2932

1.1111 1.1196

4.0287 4.0484

0.6117 0.6159

Liew et al. [40]

0.1005

0.1409

0.1014

0.1428

2.2656

1.1091

4.0197

0.6200

THT-7 QHT-11

0.0518 0.0519

0.0613 0.0614

0.0519 0.0519

0.0613 0.0614

3.0940 3.1127

0.9622 0.9735

5.0097 5.0348

0.2172 0.2179

Liew et al. [40]

0.0534

0.0629

0.0534

0.0630

3.0542

0.9617

4.9720

0.2221

THT-7 QHT-11

0.0276 0.0277

0.0351 0.0352

0.0276 0.0277

0.0351 0.0352

2.3226 2.3409

0.7009 0.7102

4.1016 4.1269

0.1437 0.1444

Liew et al. [40]

0.0284

0.0360

0.0284

0.0360

2.3136

0.7049

4.1014

0.1465

5

50 7

3

5

50 7

5

5

50 7

17

Table (3).Nondimensional stress σ x of a 3-layer symmetric cross-ply 0 ,90 , 0  laminated plate under uniform transverse load.

a/b a/h Method

CCCC

CCCS

SCSC

SSSC

CFCF

CFCS

CFSF

SCSF

1

THT-7 QHT-11

0.1648 0.1645

0.1912 0.1906

0.2327 0.2321

0.2850 0.2843

0.2495 0.2499

0.2375 0.2375

0.3901 0.3954

0.4317 0.4309

Liew et al. [40]

0.1657

0.1912

-

-

0.2490

0.2346

0.3938

-

THT-7 QHT-11

0.1698 0.1693

0.1953 0.1944

0.2246 0.2240

0.2824 0.2819

0.2492 0.2495

0.2384 0.2383

0.3786 0.3838

0.4220 0.4209

Liew et al. [40]

0.1699

0.1949

-

-

0.2488

0.2356

0.3826

-

THT-7 QHT-11

0.0139 0.0142

0.0237 0.0235

0.0146 0.0149

0.0233 0.0242

0.2538 0.2570

0.1172 0.1189

0.3929 0.4049

0.0386 0.0371

Liew et al. [40]

0.0146

0.0239

0.0153

0.0246

0.2547

0.1163

0.4020

0.0388

THT-7 QHT-11

0.0132 0.0136

0.0224 0.0222

0.0130 0.0134

0.0212 0.0222

0.2508 0.2536

0.1195 0.1214

0.3783 0.3896

0.0294 0.0278

Liew et al. [40]

0.0138

0.0225

0.0136

0.0225

0.2516

0.1186

0.3887

0.0295

THT-7 QHT-11

0.0039 0.0044

0.0086 0.0090

0.0040 0.0045

0.0083 0.0090

0.2541 0.2574

0.0516 0.0524

0.3945 0.4068

0.0025 0.0033

Liew et al. [40]

0.0044

0.0091

0.0045

0.0091

0.2549

0.0512

0.4024

0.0030

THT-7 QHT-11

0.0038 0.0043

0.0077 0.0080

0.0038 0.0043

0.0073 0.0081

0.2528 0.2561

0.0587 0.0601

0.3816 0.3942

0.0049 0.0056

Liew et al. [40]

0.0043

0.0081

0.0044

0.0081

0.2539

0.0583

0.3915

0.0053

5

50 7

3

5

50 7

5

5

50 7

In the next numeric test, a four-layered symmetrically laminated plate with stacking sequence  0 ,90 ,90 ,0  and simply supported at all the edges is analyzed for different thickness ratios  

(a/h). The shear correction factors k s 1 and k s 2 of the plate are assumed 5/6. Also, the material properties of relation (96) is considered for each layer. Due to symmetry, one quarter of plate is analyzed. By performing the analysis of this structure, the nondimensional deflection and bending moment obtained at the plate center are presented in tables (4) and (5), respectively. The results of proposed elements, is compared with those of Reddy [41], Kabir [42] and Sheikh et al. [43]. The results obtained by Reddy are based on an analytical solution, while Kabir has analyzed the plate using a three-node triangular element, and Sheikh et al. used a triangular element with 55 degrees of freedom for analysis of the plate. These tables reveal that the proposed elements have a good accuracy and rate of convergence with mesh size. In the following, a nine-layered cross-ply ( 0 ,90 ) , 0 , ( 90 , 0 )  square laminated plate is 

2

2

considered. Because of symmetry, a quarter of plate is taken in the analysis. The material properties used in this analysis are:

18

E x = 40E y , G xy = 0.6E y , G xz = 0.6E y , G yz = 0.5E y , ν xy = 0.25

(98)

The shear correction factors k s 1 and k s 2 are taken to be unity. This structure is analysed for both simply and clamped support in all edges and nondimensional results are obtained. Tables (6) and (7) demonstrate the central deflection and moment of the simply supported square plate, respectively. The obtained results of the analysis of the square plate with clamped supports are shown in tables (8) and (9). The findings are compared with published results of other researchers. It should be noted that Wang and Schweizerhof, analysed this structure by using boundary element method [27]. These tables demonstrate the high efficiency and accuracy of the proposed elements. Also, it is obvious that elements free from shear locking effect. Table (4).Nondimensional central deflection ( 10 w ) of a simply supported 4-layered cross-ply 0 , 90 ,90 , 0  square plate under uniform transverse load  

a/h 10 48 1000 10000

2×2

THT-7 4×4

8×8

2×2

QHT-11 4×4

8×8

Sheikh et al. [43]

Reddy [41]

Kabir [42]

9.626 6.860

9.958 6.905

10.060 6.928

10.199 6.944

10.119 6.947

10.102 6.948

10.250 6.987

9.69 6.92

9.66 6.91

6.745

6.789

6.795

6.788

6.796

6.797

6.798

6.79

6.76

6.745

6.789

6.795

6.788

6.796

7.797

6.798

6.79

6.76

Table (5).Nondimensional central moment ( M x ) of a simply supported 4-layered cross-ply 0 , 90 ,90 , 0  square plate under uniform transverse load  

a/h 10 48 1000 10000

2×2

THT-7 4×4

8×8

2×2

QHT-11 4×4

8×8

Sheikh et al. [43]

Reddy [41]

Kabir [42]

113.88

112.72

112.35

118.36

113.60

112.56

111.35

112.68

111.84

118.24

121.03

120.68

126.66

121.91

120.79

120.47

120.51

119.57

122.54

121.11

120.98

126.53

122.38

121.30

120.95

120.93

120.11

122.56

121.12

120.99

126.52

122.38

121.30

121.04

120.93

120.12

Table (6).Nondimensional central deflection ( 10w ) of a simply supported 9-layered cross-ply  0 ,90 

(

a/h 10 100 1000

)

2

(

, 0 , 90 , 0

2×2

THT-7 4×4

8×8

5.674 4.416

5.825 4.468

5.868 4.477

4.402

4.455

4.465

)  square plate under uniform transverse load 2

2×2

QHT-11 4×4

5.954 4.471

5.901 4.481

4.455

4.467

19

8×8

Analytical [44]

QUAD4 [45]

TRIPLT [46]

5.888 4.482

5.848 4.486

5.84 4.47

5.848 4.485

4.468

4.472

4.46

4.451

Table (7).Nondimensional central moment ( M x ) of a simply supported 9-layered cross-ply  0 ,90 

(

a/h 10 100 1000

)

2

(

, 0 , 90 , 0

)  square plate under uniform transverse load 2

2×2

QHT-11 4×4

8×8

Analytical [44]

QUAD4 [45]

TRIPLT [46]

9.071

9.612

9.186

9.091

8.424

8.840

8.422

9.637

10.060

9.756

9.662

8.879

8.880

8.811

9.593

10.045

9.755

9.668

8.885

8.880

8.418

2×2

THT-7 4×4

8×8

9.210

9.102

9.079

9.556

8.962

9.456

Table (8).Nondimensional central deflection ( 10w ) of a clamped supports 9-layered cross-ply  0 ,90 

(

a/h 10 100 1000

)

2

(

, 0 , 90 , 0

2×2

THT-7 4×4

8×8

2.167

2.298

2.331

0.884

0.956

0.870

0.931

)  square plate under uniform transverse load 2

BEM [27]

QUAD4 [45]

TRIPLT [46]

2.347

2.374

2.316

2.320

0.959

0.961

0.989

0.957

0.963

0.944

0.946

-

0.944

0.934

2×2

QHT-11 4×4

8×8

2.418

2.361

0.958

0.941

0.943

0.924

Table (9).Nondimensional central moment ( M x ) of a clamped supports 9-layered cross-ply  0 ,90 

(

a/h 10 100 1000

)

2

(

, 0 , 90 , 0

)  square plate under uniform transverse load 2

BEM [27]

QUAD4 [45]

TRIPLT [46]

5.917

5.517

5.66

5.660

7.051 7.083

6.608 -

6.61 -

6.616 6.655

2×2

THT-7 4×4

8×8

2×2

QHT-11 4×4

8×8

5.248

5.743

5.871

6.052

5.935

6.635 6.719

7.009 7.016

7.013 7.205

7.071 7.093

7.075 7.111

6.3. Angle-ply square laminated plate under uniform load At first, a nine-layered angle-ply ( 45 , −45 ) , 45 , ( −45 , 45 )  square symmetric laminated 2 2  plate with simply supported at all edges is analysed for different thickness ratios (a/h ). Similar to the previous example a quarter of the plate is taken in the analysis. The material properties in relation (98), is used in this analysis. The shear correction factors k s 1 and k s 2 are assumed unity. The nondimensional central deflection and moment are calculated using the proposed elements, and the results are compared with findings of other researchers. The nondimensional deflection and bending moment obtained at the plate center with those obtained by other investigators are presented in in tables (10) and (11), respectively. These tables indicate the high accuracy and rate of convergence of the results obtained by the proposed elements.

20

Table (10).Nondimensional central deflection ( 10w ) of a simply supported 9-layered angle-ply  45 , −45 

(

a/h 10 100 1000

)

2

(

, 45 , −45 , 45

2×2

THT-7 4×4

8×8

3.585

3.733

3.781

2.408

2.447

2.394

2.433

)  square plate under uniform transverse load 2

SQH [44]

MFE [47]

TRIPLT [46]

3.790

3.732

3.725

3.732

2.457

2.461

2.408

2.400

2.400

2.443

2.448

2.394

2.394

2.223

2×2

QHT-11 4×4

8×8

3.691

3.764

2.458

2.436

2.445

2.420

Table (11).Nondimensional central moment ( M x ) of a simply supported 9-layered angle-ply  45 , −45 

(

a/h 10 100 1000

)

2

(

, 45 , −45 , 45

2×2

THT-7 4×4

8×8

3.269

3.477

3.505

3.472

3.583

3.518

3.610

)  square plate under uniform transverse load 2

SQH [44]

MFE [47]

TRIPLT [46]

3.555

3.483

3.648

3.497

3.593

3.629

3.604

3.625

3.669

3.587

3.623

3.598

3.594

3.445

2×2

QHT-11 4×4

8×8

3.693

3.640

3.618

3.524

3.628

3.517

At the next stage, a simply supported three-layer angle-ply  β  , − β  , β   square plate is analyzed with three fiber orientation angles 15°, 30° and 45 °. Due to symmetry, a quadrant of the square plate is modeled. The material properties in relation (98), is employed for this analysis. Also, in all the case the thickness ratio (a/h ) is taken as 100. The calculated normalized central deflection and bending moments around x-axis and y-axis are listed in tables (12), (13) and (14), respectively. These tables show that the performance of the proposed elements is very good. Table (12).Nondimensional central deflection ( 10w ) of a simply supported 3-layered angle-ply  β  , − β  , β   square thin plate (a/h=0.01) under uniform transverse load  

Angle β 

15 30 45

THT-7

2×2

16×16

2×2

QHT-11 4×4 8×8

4×4

8×8

4.1876 3.6785

4.2403 3.9613

4.2666 4.0972

4.2783 4.1600

4.2452 3.8097

4.2551 4.0175

4.2731

16×16

4.1206

4.2824 4.1692

3.3150

3.6750

3.8664

3.9604

3.4561

3.7448

3.9009

3.9758

Table (13).Nondimensional central moment ( M x ) of a simply supported 3-layered angle-ply  β  , − β  , β   square thin plate (a/h=100) under uniform transverse load  

Angle β  15 30 45

THT-7 8×8

2×2

4×4

11.5126 7.7510

11.9444 7.8839

4.9295

4.7627

QHT-11 4×4 8×8

16×16

2×2

12.0645 8.0572

12.0769 8.1169

12.8328 8.4652

12.1371 7.9467

12.0799 8.0687

12.0840 8.1209

4.8348

4.8696

5.6555

4.7871

4.8336

4.8692

21

16×16

Table (14).Nondimensional central moment ( M y ) of a simply supported 3-layered angle-ply  β  , −β  , β   square thin plate (a/h=0.01) under uniform transverse load  

Angle β 15 30 45



THT-7

QHT-11 4×4 8×8

2×2

4×4

8×8

16×16

2×2

1.2315 2.8934

1.1734 2.9001

1.1636 2.9471

1.1581 2.9629

1.3316 3.1940

1.1633 2.9175

1.1510 2.9434

16×16 1.1526 2.9603

4.9296

4.7627

4.8350

4.8699

4.6549

4.7864

4.8329

4.8688

6.4. Simply supported skew laminated plate under uniform load A simply supported skew laminated plate as shown in figure (7) is analyzed, taking thickness ratio (a/h ) as 10 and 100. The plate is subjected to uniformly distributed load, and material properties in relation (96) are used. Also, the shear correction factors are taken as

k s 1 = k s 2 = 1 . The study is made for skew angle (α ) of 15°, 30° and 45° taking cross-ply  0 ,90 , 0  . It is worth emphasizing that this problem poses severe difficulties for numerical  

tests. The normalized central deflection of this plate is listed with those obtained by Kabir [42] and, Sheikh and Chakrabarti [48], in Table (15). Note that, Sheikh and Chakrabarti used a triangular element with 42 degrees of freedom for analysis of this structure. Based on these tables, it is obvious that the presented element has high accuracy even when coarse meshes are applied. Moreover, the obtained results from proposed elements are very close to the Kabir findings.



a

α

Figure (7).Simply supported skew plate with typical mesh 4×4

22

Table (15).Nondimensional central deflection ( 10w ) of a simply supported 3-layered cross-ply 0 , 90 , 0  square plate under uniform transverse load  

α 30 45 60

a/h 10 100 10 100 10 100

16×16

Sheikh [48]

Kabir [42]

7.9233

8.199

-

5.4717

5.4732

5.447

-

5.6196

5.5801

5.512

-

3.6738

3.6572

3.6731

3.628

-

2.5767

2.8058

2.6167

2.5857

2.454

2.60

1.5006

1.5183

1.5023

1.5216

1.453

1.52

4×4

THT-7 8×8

4×4

QHT-11 8×8

16×16

7.6754

7.8182

7.8838

8.2343

7.9671

5.4143

5.4428

5.5241

5.5279

5.4516

5.4896

5.5554

5.8783

3.5922

3.6137

3.6417

2.6631

2.5819

1.4774

1.4746

6.5. Simply supported rectangular laminated plate under doubly sinusoidal load As it is shown in Figure (8), a simply supported rectangular plate under doubly sinusoidal load with intensity of q = q 0 sin (π x a ) sin (π y b ) is analyzed for different aspect ratios (a/h ). Also, the meshes have been used for whole plate. Two different stacking sequences are analysed: 0 ,90 , 0  and  0 ,90 ,90 , 0  . The material properties in relation (96), is used for analysis. The shear correction factors k s 1 and k s 2 are assumed unity. For this load case, particular solution of the equation (46) is obtained as follows: q a 4b 4  πx  π y  πx   π y  w p = 0 4  A sin   sin   + B cos   cos   πQ   a   b   a   b 

(92)

Where

( B = ( 4ab D

A = b 4 D11 + 2a 2b 2 ( D12 + 2D 66 ) + a 4 D 22 3

16

+ 4a 3bD 26

)

) (93)

Q =C ×D

( D = (b

C = b 4 D11 + 4ab 3D16 + 2a 2b 2 ( D12 + 2D 66 ) + 4a 3bD 26 + a 4 D 22 4

D11 − 4ab 3D16 + 2a 2b 2 ( D12 + 2D 66 ) − 4a 3bD 26 + a 4 D 22

) )

At first, a square plate (a=b) is analysed. The normalized central deflections of three-ply  0 ,90 , 0  and four-ply  0 ,90 ,90 ,0  square laminate are presented in tables (16) and (17),    

respectively. The results of proposed elements, are compared with the analytical solution of Reddy [4] and Reis et al. [29]. In the second case, a three-layered rectangular laminated plate with aspect ratios (b=3a) and stacking sequence 0 ,90 , 0  is considered. The obtained nondimensional central deflection is listed in table (18) and the results are compared with those of Reddy [4], Reis et al. [29] and Pagano [49]. Note that, Reis et al. have analysed the plate using a boundary element method with quadratic elements. These tables reveal that the proposed elements have a good performance and give accurate solution.

23

πx q = q sin   a 0

 π y   sin     b 

Figure (8).Rectangular laminated plate under doubly sinusoidal load

Table (16).Nondimensional central deflection ( 10w ) of a simply supported 3-layered cross-ply 0 , 90 , 0  square plate (a=b) under doubly sinusoidal load  

a/h 4 10 100

8×8

2×2

QHT-11 4×4

8×8

Reddy [41]

14.904

16.758

17.272

17.198

17.373

17.429

17.760

-

5.559

6.307

6.537

6.174

6.501

6.588

6.693

6.397

4.263

4.301

4.306

4.297

4.312

4.312

4.337

4.361

2×2

THT-7 4×4

Reis et al. [29]

Table (17).Nondimensional central deflection ( 10w ) of a simply supported 4-layered cross-ply 0 , 90 ,90 , 0  square plate (a=b) under doubly sinusoidal load  

a/h 4 10 20 100

2×2

THT-7 4×4

8×8

2×2

QHT-11 4×4

8×8

Reddy [41]

Reis et al. [29]

14.361

16.178

16.686

16.781

16.833

16.854

17.100

-

5.505

6.240

6.461

6.127

6.430

6.511

6.628

6.231

4.187

4.660

4.824

4.455

4.776

4.858

4.912

-

4.254

4.299

4.306

4.297

4.312

4.312

4.337

4.368

Table (18).Nondimensional central deflection ( 10w ) of a simply supported 3-layered cross-ply 0 , 90 , 0  rectangular plate (b=3a) under doubly sinusoidal load  

a/h 4 10 20 100

2×2

THT-7 4×4

8×8

19.950

22.045

22.672

6.798

7.564

4.898 4.666

5.515 4.974

8×8

Reddy [41]

Pagano [49]

Reis et al. [29]

22.836

23.626

28.200

-

7.783

7.883

8.030

9.190

11.030

5.632 5.034

5.724 5.034

5.784 5.064

6.100 5.080

6.530 5.070

2×2

QHT-11 4×4

23.123

22.703

7.819

7.932

5.681 5.025

5.591 5.055

24

6.6. Triangular laminated plate under uniform load The analysis of an equilateral triangular laminated plate under the uniformly distributed load is considered in this test. To this end, a 4-layered angle-ply  −45 , 45 , 45 , −45  equilateral triangular laminated plate with clamped support at all edges (CCC) is analysed for different thickness ratios (a/h). Figure (9) shows this plate with typical mesh (N=6). It is worth mentioning that the shear correction factors k s 1 and k s 2 of the plate are assumed 5/6. Also, the material properties of each layer are taken as: Ex = 14 E y , Gxy = 0.5E y , Gxz = 0.5E y , Gyz = 0.385E y , ν xy = ν yz = 0.3

(99) By performing this test, the nondimensional deflection at centroid of the equilateral triangular plate are listed in table (19). The obtained results are compared with those of Hajikazemi et al. [50]. Theses researchers solved this problem by the triangular differential quadrature method (TDQM). Based on this tables, it can be deduced that the suggested elements have high accuracy and rate of convergence. It is also seen that good to excellent agreement is achieved in comparison with the results of TDQM for relatively thick plates.

(

3a 6,0

)

Figure (9).Equilateral triangular clamped plate with typical mesh (N=6) Table (19).Nondimensional deflection ( w ) at centroid of a 4-layered angle-ply  −45 , 45 ,45 , −45 

a/h 4 10 100

equilateral triangular laminated plate with clamped support at all edges THT-7 QHT-11 TDQM [50] N=6 N=12 N=18 N=6 N=12 N=18 7.3336

7.3733

7.3815

7.5326

7.4241

7.4042

7.3415

1.7581

1.7975

1.8067

1.8316

1.8186

1.8165

1.8035

0.6447

0.6610

0.6652

0.6624

0.6675

0.6721

0.5507

25

6.7. Clamped circular plate under uniformly distributed load In this test, a clamped unidirectional composite circular plate with radius a under the uniformly distributed load q, is analyzed for different aspect ratios (a/h). The material fibers are at an angle θ = 0 with respect to the global coordinates. Due to its symmetry, only a quarter of this plate is analyzed. The material properties are taken as: Ex = 5.6 × 106 , E y = 1.2 × 106 , Gxy = Gxz = G yz = 0.6 × 106 , ν xy = 0.26 (100) Also, the shear correction factors are taken as k s1 = k s 2 = 5 6 . The three different meshes that shown in figure (10), are used in this test. The normalized central deflection of the circular plate is evaluated as below:

wc = wc D qa 4

(101)

D = 3 ( D11 + D22 ) + 2 ( D12 + 2 D66 )

Where, D11 , D12 , D22 and D66 are bending rigidity coefficients of the laminate. The obtained nondimensional central deflection is presented in table (20) and the results are compared with those of Zhang and Kim (RDKQ-L24) [51], and Wilt et al. (SQUAD4) [52]. The exact solution for the thin plate (a/h=1000) is given by Lekhnitskii [53]. The accuracy and reliability of the proposed element answers are clearly illustrated by the findings. Based on this table, it is obvious that the presented element has high accuracy even when coarse meshes are applied. Moreover, the sensitivity of suggested element to the distortion is insignificant. It is important to note that no shear locking happens in the thin plate limit.

Mesh 1

Mesh 3

Mesh 2

Figure (10).Three mesh types for quadrant of a circular plate Table (20).Normalized central deflection wc for clamped support circular plate under uniformly distributed load with different aspect ratios of a/h

a/h

RDKQ-L24 [51]

Mesh 2 10 0.1244 16.67 0.1244 25 0.1244 50 0.1244 100 0.1245 1000* 0.1269 * Exact=0.1250 [53].

Mesh 3 0.1251 0.1251 0.1251 0.1251 0.1251 0.1259

SQUAD4 [52] Mesh 2 0.1355 0.1266 0.1237 0.1211 0.1193 0.1163

Mesh 3 0.1378 0.1291 0.1264 0.1247 0.1242 0.1231

Mesh 1 0.1147 0.1069 0.1043 0.1028 0.1024 0.1023

26

THT-7 Mesh 2 0.1320 0.1243 0.1217 0.1201 0.1197 0.1196

Mesh 3 0.1362 0.1282 0.1257 0.1242 0.1238 0.1237

Mesh 1 0.1130 0.1050 0.1023 0.1007 0.1003 0.1001

QHT-11 Mesh 2 Mesh 3 0.1313 0.1362 0.1234 0.1280 0.1208 0.1255 0.1192 0.1239 0.1187 0.1235 0.1186 0.1234

7. Conclusion In this study, the triangular element THT-7 and the quadrangular element QHT-11 were proposed for the analysis of anisotropic and symmetric laminated thick and thin plates. Formulations of these elements are based on the hybrid-Trefftz method. At first, the governing differential equation of such structures was derived based on FSDT theory and analytical solution of governing equation was obtained. For finite element formulation, the internal and boundary displacement fields of element were assumed to be independent. For the internal displacement field, homogenous and particular solutions of the governing differential equation are employed. The boundary field was related to the nodal degrees of freedom by interpolation functions. For deriving these functions, each edge of the element was modeled by two-node laminated Timoshenko beam element. In order to examine the accuracy and efficiency of the proposed elements, several numerical tests were conducted. As it was proven by the numerical results, the answers of the suggested elements quickly converged in all the tests and these free of shear locking for extremely thin laminates. Moreover, the present elements lead to good answers, even when the coarse meshes are utilized. In conclusion, the solution accuracies of the famous benchmark problems provide the justification of the suggested elements.

8. References [1]

G. Bao, W. Jiang, and J. Roberts, "Analytic and finite element solutions for bending and buckling of orthotropic rectangular plates," International journal of solids and structures, vol. 34, pp. 1797-1822, 1997.

[2]

E. Reissner, "The effect of transverse shear deformation on the bending of elastic plates," Journal of applied Mechanics, vol. 12, pp. 69-77, 1945.

[3]

R. D. Mindlin, "Influence of rotary inertia and shear on flexural motions of isotropic, elastic plates," J. of Appl. Mech., vol. 18, pp. 31-38, 1951.

[4]

J. N. Reddy, "A simple higher-order theory for laminated composite plates," Journal of applied mechanics, vol. 51, pp. 745-752, 1984.

[5]

P. Dey, A. Sheikh, and D. Sengupta, "A new element for the analysis of composite plates," Finite Elements in Analysis and Design, vol. 82, pp. 62-71, 2014.

[6]

Q.-H. Qin, "Trefftz finite element method and its applications," Applied Mechanics Reviews, vol. 58, pp. 316-337, 2005.

[7]

J. Jirousek and N. Leon, "A powerful finite element for plate bending," Computer Methods in Applied Mechanics and Engineering, vol. 12, pp. 77-96, 1977.

[8]

J. Jirousek and L. Guex, "The hybrid‐Trefftz finite element model and its application to plate bending," International Journal for Numerical Methods in Engineering, vol. 23, pp. 651-693, 1986. 27

[9]

J. Petrolito, "Triangular thick plate elements based on a hybrid-Trefftz approach," Computers & structures, vol. 60, pp. 883-894, 1996.

[10] Y. S. Choo, N. Choi, and B. C. Lee, "A new hybrid-Trefftz triangular and quadrilateral plate elements," Applied Mathematical Modelling, vol. 34, pp. 14-23, 2010. [11] M. Rezaiee-Pajand, M. Yaghoobi, and M. Karkon, "Hybrid Trefftz formulation for thin plate analysis," International Journal of Computational Methods, vol. 9, 2012. [12] M. Rezaiee-Pajand and M. Karkon, "Two efficient hybrid-Trefftz elements for plate bending analysis," Latin American Journal of Solids and Structures, vol. 9, pp. 43-67, 2012. [13] M. Rezaiee-Pajand and M. Karkon, "Two higher order hybrid-Trefftz elements for thin plate bending analysis," Finite Elements in Analysis and Design, vol. 85, pp. 73-86, 2014. [14] J. Jirousek, A. Wroblewski, and B. Szybiński, "A new 12 DOF quadrilateral element for analysis of thick and thin plates," International journal for numerical methods in engineering, vol. 38, pp. 2619-2638, 1995. [15] J. Jirousek and M. N'Diaye, "Solution of orthotropic plates based on p-extension of the hybrid-Trefftz finite element model," Computers & Structures, vol. 34, pp. 51-62, 1990. [16] J. Petrolito, "Analysis of orthotropic thick plates using hybrid-Trefftz elements," in Australian Conference of the Mechanics of Structures and Materials Melbourne, Australia, 2011, pp. 105-109. [17] J. Petrolito, "Vibration and stability analysis of thick orthotropic plates using hybridTrefftz elements," Applied Mathematical Modelling, vol. 38, pp. 5858-5869, 12/15/ 2014. [18] G. Shi and G. Bezine, "A general boundary integral formulation for the anisotropic plate bending problems," Journal of composite materials, vol. 22, pp. 694-716, 1988. [19] B. C. Wu and N. J. Altiero, "A new numerical method for the analysis of anisotropic thin-plate bending problems," Computer methods in applied mechanics and engineering, vol. 25, pp. 343-353, 1981. [20] C. Rajamohan and J. Raamachandran, "Bending of anisotropic plates by charge simulation method," Advances in Engineering Software, vol. 30, pp. 369-373, 1999. [21] C. Dong, S. Lo, Y. Cheung, and K. Lee, "Anisotropic thin plate bending problems by Trefftz boundary collocation method," Engineering analysis with boundary elements, vol. 28, pp. 1017-1024, 2004. 28

[22] E. Albuquerque, P. Sollero, W. Venturini, and M. Aliabadi, "Boundary element analysis of anisotropic Kirchhoff plates," International journal of solids and structures, vol. 43, pp. 4029-4046, 2006. [23] A. dos Reis, É. Lima Albuquerque, F. Luiz Torsani, L. Palermo Jr, and P. Sollero, "Computation of moments and stresses in laminated composite plates by the boundary element method," Engineering Analysis with Boundary Elements, vol. 35, pp. 105-113, 2011. [24] M. Wünsche, F. García-Sánchez, and A. Sáez, "Analysis of anisotropic Kirchhoff plates using a novel hypersingular BEM," Computational Mechanics, vol. 49, pp. 629-641, 2012. [25] K. Bhaskar and B. Kaushik, "Simple and exact series solutions for flexure of orthotropic rectangular plates with any combination of clamped and simply supported edges," Composite structures, vol. 63, pp. 63-68, 2004. [26] J. Wang and M. Huang, "Boundary element method for orthotropic thick plates," Acta Mechanica Sinica, vol. 7, pp. 258-266, 1991. [27] J. Wang and K. Schweizerhof, "Fundamental solutions and boundary integral equations of moderately thick symmetrically laminated anisotropic plates," Communications in numerical methods in engineering, vol. 12, pp. 383-394, 1996. [28] J. Sladek, V. Sladek, C. Zhang, J. Krivacek, and P. Wen, "Analysis of orthotropic thick plates by meshless local Petrov–Galerkin (MLPG) method," International journal for numerical methods in engineering, vol. 67, pp. 1830-1850, 2006. [29] A. dos Reis, É. Lima Albuquerque, and L. Palermo Júnior, "The boundary element method applied to orthotropic shear deformable plates," Engineering Analysis with Boundary Elements, vol. 37, pp. 738-746, 2013. [30] H.-T. Thai and S.-E. Kim, "Analytical solution of a two variable refined plate theory for bending analysis of orthotropic Levy-type plates," International Journal of Mechanical Sciences, vol. 54, pp. 269-276, 2012. [31] B. Pandya and T. Kant, "Flexural analysis of laminated composites using refined higher-order C° plate bending elements," Computer Methods in Applied Mechanics and Engineering, vol. 66, pp. 173-198, 1988. [32] C. Hwu, Anisotropic elastic plates: Springer, 2010. [33] Y.-T. Zhao, M.-Z. Wang, Y. Chen, and Y. Su, "Polynomial stress functions of anisotropic plane problems and their applications in hybrid finite elements," Acta Mechanica, vol. 223, pp. 493-503, 2012.

29

[34] G. Castellazzi, P. Krysl, and I. Bartoli, "A displacement-based finite element formulation for the analysis of laminated composite plates," Composite Structures, vol. 95, pp. 518-527, 2013. [35] J. N. Reddy, Mechanics of laminated composite plates and shells: theory and analysis: CRC press, 2004. [36] J. Belinha, L. Dinis, and R. N. Jorge, "Composite laminated plate analysis using the natural radial element method," Composite Structures, vol. 103, pp. 50-67, 2013. [37] M. Shahbazi, B. Boroomand, and S. Soghrati, "A mesh-free method using exponential basis functions for laminates modeled by CLPT, FSDT and TSDT–Part II: Implementation and results," Composite Structures, vol. 94, pp. 84-91, 2011. [38] S. Goswami, "A C 0 plate bending element with refined shear deformation theory for composite structures," Composite structures, vol. 72, pp. 375-382, 2006. [39] J. Belinha and L. Dinis, "Analysis of plates and laminates using the element-free Galerkin method," Computers & structures, vol. 84, pp. 1547-1559, 2006. [40] K. Liew, J.-B. Han, and Z. Xiao, "Differential quadrature method for thick symmetric cross-ply laminates with first-order shear flexibility," International Journal of Solids and Structures, vol. 33, pp. 2647-2658, 1996. [41] J. Reddy, "Energy principles and variational methods in applied mechanics. 2002," ed: Wiley, New York. [42] H. Kabir, "A shear-locking free robust isoparametric three-node triangular finite element for moderately-thick and thin arbitrarily laminated plates," Computers & structures, vol. 57, pp. 589-597, 1995. [43] A. H. Sheikh, S. Haldar, and D. Sengupta, "A high precision shear deformable element for the analysis of laminated composite plates of different shapes," Composite Structures, vol. 55, pp. 329-336, 2002. [44] A. K. Noor and M. D. Mathers, "Shear-Flexible Finite-Element Models of Laminated Composite Plates and Shells," DTIC Document1975. [45] B. Somashekar, G. Prathap, and C. R. Babu, "A field-consistent, four-noded, laminated, anisotropic plate/shell element," Computers & structures, vol. 25, pp. 345-353, 1987. [46] H. Lakshminarayana and S. S. Murthy, "A shear‐flexible triangular finite element model for laminated composite plates," International journal for numerical methods in engineering, vol. 20, pp. 591-623, 1984.

30

[47] G. Singh, P. Raveendranath, and G. Vekateswara Rao, "An accurate four‐node shear flexible composite plate element," International Journal for Numerical Methods in Engineering, vol. 47, pp. 1605-1620, 2000. [48] A. H. Sheikh and A. Chakrabarti, "A new plate bending element based on higher-order shear deformation theory for the analysis of composite plates," Finite Elements in Analysis and Design, vol. 39, pp. 883-903, 2003. [49] N. Pagano, "Exact solutions for rectangular bidirectional composites and sandwich plates," Journal of composite materials, vol. 4, pp. 20-34, 1970. [50] M. Hajikazemi, M. Sadr, and M. Ramezani-Oliaee, "Triangular Differential Quadrature Method in Bending Analysis of Triangular Symmetric Laminated Plates," Mechanics of Advanced Materials and Structures, vol. 22, pp. 305-312, 2015. [51] Y. Zhang and K. Kim, "Two simple and efficient displacement‐based quadrilateral elements for the analysis of composite laminated plates," International journal for numerical methods in engineering, vol. 61, pp. 1771-1796, 2004. [52] T. Wilt, A. Saleeb, and T. Chang, "A mixed element for laminated plates and shells," Computers & Structures, vol. 37, pp. 597-611, 1990. [53] S. G. Lekhnitskii, "Anisotropic plates," DTIC Document1968.

Appendix (A) The elastic constants Q ij( k ) of the kth layer are calculated in the following way:

(

)

(k ) Q11( k ) = Q11( k ) cos 4 θ k + 2 Q12( k ) + 2Q 66( k ) sin 2 θ k cos 2 θ k + Q 22 sin 4 θ k

(

)

(

(k ) Q12( k ) = Q11( k ) + Q 22 − 4Q 66( k ) sin 2 θ k cos 2 θ k + Q 12( k ) sin 4 θ k + cos 4 θ k

(

)

)

Q 22( k ) = Q11( k ) sin 4 θ k + 2 Q12( k ) + 2Q 66( k ) sin 2 θ k cos 2 θ k + Q 22( k ) cos 4 θ k

( ) ( ) Q ( ) = (Q ( ) − Q ( ) − 2Q ( ) ) sin θ cos θ + (Q ( ) − Q ( ) + 2Q ( ) ) sin θ cos θ Q ( ) = (Q ( ) + Q ( ) − 2Q ( ) − 2Q ( ) ) sin θ cos θ + Q ( ) ( sin θ + cos θ )

Q16( k ) = Q11( k ) − Q12( k ) − 2Q 66( k ) sin θ k cos 3 θ k + Q12( k ) − Q 22( k ) + 2Q 66( k ) sin 3 θ k cos θ k

Q

k 26

k 11

k 12

k 66

k 66

k 11

k 22

k 12

(k ) 44

(k ) 44

=Q

(

cos 2 θ k + Q

(k ) 55

3

k

k 66

k 12

k

2

k 22

2

k

k

k 66

3

k 66

k

4

k

(A.1)

4

k

k

sin 2 θ k

)

(k ) sin θ k cos θ k Q 45( k ) = Q 55( k ) − Q 44

(k ) sin 2 θ k + Q 55( k ) cos 2 θ k Q 55( k ) = Q 44

Where θk is the angle from global axis x to the principle material axis 1 (see figure (A.1)). In addition, Q ij( k ) for kth layer of laminate is obtained as: 31

(k )

Q11

E 1( ) , = k (k ) 1 −ν 12( )ν 21 k

Q 66( k ) = G12( k ) ,

(k )

Q12

ν12( k ) E 2( k ) , = k (k ) 1 −ν 12( )ν 21

Q 44( k ) = G 23( k ) ,

(k )

Q 22

E 2( ) , = k (k ) 1 −ν 12( )ν 21 k

(A.2)

Q 55( k ) = G13( k ) ,

Figure (A.1).Positive rotation principle material axes from global x-y axes

Appendix (B) The firs eleven analytical functions of homogeneous solution (55) for case ( µ1 ≠ µ2 ) are obtained as follows: n = 2:

w 1 = ( x + α1 y )2 − β12 y 2  w 2 = ( x + α1 y ) β1 y  2 2 2 w 3 = ( x + α 2 y ) − β 2 y

(B.1)

n = 3:

w 4  w 5  w 6  w 7

(B.2)

n = 4:

3

= ( x + α1 y ) − 3 ( x + α1 y ) β12 y 2 2

= −3 ( x + α1 y ) β1 y + β13 y 3 2

= −3 ( x + α 2 y ) β 2 y + β 23 y 3 = ( x + α2 y

3

)

− 3 ( x + α 2 y ) β 22 y 2

w 8 = ( x + α1 y )4 − 6 ( x + α1 y )2 β12 y 2 + β14 y 4  w 9 = −4 ( x + α1 y )3 β1 y + 4 ( x + α1 y ) β13 y 3  3 3 3 w 10 = −4 ( x + α 2 y ) β 2 y + 4 ( x + α 2 y ) β 2 y  4 2 2 2 4 4 w 11 = ( x + α 2 y ) − 6 ( x + α 2 y ) β 2 y + β 2 y

(B.3)

Also, for case ( µ1 = µ2 ) , the first eleven polynomial functions of homogeneous solution (56) (from 2, 3 and 4 order) have the following form: n = 2:

w 1 = ( x + α1 y ) 2 + β12 y 2  w 2 = ( x + α1 y ) β1 y  2 2 2 w 3 = ( x + α1 y ) − β1 y

(B.4)

32

w 4  w 5  w 6  w 7

n = 3:

3

= ( x + α1 y ) − 3 ( x + α1 y ) β12 y 2 2

= −3 ( x + α1 y ) β1 y + β13 y 3

(B.5)

2

= − ( x + α1 y ) β1 y − β13 y 3 3

= ( x + α1 y ) + ( x + α1 y ) β12 y 2

w 8 = ( x + α1 y ) 4 − 6 ( x + α1 y )2 β12 y 2 + β14 y 4  w 9 = −4 ( x + α1 y )3 β1 y + 4 ( x + α1 y ) β13 y 3  3 3 3 w 10 = −2 ( x + α1 y ) β1 y − 2 ( x + α1 y ) β1 y  4 4 4 w 11 = ( x + α1 y ) − β1 y

n = 4:

(B.6)

Appendix (C). The interpolation matrix of the boundary displacement field The entries of the interpolation matrix of the boundary displacement field have the coming appearance:

w
  N 11 N 12 N 13 N 14 N 15 
  θx  = N 21 N 22 N 23 N 24 N 25 
 N N 32 N 33 N 34 N 35 θy   31

(

)

N 12

(

) (

y ij

(

)

)

0.5 1 − s 2 + s 1 − s 2 ( −0.5 + δ )   4  1 =  2 − s −3 + s 2 − 2s 1 − s 2 δ   4

N 13 =

x ij

(

) (

(

)

)

(

)

 −0.5 1 − s 2 + s 1 − s 2 ( −0.5 + δ )   4  y ij  −0.5 1 − s 2 + s 1 − s 2 ( −0.5 + δ )  =  4 

N 15 = N 16

(C.1)

1 2 + s −3 + s 2 + 2s 1 − s 2 δ    4 x ij  0.5 1 − s 2 + s 1 − s 2 ( −0.5 + δ )  =  4 

N 11 =

N 14

w
i  
 θxi  N 16  θ
  yi  N 26    w
j N 36    θ
  xj  θ
  yj 

N 21 = N 22 = N 23 =

x ij 4l ij2

(

) (

)

(

) (

)

6 1 − s 2 

x ij y ij

(

(

(C.2)

) ( −1 + 2δ ) )

3 s 2 − 1 (1 − 2δ )   4l ij2 

1  2 (1 − s ) yij2 + −1 + s ( −2 + 3s ) + 6 1 − s 2 δ xij2  2    4lij

(

(

) ) 33

N 24 = N 25 = N 26 = N 31 = N 32 =

N 33 = N 34 =

N 35 =

N 36 =

xij

 −6 1 − s 2 4lij2 

(

) ( −1 + 2δ )

xij yij

)

4lij2

3 s 2 − 1 (1 − 2δ )   

(

1  2 (1 + s ) yij2 + −1 + s ( 2 + 3s ) + 6 1 − s 2 δ xij2  2    4lij

(

yij

(

xij yij

3 s 2 − 1 (1 − 2δ )   

6 1 − s 2 4lij2  4lij2

) )

(

) ( −1 + 2δ )

(

)

1  2 (1 − s ) xij2 + −1 + s ( −2 + 3s ) + 6 1 − s 2 δ yij2  2    4lij

(

yij

(

 −6 1 − s 2 4lij2 

x ij y ij

(

) )

) ( −1 + 2δ )

3 s 2 − 1 (1 − 2δ )   4l ij2 

1 4l ij2

(

)

 2 (1 + s ) x 2 + −1 + s ( 2 + 3s ) + 6 1 − s 2 δ y 2  ij ij   

(

(

) )

34