Enzo Cosentino
email: enzo.cosentino@airbus.com

Paul Weaver

An Enhanced Single-layer Variational Formulation for the Effect of Transverse Shear on Laminated Orthotropic Plates

Keywords: Castigliano's Theorem, Variational Approach, Mixed Formulation

A novel mixed formulation is derived by means of Reissner's variational approach based on Castigliano's principle of least work in conjunction with a Lagrange multiplier method for the calculus of variations. The governing equations present an alternative theory for modeling the important three-dimensional structural aspects of plates in a two-dimensional form. By integrating the classical Cauchy's equilibrium equations with respect to the thickness co-ordinate, and enforcing continuity of shear and normal stresses at each ply interface, condenses the effect of the thickness. A reduced system of partial differential equations of sixth order in one variable, is also proposed, which contains differential correction factors that formally modify the classical constitutive equations for composite laminates. The theory degenerates to classical composite plate analysis for thin configurations. Significant deviations from classical plate theory are observed when the thickness becomes comparable with the in-plane dimensions. A variety of case studies are presented and solutions are compared with other models available in the literature and with finite element analysis.

List of Symbols

Introduction

It is well known from experimental observations [START_REF] Reddy | Energy and Variational Methods in Applied Mechanics[END_REF] that the Classical Theory of Plates (CTP), i.e. the two-dimensional model chiefly shaped by the works undertaken by [START_REF] Kirchhoff | Uber das Gleichgewicht und die Bewegung einer elestischen Scheibe[END_REF] in the 19 th century, and then revisited by [START_REF] Love | The mathematical Theory of Elasticity[END_REF] and [START_REF] Timoshenko | Theory of Elasticity[END_REF] during the early 20 th century, underpredicts global deformations of the elastic equilibrium. The principle hypothesis on which Kirchhoff's theory is based, is that normals to the mid-plane before deformation remain straight and normal to the plane after deformation. This causes the main deficiency of the approachneglecting the effect of the transverse shear strains on the deformation of elastic two-dimensional bodies. Ignoring some deformation modes constrains the deformed, modeled structure to one degree of freedom only (namely the transverse displacement w), thus yielding overly stiff behaviour. Note, we refer to degrees of freedom as Lagrange parameters in the present work. It is also noteworthy that neglecting shear stresses leads to a contraction of the three natural (force) boundary conditions that should be satisfied along a free edge (1945), i.e. vanishing of vertical forces, bending and twisting moments.

Despite its limitations, Kirchhoff's theory is still the most widespread approach used to obtain robust and reliable prediction of the behaviour of slender plated structures. Furthermore, it underpins the foundation for composite plate analysis, the well-known Classical Laminate Analysis (CLA). The main simplification is that three-dimensional structures are treated as two-dimensional coincident with the neutral plane (gray area in Fig. 1), with a significant reduction of the total number of variables and computational effort.

Furthermore, the governing equations become easier to solve and closed-form solutions, which typically provide more intuitive and physical representation, can be developed. Nonetheless, the error induced by neglecting the effect of transverse shear stresses becomes significant in the following cases:

M A N U S C R I P T A C C E P T E D

i. Thick-sectioned configurations, for which the hypothesis of two-dimensionality becomes restrictive ii. Proximity of free edges (especially in composite plates)

iii.

Proximity of free corners

The induced error becomes more important for plates and shells made of composite layers, for which the ratio of longitudinal to transverse shear elastic modulii is relatively large compared to isotropic materials [START_REF] Reddy | Energy and Variational Methods in Applied Mechanics[END_REF]. Composite plates are, in fact, inherently subjected to transverse shear and normal stresses due to their discontinuous through-the-thickness behaviour and their global anisotropic nature. The CLA is able to accurately describe the global transverse behaviour of composite plates only within the inner domain, where the assumptions proposed by Kirchhoff are often formally recovered. However, it does not contain explanations for stress distortions and singularities that take place in thick laminates and in the proximity of free edges. Lekhnitskii proposed a first attempt to model anisotropic behaviour of composite beams [START_REF] Lekhnitskii | Strength calculation of composite beams[END_REF]. The same author proposed a series of closed form solutions and solution strategies for anisotropic plates [START_REF] Lekhnitskii | Anisotropic Plates[END_REF]. during the first decades of the 20 th century. Since then, numerous researchers have proposed a variety of approaches to model the effect of transverse shears in isotropic and composite materials. [START_REF] Timoshenko | Theory of Elasticity[END_REF] Bernoulli in the 18 th century. [START_REF] Reissner | On the theory of bending of elastic plates[END_REF] and [START_REF] Mindlin | Influence of rotatory inertia and shear on flexural motion of isotropic elastic plates[END_REF], extended the concept to the analysis of two-dimensional structures, originating the so-called First Order Shear Deformation Theory (FSDT). Other distinguished researchers, such as Vlasov, Thomas, Reissner [START_REF] Carrera | The Best on Plate/Shell Theories for Laminated Structures Analysis[END_REF] and [START_REF] Reddy | A Layerwise Shell Theory with Application to Buckling and Vibration of Cross-ply Laminated Circular Cylindrical Shells[END_REF] introduced higher order expansions of the axial and shear deformation through the plate thickness. These models provide reasonably acceptable results if applied to isotropic, orthotropic, and anisotropic single layered laminates. They are not suitable to model the behavior of thick laminated composite plates where accurate descriptions of transverse stresses are required, especially in the proximity of free edges, or other special features such as holes, cut-outs, corners and notches, where the interlaminar shears triggered by the stress field distortion are highly dependent on the real stacking sequence. The so-called free edge effect is, in fact, mainly explained by the mismatch of elastic properties, in particular the Poisson's ratio, at the interface between two adjacent plies. Hence, more accurate layer-wise field equations are required to capture such phenomena. The theories that have been developed, to account for through-the-thickness piece-wise behaviour of stresses and displacement are often referred to as Zig-Zag theories (ZZ) in the open literature. The first attempt to define a self-contained ZZ theory was possibly proposed by [START_REF] Lekhnitskii | Anisotropic Plates[END_REF]. The main limitation of the approach is that it models multi-layered beams where each layer is isotropic. [START_REF] Carrera | Historical review of Zig-Zag theories for multilayered plates and shells[END_REF] suggested that [START_REF] Ren | Bending theory of laminated plates[END_REF] later improved the work by Lekhnitskii. (1958-A and1958-B) proposed an extension of the theory developed by [START_REF] Reissner | On the theory of bending of elastic plates[END_REF] to multi-layered plates. As a general consideration, this approach is based on the analytical or numerical calculation of five Lagrange parameters, which represent the degrees of freedom of the structure in a functional sense. Indeed each degree of freedom is represented by a function that is calculated over the entire domain. Typically, the DoFs are: the two axial in-plane displacements along x and y co-ordinates (namely u x and u y ); the transverse displacement w, and the two transverse shear strains γ xz and γ yz . In uncoupled problems, i.e. when it is possible to separate the in-plane from the transverse response, the maximum number of degrees of freedom that must be determined in order to solve the transverse equilibrium problem is three. The shear strains are sometimes replaced by the two rotations of the cross section about the y and x-axes, ϕ x and ϕ y respectively. Another important class of models is symbolized by those theories that are derived by means of Reissner's variational approach. The so-called mixed variational approach-based on the variational

M A N U S C R I P T A C C E P T E D

principles developed by [START_REF] Hellinger | Die Algemeinen der Mechanik der Kontinua[END_REF] was proposed by [START_REF] Reissner | On the theory of bending of elastic plates[END_REF] and then improved by the same author [START_REF] Reissner | The effect of transverse shear deformation on the bending of elastic plates[END_REF] [START_REF] Carrera | A class of two dimensional theories for multilayered plates analysis[END_REF], the layer-wise approaches based on Reissner's principle require a large computational efforts. However, they do provide an extremely accurate description of the transverse shear stresses. The number of variables that must be computed is at least 2N+1, N being the total number of layers. The number of variables can be significantly reduced by means of a weak form of Hooke's Law [START_REF] Carrera | A class of two dimensional theories for multilayered plates analysis[END_REF], which expresses the aforementioned variables in terms of the three displacements only. Yet, the number of degrees of freedom needed is at least equal to three. Recently, other approaches have been proposed by Mittelstedt and Becker (2003[START_REF] Carrera | Historical review of Zig-Zag theories for multilayered plates and shells[END_REF]-B and 2004), which require the calculation of the three axial displacements only: u x , u y , and u z as functions of six parameters. These variables were determined by minimizing the complementary potential. Good accuracy, at relatively high computational expense, is reached, but the models are based on an arbitrary choice of functions, representing the stress field, in order to capture the decaying rate of the stress perturbation as the distance from the free edge, or free corner, increases. Also, the through-the-thickness behaviour is described by means of first order Taylor's expansion, which represents a further axiomatic assumption. Despite this limitation, the CLA was recovered within a distance from the singularity comparable to the plate thickness. The models proposed by Mittelstedt and Becker appear more suitable for accurate local computation of stress distortions in proximity of singularities, than for rapidly detecting the influence of shear deformations on the global behaviour of multi-layered panels. Recently, Tessler et al.

(2009) have developed a refined Zigzag theory based on the kinematics of FSDT. The deployment of novel piecewise linear zigzag functions provides a robust and accurate approach, at relatively low (for layer-wise theories) computational expenses (six unknown variables only must be calculated).

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT 7

There is scope to develop a novel formulation, which requires reduced computational effort when compared to the models available in the open literature. Reducing the computational effort can be obtained by means of two different strategies:

1. reduction of the number of functional degrees of freedom;

2. choice of the degrees of freedom in such a way that the associated calculation requires less time to converge.

Indeed, reducing the number of variables does not necessarily imply a consequential reduction of the calculation time. If higher order variables are used as primary variable and solutions are calculated by means of numerical approaches such as the Galerkin technique or Rayleigh-Ritz method, convergence rate could be relatively quick [START_REF] Cosentino | Prebuckling and Buckling of Unsymmetrically Laminated Composite Panels with Stringer Run-Outs[END_REF] often obtained without difficulty in terms of the primary variables. Yet, for all the secondary variables, whose computation necessitates partial or total differentiations, the convergence rate might be slower or not obtainable at all [START_REF] Reddy | Energy and Variational Methods in Applied Mechanics[END_REF]. Furthermore, the accuracy of back-calculated variables generally decreases with the order of the required differentiation [START_REF] Reddy | Energy and Variational Methods in Applied Mechanics[END_REF].. As such, systematic and careful considerations should be used when choosing the primary variables and the relations that link them to the secondary variables.

