Strong sampling surfaces formulation for laminated composite plates

Accepted Manuscript Strong sampling surfaces formulation for laminated composite plates G.M. Kulikov, S.V. Plotnikova PII: DOI: Reference:

S0263-8223(17)30628-1 http://dx.doi.org/10.1016/j.compstruct.2017.03.061 COST 8382

To appear in:

Composite Structures

Received Date: Revised Date: Accepted Date:

23 February 2017 13 March 2017 15 March 2017

Please cite this article as: Kulikov, G.M., Plotnikova, S.V., Strong sampling surfaces formulation for laminated composite plates, Composite Structures (2017), doi: http://dx.doi.org/10.1016/j.compstruct.2017.03.061

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.

Strong sampling surfaces formulation for laminated composite plates G.M. Kulikov∗, S.V. Plotnikova Laboratory of Intelligent Materials and Structures, Tambov State Technical University, Sovetskaya, 106, Tambov 392000, Russia

Abstract The paper focuses on the application of the sampling surfaces (SaS) method to exact solutions of elasticity for laminated plates. The SaS method is based on choosing the arbitrary number of SaS parallel to the middle surface in order to introduce the displacements of these surfaces as basic plate unknowns. Such choice of unknowns with the use of the Lagrange polynomials in assumed approximations of displacements through the layer thicknesses leads to a compact laminated plate formulation. The feature of the proposed approach is that all SaS are located inside the layers at Chebyshev polynomial nodes. The use of outer surfaces and interfaces is avoided that makes possible to minimize uniformly the error due to the Lagrange interpolation. Therefore, the strong SaS formulation based on direct integration of the equilibrium equations of elasticity can be applied efficiently to the obtaining of 3D exact solutions for laminated plates.

Keywords: Elasticity; Exact solution; Laminated plate; Sampling surfaces method

∗

Corresponding author. Tel.: +7 4752 630441 E-mail address: [email protected] (G.M. Kulikov)

2 1. Introduction The three-dimensional (3D) exact solutions for laminated plates have a great importance due to the fact that the validity of approximate plate theories and plate finite elements can be assessed by comparing their predictions with exact solutions. The closed-form solution of elasticity for a simply supported single-layer rectangular plate under sinusoidally distributed transverse loading was presented by Vlasov [1]. The extensions of Vlasov’s solution to laminated composite plates were done by Pagano [2, 3] and Srinivas and Rao [4, 5]. The exact solutions [2-5] were generalized to the thermal stress analysis of simply supported laminated plates in [6-8]. Pagano [9] obtained the exact solution for the cylindrical bending of laminated composite plates with general layups. The response of simply supported and clamped anisotropic laminates in cylindrical bending under temperature loading is studied analytically by Vel and Batra [10]. The antisymmetric angle-ply composite plates in the framework of 3D elasticity and thermoelasticity by using the Pagano's approach and the state space approach are considered in [11-14]. The exact solutions for anisotropic laminates in cylindrical bending and antisymmetric angle-ply plates through the use of the recently developed method of sampling surfaces (SaS) are described in [15, 16]. The more difficult problem is to find the 3D exact solutions for laminated plates made of inhomogeneous or functionally graded (FG) materials. This is due to the fact that the material properties depend on the transverse coordinate and some specific assumptions concerning their variations in the thickness direction are required [17, 18]. Therefore, the state space approach can not be applied to 3D exact solutions for FG plates directly. Nevertheless, the obtaining of 3D exact solutions for FG plates is still possible. To solve the problem, one can assume that the material properties are distributed through the thickness of the plate according to the exponential law [19-21]. The efficient analytical approach to the 3D exact

3 solutions of thermoelasticity has been proposed by Vel and Batra [22]. They studied the simply supported FG rectangular plates with the material properties presented by Taylor series expansions through the thickness coordinate. The SaS method was also applied to exact solutions for laminated rectangular plates composed of FG materials [23-26]. The feature of the SaS approach [25, 26] consists in the use of the higher-degree Lagrange polynomials in through-thickness approximations of the material properties. Therefore, it does not matter what type of the material law is adopted because only the knowledge of numerical values of the material properties on SaS is required. A new approach to closed-form solutions for isotropic and transversely isotropic FG plates is presented in [27-29]. This approach is based on a general solution of the equilibrium equations of inhomogeneous elastic media [30]. Another popular approach to 3D exact solutions is an asymptotic approach that was applied to inhomogeneous and laminated plates and shells with simply supported edges [3136]. To analyze composite plates and shells with different boundary conditions the asymptotic approach can be combined with meshless [37-39] and finite element [40, 41] methods. Note that all available approaches to 3D analytical and numerical solutions for laminated composite structures are considered and critically assessed in survey articles [42, 43]. As already said, in the literature there is a powerful tool to assess the response of laminated plates in the framework of elasticity. This is a variational SaS formulation [15] developed recently by the authors for the 3D stress analysis of laminated plates. The SaS method is based on choosing inside the nth layer I n arbitrarily located SaS parallel to the middle surface in order to introduce the displacements of these surfaces as basic plate unknowns, where I n ≥ 3 . Such choice of unknowns with the use of Lagrange polynomials of

degree I n − 1 in the through-thickness approximations of displacements of the nth layer leads to a compact form of the governing equations of the SaS laminated plate formulation.

4 The origins of the SaS concept can be traced back to contributions [44-46] in which the equally spaced SaS are utilized. The more general variational approach for homogeneous plates and shells with the SaS located at Chebyshev polynomial nodes [47] was developed later [48, 49] because of the fact that the SaS formulation with equispaced SaS does not work properly with the higher-order Lagrange interpolation. The use of only Chebyshev polynomial nodes [48, 49] improves significantly the behavior of the higher-degree Lagrange polynomials because such a choice makes possible to minimize uniformly the error due to the Lagrange interpolation. This fact gives an opportunity to obtain the displacements and stresses with a prescribed accuracy employing the sufficiently large number of SaS. It means in turn that the solutions based on the SaS concept asymptotically approach the 3D exact solutions of elasticity as the number of SaS tends to infinity. However, the implementation of the SaS variational formulation for laminated plates [15, 16, 25, 26, 50] without using the interfaces is not possible. This restriction does not allow us to have all benefits of the higherorder Lagrange interpolation with Chebyshev polynomial nodes. The present paper is intended to extend the SaS variational formulation to the strong SaS formulation for laminated plates. The latter is based on the choice of all SaS inside the layers at Chebyshev polynomial nodes and direct integration of the equilibrium equations of elasticity. The use of interfaces is avoided that allows one to minimize uniformly the error due to the higher-order Lagrange interpolation. Thus, the strong SaS formulation can be applied efficiently to the obtaining of 3D exact solutions for laminated composite plates.

