Exponential basis functions in the solution of laminated plates using a higher-order Zig–Zag theory

Exponential basis functions in the solution of laminated plates using a higher-order Zig–Zag theory

Accepted Manuscript Exponential Basis Functions in the Solution of Laminated Plates Using a HigherOrder Zig-Zag Theory F. Azhari, B. Boroomand, M. Sha...

664KB Sizes 0 Downloads 24 Views

Accepted Manuscript Exponential Basis Functions in the Solution of Laminated Plates Using a HigherOrder Zig-Zag Theory F. Azhari, B. Boroomand, M. Shahbazi PII: DOI: Reference:

S0263-8223(13)00232-8 http://dx.doi.org/10.1016/j.compstruct.2013.05.022 COST 5166

To appear in:

Composite Structures

Please cite this article as: Azhari, F., Boroomand, B., Shahbazi, M., Exponential Basis Functions in the Solution of Laminated Plates Using a Higher-Order Zig-Zag Theory, Composite Structures (2013), doi: http://dx.doi.org/ 10.1016/j.compstruct.2013.05.022

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.

Exponential Basis Functions in the Solution of Laminated Plates Using a Higher-Order Zig-Zag Theory F. Azhari a , B. Boroomand a,1, M. Shahbazi b a

b

Department of Civil Engineering, Isfahan University of Technology, Isfahan 84156-83111, Iran Department of Civil Engineering, University of British Columbia, Vancouver V6T 1Z4, Canada

Abstract The solution method presented in this paper is of Trefftz type which uses exponential basis functions (EBFs) to solve composite plate problems. Following its success in the solution of plates using well-known theories (Composite Struct. 2011;93:3112-9, 94:84-91and 2012; 94:2263– 68), here we aim to apply the method to higher order shear deformation theories. In this paper we

demonstrate the way that one can generally evaluate the EBFs for a laminated plate using a ZigZag theory. We present explicit relations for three-layer sandwich plates which have a wide range of applications in structural/mechanical engineering field. The results of our numerical experiments on the bending analysis of composites with different boundary conditions and different configurations are provided and compared with those available in the literature. For further sudies we present some results for sandwich composite plates, including those with soft cores, as new benchmarks using the Zig-Zag theory.

Keywords: mesh-free; laminated plate; shear deformation plate theories; Zig-Zag theory; exponential basis function; discrete transformation; sandwich plates

1. Introduction Nowadays it is well understood that an efficient approach for the analysis of composite laminates requires a suitable theory for assuming the through-thickness variation of the state variables as well as a suitable numerical method for the analysis and approximation of the variables throughout the area representing the mid surface. Among the theories used for the throughthickness variation of the state variables, those whose unknowns are of displacement type can be 1 Corresponding author. E-mail address: [email protected].

1

divided into three main groups; equivalent single layer models (ESL), layer-wise models (LW) and Zig-Zag theories. In the ESL models, the laminate is replaced by an equivalent single-layer anisotropic plate and the assumed displacements are distributed continuously across the thickness of the laminate. Apart from the classical plate theory (CLPT) based on Kirchhoff’s assumptions, another well-known theory in this category is the first-order shear deformation theory (FSDT) proposed by Reissner [1] and Mindlin [2]. The reader may refer to [3, 4] for a more complete literature. The theory may not predict accurate results for very thick plates. In order to increase the accuracy, higher-order shear deformation theories (HSDT) have been proposed [5-8]. Among these theories, Reddy’s third-order shear deformation theory (TSDT) is the most prevalent and efficient one [8]. More accurate models have then been introduced as the LW ones [9-15]. In these methods, a displacement field within each layer is prescribed. However, they are computationally expensive because the number of unknowns depends on the number of layers. It can be concluded that an ESL model with the characteristics similar to the LW models but with the number of unknowns independent of the number of layers may serve as an efficient theory. To this end different Zig-Zag theories have so far been introduced [16-19]. In these theories the continuity of the transverse shear stresses at the interface of each pair of layers is assured (see [20] for the history of the Zig-Zag methods). In this paper, as a sample of Zig-Zag methods, a higher-order Zig-Zag theory proposed in [19] is used to perform the analysis of laminated plates. In this theory, a displacement field with piece-wise linear variation is combined with another field with cubic variation. In the literature there are some studies focusing on the theory among which the readers may refer to studies in [21, 22]. In the realm of numerical analyses, various methods have been developed over the past decades the most well-known of which are the finite element method (FEM) and the boundary element method (BEM). The reader may refer to [23, 24] for the history of using the FEM and more general formulations for multilayered composite plates. The FEM and the BEM have a long history in the solution of problems in solid mechanics [25, 26]. However, problems regarding the expense of meshing, the connectivity of the elements and the compatibility of the state variables have led to thinking of using mesh-free methods in composite plates. Among these methods, one may refer to the element free Galerkin method (EFG) [27] used for the analysis of plates in [28, 29]. Furthermore, the studies by Ferreira et al employing radial basis functions (RBFs) for the analysis of composite laminated plates [30, 31] should be mentioned here. There is also another 2

group of mesh-free methods using a set of points just at the boundaries. The method of fundamental solutions (MFS) proposed by Kupradze and Aleksidze [32] is classified in the latter category. The BEM and the MFS are sometimes classified as the “Trefftz” type of methods [33, 34]. In this study, we shall use a mesh-free method proposed in [35] and employed for the analysis of composite laminates in [36-40]. This method is based on the use of a series of exponential basis functions (EBFs) to satisfy the governing equation (see [41-44] for its application to other engineering problems). The coefficients are found by the imposition of the boundary conditions through a collocation approach using a discrete transformation technique introduced in [35, 45, 46].Therefore, this method may also fall in the category of “Trefftz” methods and is capable of analyzing laminates with a variety of shapes and boundary conditions. In order to use the method for the analysis of composites modeled by Zig-Zag theories, we first find the corresponding EBFs. For sandwich plates with in-plane isotropy, we present the details of the procedure with explicit relations to be used as a set of library functions. In this paper, we first overview the Zig-Zag displacement field in Section 2; afterward, we elaborate on the governing equations and the boundary conditions through a variational formulation in Section 3. In section 4, we explain the procedure of finding the EBFs for a sampling plate problem. In Section 5, our numerical experiments are presented and the accuracy of the results in comparison with the available exact solutions is discussed. In Section 6, we summarize the conclusions made throughout the paper.

2. The displacement field In this paper a composite plate with N orthotropic layers and total thickness of h is analyzed based on a higher-order Zig-Zag theory [19]. The displacement field, in the kth layer ( k = 1,..., N ), may be written in a vector notation as ⎧⎪u k = u 0 ( x, y ) − z w′( x, y ) + [F( z )]k u z ( x, y ) ⎨ ⎪⎩ w = w0 ( x, y )

⎛ F11 ( z )

[F ( z )]k = ⎜

F12 ( z ) ⎞

⎟ , ⎝ F21 ( z ) F22 ( z ) ⎠k

where u k = [u, v]Tk , u 0 = [u0 , v0 ]T , w′ = [∂w0 / ∂x, ∂w0 / ∂y ]T and u z = [u z , vz ]T , so that

(u0 , v0 , w0 ) are the mid-plane displacements and (u z , vz ) are the Zig-Zag displacement terms. 3

(1)

According to the assumed displacement field, the unknown vector to be found is

[u0 , v0 , w0 , u z , vz ]

T

. The Zig-Zag part of the displacements in (1), i.e.[F( z )]k , may be expressed

as

[F( z )]k = {H1k + z 2 H 2 + z 3I 2×2 } ,

(2)

where I is an identity matrix and H1k and H 2 are obtained by the procedure given in Appendix A.

3. The governing equations and boundary conditions In this section, we proceed to find the governing equations through a variational approach (see also [24] for further details in more general cases). In a static state, the variation of the total potential of the composite plate is expressed as