The choice of primary variables and additional differential relations introduced in the present theory is therefore extensively commented upon and justified throughout the work.

The aim of the present study is to develop a flexible and robust self-contained formulation, able to reduce the computational effort by reducing the total number of variables involved, thus choosing the most appropriate set of primary variables. The present study demonstrates that an equivalent layer approach is consistent with, and can be directly derived from, a multi-layer approach under certain conditions. The theory developed is, to some extent, free from the limitations affecting the classical equivalent layer approaches, such as shear locking or membrane locking [START_REF] Reddy | Simple Finite Elements with Relaxed Continuity for Nonlinear Analysis of Plates[END_REF] . The simplifying kinematical assumptions on which CTP and CLA are based are initially retained. The through-the-thickness normal and shear stresses are then derived in an exact manner, by means of a stress-based approach, integrating the Cauchy's equations governing the elastic continuum. A mixed formulation is then derived by means of Reissner's variational approach based on the Castigliano's principle of least work in conjunction with a Lagrange multipliers method for the calculus of variations [START_REF] Reissner | The effect of transverse shear deformation on the bending of elastic plates[END_REF]. The derived field and

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT 8
boundary equations are initially expressed in terms of four functional degrees of freedom. To this end, the theory is free from the contraction of the number of boundary conditions that can be assigned, and a set of three dual boundary conditions is derived and can be fulfilled. The number of degrees of freedom can be reduced to two by means of an ad-hoc defined stress function. An additional novel boundary compatibility equation is also introduced. It is also demonstrated that, under certain simplifying conditions, the total number of degrees of freedom can be condensed to one only (the transverse displacement w), and an original transverse equilibrium equation is proposed, which contains correction terms derived from the contraction of the layer-wise, through-the-thickness properties on to the equivalent layer.

Transverse stresses derivation and formal analysis

Let us neglect the natural body forces acting upon a generic elastic continuum. The system of three Cauchy's equilibrium equations governing the elastic problem is: Let the tensors of in-plane stresses σ , in-plane strains ε , curvatures k, transverse shear stresses τ , and transverse shear strain γ be defined as follows: The following transformed in-plane and transverse stiffness matrices (namely Q and G ) are defined at each ply level [START_REF] Reddy | Energy and Variational Methods in Applied Mechanics[END_REF]:

0 = ∂ ∂ + ∂ ∂ + ∂ ∂ xz xy x z y x τ σ σ (1a) 0 = ∂ ∂ + ∂ ∂ + ∂ ∂ yz y xy z y x τ σ σ (1b) 0 = ∂ ∂ + ∂ ∂ + ∂ ∂ z yz xz z y x σ τ τ (1c)
      =       =             =             =             = yz
        =             = k k k k k k k k k k k k k k k
Q Q Q Q , Q Q Q Q Q Q Q Q Q G Q (3a, b)
In the linear regime, the following Kirchhoff's kinematical relations link membrane strains and curvatures to in-plane and transverse displacements:

x u o x ∂ ∂ = 0 ε , x v y u xy ∂ ∂ + ∂ ∂ = 0 0 0 ε , y v y ∂ ∂ = 0 0 ε (4a, b, c) 2 2 x w k x ∂ ∂ - = , 2 2 y w k x ∂ ∂ - = , y x w k xy ∂ ∂ ∂ - = 2 2
(5a, b, c) Following Classical Laminate Analysis (CLA), laminate strains and curvatures are related to in-plane forces and transverse bending moments by the constitutive equations, which are expediently expressed hereby in inverted form [START_REF] Reddy | Energy and Variational Methods in Applied Mechanics[END_REF]:

             =         M N D B B A k ε ' ' ' ' 0 T ( 6
)
where: The following constitutive equations hold true at each ply level (Fig. 2):

            =             =             =             =
k k k z z ≤ ≤ = 1 - k z , ε Q σ (8) M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT 10 k k k z z ≤ ≤ = 1 - k z , γ G τ (9)
Assuming a linear variation of the in-plane strains through the laminate thickness:

( )

) , ( ) , ( , , 0 y x z y x z y x k ε ε + = (10)
and substituting Eq. [START_REF] Cosentino | Approximate Non-linear Analysis Method for Debonding of Skin/Stringer Composite Assemblies[END_REF] in Eq. ( 6):

) ' ' ( ' ' M D N B M B N A ε + + + = T z ( 11
)
Substituting Equation [START_REF] Cosentino | Prebuckling and Buckling of Unsymmetrically Laminated Composite Panels with Stringer Run-Outs[END_REF] in Equation ( 8), the following expression for in-plane stresses is derived at the generic ply level k:

[ ] ) ' ' ( ' ' M D N B M B N A Q σ + + + = T k k z ( 12
)
Introducing the divergence ∇ and the differential operator C:

            ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ = x y y x 0 0 C ( 13
)
The governing system of equilibrium equations can be expressed in compact form at each ply level:

0 τ Cσ = ∂ ∂ + k k z (14) 0 = ∂ ∂ + ∇ k z k z σ τ ( 15
)
Substituting Eq. ( 12) in Eq. ( 14) and integrating with respect to z, a general solution for shear stresses is found as

( ) ( ) k k T k k z z a M B N A Q C M D N B Q C τ + + - + - = ' ' ' ' 2 1 2 (16)
The constants k a are calculated by enforcing the boundary conditions at each interface. Starting, for example, from the lower surface z = z 0 (in the present work, shear and normal transverse stresses are calculated starting from the lower surface of the laminate and this convention is held throughout the theory), and assuming that no shear stresses act upon the external surfaces of the laminate: [START_REF] Lekhnitskii | Anisotropic Plates[END_REF] the constants at ply level 1 are derived:

0 τ τ = = = = N z z z z 0
M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT 11 ( ) ( ) 0 2 0 1 ' ' ' ' 2 1 z z k T k M B N A Q C M D N B Q C a + - + - = (18)
The expression for shear stresses at ply level 1 is then:

( ) ( ) 1 1 2 1 1 ' ' ' ' 2 1 a z z T + + - + - = M B N A Q C M D N B Q C τ ( 19
)
The symmetry of the stress tensor in conjunction with local equilibrium satisfy the continuity of the shear stresses at each ply interface. This condition can be used to assign the boundary conditions for the subsequent ply 2:

1 1 2 1 z z z z = = = τ τ ( 20
)
Substituting Equations ( 16), [START_REF] Love | The mathematical Theory of Elasticity[END_REF], and [START_REF] Mindlin | Influence of rotatory inertia and shear on flexural motion of isotropic elastic plates[END_REF] in Equation ( 20) the following expression for 2 a is obtained:

( )( ) ( )( ) 1 1 2 2 1 1 2 1 2 ' ' ' ' 2 1 z z T M B N A Q Q C M D N B Q Q C a a + - + + - + = (21)
Equation ( 21) can be extended to the generic ply level k, and it can be straightforwardly verified that:

( )( ) ( )( ) ∑ = - - - -       + - + + - = k j j j j j T j j k z z 1 1 1 2 1 1 ' ' ' ' 2 1 M B N A Q Q M D N B Q Q C a (22)
where the transformed stiffness matrix 0 Q is conventionally made equal to the null matrix 0 .

Combining Equations ( 16) and [START_REF] Mittelstedt | A Variational finite layer technique for the investigation of thermally induced stress concentrations in composite structures[END_REF] the final expression for shear stresses at each ply level through the entire laminate thickness is derived:

( ) ( ) ( )( ) [ ] ∑ = - - + - + + - + - = k j j T j j k T k k z z z 1 2 1 1 2 ' ' 2 1 ' ' ' ' 2 1 M D N B Q Q C M B N A Q C M D N B Q C τ ( )( ) [ ] ∑ = - - + - + k j j j j z 1 1 1 ' ' M B N A Q Q C ( 23
)
Equation ( 23) can be re-written in the following more convenient form by ( ) ( )

) ( ' ' ) ( ' ' 2 1 1 2 1 2 - - - + - - + - = k k k T k k z z z z M B N A Q C M D N B Q C τ ( ) M D N B Q Q C ' ' ( ) ( 2 1 1 1 1 2 1 2 +       - + - -∑ - = - - T k j j j j j j j z z z z ( 24
)
Analyzing the expression found for the shear stresses, three components can be remarkably identified:

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT 12 i.
A local component, which is parabolic in the through-the-thickness coordinate z and that only depends upon local mechanical and geometrical properties.

ii. A local component, which is linear in the through-the-thickness coordinate z and only dependent upon local mechanical and geometrical properties.

iii.

A constant term that depends on the mechanical and geometrical properties of the laminae lying below the considered ply level.

Despite the cumbersome appearance, a physical interpretation of Equation ( 24) can be deduced. Defining the following matrices of partial in-plane and coupling stiffness (namely A and B ) :

( ) ( ) ( )dz z z z z z z z z j j k j j k k k ∫ ∑ = - + - = - - = - 0 1 1 1 1 ) ( Q Q Q A (25) ( ) ( ) ( )zdz z z z z z z z z j j k j j k k k ∫ ∑ = - + - = - - = - 0 2 1 2 1 1 2 1 2 2 1 2 1 ) ( Q Q Q B ( 26
)
Equation ( 24) can be re-written as:

( ) ( ) [ ] M D N B B M B N A A C τ ' ' ) ( ' ' ) ( + + + - = T k k k z z ( 27
)
The quantity between square brackets:

( ) ( ) M D N B B M B N A A N ' ' ) ( ' ' ) ( ) ( + + + = T k k z z z ( 28
)
represents the partial amount of normal forces carried by the part of laminate between the lower surface and the actual through-the-thickness z co-ordinate. Equation [START_REF] Reissner | On the theory of bending of elastic plates[END_REF] collapses to the well-known Jourawsky's formula for the isotropic one-dimensional case [START_REF] Jourawski | Sur la résistance d'un corps prismatique[END_REF]. Equation [START_REF] Reissner | The effect of transverse shear deformation on the bending of elastic plates[END_REF] expresses the translational in-plane equilibrium of the part of laminate above the z co-ordinate.

) ( ) ( ) ( ) ( 0 = - + ∂ ∂ + ∂ ∂ z z z N y z N x xz xz xy x τ τ (30a) 0 ) ( ) ( ) ( ) ( 0 = - + ∂ ∂ + ∂ ∂ z z z N y z N x yz yz y xy τ τ (30b)
The above system of two partial differential equations (PDE) can be readily compacted and reduced to Equation [START_REF] Reissner | On the theory of bending of elastic plates[END_REF] under the hypothesis that no shear stresses act upon the lower surface of the laminate (z = z 0 ). Let us reconsider Eq. ( 15). The transverse normal stress z σ can be obtained by means of a ply-by-ply integration with respect to z, once the shear stresses are known:

k k z z k z z z z dz k k ≤ ≤ ∇ - = - ∫ - 1 , 1 τ σ (31)
Substituting Eq. ( 16) in Eq. ( 31) and carrying out the integration, the following expressions for the transverse normal stress is obtained at the generic ply level k:

( ) ( ) k k k T k z b z z z σ k + ∇ - + ∇ + + ∇ = a M B N A Q C M D N B Q C 2 3 ' ' 2 1 ' ' 6 1 ( 32
)
Let us assume that no transverse loads act on the lower surface z = z 0 , while a transverse load per unit area p(x,y) act on the upper surface of the panel z = z N . The unknown constants b k are obtained by enforcing the required boundary conditions, i.e. the z-wise equilibrium at each ply interface

0 0 1 = =z z z σ , 1 1 1 - - - = = = k k k k z z z z z z σ σ (33a, b)
Starting from ply 1 and enforcing condition (33a) in Eq.( 32) the first constant b 1 is obtained:

( ) ( ) 0 1 2 0 1 3 0 1 1 ' ' 2 1 ' ' 6 1 z z z b T a M B N A Q C M D N B Q C ∇ + + ∇ - + ∇ - = (34)
Iterating the procedure at the generic ply level k and enforcing condition (33b):

( )( ) ( )( ) ∑ = - - - - - -       - ∇ - + - ∇ + + - ∇ = k j j j j j j j j T j j k z z z b 1 1 1 2 1 1 3 1 1 ) ( ' ' 2 1 ' ' 6 1 a a M B N A Q Q C M D N B Q Q C ( 35
)
Substituting expression (35) in Eq. ( 32) the following final expression for the transverse normal stress z σ at each ply level is derived:

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT 14 ( ) ( ) z z z k k T k z k a M B N A Q C M D N B Q C ∇ - + ∇ + + ∇ = 2 3 ' ' 2 1 ' ' 6 1 σ ( )( ) ( )( ) ∑ = - - - - - -       - ∇ - + - ∇ + + - ∇ + k j j j j j j j j T j j z z z 1 1 1 2 1 1 3 1 1 ) ( ' ' 2 1 ' ' 6 1 a a M B N A Q Q C M D N B Q Q C (36)
It is observed that the last boundary condition that can be enforced is related to the interface between ply N-1 and ply N (z = z N-1 ). Apparently, no conditions can be enforced on the value of z σ at the upper surface z = z N , which must equal the external load p. This condition is automatically fulfilled if the transverse vertical equilibrium is satisfied, as shown in Section 3.

Symmetric laminates

One major advantage of composites is their inherent ability to tailor elastic properties for maximizing performance. The prospect of selecting lay-ups with specific stacking provides a potential advantage of composite laminates against metals. However, symmetrically laminated composite plates still embody most of engineering composite structures. The possibility of modeling such configurations by means of analytical closed form solutions, due to the simplification introduced by the symmetry with respect to the mid-plane, is one explanation for their pervasive use. Quasi-isotropic and specially orthotropic lay-ups [START_REF] Ashton | Primer on Composite Materials: Analysis[END_REF], cross-ply and angle-ply composite plates remain the lay-up of choice for many composite components in aerospace, automotive and marine applications. There is, therefore, good reason for focusing the attention of the present theory on those composite plates that are symmetrical or quasisymmetrical in behaviour. For composites belonging to the class of quasi-symmetrical panels, the flexural coupling quantified by the magnitude of the components of the coupling matrix B is negligible. As a consequence, in the linear regime, the membrane and transverse structural responses can be analyzed separately, leading to a formal simplification of the mathematical model and to a significant reduction of the number of variables which must be computed simultaneously. The laminate constitutive equations [START_REF] Blaauwendraad | Timoshenko beam-column buckling[END_REF] collapse to the following:

N A ε ' 0 = , Μ = ' D k (37a, b)
De-coupling the in-plane forces N from the transverse bending moments M, and focusing on the transverse behaviour, the expressions of out-of-plane shear and normal stresses simplify:

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT 15 M D B C τ ' ) (z k k - = (38) ' ) ( M D F C z k z k ∇ = σ (39)
where

[ ] ∑ - = - - - + - = 1 1 1 1 ) ( ) ( ) ( ) ( ) ( k j j j j j k k k k z z z z z H H H H F ( 40
)
and

∫ ∑ ∫ ∫ - - - = + = = z z k j z z j k z z k k k j j dz z d d z 1 1 0 1 1 ) ( ) ( ) ( ) ( B B B H ζ ζ ζ ζ (41)
If p is the external transverse load per unit area acting upon the upper external surface z = z N , the boundary condition requires that

p z N z N = ) ( σ (42)
This condition is fulfilled only if vertical equilibrium equation is satisfied. Integrating over the total thickness with respect to z the third Cauchy's equilibrium equation (1c), the equilibrium in terms of vertical forces per unit width is derived,

0 0 0 0 = ∂ ∂ + ∂ ∂ + ∂ ∂ ∫ ∫ ∫ N N N z z z z z yz z z xz dz z dz y dz x σ τ τ (43)
Denoting V x and V y as the vertical shear resultants per unit width, by definition we have:

0 0 = - + ∂ ∂ + ∂ ∂ = = z z z z z z y x N V y V x σ σ (44)
Defining the vector V of shear force resultants:

      = y x V V V ( 45
)
and recalling the boundary conditions (33a) and (42), Equation (44) can be written as:

p - = ∇V (46)
Similarly, substituting Eqs. (33a), (38), and (41), the left hand side of Eq. ( 44) can be re-written as:

[ ] N k k k k z z z N k k k k k N k z z N k z z k k z z dz z dz = = - = = - =           - -∇ =                     -∇ = ∇ = ∇ ∑ ∑ ∫ ∑ ∫ - - σ M D H H C M D B C τ V ' ) ( ) ( ' ) ( 1 1 1 1 1 1 (47) M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT 16
Substituting Eq. ( 46) in Eq. ( 47), the boundary condition sought is obtained and Eq. ( 42) is formally recovered:

p N N z z z = = σ (48)
The identity (48), together with the fulfillment of Eq. ( 17), represents the first notable result of this theory.

The expression derived for shear and normal stresses fulfill the internal equilibrium and continuity of the stresses, intrinsically satisfying the boundary conditions as long as the transverse vertical equilibrium is satisfied. To this end, expressions for the transverse displacement (w), transverse shears (τ xz and τ yz ) and normal stress (σ z ) are derived, which are equilibrated and respect the natural boundary conditions.

Variational formulation

Under the assumption of state of plane stress throughout the domain, in order to obtain a system of governing field and boundary equations, the Castigliano's Theorem of Least Work is employed. The theorem states that, among all statically correct states of stresses, the state of stress which also satisfies the constitutive equations and the displacement boundary conditions is characterized by the condition that the variation of the following functional Π vanishes:

( )dsdz w u u dxdydz nz s ns n n T T ∫∫ ∫∫∫ + + - + = Π τ τ σ ε ) ( 2 1 γ τ σ (49)
The double integral is calculated over the thickness and the closed curve Γ defining the plate mid-plane on x-y plane. The co-ordinates n and s are, respectively, the normal and the tangent to the above curve (Fig. 3). [START_REF] Reissner | The effect of transverse shear deformation on the bending of elastic plates[END_REF], and recalling transverse equilibrium (44) expressed in terms of vertical shear resultants:

M A N U S C R I P T A C C E P T E D

0 = + ∂ ∂ + ∂ ∂ p V y V x y x (50)
The variation of expression (49) must equal zero in such a way that the equilibrium equation ( 50) is satisfied. In compliance with the rules of the calculus of variations, this is achieved by introducing a Lagrange multiplier λ(x,y) and adding an appropriate term to the variation of Castigliano's potential Π:

( )

   + + - + ∫∫ ∫∫∫ dsdz w u u dxdydz nz s ns n n T T τ τ σ δ ) ( 2 1 γ τ ε σ 0 =              + ∂ ∂ + ∂ ∂ + ∫∫ dydy p V y V x y x λ (51)
The variations are performed with respect to the three bending moments M x , M y , and M xy . Let us examine the variations of each component of the total potential separately,

0 = Π + Π + Π + Π Γ λ τ σ δ δ δ δ (52) where ∫∫ ∫ = Π N x y z z l l T dxdydz 0 0 0 2 1 ε σ σ (53a) M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT 18 ∫∫∫ = Π N x y z z l l T dxdydz 0 0 0 2 1 γ τ τ (53b) ( ) ∫ ∫ Γ Γ + + - = Π N z z nz s ns n n dsdz w u u 0 τ τ σ (53c) ∫ ∫         + ∂ ∂ + ∂ ∂ = Π x y l l y x dydy p V y V x 0 0 λ λ (53d)

a. Potential of in-plane stresses

Let us analyze the potential of internal in-plane stresses:

dxdydz z dxdydz dxdydz N k T k T T N k T k T M D Q D M ε σ ε σ ' ' 2 1 2 1 1 2 1 ∫∫∫ ∑ ∫∫∫ ∑ ∫∫∫ = = = = = Π δ δ δ δ σ ∫∫ ∫∫ ∑ ∫ =           = = - dxdy dxdy dz z T N k z z T k T T k k M D M M D Q D M ' ' ' 1 2 1 δ δ (54)
and with respect to M x , M y and M x is

∫ ∫         ∂ ∂ ∂ + ∂ ∂ ∂ + ∂ ∂ ∂ + ∂ ∂ + ∂ ∂ + ∂ ∂ = Π x y l l xy xy y xy x xy xy x y x x x M y x M y x M y x M x M x M x 0 0 2 13 2 12 2 11 2 2 13 2 2 12 2 2 11 α α α α α α δ τ     ∂ ∂ ∂ + ∂ ∂ + ∂ ∂ + ∂ ∂ +     ∂ ∂ + ∂ ∂ + ∂ ∂ x xy xy x y x x x x xy y y y x y M y x M x M x M x M M y M y M y
α α α α δ α α α M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT 19     ∂ ∂ +     ∂ ∂ + ∂ ∂ + ∂ ∂ + ∂ ∂ ∂ + ∂ ∂ ∂ + x x
∂ ∂ + ∂ ∂ + ∂ ∂ ∂ + ∂ ∂ ∂ + ∂ ∂ ∂ + ∂ ∂ + ∂ ∂ + α α α α α α α ∫ Γ         ∂ ∂ + ∂ ∂ + ∂ ∂ + ∂ ∂ + ∂ ∂ +         ∂ ∂ +
M M y M y M y M x M x M x M M y δ β β β β β β δ β     ∂ ∂ +     ∂ ∂ + ∂ ∂ + ∂ ∂ + ∂ ∂ + ∂ ∂ +     ∂ ∂ + 23 22 21 23 22 21 13 ds M M y M y M y M x M x M x xy xy y y y x y xy x y x x x         ∂ ∂ +     ∂ ∂ + ∂ ∂ + ∂ ∂ + ∂ ∂ + ∂ ∂ + δ β β β β β β

c. Potential of external contour loads

The final expression of the potential of external contour loads is :

( ) ∫ ∫ Γ Γ + + - = Π N z z nz s ns n n dsdz w u u 0 τ τ σ δ δ (57)

d. Potential of the external transverse loads

The potential of the external transverse loads is (see Appendix B for details):

∫ ∫ ∫ Γ             ∂ ∂ + ∂ ∂ + ∂ ∂ -         ∂ ∂ + ∂ ∂ ∂ + ∂ ∂ = Π x l y l xy x y x x y xy x M y n x n M x n dxdy M y M y x M x 0 0 2 2 2 2 2 2 δ λ λ δ λ δ λ δ λ δ λ δ λ ds V s M y n n y y    - ∂ ∂ + δ λ δ λ ) ( (58)

Governing equations

In order to derive the final set of governing equations, all the terms multiplying each variation must be added and equated to zero. However, a proper interpretation of the terms representing each potential is required. It is emphasized that the three mechanical primary variables chosen are the three transverse bending moments M x , M y , and M xy . However, it must be highlighted that two terms contained in Eqs. ( 57) and (58) respectively, appear as multipliers of the variation of the normal vertical shear resultant V n . These terms were not further developed for convenience. Equating their summation to zero, the following equation is obtained:

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT 20 ) ( ) ( s s w λ = (59)
Since the same result is obtained for any generic inner closed contour Γ, the Lagrange multiplier λ coincides with the transverse displacement w [START_REF] Reissner | On the theory of bending of elastic plates[END_REF]. This fundamental statement allows derivation of the final set of field and boundary equilibrium equations. Equating to zero all the aforementioned components of Eqs. ( 55), ( 56), (57), and (58), we obtain: In addition to Eqs. (60a), (60b), and (60c), the following equilibrium equation holds true:

y xy x xy xy x y x x x xy y x M y x M y x M x M x M x M D M D M D ∂ ∂ ∂ + ∂ ∂ ∂ + ∂ ∂ + ∂ ∂ + ∂ ∂ + + +
∂ ∂ - = ∂ ∂ + ∂ ∂ + ∂ ∂ + ∂ ∂ ∂ + α α α α (60a) y xy x xy xy x y x x x xy y x M y x M y x M x M x M x M D M D M D ∂ ∂ ∂ + ∂ ∂ ∂ + ∂ ∂ + ∂ ∂ + ∂ ∂ + + +
∂ ∂ - = ∂ ∂ + ∂ ∂ + ∂ ∂ + ∂ ∂ ∂ + α α α α (60b) y xy x xy xy x y x x x xy y x M y x M y x M x M x M x M D M D M D ∂ ∂ ∂ + ∂ ∂ ∂ + ∂ ∂ + ∂ ∂ + ∂ ∂ + + +
p M y M y x M x y xy x - = ∂ ∂ + ∂ ∂ ∂ + ∂ ∂ 2 2 2 2 2 2 (60d)
The essential boundary conditions are:

0 13 12 11 13 12 11 = ∂ ∂ - ∂ ∂ + ∂ ∂ + ∂ ∂ + ∂ ∂ + ∂ ∂ + ∂ ∂ x w n M y M y M y M x M x M x x xy y y y x y xy x y x x x β β β β β β (61a) 0 23 22 21 23 22 21 = ∂ ∂ - ∂ ∂ + ∂ ∂ + ∂ ∂ + ∂ ∂ + ∂ ∂ + ∂ ∂ y w n M y M y M y M x M x M x
= ∂ ∂ - ∂ ∂ - ∂ ∂ + ∂ ∂ + ∂ ∂ + ∂ ∂ + ∂ ∂ + ∂ ∂ y w n x w n M y M y M y M x M x M x
            ∂ ∂ ∂ ∂       - =             ∂ ∂ ∂ ∂
= ∂ ∂ + ∂ ∂ + ∂ ∂ + ∂ ∂ + ∂ ∂ + ∂ ∂ xy y y y x y xy x y x x x M y M y M y M x M x M x η η η η η η (66) M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT 22

Reduction of primary variables: Potential of bending moments

In order to reduce the number of primary variables, it is convenient to express Eqs. (60a), (60b), and (60c) in the following compact form:

( )

w y y x x y xy x d M α M α M α M D = ∂ ∂ + ∂ ∂ ∂ + ∂ ∂ + 2 2 2 2 2 ' ( 67
)
where the following differential operator is introduced:

T y x y x         ∂ ∂ ∂ ∂ ∂ ∂ ∂ - = 2 2 2 2 2 2 d ( 68
)
Multiplying both sides by the bending stiffness matrix D:

( )

w y y x x y xy x Dd M Dα M Dα M Dα M = ∂ ∂ + ∂ ∂ ∂ + ∂ ∂ + 2 2 2 2 2 (69a) or Dk M Dα M Dα M Dα M = ∂ ∂ + ∂ ∂ ∂ + ∂ ∂ +
x x x M y x M y x M y x M x M x M x M ∂ ∂ ∂ Λ + ∂ ∂ ∂ Λ + ∂ ∂ ∂ Λ + ∂ ∂ Λ + ∂ ∂ Λ + ∂ ∂ Λ +
∂ ∂ ∂ - ∂ ∂ - ∂ ∂ - = ∂ ∂ Λ + ∂ ∂ Λ + ∂ ∂ Λ +
M y x M y x M y x M x M x M x M ∂ ∂ ∂ Λ + ∂ ∂ ∂ Λ + ∂ ∂ ∂ Λ + ∂ ∂ Λ + ∂ ∂ Λ + ∂ ∂ Λ +
∂ ∂ ∂ - ∂ ∂ - ∂ ∂ - = ∂ ∂ Λ + ∂ ∂ Λ + ∂ ∂ Λ +
M y x M y x M y x M x M x M x M ∂ ∂ ∂ Λ + ∂ ∂ ∂ Λ + ∂ ∂ ∂ Λ + ∂ ∂ Λ + ∂ ∂ Λ + ∂ ∂ Λ +
M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT 23 A potential function Ω is now introduced, 2 2 x M x ∂ Ω ∂ = , y x M xy ∂ ∂ Ω ∂ = 2 , 2 2 y M y ∂ Ω ∂ = (72a, b, c)
By substituting Eqs. (72) in to the transverse equilibrium equation (60d), the following biharmonic equation is obtained which governs the transverse equilibrium of the panel:

p y y x x - = ∂ Ω ∂ + ∂ ∂ Ω ∂ + ∂ Ω ∂ 4 4 2 2 4 4 4 2 (73)
Differentiating Eq. (71a) twice with respect to x, Eq. (71b) twice with respect to y, and multiplying both sides of Eq.(71c) by two and differentiating twice with respect to x and y and adding the results, in view of Eq. ( 73), the following expression is obtained:

( )

y x c x c p y w D y x w D y x w D D y x w D x w D xxxxxy xxxxxx ∂ ∂ Ω ∂ + ∂ Ω ∂ + = ∂ ∂ + ∂ ∂ ∂ + ∂ ∂ ∂ + + ∂ ∂ ∂ + ∂ ∂
∂ Ω ∂ + ∂ ∂ Ω ∂ + ∂ ∂ Ω ∂ + ∂ ∂ Ω ∂ + ∂ ∂ Ω ∂ + (74)
where the following shear coefficients are introduced:

2 2 2 Λ - Λ - Λ - Λ - Λ - Λ - Λ - = (75g)
Equations ( 73) and (74) provide a complete and self-consistent mathematical representation of the elastic problem of a composite multi-layered panel. Equation (69b) establishes a novel constitutive relation between bending moments and curvatures, which involves new differential correction factors. The algebraic relations derived from the CLA, i.e. Eqs. [START_REF] Blaauwendraad | Timoshenko beam-column buckling[END_REF], are now transformed to differential equations.

Clearly, for linear problems involving transverse loads, the calculation of the internal bending moments remains the same as in the CLA. In contrast, the transverse displacements, which can be calculated once the potential Ω is known, are affected by the inclusion of the shear stresses. It is emphasized that three boundary conditions can be still fulfilled instead of two as in the CLA.

A suitable solution in terms of transverse displacements is:

( )

x w y x w n n θ sin ) , ( = (79)
Substituting Eq. (79) in Eq. ( 78), the following final expression is obtained:

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT 25 ( ) xxxxxx n n n n c D p w 2 4 11 1 θ θ + = (80)
It is evident that the solution found differs from the solution predicted by the CLA due to the presence of the coefficient c xxxxxx . The ratio between the maximum amplitudes of the transverse displacements predicted by the two models is The deviation of the present model from the CLA was calculated for five different stacking sequences. The assumed material properties are reported in Table 1. The quantity c xxxxxx /l x 2 was calculated and plotted against the ratios t/l x (Fig. 5.A) and (t/ l x ) 2 (Fig. 5.B),where t is the plate total thickness. Clearly, the differences increase with a parabolic trend as the ratio t/l x increases. Figure 5.A shows that the trends of the deviation are linear with (t/ l x ) 2 . Hence the constant c xxxxxx is proportional to the total thickness squared, according to a factor which is function of the lay up and the stacking sequence. This parabolic trend is well-known for thick isotropic plates and for thick plates of one orthotropic lamina. Therefore, results

presented in figure 5 demonstrate that a theoretical generalization to laminated composite plates is possible if shear correction factors are calculated according to Equation (75g).

For each laminate considered, the differences tend to zero as the thickness over length ratio tends to zero.

The linear relation found implies that only one extra calculation of the coefficient c xxxxxx /l x 2 is needed for each analyzed configuration. Once two points are known for each configuration, the linear relation allows calculation of the difference for any value of t/l x . As a general qualitative comment, laminates with larger bending stiffnesses show larger differences from CLA. The transverse displacement predicted by the two models were calculated in the case of an external pressure law which represents one half wave in the x-direction. Results are illustrated in Fig. 6, where the ratios w n,present /w n,CLA are plotted for each configuration against the t/l x ratio. The plate is loaded by an external in-plane load acting in the x-direction. The critical buckling loads are sought. Following the von Karman formulation for moderately large transverse displacements [START_REF] Timoshenko | Woinowsky-Krieger, S. Theory of Plates and Shells[END_REF], it is convenient to re-write the governing equations as follows:

) cos( ) sin( 4 3 2 1 x c x c x c c w ω ω + + + = (85)
where:

x xxxxxx x N c D N + = 11 2 ω (86)
The presented simplified form of the buckling equation (86) does not satisfy the complete set of three dual boundary conditions (61) and (62). The fulfillment of three dual boundary conditions requires a sixth-order equation to be solved. However, under the caveats provided by [START_REF] Blaauwendraad | Timoshenko beam-column buckling[END_REF], the use of a simplified approximate solution allows a direct comparison of the present model with the classical solution proposed by [START_REF] Timoshenko | Theory of elastic stability[END_REF]. The comparison also shows that classical accepted relations are directly derived by an approximate solution of the present theory and extendable to multi-layered composite plates. Let us assume the solution to be in the form:

) sin( x w w ω = (87)
The boundary conditions (61d) lead to the following condition:

π ω n l x = (88)
Substituting Eq. (86) in Eq (88) and choosing the minimum value of the critical load, we obtain the following expression:

        + =
Buckling loads can be readily derived by means of Equation ( 95) by performing a linear analysis on the equivalent simpler problem of a beam loaded by a sinusoidal transverse load, thus reducing significantly the calculation time.

Special two-dimensional cases. Closed form solutions a. Simply supported rectangular multi-layered plate subjected to sinusoidal transverse load

A simply supported rectangular plate is subjected to a sinusoidal transverse pressure defined by the following law:

                - = ∂ Ω ∂ + ∂ ∂ Ω ∂ + ∂ Ω ∂
∂ Ω ∂ + ∂ ∂ Ω ∂ + ∂ ∂ Ω ∂ + (98)
From Eq.( 97), a suitable solution in terms of potential Ω is obtained:

( ) ( )

sin sin y x n m mn β α Ω = Ω ,
β β α α β β α β α α (100b)
As in the one-dimensional case, the solution found differs from the solution predicted by the CLA due to the presence of the shear coefficients. The ratio between the maximum amplitudes of the transverse displacements predicted by the two models is:

Χ + Χ + Χ + Χ + Χ + + - = ( ) ( ) [ ] ( ) ( )} y x n m n m n m k k k k β α β α β α cos cos 2 2 12 31 3 23 32 3 Χ + Χ + Χ + Χ - (104)
where the components of matrix Ψ Ψ Ψ Ψ are defined in Appendix A, and the following matrix has been introduced:

(z)D' F Χ k = k (105)

b. Compression buckling of simply supported rectangular multi-layered plate

A simply supported rectangular plate is subjected to a compressive load per unit width N x . Introducing the Laplace operator ∇ , the two governing equations can be written as:

∂ Ω ∂ + ∂ ∂ Ω ∂ + ∂ ∂ Ω ∂ + (107)
The following solutions are assumed:

+ + + = (111a) ( ) ( ) ( )( ) 4
D D D D c c c c N β β α α β β α α β β α β α α α + + + + + + + + = (111b)
Equation ( 110) can be re-written in the form of : Material properties and laminate configurations and stacking sequences are provided in Tables 1 and2, respectively. Laminates are simply supported along all the edges.

s cr x CLA cr x present cr x N N N , , , , , , 1 1 1 +
Comparisons in terms of maximum deflection are presented in Tables 3 and4. Results obtained by means of different theories are indicated as

• Exact 3D Elasticity -fully three-dimensional solution developed by [START_REF] Pagano | Exact solutions for composite laminates in cylindrical bending[END_REF].

Implementation of the method and results are provided by [START_REF] Tessler | Refined Zigzag Theory for Laminated Composite and Sandwich Plates[END_REF] .

• FSDT -First order Shear Deformation Theory.

• Zigzag -Zigzag Theory developed by [START_REF] Tessler | Refined Zigzag Theory for Laminated Composite and Sandwich Plates[END_REF].