2. Kinematic description of laminated plate

Consider a laminated plate of the thickness h. Let the middle surface Ω be described by Cartesian coordinates x1 and x2 . The coordinate x3 is oriented in the thickness direction. According to the SaS concept, we choose inside the nth layer I n SaS Ω ( n )1 , Ω( n ) 2 ,..., Ω ( n ) I n

5 parallel to the middle surface. The transverse coordinates of SaS of the nth layer located at Chebyshev polynomial nodes are written as x3( n )in =

 2i − 1  1 [ n −1] 1  , ( x3 + x3[ n ] ) − hn cos π n 2 2  2I n 

(1)

where x3[ n−1] and x3[ n ] are the transverse coordinates of interfaces Ω[ n −1] and Ω[n ] depicted in Fig. 1; hn = x3[ n] − x3[ n −1] is the thickness of the nth layer; the index n identifies the belonging of any quantity to the nth layer and runs from 1 to N, where N is the number of layers; the indices in and jn , k n to be introduced later identify the belonging of any quantity to the SaS of the nth layer and runs from 1 to I n ; N SaS = ∑ I n is the total number of SaS. n

The strains of the n th layer εij(n) are given by ( n) 2εαβ = uα( n, β) + u (βn,α) ,

2εα( n3) = βα( n ) + u3( n,α) ,

( n) ε33 = β3( n) ,

(2)

where ui(n ) are the displacements of the nth layer; βi(n) are the derivatives of displacements with respect to the thickness coordinate: βi( n) = ui(,n3) .

(3)

Here and throughout this paper, Latin indices i, j , k , l range from 1 to 3, whereas Greek indices α , β range from 1 to 2; the symbol (…),i stands for the partial derivatives with respect to coordinates xi . Introduce displacements of SaS of the nth layer ui( n) in ( x1, x2 ) as basic plate unknowns by ui( n )in = ui( n ) ( x3( n) in ) .

(4)

The strains of SaS of the nth layer ε ij( n) in ( x1 , x2 ) are defined as

εij( n)in = εij( n) ( x3( n)in ).

(5)

6 The use of Eqs. (2)-(5) leads to relations between the SaS variables ( n ) in 2εαβ = uα( n, β)in + u (βn,α) in , ( n )i n 2εα( n3) in = βα( n) in + u3( n,α)in , ε33 = β3( n) in ,

(6)

where βi( n )in ( x1, x2 ) are the values of derivatives of the displacements of the nth layer with respect to the thickness coordinate at SaS: βi( n )in = βi( n ) ( x3( n )in ).

(7)

Up to this moment no assumptions concerning the displacement field have been made. We start now with the fundamental assumption of the proposed 3D laminated plate formulation. Let us assume that the displacements are distributed through the thickness of the nth layer in the following form: ui( n) = ∑ L( n) in ui( n)in ,

(8)

in

where L( n) in ( x3 ) are the Lagrange basis polynomials of degree I n − 1 corresponding to the nth layer: L( n) in =

∏

jn ≠ in

x3 − x3( n ) jn . x3( n)in − x3( n) jn

(9)

Substitution of the SaS approximation (8) in Eq. (3), results in βi( n ) = ∑ M ( n )in ui( n) in ,

(10)

in

where M ( n) jn = L(,3n) jn are the derivatives of the Lagrange basis polynomials, which are calculated at SaS as follows:

(

)

M ( n) jn x3( n)in =

(

)

x3( n) jn

M ( n) in x3( n)in = −

1 − x3( n) in

∏

k n ≠ in , jn

∑ M (n ) j (x3( n)i ). n

j n ≠ in

n

x3( n) in − x3( n ) k n for jn ≠ in , x3( n ) jn − x3( n) k n (11)

7 Using Eq. (10) and the identity M ( n) in = ∑ L( n ) jn M ( n) in ( x3( n) jn ),

(12)

jn

which follows directly from the fundamental property of the Lagrange basis polynomials

L( n) jn , one obtains βi( n) = ∑ L( n) in βi( n) in ,

(13)

in

where βi( n)in = ∑ M ( n ) jn ( x3( n )in )ui( n) jn .

(14)

jn

It is seen from Eq. (14) that the key functions βi( n )in of the SaS laminated plate formulation are represented as a linear combination of displacements of SaS of the nth layer ui( n) jn .

Remark 1. The functions βi( n)1, βi( n) 2 ,..., βi( n) I n are linearly dependent, that is, there exist

numbers α ( n)1 , α ( n) 2 ,..., α ( n) I n , which are not all zero such that

∑α (n)i βi( n)i n

n

= 0.

(15)

in

The proof of this statement can be found in paper [50] in which the weak SaS formulation for laminated plates is developed. Substituting the displacement approximations (8) and (13) into the strain-displacement relationship (2) and using Eq. (6), we arrive at the distribution of strains through the thickness of the n th layer εij( n) = ∑ L( n) in εij( n) in . in

(16)

8 As can be seen, the through-thickness strain distribution (16) is similar to the SaS displacement distributions (8) and (13). This is important because the strains of the nth layer explicitly depend on the SaS strains of the same layer.

3. Strong SaS formulation for laminated plate

For simplicity, we consider the case of linear elastic materials, which are described by (n ) ( n) σ ij( n) = Cijkl εkl ,

(17)

(n ) where Cijkl are the elastic constants of the nth layer. Here and in the following developments,

the summation on repeated Latin indices is implied. The equilibrium equations of the laminated plate can be written as σ ij( n, )j = 0,