δΠ = δ U + δ V = 0 ,

(3)

in which h 2 h − 2

δ U = ∫ {∫ [σ xxδε xx + σ yyδε yy + 2σ xyδε xy + 2σ xzδε xz + 2σ yzδε yz ]dz}d Ω0 , Ω0

(4)

In the above relations, Ω 0 denotes the domain of the mid-plane defined in x-y. Also σ and ε represent the components of the stress and strain tensors, respectively. The second term in (3),

δ V , pertains to the variation of the potential of the applied loads (we shall refer to this part when boundary conditions are to be derived). By substituting the strain-displacement relations in (4) as well as the constitutive relations and performing integral by parts, a variational expression in terms of δ u0 , δ v0 , δ w0 , δ u z and δ vz is resulted. The expressions conjoint with the displacement variations represent the governing differential equation as Lu = q

in Ω0

(5)

where L is an operator matrix, not elaborated here for the sake of brevity, u = [u0 , v0 , w0 , u z , vz ]

T

is the unknown vector to be found, and q is a vector containing the loading components defined 4

as q = [ 0, 0, qz , 0, 0] . It is important to note that when L is operated on a displacement field of a T

plate with symmetric layer sequence, some of its elements vanish. Therefore the operator matrix may be condensed to Lsub (3×3) as Lsubu = q ,

(6)

where u is now defined as [ w0 , u z , vz ] , and q is the corresponding loading vector defined T

as [ qz , 0, 0] . T

3.1 Boundary conditions

For derivation of the boundary condition, we focus on δ V in (3) with a general expression as

δ V = δ VΩ + δ V∂Ω . 0

(7)

0

In the above relation ∂Ω0 represents the boundary of Ω0 . Also δ VΩ0 denotes part of δ V pertaining to the transverse load while δ V∂Ω0 denotes part of δ V affected by the boundary conditions. Regardless of the type of the boundary conditions, δ VΩ0 may be written as

δ V∂Ω = −∫ 0

s ≡∂Ω0

h

{∫ 2h [σˆ nnδ un + σˆ nsδ us + σˆ nzδ w]dz}ds , −

(8)

2

in which Neumann and Dirichlet conditions appear as a set of conjugate pairs. In the above relation σˆ nn , σˆ ns and σˆ nz are the stress components on ∂Ω0 in the direction of n and s, normal and tangent to the boundary. For Neumann conditions the stresses are prescribed but at this stage there is no need to distinguish them from Dirichlet ones. In (8) δ un and δ us denote the variation of the components of u in (1) evaluated in the directions of n and s. By defining

δ u n = [δ un δ us ] , δ u = [δ u δ v ] , T

T

(9)

one may write

δ u n = nδ u ,

⎡ nx n=⎢ ⎣−ny

ny ⎤ , nx ⎥⎦

(10)

5

where nx and n y are the components of the unit vector normal to the boundary. In view of (1), for kth layer we have

δ u kn = δ u 0 n − z δ w′n + [F( z )]k δ u nz ,

δ w′ = [∂ (δ w0 ) / ∂n, ∂ (δ w0 ) / ∂s ]T .

(11)

By using (11) in (8), and evaluating the integrals in z direction, we arrive at δ V∂Ω = − ∫ 0

s ≡∂Ω0

{δ u

0n

}

N nn + δ u0 s N ns + δ u zn M nn* + δ u zs M ns* + [ ∂ (δ w0 ) / ∂n ] M nn − [ ∂ (δ w0 ) / ∂s ] M ns + δ w0 Qn ds

(12) The signs in the above relation depend on the positive direction of each component of the conjugate pairs.

Free boundary conditions

By considering δ un and δ us as arbitrary values in (9), the variation of displacement components

δ u0 n , δ u0 s , δ w0 , δ u zn , δ u zs , ∂ (δ w0 ) / ∂n and ∂(δ w0 ) / ∂s take on arbitrary values (see (12)). The boundary conditions needed for free edges are then obtained as * * N nn = 0 , N ns = 0 , M nn = 0 , M ns = 0 , M nn = 0 , M ns = 0 , Qn = 0 .

(13)

However, in view of (12), and similar to the classical plate theory (CLPT), one may evaluate the sixth term in (12) by the use of integration by parts and combine it with the seventh term to define a quantity as Vn = Qn +

∂M ns . ∂s

(14)

In that case the conditions at free edges will be as * * N nn = 0 , N ns = 0 , M nn = 0 , M ns = 0 , M nn = 0 , Vn = 0 .

6

(15)

Simply supported conditions In this case the variation of displacement components δ u0 n , δ u zn , and ∂ (δ w0 ) / ∂n may take on arbitrary values. This leads to the following mixed conditions * N nn = 0 , u0 s = 0 , w0 = 0 , M nn = 0 , u zs = 0 , M nn = 0 .

(16)

Clamped conditions The following conditions are intuitively considered for clamped edges u0 n = 0 , u0 s = 0 , w = 0 , u zn = 0 , u zs = 0 , ∂w0 / ∂n = 0 , ( ∂w0 / ∂s = 0 ).

(17)

Note that by prescribing w0 as the third conditions, we have ∂w0 / ∂s = 0 and thus the last condition in (17) is automatically satisfied. It can be shown that a paradox exists in the formulation when using the above conditions (see [39] for similar discussion in other shear deformation theories). However, the discussion is beyond the scope of this paper.

4. The solution method In this study we use a Trefftz method, the full account of which has been presented for CLPT, FSDT and TSDT in [37, 38]. We split the solution into two parts, i.e. the homogeneous solution denoted by subscript “h”, and the particular solution denoted by subscript “p”, so that u = u h + u p . Therefore (see [37, 38] for more details)

Lu h = 0 ,

Lu p = q in Ω0

(18)

ˆ  = u Bu h Bh