• Refined Zigzag -Transverse shear stresses obtained by integration of the equilibrium equations (1) using the in-plane stresses derived by [START_REF] Tessler | Refined Zigzag Theory for Laminated Composite and Sandwich Plates[END_REF].

• CLA -Classical Laminate Analysis.

• FEM -Results extracted from ad-hoc three-dimensional finite element models.

In Table 3, comparisons of the maximum displacements calculated by means of a variety of numerical solutions are presented, for a span-to-thickness ratio equal to 5. Results show that the present theory is in very good agreement with solutions calculated by [START_REF] Tessler | Refined Zigzag Theory for Laminated Composite and Sandwich Plates[END_REF] and [START_REF] Pagano | Exact solutions for composite laminates in cylindrical bending[END_REF] for configurations B and C. Some discrepancies were found in comparing results for configurations A and D, differences of the order of 15% and 20% respectively. For both configurations A and D, the maximum deflections calculated by means of the present approach appear to be slightly over-estimated. To understand the cause of such a discrepancy, the normalized transverse shear stresses τ xz (0,l x /2,z) (normalization factor equal to 2h/p 0 l x ) are plotted against the normalized through-the-thickness coordinate in Fig. 9. Clearly, the present model calculates the shear stresses within the sandwich core with good accuracy, but under- Good agreement is also found in the case of configuration D.

To provide an overview on the range of applicability of the present model, a sensitivity study was undertaken on laminate configuration B. Several span-to-thickness ratios were analyzed and normalized maximum transverse displacement were calculated and compared to results calculated by [START_REF] Tessler | Refined Zigzag Theory for Laminated Composite and Sandwich Plates[END_REF]. Results are presented in Table 4. It is noted that the FSDT generally yields over stiff results, while the present theory provide very accurate results for values of the span-to-thickness ratio equal to 10. For relatively thicker configurations , i.e. span-to-thickness ratio equal to 4, the error induced increases up to 30%.

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT

36

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT

41

A three-dimensional high-fidelity finite element model was created with ABAQUS. The domain was regularly discretized using three elements per each layer and 50 subdivision along each span direction, for a total of 60,000 quadratic hexahedron -type C3D20R elements (ABAQUS User manual, 2004). In consideration of potential numerical problems, such as shear locking or hourglassing, quadratic elements with reduced integration were used to counter the former, whilst mesh convergence studies and comparisons with other methods [START_REF] Tessler | Refined Zigzag Theory for Laminated Composite and Sandwich Plates[END_REF] provided confidence that the latter was suppressed.

As a first case-study, and for all configurations, buckling loads were calculated for a compressive longitudinal load N x . To assess the effect of the span-to-thickness ratio l x /2h, calculations were performed for five different ratios (l x /2h) equal to 500, 50, 20, 10 and 5 respectively). Results are presented in Table 6.

As expected, the normalized buckling loads calculated by CLA are constant and do not depend upon the span-to-thickness ratio. For span-to-thickness ratios that are larger than 10, the difference between the present theory and FEM are within a 10%-15% range. Significant discrepancies (larger than 15%) are found as soon as l x /2h=5. However, within the range of reasonably thick laminates, the present analysis shows good accuracy and the inclusion of through-the-thickness normal stresses in the variational statement (49) does not appear to be required. To assess the convergence rate of the solution, a preliminary convergence study was performed on configuration L1 for a series of values of M and N. For the sake of simplicity it was assumed M=N in all the cases analyzed. Convergence of transverse displacements is obtained at relatively small computational expenses The maximum transverse displacement appears to converge when more than three base functions are used per each span co-ordinate (Fig. 16). In contrast, more base functions are needed to guarantee convergence of higher order variables. The case of through-the-thickness normal stresses was assessed in configuration L1, subjected to a uniform pressure and with a span-to-thickness ratio equal to 1. Results are reported in Fig. 17, where the normalized stresses are plotted against the normalized thickness co-ordinate for values of M and N varying from one to thirteen.

M A N U S C R I P T A C C E P T E D

To fulfill the natural boundary condition at the upper interface, the value of σ z /p 0 must equal 1, i.e. σ z in equilibrium with the external pressure. Notably, the fulfillment of such a boundary condition is obtained when 13 or more base functions are used. However, it must be emphasized that the convergence rate is, in this case, directly influenced by the degree of accuracy and convergence rate associated with the discretization of the external load in generalized Fourier series obtained by the Galerkin technique employed. More comparison between the present theory and results obtained by means of FEM are presented in Figs. 18 to 21, where the case of uniform pressure is analyzed and normalized normal through-the-thickness stresses, calculated in the center of the laminate (x=l x /2, y=l y /2) are plotted against the normalized thickness co-ordinate. The values M=N=13 were chosen to perform the simulations. Comparisons between the calculated transverse displacements are presented for all the configurations in Table 7. A uniform pressure was assumed to load the top surface of the panel and two span-to-thickness ratios, five and ten respectively, were chosen.