(18)

where σ ij(n) are the stresses of the nth layer. The boundary conditions on the bottom and top surfaces are defined as ui(1) (−h / 2) = wi− or σ i(31) (− h / 2) = pi− ,

(19)

ui( N ) (h / 2) = wi+ or σ i(3N ) (h / 2) = pi+ ,

(20)

where wi− ( x1, x2 ) and wi+ ( x1, x2 ) are the prescribed displacements; pi− ( x1, x2 ) and pi+ ( x1, x2 ) are the external mechanical loads acting on the bottom and top surfaces. The continuity conditions at interfaces are ui( m) ( x3[ m] ) = ui( m +1) ( x3[ m ] ),

(21)

σ i(3m) ( x3[ m] ) = σi(3m +1) ( x3[ m] ),

(22)

where the index m identifies the belonging of any quantity to the interface Ω[m] and runs from 1 to N − 1 .

9 Following the SaS technique developed, we introduce stresses at SaS of the nth layer σ ij( n) in ( x1 , x2 ) by

σ ij( n)in = σ ij( n) ( x3( n)in ).

(23)

According to Eqs. (5), (17) and (23) the constitutive equations can be expressed in terms of SaS variables as ( n ) ( n ) in σ ij( n) in = Cijkl εkl .

(24)

The use of Eqs. (16), (17) and (24) yields the distribution of stresses through the thickness of the nth layer σ ij( n) = ∑ L( n )in σ ij( n) in .

(25)

in

Satisfying the equilibrium equations (18) at inner points x3( n ) mn inside the layers, the following differential equations are obtained: σ i(1n,1) ( x3( n) mn ) + σ i(2n,)2 ( x3( n ) mn ) + σ i(3n,)3 ( x3( n) mn ) = 0

(mn = 2, 3,..., I n − 1)

(26)

that can be written by using Eqs. (23) and (25) as follows: σ i(1n,1) mn + σ i(2n,)2mn + ∑ M ( n) in ( x3( n ) mn )σ i(3n) in = 0

( mn = 2, 3,..., I n − 1).

(27)

in

Next, we satisfy the boundary conditions (19) and (20) that results in

∑ L(1)i (−h / 2)ui(1)i 1

1

= wi− or

∑ L(1)i (−h / 2)σi(31)i

i1

1

1

= pi− ,

(28)

i1

∑ L( N )i

N

(h / 2)ui( N ) iN = wi+ or

iN

∑ L( N )i

N

(h / 2)σ i(3N ) i N = pi+ ,

(29)

iN

and the continuity conditions (21) and (22):

∑ L(m )i

m

( x3[ m] )ui( m )im =

im

∑ L(m )i im

∑ L( m +1)i

m +1

( x3[ m] )ui( m +1)im +1

( m = 1, 2,..., N − 1),

(30)

( x3[ m ] )σ i(3m +1) im +1

( m = 1, 2,..., N − 1).

(31)

i m +1

m

( x3[ m] )σ i(3m)im =

∑ L(m +1)i i m +1

m +1

10 Thus, the proposed formulation deals with 3( I1 + I 2 + ... + I N ) governing equations (27)(31) for finding the same number of SaS displacements ui( n )in . These differential and algebraic equations have to be solved to describe the response of the laminated composite plate with different boundary conditions. For this purpose, the differential quadrature method [51, 52] could be applied effectively. Here, however, we restrict ourselves to the 3D exact solution for the simply supported rectangular plate to assess the potential of the strong SaS formulation developed. The general boundary conditions will be discussed in our future developments.

4. Analytical solution for simply supported laminated rectangular plate

In this section, we study the simply supported laminated rectangular plate subjected to the sinusoidally distributed transverse load + p3+ = p30 sin r x1 sin s x2 , r = rπ / a, s = sπ / b,

(32)

where a and b are the length and width of the plate; r and s are the half-wave numbers in x1 and x2 directions. The boundary conditions are written as (n) σ11 = u 2( n) = u3( n) = 0 at x1 = 0 and x1 = a , (n) σ 22 = u1( n ) = u3( n) = 0 at x2 = 0 and x2 = b ,

(33)

To satisfy the boundary conditions (33), we seek the analytical solution in the following form: ( n )i n u1( n)in = u10 cos r x1 sin s x2 ,

( n )i n u2( n) in = u20 sin r x1 cos s x2 ,

( n )i n u3( n) in = u30 sin r x1 sin s x2 .

(34)

The use of Eqs. (6), (24) and (34) leads to ( n ) in ( n ) in ( n ) in ( n )i n ( n ) in ( n ) in ( n ) in ( n )i n ( ε11 , σ11 ) = ( ε110 , σ110 ) sin r x1 sin s x2 , (ε 22 , σ 22 ) = ( ε220 , σ 220 ) sin r x1 sin s x2 , ( n ) in ( n ) in ( n ) in ( n )i n ( n ) in ( n )i n ( n )i n ( n ) in ( ε12 , σ12 ) = ( ε120 , σ120 ) cos r x1 cos s x2 , ( ε13 , σ13 ) = ( ε130 , σ130 ) cos r x1 sin s x2 ,

11 ( n ) in ( n ) in ( n ) in ( n ) in ( n )i n ( n ) in ( n )i n ( n ) in ( ε23 , σ 23 ) = ( ε230 , σ 230 ) sin r x1 cos s x2 , ( ε33 , σ 33 ) = ( ε330 , σ 330 ) sin r x1 sin s x2 ,

(35) where ( n ) in ( n ) in ( n )i n ( n ) in ( n ) in ( n )i n ( n )i n ε110 = − r u10 , ε220 = − s u20 , 2ε120 = s u10 + r u 20 , ( n )i n ( n ) in ( n ) in ( n ) in ( n )i n ( n ) in ( n ) in ( n ) in 2ε130 = r u30 + β10 , 2ε230 = s u30 + β20 , ε330 = β30 ,

βi(0n)in = ∑ M ( n ) jn ( x3( n )in )ui(0n) jn .

(36)

jn