ˆ ) on ∂Ω 0 , (u Bh = u B − Bu p

(19)

In the above relations, u B is a vector composed of the prescribed boundary conditions and Bˆ is the boundary operator matrix.

7

The particular solution u p is found in manner similar to the procedure explained in [35, 37, 38]. Here we elaborate on the homogenous part u h which is assumed as a summation of a series of exponential basis functions (EBFs) with a generic expression as below u h ( x, y, α , β ) = h (α , β ) eα x + β y ,

{

where, h(α , β ) = hu

hv

hw

∀(x, y) ∈ Ω0 and (α , β ) ∈ » 2 ,

hu z

h vz

}

T

(20)

. Note that we have started from a displacement field

with a complete set of elements as used in (5). By inserting u h in the homogeneous equation (18) a system of equations in terms of α and β is resulted as Λ (α , β ) h ( α , β ) = 0 .

(21)

In order to obtain a non-trivial solution for the above equation, the determinant of the matrix Λ (α , β ) must be set to zero. In that case h(α ,β ) is a vector in the null-space of matrix Λ (α , β ) .

Considering a numerical scheme (see [35, 37] for more explanations) the homogeneous part of the solutions is written as below  h = ∑ ci h (αi , βi ) eαi x + βi y . u

(22)

i

with ci being unknown coefficients to be determined from the information at the boundaries.

4.1 EBFs for Sandwich plates with in-plane isotropy:

In this section we present the essential and explicit relations for the evaluation of the EBFs for a class of sandwich plates. The schematic presentation of the layer sequence and the properties are shown in Figure 1. The system of equations for such a sandwich plate is in the form of (6).

Please see the figures Fig. 1. The sandwich plate with in-plane isotropy.

8

Introduction of (20) in (6) results in the following relation ⎛ D11α i4 + (2 D12 + 4 D66 )α i2 β i2 + D22 β i4 ⎜ − E11Z α i3 − ( E12Z + 2 E66Z )α i β i2 ⎜ ⎜ −( E12Z + 2 E66Z )α i2 β i − E22Z β i3 ⎝

− E11Z α i3 − ( E12Z + 2 E66Z )α i β i2

−( E12Z + 2 E66Z )α i2 β i − E22Z β i3 ⎞ ⎧ hiw ⎫

E α + E β − J 55

( E + E )α i β i

( E + E )α i β i

E α i2 + E22E β i2 − J 44

E 11

2 i

E 66

E 12

2 i

E 12

E 66

Z 66

E 66

⎧0 ⎫ ⎟ ⎪ uz ⎪ ⎪ ⎪ . (23) ⎟ ⎨hi ⎬ = ⎨0 ⎬ ⎟ ⎪ h vz ⎪ ⎪0 ⎪ ⎠⎩ i ⎭ ⎩ ⎭

In the above relation N

( Dijkl , E , E , J 3i 3 j ) = ∑ ∫ Z ijkl

E ijkl

m =1

zm+1

zm

(m) ijkl

Q

2 ⎛ 2 ⎞ 2 ⎛ ∂ ( F11 ) m ⎞ ⎜⎜ z , ( ( F11 ) m z ) , ( F11 ) m , ⎜ ⎟ ⎟⎟ dz ∂ z ⎝ ⎠ ⎝ ⎠

(24)

where zm and zm +1 represent the lower and upper z coordinate of the mth lamina, respectively. (m) Also Qijkl and Q3(im3)j denote the in-plane and transverse transformed plane stress stiffness

components, respectively, noting that i, j, k, l= 1, 2 .The components of the fourth-order tensor (m) (m) Q ( m ) can be distinguished by a double-indicial notation for instance Q1112 , A2211 ≡ A21 , ≡ Q16

A3131 ≡ A55 etc. The simplified forms of such quantities, for a sandwich plate shown in Figure 1, are given in Appendix B. The characteristic equation in this case is as D11α i4 + (2 D12 + 4 D66 )α i2 β i2 + D22 β i4

− E11Z α i3 − ( E12Z + 2 E66Z )α i β i2

−( E12Z + 2 E66Z )α i2 β i − E22Z β i3

− E11Z α i3 − ( E12Z + 2 E66Z )α i β i2

E11Eα i2 + E66E β i2 − J 55

( E12E + E66Z )α i β i

−( E + 2 E )α β i − E β

( E + E )α i β i

E α + E β − J 44

Ψ (α i , β i ) =

Z 12

Z 66

2 i

Z 22

3 i

E 12

E 66

E 66

2 i

E 22

2 i

(25) which yields to Ψ (α i , βi ) = (α i2 + βi2 ) 2 (C1 (α i2 + βi2 ) 2 + C2 (α i2 + βi2 ) + C3 ) = 0 ,

(26)

where C1 =

1 1 h 4C3′ ( 2 C1′ − 336C2′ ) , 7 2 λ 150528000 A Bλ

(27)

3(C4′ ) 2 C5′ 1 2 20 2 ′ ′ ′ ′ ′ C2 = − h C4 ( 2 C1 − 672C2 − 2 2 C3C5 ) , C3 = − , 1600 A5 B 3 ABλ 2 λ 1792000 A6 Bλ B λ AB while

9

(28)

=0

A = 1 +ν 2 , A = ν 2 − 1 , B = 1 + ν 1 , B = ν 1 − 1 , λ = E1 / E2 ,

(29)

and ⎛ 8 C1′ = 5 A ⎜ − (1260 B 2δ 2 (1 + δ ) 2 + 420 ABδ (1 + δ ) 2 (5 + (δ − 1)δ )λ ⎝ A + A2 (1088 + 7δ (244 + δ (96 + 5δ (10 + 8δ + δ 3 ))))λ 2 )



(30)

4 Aδ 3λ (420 B 2 (1 + δ ) 2 + 168 ABδ 2 (1 + δ )λ + 17 A2δ 4 λ 2 ) ⎞⎛ 8 + 12δ + 6δ 2 δ 3λ ⎞ − ⎟⎜ − ⎟, BB AA BB ⎠ ⎠⎝ 2

⎛ 30 Bδ (1 + δ ) 2 + 32 Aλ + 5 Aδλ (10 + δ (4 + δ + δ 2 )) Aδ 3 (5 B (1 + δ ) + Aδ 2 λ ) ⎞ C2′ = ⎜ + ⎟ , λA BB ⎝ ⎠

C3′ = 2520 B 3δ 2 (1 + δ ) 2 + 420 AB 2 λδ (1 + δ ) 2 (10 − 2δ + 3δ 2 ) + 2 A2 Bλ 2 (1088 + 1708δ +672δ 2 + 350δ 3 + 280δ 4 + 84δ 5 + 119δ 6 ) + 17 A3δ 7 λ 3 ,

(31)

(32)

C4′ = 30 B 2δ (1 + δ ) 2 + 2 ABλ (16 + 25δ + 10δ 2 + 5δ 3 + 5δ 4 ) + A2δ 5λ 2 ,

(33)

C5′ = AAδ 3λ + 2 BB (4 + 6δ + 3δ 2 ) ,

(34)

Note that the height of the top and bottom layers is taken as h while the height of the mid layer is δ h (see Figure 1) and thus the height of the composite is

h = (δ + 2)h .

(35)

By evaluating βi in terms of α i in (26), the following roots are evaluated

βi = ±iα i (double roots for each sign), βi = ±

−C2 + C22 − 4C1C3 − 2C1α i2 2C1

, βi = ±

h (αi , βi ) = {1, 0, 0 } , T

−C2 − C22 − 4C1C3 − 2C1α i2 2C1

(36)

(37)

For each βi in (37), after introducing the plate properties in (27) to (34), h (αi , βi ) is obtained by finding the null-space of the coefficient matrix in (23).

10

For the missing EBFs in (36), we may assume another set of EBFs, as the product of polynomial terms and the former exponential functions (see [35, 37, 42] for more details). The final expression is obtained as 0 ⎛ ⎧1 ⎫ ⎡⎧1 ⎫ ⎧1 ⎫ ⎧ ⎫⎤ ⎧1 ⎫ ⎜ 1 ⎪ ⎪ α i ( x + iy ) 2 ⎪ ⎪ α i ( x − iy ) 3 ⎢ ⎪ ⎪ ⎪ ⎪ ⎪ 2 ⎪ ⎥ α ( x + iy )  h = ∑ ⎜ ci ⎨ 0 ⎬ e u + ci ⎨ 0 ⎬ e + ci ⎢ ⎨ 0 ⎬ x + ⎨ 0 ⎬ y + ⎨Cα i (1 + i ) ⎬⎥ e i i ⎜ ⎪0 ,⎪ ⎪0 ,⎪ ⎪ 2 ⎪⎥ ⎢ ⎪⎩0 ,⎪⎭ ⎪⎩0 ,⎪⎭ ⎩ ⎭ ⎩Cα i (i − 1) ⎭⎦ ⎣ ⎝ ⎩ ⎭

+ ci4

0 ⎡⎧1 ⎫ ⎧1 ⎫ ⎧ ⎫⎤ ⎢⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎥ αi ( x −iy ) 2 + ⎢ ⎨ 0 ⎬ x + ⎨ 0 ⎬ y + ⎨ Cα i (1 − i ) ⎬⎥ e ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ 2 ⎢ ⎩0 ,⎭ ⎩0,⎭ ⎥ ⎩−Cα i (i + 1) ⎭⎦ ⎣ αi x −

+ci5h5e

− C2 + C22 − 4C1C3 − 2C1αi2 y 2C1

αi x −

− C2 − C22 − 4 C1C3 − 2 C1αi2

+c h7 e 7 i

2 C1

αi x +

+ ci6h 6e y

αi x +

+c h e 8 i 8

− C2 + C22 − 4C1C3 − 2C1αi2 y 2C1

− C2 − C22 − 4 C1C3 − 2 C1α i2 2 C1

y

+

⎞ ⎟ ⎟, ⎟ ⎠

(38)

in which, 2(δ λ AA + 32λ ABB + 20δ BB (3 B + λ A) + 5δ λ AB (1 + δ )( B + A ) + 10δ BB (3 B + 5λ A) + 30δ B B )) 5