M A N U S C R I P T A C C E P T E D

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT 47 analyzed.
Notably, the error associated with the results provided by CLA is of the order of 10% for span-tothickness ratios larger than five.

Conclusions

A novel theory was developed to assess the effect of transverse shear stresses on bending and buckling of reasonably thick composite laminates and sandwich plates. The contribution of transverse stresses was included in the variational statement governing the elastic equilibrium of the body and a novel set of equilibrium equations and boundary conditions was derived. To reduce the final number of unknown variables, a stress function was introduced. Simple problems were solved and a variety of closed form solutions were derived for simple one-dimensional and two-dimensional cases. The continuity of the shear stresses at each lamina interface, which is one of the most common weaknesses of other theories, is intrinsically satisfied by the present approach. Another important attribute of the present theory is that the total number of unknown variables needed to fully determine the stress field throughout the domain is reduced to two only, rendering the present approach extremely efficient in terms of computational expenses. Furthermore, the theory is particularly suitable for meshless approaches, which represent a robust alternative to FEM. An intrinsically three-dimensional problem can be reduced to an equivalent twodimensional single layer approach by appropriately condensing the layer-wise properties. The present theory is, of course, less appropriate when very thick configurations are analyzed. However, the range of validity is sufficiently large to include a significant variety of structural problems, especially in the aerospace industry.

Results were found to be in very good agreement with other more refined theories and with high-fidelity FEM models, suggesting that the present approach can be an alternative, within its range of validity to other more accurate, but less computationally efficient, approaches.

APPENDIX A Potential of transverse shear stresses

In order to carry out the variations, it is convenient to express the potential as an explicit function of the transverse shear stresses τ xz and τ yz . Defining the following matrix:

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT 48 ' ) ( D B Ψ z k k - = (A.1)
and considering Eq. ( 38), the following explicit expressions are derived for shear stresses as functions of the transverse bending moments,

xy y x xy y x xz M y Ψ M y Ψ M y Ψ M x Ψ M x Ψ M x Ψ k k k k k k k ∂ ∂ + ∂ ∂ + ∂ ∂ + ∂ ∂ + ∂ ∂ + ∂ ∂ =
Ψ M x Ψ M x Ψ M y Ψ M y Ψ M y Ψ k k k k k k k ∂ ∂ + ∂ ∂ + ∂ ∂ + ∂ ∂ + ∂ ∂ + ∂ ∂ =
∫ ∫ ∑ ∫ ∫∫∫           = = Π = - - y y k k l l N k z z k T T dxdy dz dxdydz 0 0 1 1 1 2 1 2 1 τ G τ γ τ τ (A.4)
The variation of the potential is

∫ ∫ ∑ ∫ ∫ ∫ ∑ ∫           +           = Π = = - -
∫ ∫ ∫ ∑                 ∂ ∂ Ψ + ∂ ∂ Ψ + ∂ ∂ Ψ + ∂ ∂ Ψ + ∂ ∂ Ψ + ∂ ∂ Ψ = Π - = x y k k k k k k k k k l l z z xy y x xy y x N k xz xz M y M y M y M x M x M x G 0 0 33 32 31 13 12 11 1 1 1 δ dxdy dz M y M y M y M x M x M x xy y x xy y x k k k k k k                 ∂ ∂ Ψ + ∂ ∂ Ψ + ∂ ∂ Ψ + ∂ ∂ Ψ + ∂ ∂ Ψ + ∂ ∂ Ψ × δ δ δ δ δ δ
∫∫ ∫ ∑                  ∂ ∂ Ψ + ∂ ∂ Ψ + ∂ ∂ Ψ + ∂ ∂ Ψ + ∂ ∂ Ψ + ∂ ∂ Ψ = Π - = x y k k k k k k k k k l l z z xy y x xy y x N k yz yz M x M x M x M y M y M y G 0 0 33 32 31 23 22 21 1 1 1 δ dxdy dz M x M x M x M y M y M y xy y x xy y x k k k k k k                  ∂ ∂ Ψ + ∂ ∂ Ψ + ∂ ∂ Ψ + ∂ ∂ Ψ + ∂ ∂ Ψ + ∂ ∂ Ψ × δ δ δ δ δ δ
∫ ∫ ∑ ∫ =              ∂ ∂ ∂ Ψ + ∂ ∂ ∂ Ψ + ∂ ∂ Ψ + ∂ ∂ Ψ + ∂ ∂ Ψ Ψ - = Π - x y k k k k k k k k k xz l l N k z z y x xy y x xz M y x M y x M x M x M x G 0 0 1 2 32 2 31 2 2 13 2 2 12 2 2 11 11 1 τ δ     ∂ ∂ Ψ + ∂ ∂ Ψ + ∂ ∂ ∂ Ψ + ∂ ∂ ∂ Ψ + ∂ ∂ ∂ Ψ Ψ +     ∂ ∂ ∂ Ψ + y x xy y x xz xy M y M y M y x M y x M y x G M y x k k k k k k k k 2 2
∫ ∑ ∫ Γ =             ∂ ∂ Ψ + ∂ ∂ Ψ + ∂ ∂ Ψ + ∂ ∂ Ψ Ψ +              ∂ ∂ Ψ + - N k z z x xy y x xz x xy k k k k k k k k k M y M x M x M x G M dxdy dz M y 1 31 13 12 11 11 2 2 33 1 δ     ∂ ∂ Ψ + ∂ ∂ Ψ + ∂ ∂ Ψ + ∂ ∂ Ψ + ∂ ∂ Ψ Ψ +     ∂ ∂ Ψ + ∂ ∂ Ψ + y x xy y x xz x xy y M y M y M x M x M x G n M y M y k k k k k k k k k
∫∫ ∑ ∫ =              ∂ ∂ ∂ Ψ + ∂ ∂ Ψ + ∂ ∂ Ψ + ∂ ∂ Ψ Ψ -              ∂ ∂ Ψ + - x y k k k k k k k k k l l N k z z x xy y x xz x y xy M y x M x M x M x G M ds n M y 0 0 1 2 31 2 2 13 2 2 12 2 2 11 12 33 1 δ     ∂ ∂ Ψ + ∂ ∂ ∂ Ψ + ∂ ∂ ∂ Ψ + ∂ ∂ ∂ Ψ Ψ +     ∂ ∂ ∂ Ψ + ∂ ∂ ∂ Ψ +
    ∂ ∂ Ψ + ∂ ∂ Ψ + ∂ ∂ Ψ + ∂ ∂ Ψ Ψ +     ∂ ∂ Ψ + ∂ ∂ Ψ + ∂ ∂ Ψ + x xy y x xz x xy y x M y M x M x M x G n M y M y M y k k k k k k k k k 31 13
∫∫ ∑ ∫ =              ∂ ∂ Ψ + ∂ ∂ Ψ + ∂ ∂ Ψ Ψ -              ∂ ∂ Ψ + ∂ ∂ Ψ + - x y k k k k k k k k k l l N k z z xy y x xz y y xy y M x M x M x G M ds n M y M y 0 0 1 2 2 13 2 2 12 2 2 11 13 33 32 1 δ     ∂ ∂ ∂ Ψ + ∂ ∂ ∂ Ψ + ∂ ∂ ∂ Ψ Ψ +     ∂ ∂ ∂ Ψ + ∂ ∂ ∂ Ψ + ∂ ∂ ∂ Ψ + xy y x xz xy y x M y x M y x M y x G M y x M y x M y x k k k k k k k k
∫ ∑ ∫ Γ =             ∂ ∂ Ψ + ∂ ∂ Ψ Ψ +              ∂ ∂ Ψ + ∂ ∂ Ψ + ∂ ∂ Ψ + - N k z z y x xz x xy y x k k k k k k k k k M x M x G M dxdy dz M y M y M y
  ∂ ∂ Ψ + ∂ ∂ Ψ + ∂ ∂ Ψ Ψ +     ∂ ∂ Ψ + ∂ ∂ Ψ + ∂ ∂ Ψ + ∂ ∂ Ψ + xy y x xz x xy y x xy M x M x M x G n M y M y M y M x k k k k k k k k k
∫ ∫ ∑ ∫ =              ∂ ∂ ∂ Ψ + ∂ ∂ ∂ Ψ + ∂ ∂ Ψ + ∂ ∂ Ψ + ∂ ∂ Ψ Ψ - = Π - x y k k k k k k k k k
    ∂ ∂ Ψ + ∂ ∂ Ψ + ∂ ∂ ∂ Ψ + ∂ ∂ ∂ Ψ + ∂ ∂ ∂ Ψ Ψ +     ∂ ∂ ∂ Ψ + y x xy y x yz xy M x M x M y x M y x M y x G M y x k k k k k k k k