Considering a special case of the orthotropic material with (n ) ( n) ( n) ( n) (n ) (n) ( n) ( n) ( n) (n) (n ) (n ) C1112 = C1113 = C1123 = C2212 = C2213 = C2223 = C1213 = C1223 = C1233 = C1323 = C1333 = C2333 = 0,

one can write the constitutive equations in a compact form as ( n ) ( n ) in σ ij( n0)in = Cijkl ε kl 0 .

(37)

Substituting (34) and (35) in Eqs. (27)-(31), we arrive at the system of linear equations ( n)mn ( n ) mn ( n )i n r σ110 − s σ120 + ∑ M ( n )in ( x3( n) mn )σ130 =0

( mn = 2, 3,..., I n − 1),

in

( n ) mn (n )mn ( n ) in − r σ120 + s σ 220 + ∑ M ( n)in ( x3( n) mn )σ 230 =0

(mn = 2, 3,..., I n − 1),

in

( n ) mn ( n ) mn ( n ) in − r σ130 − s σ 230 + ∑ M ( n) in ( x3( n ) mn )σ 330 =0

( mn = 2, 3,..., I n − 1),

(38)

in

∑ L(1)i (−h / 2)σi(301)i 1

1

= 0,

(39)

i1

∑ L( N )i

N

N )i N (h / 2)σ i(30 = pi+0 ,

+ + p10 = p20 = 0,

(40)

iN

∑ L(m )i

m

( x3[ m] )ui(0m )im =

im

∑ L(m )i im

∑ L( m +1)i

m +1

( x3[ m] )ui(0m +1)im +1

(m = 1, 2,..., N − 1),

(41)

m +1) im +1 ( x3[ m] )σ i(30

( m = 1, 2,..., N − 1)

(42)

im +1

m

m )i m ( x3[ m] )σ i(30 =

∑ L( m +1)i im + 1

m +1

12 of order 3( I1 + I 2 + ... + I N ) . Therefore, the displacements of SaS of the nth layer ui(0n)in can be found readily with the Symbolic Math Toolbox, which incorporates symbolic computations into the numeric environment of MATLAB. For the solution of the system of algebraic equations (38)-(42) the function solve is utilized. This makes it possible to obtain the analytical solution for laminated orthotropic rectangular plates in the framework of the SaS formulation.

5. Numerical examples

5.1. Single-layer isotropic rectangular plate Consider a simply supported single-layer isotropic rectangular plate subjected to a sinusoidally distributed transverse load p3+ = p0 sin

πx1 πx sin 2 . a b

(43)

There is a closed-form solution of the problem obtained by Vlasov [1] as early as 1957. The Vlasov's solution is presented in Appendix A. To analyze the results efficiently, we introduce the dimensionless variables at crucial points as functions of the thickness coordinate: u1 = 100 Eh 2u1( 0, b / 2, z ) / a 3 p0 , u3 = 100 Eh3u3 ( a / 2, b / 2, z ) / a 4 p0 ,

σ11 = 10h 2σ11 ( a / 2, b / 2, z ) / a 2 p0 , σ12 = 10h 2σ12 (0, 0, z ) / a 2 p0 , σ13 = 10hσ13 (0, b / 2, z ) / ap0 , σ 33 = σ 33 ( a / 2, b / 2, z ) / p0 , z = x3 / h,

(44)

where p0 = 1 . The material properties and plate dimensions are E = 107 , ν = 0.3 and

a = b = 1. Table 1 lists the results of the convergence study due to increasing the number of SaS. A comparison with Vlasov's exact solution [1] is also given. As it turned out, the proposed strong SaS formulation provides from 15 to 16 right digits for displacements and stresses

13 choosing 15 SaS for the thick plate with a / h = 4 . However, the use of the more number of SaS does not yield a better accuracy. This is due to the fact that the authors restrict themselves to default MATLAB precision. Note also that we did not discover the divergence of the developed algorithm of symbolic computations even taking 55 SaS. This is illustrated in Fig. 2 by using the theoretical and computational errors log( δ3 ) and log( ∆3 ) with δ3 = 100 Eh3δ3 / a 4 p0 ,

∆3 = u3 (0.5) − u3 exact ( 0.5) ,

(45)

where δ3 is the theoretical error bound defined in Appendix B. The computational error is evaluated additionally with 64-digit precision.

5.2. Sandwich rectangular plate Here, we study a simply supported sandwich rectangular plate with different boundary conditions on the top surface, namely, transverse loading (43) and imposed transverse deformation defined according to Eq. (29) as w3+ = w0 sin

πx1 πx sin 2 . a b

(46)

The face sheets are made of the graphite-epoxy composite with the material properties E L = 25ET , GLT = 0.5ET , GTT = 0.2 ET , ET = 106 , ν LT = ν TT = 0.25 . The fiber directions

coincide with x1 -direction in both face sheets. The material properties of the core are taken to be E1 = E2 = 4 × 10 4 , E3 = 5 × 105 , G13 = G23 = 6 × 104 , G12 = 1.6 × 10 4 , ν12 = ν31 = ν32 = 0.25 . The thicknesses of face sheets and a core are h1 = h3 = 0.1h and h2 = 0.8h . In the case of transverse loading (43), we utilize the Pagano's dimensionless variables as functions of the thickness coordinate: u1 = 10 ET h 2u1 (0, b / 2, z ) / a 3 p0 , u3 = 100ET h3u3 ( a / 2, b / 2, z ) / a 4 p0 , σ11 = h 2σ11 ( a / 2, b / 2, z ) / a 2 p0 , σ12 = 10h 2 σ12 (0, 0, z ) / a 2 p0 ,

14

σ13 = 10hσ13 (0, b / 2, z ) / ap0 , σ 23 = 10hσ 23 (a / 2, 0, z ) / ap0 , σ 33 = σ33 (a / 2, b / 2, z ) / p0 , z = x3 / h ,

(47)

where p0 = 1 and a = b = 1 , whereas the case of imposed transverse deformation (46) is characterized by the following dimensionless variables:

u1 = 100u1 (0, b / 2, z ) / w0 , u3 = u3 (a / 2, b / 2, z ) / w0 ,

σ11 = 10aσ11 (a / 2, b / 2, z ) / ET w0 , σ12 = 10aσ12 (0, 0, z ) / ET w0 ,

σ13 = 10a 2σ13 (0, b / 2, z ) / hET w0 , σ 23 = 10a 2σ 23 (a / 2, 0, z ) / hET w0 , σ 33 = 10a 3σ 33 ( a / 2, b / 2, z ) / h 2 ET w0 , z = x3 / h,

(48)

where w0 = 1 . Tables 2-5 list the results of the convergence study due to increasing the number of SaS I n inside each layer for two values of the slenderness ratio a/h. The derived analytical solution for the sandwich plate under transverse loading is compared with Pagano's exact solution [3]. As can be seen, the algorithm of symbolic computations provides already 15 right digits for all basic variables since 13 and 11 SaS for thick and moderately thick plates, respectively, through default MATLAB precision. Figs. 3 and 4 show the through-thickness distributions of displacements and stresses for different values of the slenderness ratio a / h taking seven SaS for each layer. These results demonstrate convincingly the high potential of the strong SaS formulation because the boundary conditions on bottom and top surfaces and the continuity conditions at interfaces for the transverse stresses are satisfied exactly. Note that in a variational SaS formulation [15, 16] the boundary conditions on outer surfaces and the continuity conditions at interfaces are satisfied approximately.

5.3. Three-layer cross-ply rectangular plate Consider next a symmetric three-layer rectangular plate composed of the graphite-epoxy

15 composite with the stacking sequence [0/90/0], that is, the fiber directions coincide with x1 direction in outer layers and x2 -direction in a central layer. The material properties of the graphite-epoxy are given in a previous section. The plate dimensions are h1 = h2 = h3 = h / 3 , a = 1 and b = 3 . The data listed in Tables 6 and 7 show the results of the convergence study due to increasing the number of SaS for the moderately thick plate with a / h = 10 for both types of boundary conditions (43) and (46). A comparison with Pagano's exact solution [3] is also given. Figs. 5 and 6 present the distributions of displacements and stresses in the thickness direction for different values of the slenderness ratio by choosing seven SaS for each layer. These results demonstrate again the high potential of the strong SaS formulation for laminated composite plates.

6. Conclusions

An efficient strong SaS formulation based on a new concept of SaS located at Chebyshev polynomial nodes throughout the layers and direct integration of the equilibrium equations of elasticity for the 3D stress analysis of laminated plates has been proposed. The use of only Chebyshev polynomial nodes makes it possible to minimize uniformly the error due to Lagrange interpolation. Therefore, the developed strong SaS formulation gives an opportunity to obtain the analytical solutions for the laminated composite plates with a prescribed accuracy, which asymptotically approach the exact solutions of elasticity as the number of SaS tends to infinity.

Acknowledgements

This work was supported by the Russian Science Foundation (Grant No. 15-19-30002) and the Russian Ministry of Education and Science (Grant No. 9.4914.2017).

16 Appendix A. Vlasov's solution for homogeneous rectangular plate

The Vlasov's closed-form solution [1] for a simply supported isotropic rectangular plate subjected to transverse loading (43) can be expressed as u1 = U ( x3 ) cos

πx1 πx πx πx πx πx sin 2 , u 2 = V ( x3 ) sin 1 cos 2 , u3 = W ( x3 ) sin 1 sin 2 , a b a b a b

(A1)

U = ( A1 + A2 x3 ) exp( λx3 ) + ( A3 + A4 x3 ) exp(− λx3 ), V = ( B1 + B2 x3 ) exp( λx3 ) + ( B3 + B4 x3 ) exp( − λx3 ), W = (C1 + C2 x3 ) exp( λx3 ) + (C3 + C4 x3 ) exp( − λx3 ),

(A2)

where λ=

πc , ab

A2 =

c = a 2 + b2 ,

b C2 , c

B2 =

a C2 , c

b a A4 = − C4 , B4 = − C4 , c c

C1 =

b a (3 − 4ν) ab A1 + B1 − C2 , c c πc

b a (3 − 4ν)ab C3 = − A3 − B3 + C4 . c c πc

(A3)

Six coefficients A1, A3, B1, B3, C2 and C4 are found from boundary conditions σ13 (−h / 2) = σ 23 (−h / 2) = σ 33 (−h / 2) = 0, σ13 ( h / 2) = σ 23 ( h / 2) = 0,

σ 33 ( h / 2) = p3+ .

(A4)

Appendix B. Error estimate of Lagrange interpolation

Let U L −1 ( x3 ) , VL −1 ( x3 ) and WL −1 ( x3 ) be the Lagrange polynomials of degree L − 1 , which interpolate, respectively, functions U ( x3 ) , V ( x3 ) and W ( x3 ) from Eq. (A2) at Chebyshev nodes tl ∈ (− h / 2, h / 2) , where l = 1, 2,..., L . These nodes are the roots of the Chebyshev polynomial of the first kind of degree L given by 1  2l − 1  tl = − h cos π . 2  2L 

(B1)

17 The Chebyshev nodes are important because their choice makes possible to minimize uniformly the error of the Lagrange interpolation due to a remainder estimate [47]: U − U L −1 ≤

where

f

f =

hL

2

2 L −1

L!

U (L) ,

(B2)

is the norm defined as max

t∈[− h / 2, h / 2 ]

f (t ) .

(B3)

Introduce next functions F1 = exp( λx3 ), F2 = x3 exp( λx3 ), F3 = exp(− λx3 ), F4 = x3 exp(− λx3 ).

(B4)

The norms of Lth order derivatives of these functions are F1( L ) = F3( L ) = λL exp( λh / 2),

F2( L ) = F4( L ) = λL −1 ( L + λh / 2) exp( λh / 2).

(B5)

The use of relations (A2), (B2), (B4) and (B5) yields an improved remainder estimate: U − U L −1 ≤

δ1 =

hL 2

2 L −1

4

L!

∑ Aq Fq( L)

≤ δ1,

q =1

h L λL −1 exp( λh / 2) (λ( A1 + A3 ) + ( L + λh / 2)( A2 + A4 )). 22 L −1 L!

(B6)

The similar error estimates can be obtained for the solutions V ( x3 ) and W ( x3 ) . Here, for conciseness we present the latter estimate: W − WL −1 ≤ δ3 ,

δ3 =

h L λL −1 exp( λh / 2) (λ( C1 + C3 ) + ( L + λh / 2)( C2 + C4 )). 22 L −1 L!

(B7)

18 References

[1] Vlasov BF. On the bending of rectangular thick plate. Vestnik Moskovskogo Universiteta. Matematika, Mekhanika, Astronomia 1957;2:25–31 (in Russian). [2] Pagano NJ. Exact solutions for composite laminates in cylindrical bending. J Compos Mater 1969;3:398–411. [3] Pagano NJ. Exact solutions for rectangular bidirectional composites and sandwich plates. J Compos Mater 1970;4:20–34. [4] Srinivas S, Rao AK, Rao CVJ. Flexure of simply supported thick homogeneous and laminated rectangular plates. ZAMM 1969;49:449–458. [5] Srinivas S, Rao AK. Bending, vibration and buckling of simply supported thick orthotropic rectangular plates and laminates. Int J Solids Struct 1970;6:1463–1481. [6] Murakami H. Assessment of plate theories for treating the thermomechanical response of layered plates. Compos Eng 1993;3:137–149. [7] Tungikar VB, Rao KM. Three dimensional exact solution of thermal stresses in rectangular composite laminate. Compos Struct 1994;27:419–430. [8] Bhaskar K, Varadan TK, Ali JSM. Thermoelastic solutions for orthotropic and anisotropic composite laminates. Compos Part B 1996;27:415–420. [9] Pagano NJ. Influence of shear coupling in cylindrical bending of anisotropic laminates. J Compos Mater 1970;4:330–343. [10] Vel SS, Batra RC. Generalized plane strain thermoelastic deformation of laminated anisotropic thick plates. Int J Solids Struct 2001;38:1395–1414. [11] Noor AK, Burton WS. Three-dimensional solutions for antisymmetrically laminated anisotropic plates. J Appl Mech 1990;57:182–188. [12] Savoia M, Reddy JN. A variational approach to three-dimensional elasticity solutions of laminated composite plates. J Appl Mech 1992;59:S166–S175.

19 [13] Savoia M, Reddy JN. Three-dimensional thermal analysis of laminated composite plates. Int J Solids Struct 1995;32:593–608. [14] Loredo A. Exact 3D solution for static and damped harmonic response of simply supported general laminates. Compos Struct 2014;108:625–634. [15] Kulikov GM, Plotnikova SV. Exact 3D stress analysis of laminated composite plates by sampling surfaces method. Compos Struct 2012;94:3654–3663. [16] Kulikov GM, Plotnikova SV. Three-dimensional thermal stress analysis of laminated composite plates with general layups by a sampling surfaces method. Eur J Mech A/Solids 2015;49:214–226. [17] Reddy JN. Mechanics of laminated composite plates and shells: Theory and analysis. 2nd ed. Boca Raton: CRC Press; 2004. [18] Jha DK, Kant T, Singh RK. A critical review of recent research on functionally graded plates. Compos Struct 2013;96:833–849. [19] Zhong Z, Shang ET. Three-dimensional exact analysis of a simply supported functionally gradient piezoelectric plate. Int J Solids Struct 2003;40:5335–5352. [20] Pan E, Han F. Exact solution for functionally graded and layered magneto-electroelastic plates. Int J Eng Sci 2005;43:321–339. [21] Alibeigloo A. Exact solution for thermo-elastic response of functionally graded rectangular plates. Compos Struct 2010;92:113–121. [22] Vel SS, Batra RC. Exact solution for thermoelastic deformations of functionally graded thick rectangular plates. AIAA J 2002;40:1421–1433. [23] Kulikov GM, Plotnikova SV. A new approach to three-dimensional exact solutions for functionally graded piezoelectric laminated plates. Compos Struct 2013;106:33–46. [24] Darabi M, Ganesan R. Exact 3-D stress and stiffness analysis of functionally graded sandwich plates using sampling surfaces method. In: Proceedings of the ASME

20 international mechanical engineering congress and exposition: mechanics of solids, structures and fluids, vol. 9. p. 1–11. V009T12A065. [25] Kulikov GM, Plotnikova SV. Three-dimensional exact analysis of functionally graded laminated composite plates. Adv Struct Mater 2015;45:223–241. [26] Kulikov GM, Plotnikova SV. A sampling surfaces method and its implementation for 3D thermal stress analysis of functionally graded plates. Compos Struct 2015;120:315– 325. [27] Kashtalyan M. Three-dimensional elasticity solution for bending of functionally graded rectangular plates. Eur J Mech A/Solids 2004;23:853–864. [28] Kashtalyan M, Menshykova M. Three-dimensional elasticity solution for sandwich panels with a functionally graded core. Compos Struct 2009;87:36–43. [29] Woodward B, Kashtalyan M. Three-dimensional elasticity solution for bending of transversely isotropic functionally graded plates. Eur J Mech A/Solids 2011;30:705– 718. [30] Plevako VP. On the theory of elasticity of inhomogeneous media. J Appl Math Mech 1971;35:806–813. [31] Wang YM, Tarn JQ. A three-dimensional analysis of anisotropic inhomogeneous and laminated plates. Int J Solids Struct 1994;31:497–515. [32] Tarn JQ, Yen CB. A three-dimensional asymptotic analysis of anisotropic inhomogeneous and laminated cylindrical shells. Acta Mech 1995;113:137–153. [33] Wu CP, Tarn JQ, Chi SM. Three-dimensional analysis of doubly curved laminated shells. J Eng Mech 1996;122:391–401. [34] Cheng ZQ, Batra RC. Three-dimensional asymptotic analysis of multiple-electroded piezoelectric laminates. AIAA J 2000;38:317–324.

21 [35] Reddy JN, Cheng ZQ. Three-dimensional thermomechanical deformations of functionally graded rectangular plates. Eur J Mech A/Solids 2001;20:841–855. [36] Vetyukov Yu, Kuzin A, Krommer M. Asymptotic splitting in the three-dimensional problem of elasticity for non-homogeneous piezoelectric plates. Int J Solids Struct 2011;48:12–23. [37] Wu CP, Lo JY. Three-dimensional elasticity solutions of laminated annular spherical shells. J Eng Mech 2000;126:882–885. [38] Wu CP, Hung YC, Lo JY. A refined asymptotic theory of laminated circular conical shells. Eur J Mech A/Solids 2002;21:281–300. [39] Wu CP, Jiang RY. An asymptotic meshless method for sandwich functionally graded circular hollow cylinders with various boundary conditions. J Sandwich Struct Mater 2015;17:469–510. [40] Tarn JQ, Wang YB, Wang YM. Three-dimensional asymptotic finite element method for anisotropic inhomogeneous and laminated plates. Int J Solids Struct 1996;33:1939– 1960. [41] Tarn JQ, Wang YB. A refined asymptotic theory and computational model for multilayered composite plates. Comput Methods Appl Mech Eng 1997;145:167–184. [42] Wu CP, Chiu KH, Wang YM. A review on the three-dimensional analytical approaches of multilayered and functionally graded piezoelectric plates and shells. Comput Mater Continua 2008;8:93–132. [43] Wu CP, Liu YC. A review of semi-analytical numerical methods for laminated composite and multilayered functionally graded elastic/piezoelectric plates and shells. Compos Struct 2016;147:1–15. [44] Kulikov GM. Refined global approximation theory of multilayered plates and shells. J Eng Mech 2001;127:119–125.

22 [45] Kulikov GM, Carrera E. Finite deformation higher-order shell models and rigid-body motions. Int J Solids Struct 2008;45:3153–3172. [46] Kulikov GM, Plotnikova SV. On the use of a new concept of sampling surfaces in a shell theory. Adv Struct Mater 2011;15:715–726. [47] Bakhvalov NS. Numerical methods: Analysis, algebra, ordinary differential equations. Moscow: MIR Publishers; 1977. [48] Kulikov GM, Plotnikova SV. On the use of sampling surfaces method for solution of 3D elasticity problems for thick shells. ZAMM 2012;92:910–920. [49] Kulikov GM, Plotnikova SV. Solution of three-dimensional problems for thick elastic shells by the method of reference surfaces. Mech Solids 2014;49:403–412. [50] Kulikov GM, Plotnikova SV. Hybrid-mixed ANS finite elements for stress analysis of laminated composite structures: Sampling surfaces plate formulation. Comput Methods Appl Mech Eng 2016;303:374–399. [51] Shu C. Differential quadrature and its application in engineering. Berlin: Springer; 2000. [52] Tornabene F, Fantuzzi N, Ubertini F, Viola E. Strong formulation finite element method based on differential quadrature: A survey. Appl Mech Rev 2015;67. p. 1–55. 020801.

23

Fig. 1. Geometry of the laminated plate.

24

Fig. 2. Theoretical log(δ3 ) and computational log( ∆3 ) errors of the Lagrange interpolation vs. the number of SaS I1 for a single-layer square plate: 16-digit precision (
) and 64-digit precision (□).

25

Fig. 3. Through-thickness distributions of transverse shear stresses for a sandwich square plate under transverse loading for I1 = I 2 = I 3 = 7 : SaS formulation (––) and Pagano's solution [3] (
).