C=

2

2

2

3

3

2

3 AB (δ λ A + 32λ AB + 10δ λ AB + 10δ B (3 B + λ A)(2 + δ ) + 10δ B (3 B + 5λ A)) 5

2

2

4

2

(39)

In (38), h5 to h8 are easily evaluated by substituting (37) in (23) and searching for the nullspace of the matrix in each case.

4.2 Imposition of boundary conditions

Having found the EBFs, as the next step we find the unknown coefficients ci in (22). By satisfying the boundary conditions through a collocation method over a set of points on the boundary ∂Ω0 , the coefficients ci are found. Therefore, considering (22), the coefficients ci must be determined so that Ne

1

i =1

i

∑s

ci v i = u B h ,

(40)

11

where N e is the number of EBFs in the series which is not necessarily equal to the number of points, M, and u B h is a vector containing the modified prescribed boundary conditions at each boundary point. Also in (40), vi is the contribution of the ith EBF to the boundary conditions. In (40) si is the norm of vector vi (see [35] or [37] for other forms of the normalization required for Neumann boundary conditions). The formulation is similar to that explained in [35, 37] and thus we avoid repeating here for the sake of conciseness. It may suffice to mention that in (40) the coefficients ci are found as ci =

1 T vi R h uB h , si

(41)

in which R h plays the role of a projection matrix. By substitution of (41) in (40) one may evaluate the projection matrix as ⎛ Ne 1 ⎞ T ⎜⎜ ∑ 2 v i v i ⎟⎟ R h u B h = u B h , ⎝ i =1 si ⎠

+

⎛ Ne 1 ⎞ R h = ⎜ ∑ 2 v i vTi ⎟ . ⎜ s ⎟ ⎝ i =1 i ⎠

(42)

In the above relations (.) + denotes the pseudo-inverse of the matrix. The final solution is as N ⎫⎪ ⎞ ⎪⎧⎛ e 1 u h = Re ⎨⎜ ∑ h (αi , βi ) eαi x + βi y vTi ⎟ R h u Bh ⎬ , ⎠ ⎩⎪⎝ i =1 si ⎭⎪

(43)

where, Re{.} stands for the real part of the complex value. In [35] two strategies are presented for choosing α and β for the homogenous solution. In this paper we follow strategy II (see Appendix C) for the numerical experiments.

5. Numerical experiments Here, the bending solution is presented for several isotropic/orthotropic laminates with different boundary conditions shapes and loadings. In all examples, the material properties of the orthotropic laminates are E2 = 10GPa , E1 = 25 E2 , G12 = G13 = 0.5 E2 , G23 = 0.2 E2 , ν 12 = ν 23 = ν 13 = 0.25 ,

(44)

and the material properties used for isotropic laminates are as follows E = 30GPa , G = E / 2(1 + ν ) , D = Eh3 / 12(1 −ν 2 ) , ν = 0.3 .

12

(45)

5.1. Square plates Square plates with side length of a and total thickness of h are analyzed. The coordinate origin is placed at the center of the plates and in the mid-plane. The thicknesses of all layers are the same. In addition, the following non-dimensional deflection and bending stresses are defined to present the results w= (k ) σ xy

100h3 E2 q0 a 4 =

h2 q0 a 2

(k ) w(0, 0), σ xx =

h2 q0 a

σ xy (a / 2, a / 2, z ),

(k ) σ (0, 0, z ), σ yy = 2 xx

(k ) σ yz

h2 q0 a 2

σ yy (0, 0, z ),

(46)

h h (k ) = σ yz (0, −a / 2, z ), σ xz = σ xz (−a / 2, 0, z ) , q0 a q0 a

where the superscript k denotes the layer number. As a reference solution for the simply supported composites, we find a Navier type of solution by assuming

u0 = ∑∑ u0mn , v0 = ∑∑ v0mn , w0 = ∑∑ w0mn , u z = ∑∑ u zmn , vz = ∑∑ vzmn , m

n

m

n

m

n

m

n