∫ ∑ ∫ Γ =              ∂ ∂ Ψ + ∂ ∂ Ψ + ∂ ∂ Ψ + ∂ ∂ Ψ Ψ +              ∂ ∂ Ψ + - N k z z x xy y x xz x xy k k k k k k k k k M x M y M y M y G M dxdy dz M x
M x M x M y M y M y G n M x M x k k k k k k k k k ∂ ∂ Ψ +     ∂ ∂ Ψ + ∂ ∂ Ψ + ∂ ∂ Ψ + ∂ ∂ Ψ Ψ +    ∂ ∂ Ψ + ∂ ∂ Ψ +
∫∫ ∑ ∫ =              ∂ ∂ ∂ Ψ + ∂ ∂ Ψ + ∂ ∂ Ψ + ∂ ∂ Ψ Ψ -          ∂ Ψ + - x y k k k k k k k k k l l N k z z x xy y x yz x x xy M y x M y M y M y G M ds n M x 0 0 1 2 31 2 2 23 2 2 22 2 2 21 22 33 1 δ     ∂ ∂ Ψ + ∂ ∂ ∂ Ψ + ∂ ∂ ∂ Ψ + ∂ ∂ ∂ Ψ Ψ +     ∂ ∂ ∂ Ψ + ∂ ∂ ∂ Ψ + x xy y x yz xy y M x M y x M y x M y x G M y x M y x k k k k k k k k
∫ ∑ ∫ Γ =              ∂ ∂ Ψ + ∂ ∂ Ψ + ∂ ∂ Ψ Ψ +              ∂ ∂ Ψ + ∂ ∂ Ψ + - N k z z xy y x xz y xy y k k k k k k k k k M y M y M y G M dxdy dz M x M x 1 23 22 21 22 2 2 33 2 2 32 1 δ     ∂ ∂ Ψ + ∂ ∂ Ψ + ∂ ∂ Ψ + ∂ ∂ Ψ Ψ +    ∂ ∂ Ψ + ∂ ∂ Ψ + ∂ ∂ Ψ + x xy y x xz y xy y x M x M y M y M y G n M x M x M x k k k k k k k k k 31 23 22 21 32 33 32 31 ∫ ∫ ∑ ∫ =              ∂ ∂ Ψ + ∂ ∂ Ψ + ∂ ∂ Ψ Ψ -          ∂ ∂ Ψ + ∂ ∂ Ψ + - x y k k k k k k k k k l l N
δ     ∂ ∂ ∂ Ψ + ∂ ∂ ∂ Ψ + ∂ ∂ ∂ Ψ Ψ +     ∂ ∂ ∂ Ψ + ∂ ∂ ∂ Ψ + ∂ ∂ ∂ Ψ + xy y x yz xy y x M y x M y x M y x G M y x M y x M y x k k k k k k k k
∫ ∑ ∫ Γ =              ∂ ∂ Ψ + ∂ ∂ Ψ Ψ +              ∂ ∂ Ψ + ∂ ∂ Ψ + ∂ ∂ Ψ + - N k z z y x xz xy xy y x k k k k k k k k k M y M y G M dxdy dz M x M x M x 1 22 21 23 2 2 33 2 2 32 2 2 31 1 δ     ∂ ∂ Ψ + ∂ ∂ Ψ + ∂ ∂ Ψ Ψ +    ∂ ∂ Ψ + ∂ ∂ Ψ + ∂ ∂ Ψ + ∂ ∂ Ψ + xy y x xz y xy y x xy M y M y M y G n M x M x M x M y k k k k k k k k k 23 22 21 33 33 32 31 23 xy x xy y x M ds n M x M x M x k k k δ          ∂ ∂ Ψ + ∂ ∂ Ψ + ∂ ∂ Ψ + 33 32 31 (A.8.b)
Substituting expressions (A.8) in Eq. (A.5), performing the integrations with respect to z, and factoring out with respect to the three virtual variations, namely δM x , δM y , and δM xy , the following field and boundary integrals are obtained: Generally, exact solutions can be found for few particular cases only, and the use of approximate solutions in conjunction with numerical methods to minimize the induced error is required. Therefore Eq. (C.5) represents a supplementary compatibility condition over the boundary of the domain.

∫ ∫         ∂ ∂ ∂ + ∂ ∂ ∂ + ∂ ∂ ∂ + ∂ ∂ + ∂ ∂ + ∂ ∂ = Π
M M y M y M y M y x M y x M y x δ α α α α α α     ∂ ∂ + ∂ ∂ + ∂ ∂ + ∂ ∂ ∂ + ∂ ∂ ∂ + ∂ ∂ ∂ +
M y M y M x M x M x dxdy M M y M y ∂ ∂ + ∂ ∂ + ∂ ∂ + ∂ ∂ +       ∂ ∂ +         ∂ ∂ + ∂ ∂ + ∫ Γ 12
M M y M y M y M x M x M x M M y δ β β β β β β δ β     ∂ ∂ +     ∂ ∂ + ∂ ∂ + ∂ ∂ + ∂ ∂ + ∂ ∂ +     ∂ ∂ +
       ∂ ∂ +     ∂ ∂ + ∂ ∂ + ∂ ∂ + ∂ + ∂ ∂ + δ β β β β β β
∑ ∫ =                 Ψ Ψ + Ψ Ψ - = - N k z z yz xz x dz G G k k k k k k k k 1 33 31 13 11 13 1 α (A.11.c) ∑ ∫ =                 Ψ + Ψ - = - N k z z yz xz x dz G G k k k k k k 1 2 32 2 12 22 1 α (A.11.d) ∑ ∫ =                 Ψ Ψ + Ψ Ψ - = - N k z z yz xz x dz G G k k k k k k k k 1 33 32 13 12 23 1 α (A.11.e) ∑ ∫ =                 Ψ + Ψ - = - N k z z yz xz x dz G G k k k k k k 1 2 33 2 13 33 1 α (A.11.f) x x 12 21 α α = , x x 13 31 α α = , x x 23 32 α α = (A.11.g, h, i) ∑ ∫ =                 Ψ Ψ + Ψ Ψ - = - N k z z yz xz xy dz G G k k k k k k k k
∑ ∫ =                 Ψ Ψ + Ψ Ψ + Ψ Ψ + Ψ Ψ - = - N k z z yz xz xy dz G G k k k k k k k k k k k k 1 22
∑ ∫ =                 Ψ Ψ + Ψ Ψ + Ψ Ψ + Ψ Ψ - = - N k z z yz xz xy dz G G k k k k k k k k k k k k 1 23
∑ ∫ =                 Ψ Ψ + Ψ Ψ - = - N k z z yz xz xy dz G G k k k k k k k k 1 32 22 32 12 22 1 2 α (A.12.d) ∑ ∫ =                 Ψ Ψ + Ψ Ψ + Ψ Ψ + Ψ Ψ - = - N k z z yz xz xy dz G G k k k k k k k k k k k k 1 23
∑ ∫ =                 Ψ Ψ + Ψ Ψ - = - N k z z yz xz xy dz G G k k k k k k k k 1 33
∑ ∫ =                 Ψ + Ψ - = - N k z z yz xz y dz G G k k k k k k 1 2 21 2 31 11 1 α (A.13.a) ∑ ∫ =                 Ψ Ψ + Ψ Ψ - = - N k z z yz xz y dz G G k k k k k k k k 1 22 21 32 31 12 1 α (A.13.b) ∑ ∫ =                 Ψ Ψ + Ψ Ψ - = - N k z z yz xz y dz G G k k k k k k k k
∑ ∫ =                 Ψ Ψ + Ψ Ψ - = - N k z z yz xz y dz G G k k k k k k k k 1 23 22 33 32 23 1 α (A.13.e) ∑ ∫ =                 Ψ + Ψ - = - N k z z yz xz y dz G G k k k k k k 1 2
α α = (A.13.g, h, i) ∑ ∫ =                 Ψ Ψ + Ψ + Ψ Ψ + Ψ = - N k z z yz y x xz y x x dz G n n G n n k k k k k k k k k k
∑ ∫ =                 Ψ Ψ + Ψ Ψ + Ψ Ψ + Ψ Ψ = - N k z z yz y x xz y x x dz G n n G n n k k k k k k k k k k k k 1 21
∑ ∫ =                 Ψ Ψ + Ψ Ψ + Ψ Ψ + Ψ Ψ = - N k z z yz y x xz y x x dz G n n G n n k k k k k k k k k k k k

  and shear stresses at ply level k k k b a , = constants of integration of shear and normal stresses at ply level k C = differential operator defined in equation (13) in equation (105)

Figure 1 .

 1 Figure 1. Two-dimensional elastic structure

  introduced a second variable, namely the shear the study of isotropic beams, improving the classical model developed by Euler and

Figure 2 .

 2 Figure 2. Ply levels and interfaces numbering

  (or static) moment of the stiffness of the partial cross-section lying between the lower surface and the current z co-ordinate. The same result can be formally derived considering the first two Cauchy's equilibrium equations. Equations (1a) and (1b) are valid throughout the inner domain with the exception of the boundary. Integrating indefinitely both equations with respect to z we can re-write Eqs. (29a) and (29b) as follows: 0

Figure 3 .

 3 Figure 3. Boundary tangential (s) and normal (n) local co-ordinate systemFollowing Reissner's variational principle[START_REF] Reissner | The effect of transverse shear deformation on the bending of elastic plates[END_REF], and recalling transverse equilibrium (44)

  of transverse shear stressesThe field and boundary integrals that express the potential of transverse shear stresses are (see Appendix A for details)

  Fig.3):

  the final set of essential boundary conditions can be expressed in terms of the four primary variables M x , M y , M xy , and w: (see Appendix C for details) addition, the following boundary compatibility condition must be satisfied:

  b, c) and expanding (69), we obtain the following set of governing equations,

7 .

 7 Special one-dimensional cases. Closed form solutions a. Simply supported multi-layered panel subjected to transverse loadLet us consider a rectangular multi-layered panel as illustrated in Fig.4. The panel is simply supported along the two opposite edges x = 0 and x = l x . The remaining two edges y = 0 and y = l y are free.

Figure 4 .

 4 Figure 4. Multi-layered composite plate subjected to transverse loading

Figure 5 .

 5 Figure 5. Difference between proposed model and CLA.

Figure 6 .

 6 Figure 6. Deviation from CLA. Compared maximum transverse displacements

Figure 7 .

 7 Figure 7. Multi-layered composite plate subjected to axial in-plane load

  the transverse equilibrium equation (82a) in to the constitutive equation (82b): (82a) twice with respect to x and substituting the result in Eq. (83) the following buckling expression is derived:

Figure 8 .

 8 Figure 8. Deviation from CLA. Compared critical buckling loads

  assumption of symmetric lay-up and neglecting, as a first approximation, the terms D 16 and D 26 Eqs. (73) and (74) become:

  b) Substituting Eqs. (108) in Eq. (106):

  following relation between buckling loads and transverse displacements holds true also in the two-dimensional case: To validate the accuracy of the present model, solutions for a variety of examples are considered and compared to finite element analysis (FEA) and to other solutions available in the open literature. In particular, the refined Zigzag Theory recently developed by Tessler et al. (2009), the three-dimensional elasticity solution by Pagano (1969) and a first order shear deformation theory (FSDT) were compare to the present analysis. A first set of simulations was performed on sandwich square laminates defined in a domain identified by [ ] are subjected to a sinusoidal transverse pressure p defined by p=p 0 sin(πx/a)sin(πy/a). Examples include three variations of a three-layer square uniaxial sandwich M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT 34 laminate, indicated as A, B, and C respectively, having uniaxial carbon-epoxy external sheets and polyvinyl chloride (PVC) core (modeled as an isotropic material). Furthermore, a sandwich laminate made of titanium external sheets and titanium honeycomb core, indicated herein as laminate D, is simulated.

  estimates the stresses arising within the face sheets. The shear stress resultants, i.e. the integral over the thickness of the shear stress distribution, are the same for the three models, i.e. present, Tessler et al. (2009), and Pagano (1969). Therefore, the amount of shear carried by the face sheet predicted by the present theory is smaller than the one predicted by Pagano and Tessler et al. Consequently, the present M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT35 theory predicts a larger amount of deformation energy concentrated within the core, which is considerably less stiff than the face sheets. Since the face sheets represent 20% of the total laminate thickness, their influence on the global behaviour is significant, therefore predictions obtained by the present theory are to some extent more conservative (in terms of transverse displacements). Furthermore, the shear stresses distribution calculated by means of Pagano's exact solution suggests that state of stress and state of strain are fully three-dimensional within the face sheets. Therefore the hypothesis of state of plane stress on which the present theory is based, does not hold true. For a more accurate description of the real state of stress throughout the domain, the work performed by the through-the-thickness normal stress σ z should be included in the variational statement (49). This is beyond the scope of the present study, which focuses on the behaviour of reasonably thick composite laminates, thus neglecting the influence of σ z .Further comparisons in terms of normalized transverse shear stresses are presented in Figs.10 to 15.Notably, as the face sheet thickness decreases (configurations B and C), the present theory captures adequately the shear stress distribution in the proximity of the upper and lower surfaces of the laminates.

Fig. 9 .Fig. 10 .Fig. 11 .Fig. 12 .Fig. 13 .Fig. 14 .Fig. 15 .

 9101112131415 Fig.9. Normalized shear stresses for laminate A. Normalization factor: 2h/p 0 l x .

  case study, the plate was loaded with a uniform transverse pressure p 0 . Normalized transverse maximum displacements (normalization factor equal to 10 2 D 11 /p o l x 4 ) and non-dimensional through-thethickness normal stress (σ z /p 0 ) were calculated and results were compared to FEM. Solutions are calculated by linear superposition of basic solutions (99a) and (100a): where the amplitudes Ω mn and w mn are given by Eqs (99b) and (100b), and M and N are reasonably large integers Cosentino and Weaver (2009). Coefficients p mn , i.e. the components of the transverse pressure along the sinusoidal base functions, are calculated by means of a Galerkin technique Cosentino and Weaver (2008).

Fig. 16 .

 16 Fig.16. Comparison of maximum transverse displacements.

Fig. 17 .

 17 Fig.17. Comparison of through-the-thickness stresses -σ σ σ σ z (l x /2,l y /2,z) for configuration L1.

Fig. 18 .Fig. 19 .Fig. 20 .Fig. 21 .

 18192021 Fig.18. Comparison of through-the-thickness stresses -σ σ σ σ z (l x /2,l y /2,z)/p 0 for configuration L2.

  substituting the constitutive equation (A.3) and Eqs. (57) in Eq. (53b):

5 )

 5 For convenience, the two contributions Π xz and Π yz due to τ xz and τ yz respectively, are analyzed individually. Defining the following 6a, b) and substituting Eqs. (A.2) in Eqs. (A.6) we obtain:

  and conveniently re-arranging

  10.a, b) The matrices x α , y α , xy α , x β , and y β were introduced. Their components are hereby defined as

2 ) 5 )

 25 Referring to Eq. (53d), and performing the variations, we obtain:Projecting the shear resultants in the boundary integral of Eq. (B.2): example, Reissner [8], the following vertical equilibrium equations are introduced: (B.4) in the domain integral of Eq. (B) and carrying out the variations: Substituting Eq. (63) in Eqs. (61a), (61b), and (61c) gives 4.d, e, f) Substituting Eqs. (C.2) in Eq. (C.1.c), the following expression is obtained: is identically satisfied if exact solutions are found in terms of transverse bending moments.

  to solve the problem of stress concentration around holes in transversely bent panels. Recently, Shimpi et al. (2007) derived two novel formulations with only two variables, which provide fast and accurate results for moderately thick isotropic plates and includes classical laminate analysis as a degenerate case. However, it requires ad-hoc calculated shear correction factors for transverse shear stresses in multilayered composite panels. The governing field equations are derived from Castigliano's Theorem of Least Work. A variety of conditions can be assigned to kinematical and mechanical variables contemporarily. As thoroughly summarized by

Table 1

 1

	. Material properties

32

  Transverse shears τ xz and τ yz , and transverse normal stress σ z are derived by substituting Eqs. (99) in Eqs.

	(38) and (39) respectively:
	( ) mn k xz τ	=	-	α	m 4	+	p m mn β α 2 2	n 2	+	β	n 4	[ { α	m 3	k 11 Ψ	+	α	m	β	n 2	( Ψ 12	k	+	33 Ψ	k	) ] ( ) ( ) y x n m β α sin cos
									+	[ β	n 3	32 Ψ	k	+	α	m 2	β	n	( Ψ 13	k	+	k 31 Ψ	) ] ( ) ( )} y x n m β α cos sin	(102)
	( ) mn k yz τ	=	-	α	m 4	+	p m mn β α 2 2	n 2	+	β	n 4	[ { α	m 3	k 31 Ψ	+	α	m	β	n 2	( Ψ	23	k	+	32 Ψ	k	) ] ( ) ( ) y x n m β α sin cos
									+	[ β	n 3	22 Ψ	k	+	α	m 2	β	n	( Ψ 21 k	+	33 Ψ	k	) ] ( ) ( )} y x n m β α cos sin	(103)
	( ) mn k z σ	α	m 4	p m mn β α 2 2	n 2		β	n 4	[ { α	m 4		k 31			α	m 2	β	n 2	(	12	k	33	k	k 31	)	β	n 4	22	k	] ( ) ( ) y x n m β α sin sin

Table 1 . Material properties

 1

		Orthotropic materials				Isotropic materials
	Properties Carbon-epoxy (C) Titanium core (H) PVC core (P) Titanium lamina (T)
	E 1 , GPa	1.579×10 2	1.915×10 -1
	E 2 , GPa	9.584	1.915×10 -1	1.040×10 -1	1.041×10 2
	E 3 , GPa	9.584	1.915
	G 12 , GPa	5.930	4.227×10 -8
	G 13 , GPa	5.930	5.651×10 -1	G	=	2	( ) ν + 1 E	G	=	2	( ) ν + 1 E
	G 23 , GPa	3.227	1.248
	ν 12	0.32	0.658×10 -2
	ν 13	0.32	0.643×10 -6			0.3			0.31
	ν 23	0.49	0.643×10 -6

Table 2 . Configurations and stacking sequences

 2

	Laminate	Normalized ply	Laminate	Stacking
	configuration	thicknesses (h k /2h)	materials	sequence
	A	(0.1/0.8/0.1)	(C/P/C)	(0/0/0)
	B	(0.025/0.95/0.025)	(C/P/C)	(0/0/0)
	C	(0.0025/0.995/0.0025)	(C/P/C)	(0/0/0)
	D	(0.1/0.8/0.1)	(T/H/T)	(0/0/0)

Table 3 . Normalized maximum (central) deflections ( Normalization factor: 10 2 D 11 /p o l x 4 )

 3

	Laminate	3-D Elasticity	FSDT	Zigzag	Present
	configuration	(Pagano)		(Tessler et al.)
	A	29.761	2.731	29.769	35.665
	B	11.645	2.819	11.693	11.303
	C	2.080	1.728	2.103	2.001
	D	1.331	0.389	1.333	1.672

Table 4 . Normalized maximum (central) deflections ( Normalization factor: 10 2 D 11 /p o l x 4 )

 4 To further validate the present theory, a variety of solutions were derived for classical multilayered 8-layer composite laminates. A simply supported square plate was studied (l x =5mm). Five different quasi-isotropic stacking sequences were analyzed. Lamina properties are the same as the properties reported in Table1for the carbon-epoxy unidirectional lamina. Configurations are indicated in Table5.

	Laminate	Span-to-	3-D Elasticity	FSDT	Zigzag	Present
	configuration	thickness ratio	(Pagano)		(Tessler et al.)
	A	4	42.42	3.739	42.124	55.23
	A	10	9.734	1.321	9.738	9.981
	A	20	3.487	0.948	3.489	3.664
	A	50	1.303	0.841	1.305	1.367
	A	100	0.945	0.826	0.945	0.973
	A	200	0.852	0.822	0.852	0.882

Table 5 . Laminate configurations

 5

	Laminate	Stacking
	configuration	sequence
	L1	[0/90/45/-45]s
	L2	[0/45/-45/90]s
	L3	[90/45/-45/0]s
	L4	[0/90] 2 s
	L5	[45/-45] 2 s

Table 6 . Normalized buckling loads. Normalization factor: l x 2 /D 11

 6

	Laminate configuration	Span-to-thickness ratio, lx/2h	CLA	Present	FEM
		500	17.6	17.5	17.6
		50	17.6	17.3	17.0
	L1	20	17.6	17.1	15.1
		10	17.6	16.2	14.4
		5	17.6	13.4	10.1
		500	18.5	18.5	18.5
		50	18.5	18.3	17.7
	L2	20	18.5	18.1	16.7
		10	18.5	17.3	14.1
		5	18.5	14.8	11.4
		500	69.2	68.9	69.2
		50	69.2	68.5	68.2
	L3	20	69.2	67.8	66.2
		10	69.2	64.0	57.9
		5	69.2	52.2	43.7
	L4	500	18.4	18.3	18.4
		50	18.4	18.4	17.4

Table 7 . Comparison of maximum transverse displacements for panels under uniform pressure.

 7

	Maximum transverse displacement, µm

Results are very accurate in terms of through-the-thickness normal stresses. The calculation of transverse displacements are less accurate but still acceptable if considering that the case of a very thick laminate is