26

Fig. 4. Through-thickness distributions of displacements and stresses for a sandwich square plate under imposed transverse deformation for I1 = I 2 = I 3 = 7 .

27

Fig. 5. Through-thickness distributions of transverse shear stresses for a symmetric three-ply rectangular plate under transverse loading for I1 = I 2 = I 3 = 7 : SaS formulation (––) and Pagano's solution [3] (
).

28

Fig. 6. Through-thickness distributions of displacements and stresses for a symmetric threeply rectangular plate under imposed transverse deformation for I1 = I 2 = I 3 = 7 .

Table 1 Results of the convergence study for a single-layer isotropic plate with a / h = 4 .

I1

u1 (0.5)

u3 (0.5)

σ11 (0.5)

σ12 (0.5)

σ13 (0)

σ 33 (0)

5

-4.325773721611378

3.668142517752232

2.209259849272290

-1.045370688069702

2.374037736473570

0.4969261935948460

7

-4.249048061124980

3.623967815103729

2.174825453368604

-1.026829090275407

2.361932373774121

0.4981277788146725

9

-4.248634179691244

3.623728522595966

2.174639703815501

-1.026729071285270

2.361901499297424

0.4981386504489832

11

-4.248633189690484

3.623727947100229

2.174639259504200

-1.026728832040723