m

)47(

n

where m and n denote the mode numbers in x and y directions, respectively, and u0mn = uomn cos(

mπ x0 nπ y0 mπ x0 nπ y0 ) sin( ) , v0mn = v0mn sin( ) cos( ), a a a a

(48)

w0mn = w0mn sin(

mπ x0 nπ y0 mπ x0 nπ y0 ) sin( ) , u zmn = u zmn cos( ) sin( ), a a a a

(49)

vzmn = vzmn sin(

mπ x0 nπ y0 ) cos( ) a a

.

(50)

while x0 = x − a / 2 and y0 = y − a / 2 . In the above relations the unknown amplitudes of the displacements uomn , v0mn , w0mn , u zmn and vzmn are determined such that the following equation is satisfied.

Lu mn = q mn with q mn = ⎡⎣0 0 q mn

)51( T

0 0 ⎤⎦ and u mn = ⎡⎣u0mn

v0mn

13

w0mn

u zmn

vzmn ⎤⎦

T

where

q mn = q mn sin(

mπ x0 nπ y0 ) sin( ), a a

(52)

with qmn =

4 a2

a a

∫∫q

z

sin(

0 0

mπ x0 mπ x0 ) sin( )dx0 dy0 . a a

(53)

The rest of the procedure is straightforward. In order to use the method proposed in Section 4, an appropriate set of parameters α (or β ) is needed. For generating the EBFs for the homogenous solution we have followed one of the strategies in [35], i.e. strategy II explained in Appendix C (using γ = 2π and M = N = 3 ). This leads to 432 and 288 EBFs for anti-symmetric and symmetric plates, respectively. As shown in Figure 2, in this problem the boundaries of the plates are discretized with 100 points (25 evenly spaced points in each side) and a 20 × 20 grid of points is used to discretize the domain of the plates for the particular solution. In Table 1, the dimensionless center deflection of simply supported (four edges) and clamped (four edges) anti-symmetric cross-ply square laminates with layer orientation of ( 0° ,90° , 0° ,90° ) subjected to a uniformly distributed transverse load is presented for different span to thickness ratios. The results are compared with those given in [21] based on the model used in [19] and the spline finite strip method. For simply supported plates, the results of Navier solution are also reported. As is observed the results of the proposed method are in excellent agreement with those of Navier solution.

Please see the figures

(a)

(b)

(c)

Fig. 2. The square plate in Example 5.1: (a) plate geometry; (b) boundary points; (c) domain points used for the particular solution.

To give more insight to the model and the solution method employed, we compare the results obtained by the proposed method with those in [3] obtained by various methods including a Zig14

Zag model proposed in [13]. Table 2 contains the non-dimensional transverse shear stress for orthotropic three- and nine-layer square laminates subjected to a sinusoidal loading. Again simply supports (four edges) are assumed for the laminates. The layer sequence is symmetric with respect to the mid-plane and fiber orientations alternate between 0° and 90° as (0°,90°,0°) for the three-layer and (0°,90°,0°,90°,0°,90°,0°,90°,0°) for the nine-layer laminate. The total thickness of the 0° layers is equal to that of 90° layers and the layers with similar orientations have similar thicknesses. The results for various side-to-thickness ratios are also compared with the 3D solutions provided in [47]. The acronyms ESL-1, ESL-2 and ESL-3 stand for, respectively, results obtained from equivalent single-layer models using linear, parabolic, and cubic displacements through thickness direction defined in [13] and used in [3]. The ESL models reported in [3] in fact use a Zig-Zag variation of through-thickness displacements. Likewise, the acronyms LW-1, LW-2 and LW-3 refer to, respectively, layer-wise models using linear, parabolic, and cubic displacements and displacements and stresses through the thickness of the layers (see [48, 49]). We have reported the results in [3] obtained with no post processing, e.g. no additional integrals for evaluating shear stresses, since in the presented method we just find the stresses from the constitutive equations. Moreover we have included the results of Navier solutions in the table. As is observed the results of the proposed method are in excellent agreement with those of Navier solution. For the three-layer composites, N = 3 , the results of the Zig-Zag model used are comparable with those of ESL-1 and LW-1models but are inferior to those obtained from other models. However, for nine-layer composites, N = 9 , the results of the Zig-Zag model used in this paper are superior to those of ESL models (and even superior to LW1 model). However, still in this latter case, the results are inferior to those obtained by LW-2 and LW-3. The reader may note that the Zig-Zag model used in this paper may be categorized in the ESL models. Therefore it is not surprising that the results are less accurate than those obtained by LW models. Nevertheless, categorized as an ESL model, the model used in this paper is still capable of producing reasonably accurate results in comparison with other ESL models as is seen in Table 2. Considering the computational time for LW models the ESL models, the model used in this paper may be considered as a good and efficient one for the analysis of composite plates (see also [3] for similar discussion). Table 2 shows that this effect becomes prominent when relatively high number of layers are to be used in a composite plate.

15

5.2. Skew plates In this section isotropic skew plates with total thickness of h, the inclination angle of ψ and length of a = 20m for all four sides, subjected to a uniform transverse load with q0 = 50 kPa intensity are analyzed. For the particular solution, the domain of the plates is

discretized by a grid of points so that the horizontal/vertical distance between adjacent points is 1m. The boundaries are discretized by 120 points (30 evenly spaced points in each side). The geometry of the plates, the coordinate origin, the sequence of the boundary edges and the distribution of the points on the boundary and the domain are displayed in Figure 3. For generating the EBFs of the homogenous solution we have adopted similar parameters as used in Example 5.1, yielding to total number of 288 bases.

Please see the figures

(a)

(b)

(c)

Fig. 3. The skew plates in Example 5.2: (a) plate geometry; (b) boundary points; (c) domain points used for the particular solution.

Table 3 shows the non-dimensional deflection and moment resultants at the centroid of simply supported (SSSS) isotropic skew plates with different inclination angles. In this table the following normalized quantities are used for simplicity w=

qa4 qa2 qa2 w, M max = M max , M min = M min , 1600D 40 40

in which D =

(54)

Eh3 . 12(1 −ν 2 )

As is observed, the results of the present model are in good agreement with those in reference [50, 51]. Note that the results reported in reference [51] are in fact obtained by a finite element (FE) analysis. Nevertheless, in order to provide a sound basis for the comparison we have performed a 3D finite element analysis using solid elements. To this end, a number of 3200 ( 40 × 40 × 2 ) quadratic elements are employed for the analysis. As expected, the results obtained by the FE model are in good agreement with those reported in [51] for all values of ψ . Taking

16

the FE analysis as a basis for the comparisons, it can be seen that for ψ = 75o (or even for

ψ ≥ 45o ) the results of the present model are superior to those of the reference [50]. However some discrepancies are seen between the results of the present model and those obtained by the FE analyses when acute angles are used for ψ . Interesting to note is that although a first order shear deformation theory is used in [50] and [51], some discrepancies are seen between the results reported in these references. The reason may be traced in the fact that when acute angles are chosen for ψ , two strong singular points appear at corners with obtuse angles. It is a priori known that methods such as the finite element method are prone to the pollution error (see [52] for instance). However, the solution method presented in this paper, and also that in [50], may also be affected by such singularities. The reason lies in the fact that the EBFs are in fact smooth functions and do not generally exhibit singular behavior. To alleviate the effect one may add some functions exhibiting singular behavior at the corners of the domain. It may be worthwhile to mention that with appropriate arrangement of EBFs one may define special forms of nearlysingular functions (see [46] for instance), however, the discussion is beyond the scope of this paper. 5.3. Sandwich plates Sandwich plates are of special interest in many industries. The early history of the theories used for their solution can be found in [53]. Most important application of the sandwich plates is in the form of plates with soft cores. The latest history of the theories used may be found in references [54-56 ]. The application of the theory proposed in [19] to sandwich plates, including those with soft core, may be found in [57, 58]. In this paper we present further results using the closed form relations given in Section 4.1. The results are reported for two sets of benchmark problems. The first set consists SSSS square sandwich plates with isotropic Aluminum core ( E = 70 GPa, ν = 0.3 with thickness of h / 2 ) and isotropic Ceramic faces ( E = 151GPa, ν = 0.3 with thickness of h / 4 ) subjected to a sinusoidal transverse load with q0 intensity. The domain of the plates is discretized with a

20 × 20 grid of points (for the particular solution), the boundaries are discretized with 100 equally spaced points and the coordinate origin is placed at the center of the plates in the midplane. Table 4 presents the non-dimensional center deflection and stresses of the afore17

mentioned sandwich plates with different side-to-thickness ratios. In this table, the following quantities are used. w= (k ) σ xy

100h3 Ec q0 a 4 =

10h 2 q0 a 2

(k ) w(0, 0), σ xx =

10h 2

10h 2

q0 a

q0 a 2

σ xy (a / 2, a / 2, z ),

(k ) σ (0, 0, z ), σ yy = 2 xx

(k ) σ yz

σ yy (0, 0, z ),

(55)

h h (k ) = σ yz (0, −a / 2, z ), σ xz = σ xz (−a / 2, 0, z ), q0 a q0 a

where Ec is the core material modulus of elasticity. The reported results are to be used as benchmark values for further studies. In order to give some insight to the validity of the results, those obtained by other theories for the center deflection of the plates are given in Table 5. We have included the results obtained by a FE analysis using 3200 ( 20 × 20 × 8 ) quadratic elements as a numerical 3D analysis. Table 5 also contains the results evaluated by a Navier type of solutions including those obtained by the relations given in (47) to (53). Again taking the 3D analysis as a basis for the comparisons, it can be seen that the results obtained in this paper (either by the EBFs or the Navier solution) are superior to those obtained by other methods especially when the plate is moderately thick. However, still some discrepancies are seen between the results of this paper and those obtained by the 3D model for moderately thick plates. Acknowledging the fact that the Zig-Zag model used in this paper falls in the category of ESL models, the discrepancies are considered acceptable in an engineering standpoint. The readers may note that for more accurate solutions one may need to use layerwise theories which are much more expensive in a computational point of view (see [55, 56]). The second set of problems includes SSSS plates with soft cores. Plates with isotropic Calcium alginate core ( E = 217.789 MPa, ν = 0.3 with thickness of 5h / 6 ) and isotropic Aluminum faces ( E = 68970 MPa, ν = 0.3 with thickness of h /12 ) subjected to a sinusoidal transverse load with

q0 intensity are analyzed. Table 6 shows the center deflection of the plates obtained by a 3D finite element analysis, the present method using EBFs, and a series of Navier solutions corresponding to various theories including the one used in this paper. The FE model is similar to the previous one and consists of eight layers of elements two of which are used for the two faces. As is seen in Table 6, although all models can accurately predict the deflections when the plate is thin, some discrepancies are observed between the results of the 3D model and all other ones when the plates become moderately thick. Surprisingly, the results of FSDT and TSDT 18

models are close to the ones obtained for thin plates even when the plates is moderately thick. This effect is not clearly seen in Table 5, however due to the presence of the soft core the effect becomes noticeable in Table 6. The table clearly shows that the results of the Zig-Zag model used in this paper are the most accurate ones though some discrepancies are still seen between these results and those of the 3D model. Considering the cost of the computation, which is similar to FSDT and TSDT, the Zig-Zag model used in this paper is of interest by engineers.

6. Conclusions In this paper, we have presented a mesh-free method based on the use of exponential basis functions (EBFs) to solve composite problems by employing a Zig-Zag theory. The results of our numerical experiments have been given for composite plates with different boundary conditions and different shapes including rectangular and skew configurations. The results have been compared with those of known references. The method has shown the capability of yielding results with excellent accuracy.

Acknowledgement The author wish to thank Mr. Zibasokhan for his valuable help in the finite element modeling of the composite plates.

Appendix A

In this appendix the procedure of evaluating vectors H1k and H 2 in (2) is explained. According to Cho and Parmerter’s theory, the Zig-Zag part of the displacement field may be considered as u kzig = [F ( z )]k u z ( x, y ) . To derive the elements of [F( z )]k we start from the following Zig-Zag

displacement u kzig = N k u k + z 2 a + z 3b ,

(A-1)

where

19

⎡ N1 Nk = ⎢ ⎣0

0

N1

N1

0

⎧ uk ⎫ ⎪v ⎪ 0⎤ ⎧ a1 ⎫ ⎧ b1 ⎫ 1 1 ⎪ k ⎪ a = , b = , u = N (1 ), N (1 ) = − ξ = + ξ , , ⎨ ⎬ ⎨ ⎬ k ⎨ ⎬. 1 2 ⎥ N1 ⎦ 2 2 ⎩a2 ⎭ ⎩b2 ⎭ ⎪uk +1 ⎪ ⎪⎩ vk +1 ⎪⎭

(A-2)

We note that u k , a and b do not depend on z but can be considered as functions of x and y . Moreover uk and vk are the bottom unknowns of the in-plane displacements of the kth lamina and

ξ represents the normalized z coordinate of the laminate which varies from -1 to 1 in each lamina. Considering (1), the transverse shear stress vector of the kth lamina is as ˆ k L′u k , with τ k = ⎡σ τ =Q zig ⎣ xz k

ˆ k = ⎡Q55 σ yz ⎤⎦ k , Q ⎢Q ⎣ 45 T

Q45 ⎤

k

0 ⎤ ⎡ ∂ / ∂z while L′ = ⎢ , ⎥ ∂ / ∂z ⎥⎦ Q44 ⎦ ⎣ 0

(A-3)

where, Qij ( i or j = 4,5 ) denotes the transformed plane stress stiffness components. By introducing (A-1) in (A-3) one may rewrite the above relation as ˆ k (B k u k + B′a + B′′b ) , with B k = L′( N k ) , B′ = L′( z 2 ) , B′′ = L′( z 3 ) . τk = Q

(A-4)

In order to define [F( z )]k in (1), the shear stress equilibrium is satisfied between the layers as below τ1 |ξ =−1 = 0, τ k |ξ =1 − τ k +1 |ξ =−1 = 0, τ N |ξ =1 = 0 ,

(A-5)

for k = 1,..., N − 1 where N is the total number of layers. We consider u kzig= J |z =0 = 0 in which u kzig = [F ( z )]k u z ( x, y ) denotes the Zig-Zag part of the displacements and J in the superscript denotes the ordinal number of the layer in which z J × z J +1 ≤ 0 , noting that zk represents the z coordinate of the bottom of the kth lamina. After satisfaction of the equilibrium equations (A-5), the following system of equations is resulted

20

⎡⎣ A ( 2 N + 4)×( 2 N + 4)

⎧u ⎫ ⎪ ⎪  D( 2 N + 4)×2 ⎤⎦ ⎨ a ⎬ = 0 , with u = {u1 ⎪b ⎪ ⎩ ⎭

v1

... u N +1

vN +1} , T

(A-6)

The system of equations of (A-6) can be rewritten as ⎧u ⎫ ⎧u ⎫ A ⎨ ⎬ = − Db ⇒ ⎨ ⎬ = − A -1Db . ⎩a ⎭ ⎩a ⎭

(A-7)

By splitting the above equations we may have u k = −( A −1 ) k4×( 2 N + 4) Db ,

(A-8)

a = −( A −1 ) 2a×( 2 N + 4) Db .

(A-9)

In (A-8) the array ( A -1 ) k4×(2 N + 4) contains the four rows of matrix A −1 which are multiplied by the vector Db to form the vector u k (see (A-2)). Also ( A -1 ) a2×(2 N + 4) contains the two rows of the matrix A −1 which are multiplied by the vector Db to form the vector a. By substituting (A-8) and (A-9) in (A-1) the following relation is resulted u kzig = ( H1k + z 2 H 2 + z 3I 2×2 )b ,

(A-10)

where −1 a H1k = − N k ( A −1 ) k4×( 2 N + 4) D , H 2 = − ( A ) 2×( 2 N + 4) D .

(A-11)

Since b can still be considered as a function of x and y , we rename it as b = u z ( x, y ) .

(A-12)

In view of (A-10) we write u kzig = ( H1k + z 2 H 2 + z 3I 2×2 )u z ( x, y ) .

(A-13)

This means that in the present Zig-Zag theory one may write [F( z )]k as the form given in (2). For a given layer sequence, the elements of H1k and H 2 are derived and stored as library functions usable for all composites with similar layer sequence.

21

Appendix B

In this appendix, the simplified elements of matrices D, E Z , E E and J , used in (23) and defined in (24), for a sandwich plate shown in Figure 1 are given as the following relations. The elements of matrix D : A

D11 = D1′ + D2′ , D12 = ν 2 D1′ + ν 1 D2′ , D22 = D11 , D66 = −(

2

B

D1′ +

D2′ ) , D16 = D26 = 0 ,

2

(B-1)

where D1′ = −

Eh 3 (4 + 6δ + 3δ 2 ) 6 AA

, D2′ = −

Eh 3δ 3λ 12 BB

.

(B-2)

See (29) for definition of A , A , B , B and λ . The elements of matrix E Z : Z E11Z = D3′ + D4′ , E12Z = ν 2 D3′ + ν 1 D4′ , E22Z = E11Z , E66 = −(

A 2

D3′ +

B 2

D4′ ) , E16 = E26 = 0 , Z

Z

(B-3)

where D3′ =

Eh 5 (30δ (1 + δ ) 2 B + (32 + 5δ (10 + δ (4 + δ + δ 2 )))λ A) 20λ A2 A

, D4′ =

Eh 5δ 3 (5 B (1 + δ ) + δ 2 λ A) 20 ABB

.

(B-4)

The elements of matrix E E : E E11E = D5′ + D6′ , E12Z = ν 2 D5′ + ν 1 D6′ , E22E = E11E , E66 = −(

A 2

D5′ +

B 2

D6′ ) , E16 = E26 = 0 , E

E

(B-5)

where D5′ = − Eh 7 (1260 B 2δ 2 (1 + δ ) 2 + 420 ABδ (1 + δ ) 2 (5 − δ + δ 2 )λ + A2 (1088 + 1708δ + 672δ 2 +350δ 3 + 280δ 4 + 35δ 6 )λ 2 ) / 280λ 2 A3 A, D6′ = −

Eh 7δ 3 (17 A2δ 4 λ 2 + 168 ABδ 2 (1 + δ )λ + 420 B 2 (1 + δ ) 2 ) 560λ A2 BB

(B-6)

.

The elements of matrix J : J 44 = D7′ + D8′ , J 55 = J 44 , J 45 = 0 ,

(B-7)

where D7′ =

3Eh 5 (16 + 25δ + 10δ 2 ) 10 A

, D8′ =

3Eh 5δ (30 B 2 (1 + δ ) 2 + 10 ABδ 2 (1 + δ )λ + A2δ 4 λ 2 ) 20λ BA2

22

.

(B-8)

Appendix C

The heuristic strategy for choosing α and β proposed in [35], i.e. strategy II, is as follows

α =±

⎞ mγ ⎛ k ⎜ + i ⎟ , m = 1,..., M , k = 1,..., N , i = −1 , L ⎝N ⎠

(C-1)

for β = ± i α . In a similar manner, when α = ± i β we select

β =±

⎞ mγ ⎛ k ⎜ + 2 i ⎟ , m = 1,..., M , k = 1,..., N , L ⎝N ⎠

(C-2)

In the above relations, we choose (M , N ) ∈ N 2 , γ ∈ R and L is a characteristic length. The following bounds are found to be appropriate for many cases

5.6 ≤ γ ≤ 7.2 (say γ = 2π ), L = 1.6 max(L x , L y ) , M min = 4 , N min = 2 , N max = 8

(C-3)

where Lx and Ly are the dimensions of the rectangle that circumscribes the domain.

References [1] Reissner, E. The effect of the transverse shear deformation on the bending of elastic plates. ASME J Appl Mech 1945;12:69-77.

[2] Mindlin RD. Influence of rotatory inertia and shear on flexural motions of isotropic, elastic plates. J Appl Mech 1951;18:31-8. [3] Carrera E. Single- vs multilayer plate modelings on the basis of Reissner’s mixed theorem. AIAA 2000;38:342-52.

[4] Carrera E. Developments, ideas and evaluations based upon the Reissner's Mixed variational theorem in the modeling of multilayered plates and shells. App Mech Rev 2001;54:301-29.

[5] Whitney JM, Sun CT. A refined theory for laminated, anisotropic, cylindrical shells. J Appl Mech 1974;41:471-476.

[6] Gautham BP, Ganesan N. Free vibration analysis of thick spherical shells. Comput Struct 1992;45:307-313.

23

[7] Levinson M. An accurate, simple theory for statics and dynamics of elastic plates. Mech Res Commun 1980;7:343-350. [8] Reddy JN. A simple higher-order theory for laminated composite plates. ASME J Appl Mech 1984;51:745-52. [9] Hsu T, Wang, JT. A theory of laminated cylindrical shells consisting of layers of orthotropic laminate. AIAA 1970;8:2141-2146. [10] Srinivas, S. A refined analysis of composite laminates. J Sound Vib 1973; 30: 495-507. [11] Robbins DHJ, Reddy JN. Modeling of thick composites using a layer-wise theory. Int J Numer Methods Eng 1993;36:655-677. [12] Reddy JN. Generalization of two-dimensional theories of laminated composite plates. Commun Appl Numer Methods 1987;3:173-180. [13] Toledano A. and Murakami, H. A composite plate theory for arbitrary laminate configuration. J Appl Mech 1987;54:181-9. [14] Lu X, Liu D. An interlaminar shear stress continuity theory for both thin and thick composite laminates. J Appl Mech 1992;59:502-9. [15] Lee CY, Liu D. An interlaminar stress continuity theory for laminated composite analysis. Comput Struct 1992;42:69-78. [16] Di Sciuva M. Development of an anisotropic multilayered shear deformable rectangular plate element. Comput Struct 1985;21:789-96. [17] Di Sciuva M. An improved shear-deformation theory for moderately thick multilayered anisotropic shells and plates. J Appl Mech 1987;54:589- 96. [18] Murakami H. Laminated composite theory with improved in-plane respons. J Appl Mech 1986;53:661-6. [19] Cho M, Parmerter RR. Efficient higher order composite plate theory for general lamination configurations. AIAA 1993;31(7):1299–306. [20] Carrera E. Historical review of Zig-Zag theories for multilayered plates and shells. App Mech Rev 2003;56:287-308. [21] Akhras G, Li W. Spline finite strip analysis of composite plates based on higher-order zigzag composite plate theory. Compos Struct 2007;78:112–118.

24

[22] Akhras G, Li W. Stability and free vibration analysis of thick piezoelectric composite plates using spline finite strip method. Int J Mech Sciences 2011;53:575–584. [23] Carrera E. Theories and finite elements for multilayered, anisotropic, composite plates and shells. Archives of Computational Methods in Engineering 2002;9:87–140. [24] Carrera E, Ciuffreda A. A unified formulation to assess theories of multilayered plates for various bending problems. Compos Struct 2005;69:271-93. [25] Zienkiewicz OC, Taylor RL. The finite element method, (5th edn), Butterworth-Heineman, 2000. [26] Brebbia CA, Dominguez J. Boundary Elements an Introductory Course (2nd edn). WIT Press, Computational Mechanics Publications: Southampton, 1998. [27] Belytschko T, Lu YY, Gu L. Element-free Galerkin methods. Int J Numer Methods Eng 1994;37:229-56. [28] Donning BM, Liu WK. Meshless methods for shear-deformable beams and plates. Comput Meth Appl Mech Eng 1998;152:47-71. [29] Belinha J, Dinis L. Analysis of plates and laminates using the element-free Galerkin method. Comput Struct 2006;84:1547-59. [30] Ferreira AJM. A formulation of multiquadric radial basis function method for the analysis of laminated composite plates. Compos Struct 2003;59:385–92. [31] Ferreira AJM, Roque CMC, Martins PALS. Analysis of composite plates using higher order shear deformation theory and a finite point formulation based on the multiquadric radial basis function method. Composites Part B 2003;34:627–36. [32] Kupradze VD, Aleksidze MA. The method of functional equations for the approximate solution of certain boundary value problems. USSR Comput Math Math Phys 1964;4:82–126. [33] Cho HA, Golberg MA, Muleshkov AS, Li X. Trefftz methods for time dependent partial differential equations. Comput Mater Continua 2004; 1:1–37. [34] Dong CY, Lo SH, Cheung YK, Lee KY. Anisotropic thin plate bending problem by Trefftz boundary collocation method. Eng Anal Boundary Elem 2004;28:1017–24. [35] Boroomand B, Soghrati S, Movahedian B. Exponential basis functions in solution of static and time harmonic elastic problems in a meshless style. Int J Numer Methods Eng 2010;81:971-1018. [36] Shamsaei B, Boroomand B. Exponential basis functions in solution of laminated structures. Compos Struct 2011;93:2010-19. 25

[37] Shahbazi M, Boroomand B, Soghrati S. A mesh-free method using exponential basis functions for laminates modeled by CLPT, FSDT and TSDT – Part I: Formulation. Compos Struct 2011;93:31129. [38] Shahbazi M, Boroomand B, Soghrati S. A mesh-free method using exponential basis functions for laminates modeled by CLPT, FSDT and TSDT – Part II: Implementation and results. Compos Struct 2011;94:84-91. [39] Shahbazi M, Boroomand B, Soghrati S. On using exponential basis functions for laminates Modeled by CLPT, FSDT and TSDT; Further tests and results. Compos Struct 2012;94: 2263–68. [40] Boroomand B, Azhari F, Shahbazi M. On definition of clamped conditions in TSDT and FSDT; The use of exponential basis functions in solution of laminated composites. Compos Struct 2013; 97 129–135. [41] Zandi SM, Boroomand B, Soghrati S. Exponential basis functions in solution of incompressible fluid problems with moving free surfaces. J. Comput. Physics 2012; 231: 505–527. [42] Zandi SM, Boroomand B, Soghrati S. Exponential basis functions in solution of problems with fully incompressible materials: A mesh-free method. J. Comput. Physics 2012; 231: 7255–7273. [43] Hashemi SH, Boroomand B, Movahedian B. Exponential basis functions in space and time; A meshless method for 2D time dependent problems. J. Comput. Physics 2013; 241: 526–545. [44] Movahedian B, Boroomand B, Soghrati S. A Trefftz method in space and time using exponential basis functions; application to direct and inverse heat conduction problems. Engg Anal Boundary Elem (in press) http://dx.doi.org/10.1016/j.enganabound.2013.03.001. [45] Boroomand B, Mossaiby F. Generalization of robustness test procedure for error estimators. Part I: formulation for patches near kinked boundaries, Int J Numer Methods Eng 2005; 64: 427–460. [46] Boroomand B, Mossaiby F. Dynamic solution of unbounded domains using finite element method: Discrete Green’s functions in frequency domain. Int J Numer Methods Eng 2006;67:1491–530. [47] Pagano NJ, Hatfield S J. Elastic behavior of multilayered bidirectional composites. AIAA 1972;10: 931–33. [48] Carrera E. Evaluation of layerwise mixed theories for laminated plates analysis. AIAA 1998; 36: 830–9.

26

[49] Carrera E. Layer-wise mixed models for accurate vibration analysis of multilayered plates,” J Applied Mech 1998; 65: 820–8. [50] Liu FL, Liew KM. Differential cubature method for static solutions of arbitrarily shaped thick plates. Int J Solids Struct 1998;35:3655-3674. [51] Sengupta, D. Performance study of a simple finite element in the analysis of skew rhombic plates. Comput Struct 1995;54:1173-1182. [52] Babuška I, Li L. The problem of plate modeling: Theoretical and computational results. Comput Meth Appl Mech Eng 1992;100: 249-273. [53] Pandya BN, Kant T. Higher-order shear deformable theories for flexure of sandwich plates – finite element evaluations. Int J Solids Structures 1988; 24: 1267-1286. [54] Liu Q, Zhao Y. Natural frequency analysis of a sandwich panel with soft core based on a refined shear deformation model. Compos Struct 2006; 72: 364–374. [55] Moreira RAS, Rodrigues JD. A layerwise model for thin soft core sandwich plates. Comput Struct 2006; 84: 1256–1263. [56] Moreira RAS, Rodrigues JD. Static and dynamic analysis of soft core sandwich panels with throughthickness deformation. Compos Struct 2010; 92: 201–215. [57] Chakrabarti A, Sheikh AH. A new triangular element to model inter-laminar shear stress continuous plate theory. Int J. Numer Meth Engng 2004; 60:1237–1257. [58] Pandit MK, Sheikh AH, Singh BN. An improved higher order zigzag theory for the static analysis of laminated sandwich plate with soft core. Finite Elem Anal Design 2008; 44: 602 – 610.

Caption of Figures Fig. 1. The sandwich plate with in-plane isotropy. Fig. 2. The square plate in Example 5.1: (a) plate geometry; (b) boundary points; (c) domain points used for the particular solution. Fig. 3. The skew plates in Example 5.2: (a) plate geometry; (b) boundary points; (c) domain points used for the particular solution. 27

Tables Table 1. Non-dimensional center deflections w of anti-symmetric cross-ply (0o, 90o, 0o, 90o) laminated plates under a uniform transverse load based on the present higher-order Zig-Zag deformation theory.

`

Four edges simply

Method

4

10

20

100

Four edges clamped

supported

Present (EBFs)

2.9786

1.6842

Spline Finite Strip [21]

2.9788

1.6733

Navier

2.9786

---

Present (EBFs)

1.1706

0.4948

Spline Finite Strip [21]

1.1705

0.4872

Navier

1.1706

---

Present (EBFs)

0.8992

0.2625

Spline Finite Strip [21]

0.8993

0.2599

Navier

0.8992

---

Present (EBFs)

0.8122

0.1798

Spline Finite Strip [21]

0.8122

0.1797

Navier

0.8122

---

Table 2. Non-dimensional transverse shear stresses

σ xz (0)

of simply supported (SSSS conditions) symmetric

three- and nine-layer plates under a sinusoidally distributed transverse load.

N=3

Theory

N=9

Method

a/h =2

4

10

100

a/h =2

4

10

100

EBFs

0.1630

0.2298

0.3138

0.3519

0.2092

0.2357

0.2557

0.2665

Navier

0.1631

0.2298

0.3138

0.3519

0.2092

0.2357

0.2557

0.2665

ESL-1

Ref. [3]

0.159

0.228

0.315

0.355

0.146

0.143

0.155

0.165

ESL-2

Ref. [3]

0.147

0.217

0.303

0.342

0.125

0.138

0.156

0.167

ESL-3

Ref. [3]

0.157

0.225

0.310

0.350

0.197

0.232

0.259

0.273

LW-1

Ref. [3]

0.167

0.229

0.302

0.336

0.193

0.220

0.245

0.257

LW-2

Ref. [3]

0.155

0.221

0.302

0.339

0.204

0.223

0.247

0.259

LW-3

Ref. [3]

0.154

0.219

0.301

0.339

0.204

0.223

0.247

0.259

3D

Ref. [47]

0.153

0.219

0.301

0.339

0.204

0.223

0.247

0.259

Zig-Zag

28

Table 3. Non-dimensional deflections and moment resultants at the centroid point of simply supported SSSS uniformly loaded equilateral skew plates with isotropic material

(h / a = 0.01) .

ψ

Method

w

M max

M min

75o

Present-Zigzag FSDT [50] FSDT [51] 3D FEM

5.8114 5.7341 5.8468 5.8636

1.9162 1.9061 1.9241 ---

1.7026 1.6741 1.7097 ---

60o

Present-Zigzag FSDT [50] FSDT [51] 3D FEM

3.8814 3.7591 4.1123 4.1377

1.6355 1.6492 1.7075 ---

1.2755 1.2418 1.3391 ---

Present-Zigzag FSDT [50] FSDT [51] 3D FEM

1.7114 1.8528 2.1330 2.1322

1.1086 1.1988 1.2995 ---

0.7386 0.7431 0.8866 ---

Present-Zigzag FSDT [50] FSDT [51] 3D FEM

0.4675 0.5449 0.6690 0.6597

0.5940 0.6799 0.7734 ---

0.3308 0.3321 0.4481 ---

45

o

30o

Table 4. Non-dimensional deflections and stresses of SSSS sandwich plates with Aluminum core and Ceramic face layers under a sinusoidally distributed transverse load obtained by the presented method.

a/h

w

(3) σ xx (h / 2)

(2) σ yy ( h / 6)

(1) σ xy (−h / 2)

(2) σ yz (0)

(2) σ xz (h / 6)

5

1.92452

0.32989

0.07040

0.17803

0.21974

0.21892

10

1.52645

0.47922

0.07894

0.25806

0.22109

0.22027

20

1.42612

0.51687

0.08109

0.27832

0.22145

0.22063

50

1.39797

0.52743

0.08170

0.28400

0.22155

0.22073

100

1.39394

0.52894

0.081783

0.28481

0.22157

0.22074

29

Table 5. Non-dimensional deflection w of SSSS sandwich plates with Aluminum core and Ceramic face layers under a sinusoidally distributed transverse load.

Navier solutions

a/h

FEA (3D)

Present

Zig-Zag*

CLPT

FSDT

TSDT

10

1.67945

1.52645

1.52645

1.39260

1.49273

1.40095

20

1.49693

1.42612

1.42612

1.39260

1.41764

1.39469

50

1.42346

1.39797

1.39797

1.39260

1.39661

1.39294

100

1.40320

1.39394

1.39394

1.39260

1.39360

1.39269

1000

1.38723

1.39262

1.39262

1.39260

1.39261

1.39260

* The results are obtained by Navier solution given in (47) to (53).

Table 6. Non-dimensional deflection w of SSSS sandwich plates with Calcium alginate core and Aluminum face layers under a sinusoidally distributed transverse load.

Navier solutions

a/h

FEA (3D)

Present

Zig-Zag*

CLPT

FSDT

TSDT

10

0.1665565

0.148699

0.148699

0.020915

0.023864

0.021185

20

0.0599519

0.053291

0.053291

0.020915

0.021652

0.020983

50

0.0292966

0.026114

0.026114

0.020915

0.021033

0.020926

100

0.0238882

0.022216

0.022216

0.020915

0.020845

0.020918

1000

0.0209558

0.020928

0.020928

0.020915

0.020916

0.020915

* The results are obtained by Navier solution given in (47) to (53).

30

Figure 1

Figure 2-a

Figure 2-b

Figure 2-c

Figure 3-a

Figure 3-b

Figure 3-c