2.361901343377887

0.4981386741941986

13

-4.248633188381382

3.623727946335935

2.174639258916675

-1.026728831724364

2.361901343235282

0.4981386742328633

15

-4.248633188380300

3.623727946335301

2.174639258916189

-1.026728831724102

2.361901343235109

0.4981386742328927

17

-4.248633188380298

3.623727946335302

2.174639258916189

-1.026728831724102

2.361901343235110

0.4981386742328930

-4.248633188380297

3.623727946335301

2.174639258916189

-1.026728831724102

2.361901343235109

0.4981386742328930

a

Vlasov a

The exact results have been obtained by using Eqs. (A1)-(A3)

Table 2 Results of the convergence study for a sandwich plate with a / h = 4 under transverse loading.

In

u1 (0.5)

u 3 ( 0)

σ11 (0.5)

σ12 (0.5)

σ13 (0)

σ 23 (0)

3

-2.012411861508411

7.925643598036384

1.665145792457571

-1.534250874975219

2.475763108985881

1.054447221538144

5

-1.878862849186595

7.597131275954131

1.556165159317271

-1.436892394047969

2.386981422359567

1.071914037778667

7

-1.878477083519309

7.596236363941784

1.555849800904065

-1.436599990706936

2.386681307396541

1.071914921686850

9

-1.878476690136775

7.596235469230204

1.555849478729198

-1.436599680764365

2.386681061835270

1.071914999074346

11

-1.878476689941961

7.596235468793107

1.555849478569353

-1.436599680604989

2.386681061691722

1.071914999095251

13

-1.878476689941904

7.596235468792981

1.555849478569307

-1.436599680604942

2.386681061691688

1.071914999095267

15

-1.878476689941904

7.596235468792979

1.555849478569306

-1.436599680604941

2.386681061691688

1.071914999095267

1.556

-1.437

2.39

1.072

Pagano

Table 3 Results of the convergence study for a sandwich plate with a / h = 10 under transverse loading.

In

u1 (0.5)

u3 ( 0)

σ11 (0.5)

σ12 (0.5)

σ13 (0)

σ 23 (0)

3

-1.461521472302024

2.237672111578537

1.178581690589604

-0.7222543810461195

3.033238359574364

0.5023250263118528

5

-1.429879352761541

2.200441783659711

1.153135468244837

-0.7066651642153221

2.997863094070032

0.5269276974680131

7

-1.429862290670034

2.200422528621556

1.153121744401966

-0.7066567041496520

2.997840370444016

0.5269358928828287

9

-1.429862287655427

2.200422525294808

1.153121741976948

-0.7066567026503436

2.997840367265042

0.5269358952491626

11

-1.429862287655172

2.200422525294530

1.153121741976742

-0.7066567026502162

2.997840367264738

0.5269358952493245

13

-1.429862287655172

2.200422525294531

1.153121741976742

-0.7066567026502160

2.997840367264738

0.5269358952493254

1.153

-0.707

3.00

0.527

Pagano

Table 4 Results of the convergence study for a sandwich plate with a / h = 4 under imposed transverse deformation.

In

u1 (0.5)

u 3 ( 0)

σ11 (0.5)

σ13 (0)

σ 23 (0)

σ 33 (0.5)

3

-6.179780768982413

0.9735331186577605

51.13384661760139

7.602655073797430

3.238031332554400

3.070833007489000

5

-6.010787913556095

0.9721783550294466

49.78425506242155

7.636341891380344

3.429259380651248

3.199162682980458

7

-6.010231434914991

0.9721734479948576

49.77977886150624

7.636242754670397

3.429616903181413

3.199523426527451

9

-6.010230848889353

0.9721734422857414

49.77977412435065

7.636242823570991

3.429617534596478

3.199523784589383

11

-6.010230848592807

0.9721734422826567

49.77977412194280

7.636242823526876

3.429617534849826

3.199523784763336

13

-6.010230848592719

0.9721734422826557

49.77977412194206

7.636242823526886

3.429617534849931

3.199523784763386

15

-6.010230848592720

0.9721734422826557

49.77977412194208

7.636242823526886

3.429617534849932

3.199523784763387

Table 5 Results of the convergence study for a sandwich plate with a / h = 10 under imposed transverse deformation.

In

u1 (0.5)

u 3 ( 0)

σ11 (0.5)

σ13 (0)

σ 23 (0)

σ 33 (0.5)

3

-6.523431649060901

0.9987743149541511

52.60543376954998

13.53871529702984

2.242103887530790

4.463452486117721

5

-6.489776958651818

0.9987119793617310

52.33722675669363

13.60636671584480

2.391560674891188

4.538688488730214

7

-6.489756075144884

0.9987119435337858

52.33705997384936

13.60638215536610

2.391618713467690

4.538728042197533

9

-6.489756071231461

0.9987119435272305

52.33705994256759

13.60638216141934

2.391618727807949

4.538728049029698

11

-6.489756071231114

0.9987119435272299

52.33705994256481

13.60638216141967

2.391618727808986

4.538728049030269

13

-6.489756071231113

0.9987119435272299

52.33705994256482

13.60638216141967

2.391618727808988

4.538728049030267

Table 6 Results of the convergence study for a symmetric three-ply plate with a / h = 10 under transverse loading.

In

u1 (0.5)

u3 ( 0)

σ11 (0.5)

σ12 (0.5)

3

-1.0662550311327100

1.0313639907775190

0.8440235265935225

-0.1370524596048226 4.729378164244301

0.4998653788171814

5

-0.9172023648951146

0.9194154336410869

0.7265073857024849

-0.1198484666720689 4.203965814463195

0.5002525718599943

7

-0.9165308368312481

0.9189063085004530

0.7259779581682628

-0.1197720821035059 4.201372555647402

0.5002521533924405

9

-0.9165297276388366

0.9189054637672304

0.7259770836890195

-0.1197719559046020 4.201368155991792

0.5002521523065288

11

-0.9165297267579845

0.9189054630948098

0.7259770829945604

-0.1197719558042799 4.201368152451253

0.5002521523056154

13

-0.9165297267575746

0.9189054630944964

0.7259770829942374

-0.1197719558042332 4.201368152449593

0.5002521523056145

15

-0.9165297267575749

0.9189054630944966

0.7259770829942376

-0.1197719558042332 4.201368152449592

0.5002521523056154

0.919

0.726

-0.120

Pagano

σ13 (0)

4.20

σ 33 (0)

Table 7 Results of the convergence study for a symmetric three-ply plate with a / h = 10 under imposed transverse deformation.

In

u1 (0.5)

u 3 ( 0)

σ11 (0.5)

σ12 (0.5)

σ13 (0)

σ 33 (0.5)

3

-10.32553460253789

0.9987652356745952

81.73461202747450

-1.327205849157104

45.79894721013859

9.683925797347765

5

-9.957305612327193

0.9981331064820719

78.87088330697611

-1.301095434876105

45.63897128980135

10.85617088816147

7

-9.955495990703302

0.9981298721750038

78.85681923011865

-1.300982394986435

45.63594147885691

10.86215061254541

9

-9.955493040926806

0.9981298668153145

78.85679631099092

-1.300982213177405

45.63593539633644

10.86216053959771

11

-9.955493038601379

0.9981298668110505

78.85679629292525

-1.300982213034143

45.63593539107826

10.86216054749983

13

-9.955493038600306

0.9981298668110485

78.85679629291690

-1.300982213034077

45.63593539107570

10.86216054750351

15

-9.955493038600302

0.9981298668110485

78.85679629291691

-1.300982213034077

45.63593539107568

10.86216054750351