An accurate zigzag theory for bending and buckling analysis of thick laminated sandwich plates with soft core

Abstract

An accurate and computationally attractive zigzag theory is developed for bending and buckling analysis of thick laminated soft core sandwich plates. The kinematic assumptions of the proposed zigzag theory are obtained by superimposing a nonlinear zigzag function on the first-order shear deformation theory. In order to obtain the accurate transverse shear stresses, a preprocessing approach based on the three-dimensional equilibrium equations and the Reissner mixed variational theorem is used. It is significant that the second-order derivatives of in-plane displacement variables have been removed from the transverse shear stresses, such that the finite element implementation is greatly simplified. Thus, based on the proposed zigzag model, a computationally efficient four-node C⁰ quadrilateral plate element with linear interpolation function is proposed for bending and buckling analysis of soft core sandwich plates. The advantage of the present formulation is that no post-processing approach is needed to calculate the transverse shear stresses while maintaining the computational accuracy of a linear plate element. Moreover, the accurate transverse shear stresses can be involved in the strain energy which can actively improve the accuracy of critical loads. Performance of the proposed model is assessed by comparing with several benchmark solutions. Agreement between the present results and the reference solutions is very good, and the proposed model only includes the seven displacement variables which can demonstrate the accuracy and effectiveness of the proposed model.

Keywords

Soft core sandwich plates four-node quadrilateral element mixed zigzag model Reissner mixed variational theorem

Introduction

Laminated composite structures have been increasingly used in weight-sensitive engineering fields (civil, automotive, aerospace, etc.) due to their high stiffness and strength-to-weight ratio. Sandwich structures are constituted by soft cores and stiff face sheets in the form of composite laminates. Owing to the material characteristics which are quite different between the core and face layers, the effect of shear deformation is quite significant and it becomes more complex in sandwich structures. For reliable and safe designs, it is necessary to well understand the structural behaviors of sandwich structures. Thus, adequate analysis models must be established to accurately predict the responses of sandwich structures.

Since the transverse variation of the material characteristics is very large in the core and the face sheet and the thickness of sandwich plates is often greater than that of laminated composite plates, the traditional global theories (equivalent single layer theories) such as the classical lamination theory (CLT),¹ the first-order shear deformation theory (FSDT),² and the higher order shear deformation theories (HSDT)^3–6 may not provide accurate results in the cases wherein the core is made up of very soft materials or the number of the layers is large. To overcome these drawbacks, layerwise theories^7–12 have been proposed to model the local distributions of the displacement components. Layerwise theories can produce accurate displacements and stresses, while huge computational efforts are required as the number of the displacement variables grows with the number of the layers. In order to retain the merits of equivalent single layer and layerwise models, a number of zigzag models have been proposed. By introducing a linear zigzag function to the first-order shear deformation model, Murakami¹³ developed a zigzag model that can model the zig-zag distribution of in-plane displacements through the thickness direction. The performance of the Murakami zigzag function (MZZF)¹³ in 2D plate and shell models has been evaluated in detail.¹⁴ Additionally, Carrera¹⁴ indicated that the higher order model including the MZZF is more efficient than the introduction of higher order polynomials. Nevertheless, such zig-zag models^13,14 violate the continuity conditions of transverse shear stresses at interfaces. In view of this situation, the higher order zig-zag models as well as the higher order global-local models with interlaminar continuous transverse stresses were proposed. In view of this situation, Di Sciuva¹⁵ developed an equivalent single-layer zigzag model with continuous transverse shear stresses. Cho and Parmerter^16,17 proposed an efficient zigzag theory satisfying the continuity of transverse shear stresses at the interfaces between layers as well as the shear-free conditions at the surfaces. Although the continuity conditions of transverse shear stresses have been satisfied, the transverse shear stresses computed directly from the constitutive relations are still less accurate.^15–17 Using the double superposition hypothesis approach, Li and Liu¹⁸ proposed a third-order global-local model satisfying the continuity conditions of transverse shear stresses at interfaces that can accurately calculate the transverse shear stresses of laminated structures directly from constitutive equations. Based on the third-order global-local model, Wu et al.¹⁹ and Wu and Chen²⁰ developed the refined quadrilateral and triangular plate elements to investigate the bending and buckling behaviors of laminated plates. Based on a fifth-order global-local model, Wu and Ren²¹ developed a refined triangular plate element to investigate hygrothermal behaviors of angel-ply composite plates. Recently, Wu et al.^22,23 suggested a new global higher order zigzag model derived from the HW variational theorem for interlaminar stresses analysis of multilayered composite structures. By accounting for transverse normal deformation, Nath and Mishra²⁴ developed an enhanced zigzag theory to study bending behaviors of composite plates. Reviewing the above-mentioned literature,^15–24 it is found that the first-order derivatives of transverse displacement are included in the displacement fields of the zigzag and global-local models. Therefore, C¹-continuous interpolation functions are required to approximate the transverse displacement due to the presence of its second-order derivatives in the strain energy. Aiming to improve this situation, Wu and Chen²⁵ developed a C⁰-type global-local higher order theory. One major benefit arising from the analytical form of the C⁰-type global-local model²⁵ is its ideal suitability to finite element modeling, where the kinematic-variable approximations need not exceed C⁰ continuity. The implication of this attribute is that computationally efficient finite elements can be developed and used in general-purpose finite element codes. Based on the global-local theory,²⁵ Wu et al.^26–28 proposed a set of C⁰ plate element for bending, buckling and vibration analysis of laminated composite and sandwich plates. By considering transverse normal strain, an improved C⁰-type global-local theory²⁹ was proposed for hygrothermal analysis of laminated composite plates. Shariyat et al.³⁰ suggested a global-local theory with a new post-processing technique for stress analysis of asymmetric orthotropic sandwich plates with single/dual cores. Ren et al.³¹ developed a zigzag model in which the first-order derivatives of transverse displacement have been taken out from the displacement fields. In terms of a C⁰ inverse trigonometric zigzag theory, Sahoo and Sinbh³² investigated the buckling behaviors of laminated composite and sandwich plates. By accounting for transverse normal deformation, an improved C⁰ zigzag plate theory^33,34 was developed for bending and buckling analysis of soft core sandwich plates. However, for the multilayered composite and sandwich plates, the stress smoothing technique of integrating 3D equilibrium equations must be used to obtain accurate transverse shear stresses in these models.^25–34 This approach means that the first-order derivatives of in-plane stresses will be involved in the 3D equilibrium equations. To accurately calculate transverse shear stresses by means of the equilibrium equation stress smoothing technique, the higher order elements are usually developed for the analysis of multilayered composite and sandwich plates.^{26,30–31,33} Nevertheless, such elements contain the second-order derivatives of in-plane displacements as nodal degrees of freedom, which will decrease the computational efficiency.

By assuming two independent fields for transverse stresses and displacement variables, Murakami³⁵ was the first to apply the RMVT³⁶ to laminated structures. In review article,³⁷ Carrera indicated that the RMVT remains a valuable tool to analyze multilayered plates. In terms of the RMVT, Messina³⁸ developed a mixed two-dimensional model for free vibration of laminated plates. Kim and Cho³⁹ proposed an enhanced first-order theory based on mixed formulation for static analysis of laminated and sandwich plates. Lezgy-Nazargah⁴⁰ developed a high-performance finite element model in conjunction with the parametrized mixed variational principle for bending and vibration analysis of thick plates. Recently, an efficient mixed finite element model⁴¹ is proposed for static and free vibration analysis of FGM plates rested on elastic foundations. Lezgy-Nazargah and Salahshuran⁴² investigated the bending and vibration behaviors of multilayered composite plate based on a new mixed-field global-local theory. Tessler et al.^43,44 developed a mixed-form refined zigzag theory for the linearly elastic analysis of composite beams and plates. Based on the mixed-form zigzag model,⁴³ Groh and Tessler⁴⁵ proposed the three-node and two-node shear locking-free beam elements to analyze sandwich beam and laminated composite with delaminations. However, Wu and Chen²¹ indicated that the accuracy of the mixed-form zigzag model⁴⁴ decreases as the number of layers is increased, so the mixed-form zigzag model is inadequate for the analysis of multilayered composite structures.

In this paper, an accurate zigzag theory for thick laminated soft core sandwich plates is proposed. In the proposed theory, the displacement field is assumed to be a combination of the first-order shear deformation theory and a nonlinear zigzag function. Additionally, by applying 3D elasticity equations, a set of accurate transverse shear stresses is obtained. With the aid of the strain-compatibility variational statement of the RMVT, the higher-order derivatives of in-plane displacement variables have been taken out from the transverse shear stress fields. The major benefit of this mixed-field formulation is that no post-processing approach is required to accurately determine the transverse shear stresses, while the low-order plate element can be easily developed. Moreover, the accurate transverse shear stresses can be included in the strain energy which can actively improve the accuracy of critical loads. Compared to the higher order elements, the C⁰-type four-node quadrilateral element is computationally more efficient and economical to implement in general-purpose finite element codes. Thus, a simple four-node C⁰ quadrilateral plate element is developed to accurately predict displacements, stresses and critical loads of laminated soft core sandwich plates. Performance of the proposed zigzag model is numerically assessed by analyzing laminated soft core sandwich plates with different lamination sequences, material characteristics, boundary conditions and span-to-thickness ratios. By comparing with the 3D elasticity solutions and layerwise solutions, the numerical results shows that the proposed zigzag model can accurately predict the transverse shear stresses and the critical loads for thick and highly heterogeneous sandwich plates.

Kinematic assumptions

Displacement field

Consider a plate of thickness h consisting of N perfectly bonded orthotropic plies, as shown in Figure 1. The plate is studied in the Cartesian coordinate system (x, y, z) where the thickness coordinate z ranges from −h/2 to h/2. By introducing the nonlinear zigzag functions into the first-order displacement model, the displacement field of the zigzag model in the present work can be written as

\begin{matrix} u^{k} (x, y, z) = u_{0} (x, y) + {zu}_{1} (x, y) + Φ^{k} (z) u_{z} (x, y) \\ v^{k} (x, y, z) = v_{0} (x, y) + {zv}_{1} (x, y) + Ψ^{k} (z) v_{z} (x, y) \\ w^{k} (x, y, z) = w_{0} (x, y) \end{matrix}

(1)

where the superscript k denotes the k-th ply of the laminated plate, where u₀ and v₀ are the in-plane uniform displacements, u₁ and v₁ are the bending rotations, w₀ is the transverse deflection and u_z and v_z are the amplitude functions.

The zigzag functions $Φ^{k} (z)$ and $ψ^{k} (z)$ can be explicitly written as

\begin{matrix} Φ^{k} (z) = φ^{k} (z) + z^{2} φ^{k} (z) \\ Ψ^{k} (z) = ψ^{k} (z) + z^{2} ψ^{k} (z) \end{matrix}

(2)

where

φ^{k} (z)

ψ^{k} (z)

are the linear zigzag functions proposed by Tessler et al.,⁴⁶ which can be expressed as

\begin{matrix} \begin{matrix} \begin{matrix} φ^{k} (z) = \frac{1}{2} (1 - ζ_{k}) {\bar{u}}^{k} + \frac{1}{2} (1 + ζ_{k}) {\bar{u}}^{k + 1} \\ ψ^{k} (z) = \frac{1}{2} (1 - ζ_{k}) {\bar{v}}^{k} + \frac{1}{2} (1 + ζ_{k}) {\bar{v}}^{k + 1} \end{matrix} & (k = 1, 2, 3, \dots, N) \end{matrix} \end{matrix}

(3)

where

{\bar{u}}^{k}, {\bar{v}}^{k} (k = 1, 2, 3, \dots, N)

are the interfacial axial displacements.

ζ_{k} \in [- 1, 1]

ζ_{k} = a_{k} z - b_{k}

a_{k} = 2 / (z_{k + 1} - z_{k})

b_{k} = (z_{k + 1} + z_{k}) / (z_{k + 1} - z_{k})

, and z_k is defined in Figure 1. The origin of the coordinates (

z = 0

) is taken at the mid-plane of the sandwich plate.

Consistent with the kinematic assumptions (equation (1)) and by using the linear strain-displacement relations, the transverse shear strains read as

\begin{matrix} γ_{xz}^{k} = \frac{\partial w_{0}}{\partial x} + u_{1} + \frac{\partial Φ^{k} (z)}{\partial z} u_{z} \\ γ_{yz}^{k} = \frac{\partial w_{0}}{\partial y} + v_{1} + \frac{\partial Ψ^{k} (z)}{\partial z} v_{z} \end{matrix}

(4)

Following the RZT strategy,⁴⁶ the strain measures are introduced, which are defined as

\begin{matrix} η_{x} = \frac{\partial w_{0}}{\partial x} + u_{1} - u_{z} \\ η_{y} = \frac{\partial w_{0}}{\partial y} + v_{1} - v_{z} \end{matrix}

(5)

The transverse shear strains can be rewritten as

\begin{matrix} γ_{xz}^{k} = η_{x} + [1 + \frac{\partial Φ^{k} (z)}{\partial z}] u_{z} \\ γ_{yz}^{k} = η_{y} + [1 + \frac{\partial Ψ^{k} (z)}{\partial z}] v_{z} \end{matrix}

(6)

The transverse shear stresses for the kth ply can be given by

\begin{matrix} τ_{xz}^{k} = Q_{44}^{k} η_{x} + Q_{44}^{k} [1 + \frac{\partial Φ^{k} (z)}{\partial z}] u_{z} = Q_{44}^{k} η_{x} + {\bar{τ}}_{xz}^{k} \\ τ_{yz}^{k} = Q_{55}^{k} η_{y} + Q_{55}^{k} [1 + \frac{\partial Ψ^{k} (z)}{\partial z}] v_{z} = Q_{55}^{k} η_{y} + {\bar{τ}}_{yz}^{k} \end{matrix}

(7)

It can be seen from equation (7) that the transverse shear stresses are composed by two contributions: the continuity conditions are enforced only on the zigzag-dependent contribution, ${\bar{τ}}_{xz}^{k}$ and ${\bar{τ}}_{yz}^{k}$ , that is, on the transverse shear stresses obtained by vanishing the shear measures $η_{x}$ and $η_{y}$ . The constraint reads as

\begin{matrix} {\bar{τ}}_{xz}^{k} (x, y, z_{k + 1}) = {\bar{τ}}_{xz}^{k + 1} (x, y, z_{k + 1}) \\ {\bar{τ}}_{yz}^{k} (x, y, z_{k + 1}) = {\bar{τ}}_{yz}^{k + 1} (x, y, z_{k + 1}) \end{matrix}

(8)

Equation (8) provides 2(N−1) conditions for the interfacial axial displacements. The other four conditions derive from the zero-value constraint on the outer plate surfaces, that is

\begin{matrix} {\bar{τ}}_{xz}^{k} (x, y, \pm h / 2) = 0 \\ {\bar{τ}}_{yz}^{k} (x, y, \pm h / 2) = 0 \end{matrix}

(9)

Finally, equation (8) along with equation (9) constitute a set of 2(N + 1) conditions that allow to obtain the interfacial axial displacements in equation (3).

For linear elasticity, the strain of the present model can be written as

{\begin{matrix} ɛ_{x} \\ ɛ_{y} \\ γ_{xy} \\ γ_{xz} \\ γ_{yz} \end{matrix}}^{k} = [\begin{matrix} B_{IX} & B_{IY} & B_{IXY} & B_{TXZ} & B_{TYZ} \end{matrix}] {\begin{matrix} E_{IX} \\ E_{IY} \\ E_{IXY} \\ E_{TXZ} \\ E_{TYZ} \end{matrix}}

(10)

where

\begin{matrix} B_{IX} = {\begin{matrix} 1 & z & Φ^{k} \end{matrix}} \\ B_{IY} = {\begin{matrix} 1 & z & Ψ^{k} \end{matrix}} \\ B_{IXY} = {\begin{matrix} 1 & z & Ψ^{k} & 1 & z & Φ^{k} \end{matrix}} \\ B_{TXZ} = {\begin{matrix} 1 & \frac{\partial Φ^{k}}{\partial z} & 1 \end{matrix}} \\ B_{TYZ} = {\begin{matrix} 1 & \frac{\partial Ψ^{k}}{\partial z} & 1 \end{matrix}} \end{matrix}

(11)

\begin{matrix} E_{IX} = {\begin{matrix} \frac{\partial u_{0}}{\partial x} & \frac{\partial u_{1}}{\partial x} & \frac{\partial u_{z}}{\partial x} \end{matrix}}^{T} \\ E_{IY} = {\begin{matrix} \frac{\partial v_{0}}{\partial y} & \frac{\partial v_{1}}{\partial y} & \frac{\partial v_{z}}{\partial y} \end{matrix}}^{T} \\ E_{IXY} = {\begin{matrix} \frac{\partial v_{0}}{\partial x} & \frac{\partial v_{1}}{\partial x} & \frac{\partial v_{z}}{\partial x} & \frac{\partial u_{0}}{\partial y} & \frac{\partial u_{1}}{\partial y} & \frac{\partial u_{z}}{\partial y} \end{matrix}}^{T} \\ E_{TXZ} = {\begin{matrix} u_{1} & u_{z} & \frac{\partial w_{0}}{\partial x} \end{matrix}}^{T} \\ E_{TYZ} = {\begin{matrix} v_{1} & v_{z} & \frac{\partial w_{0}}{\partial y} \end{matrix}}^{T} \end{matrix}

(12)

The constitutive equations for the kth lamina of cross-ply composite plate can be expressed as

{\begin{matrix} σ_{x} \\ σ_{y} \\ τ_{xy} \\ τ_{xz} \\ τ_{yz} \end{matrix}}^{k} = [\begin{matrix} Q_{11} & Q_{12} & 0 & 0 & 0 \\ Q_{21} & Q_{22} & 0 & 0 & 0 \\ 0 & 0 & Q_{33} & 0 & 0 \\ 0 & 0 & 0 & Q_{44} & 0 \\ 0 & 0 & 0 & 0 & Q_{55} \end{matrix}]^{k} {\begin{matrix} ɛ_{x} \\ ɛ_{y} \\ γ_{xy} \\ γ_{xz} \\ γ_{yz} \end{matrix}}^{k}

(13)

where

Q_{ij}^{k}

is the kth ply transformed elastic stiffness coefficients referring to the (x, y, z) coordinate system.

Improved transverse shear stresses

The inherent drawback of the proposed model is that the transverse shear stresses obtained from the constitutive equations are too poor to be useful. The 3D equilibrium equations are commonly used in an attempt to derive improved layer interface-continuous transverse shear stresses. Moreover, Wu and Chen⁴⁷ showed that the accurate prediction of critical buckling loads of laminated composite structures needs the models which can describe the interlaminar continuous transverse stresses. Thus, the 3D equilibrium equations will be employed to derive the improved transverse shear stresses of sandwich plates.

Discarding the body forces, the 3D equilibrium equations can be given by

\begin{matrix} \frac{\partial σ_{x}^{k}}{\partial x} + \frac{\partial τ_{xy}^{k}}{\partial y} + \frac{\partial τ_{xz}^{k}}{\partial z} = 0 \\ \frac{\partial σ_{y}^{k}}{\partial y} + \frac{\partial τ_{xy}^{k}}{\partial x} + \frac{\partial τ_{yz}^{k}}{\partial z} = 0 \end{matrix}

(14)

By integrating equation (14) across the z-coordinate direction and enforcing the transverse shear traction conditions at the bottom surface, the transverse shear stresses can be expressed as

\begin{matrix} {\tilde{τ}}_{xz}^{k} (x, y, z) = - \int_{- h / 2}^{z} [\frac{\partial σ_{x}^{k}}{\partial x} + \frac{\partial τ_{xy}^{k}}{\partial y}] d z \\ {\tilde{τ}}_{yz}^{k} (x, y, z) = - \int_{- h / 2}^{z} [\frac{\partial σ_{y}^{k}}{\partial y} + \frac{\partial τ_{xy}^{k}}{\partial x}] d z \end{matrix}

(15)

Substituting equations (10) to (13) into equation (15) may result in a number of second-order partial derivatives of displacement parameters. For the cross-ply laminated sandwich plates, making a compromise between accuracy and efficiency, a simplified expression of transverse shear stresses proposed by Iurlaro et al.⁴⁴ will be adopted in equation (15). The simplified transverse shear stresses are given by

\begin{matrix} {\tilde{τ}}_{xz}^{k} (x, y, z) = - C_{x} (z) U_{x} \\ {\tilde{τ}}_{yz}^{k} (x, y, z) = - C_{y} (z) U_{y} \end{matrix}

(16)

where

\begin{matrix} C_{x} (z) = \int_{- h / 2}^{z} Q_{11}^{k} [\begin{matrix} 1 & z & Φ^{k} \end{matrix}] d z \\ C_{y} (z) = \int_{- h / 2}^{z} Q_{22}^{k} [\begin{matrix} 1 & z & Ψ^{k} \end{matrix}] d z \end{matrix}

(17)

\begin{matrix} U_{x} = \frac{\partial^{2}}{\partial x^{2}} [\begin{matrix} u_{0} & u_{1} & u_{z} \end{matrix}]^{T} \\ U_{y} = \frac{\partial^{2}}{\partial y^{2}} [\begin{matrix} v_{0} & v_{1} & v_{z} \end{matrix}]^{T} \end{matrix}

(18)

By employing the top surface (z = h/2) transverse shear free conditions, the second-order partial derivatives of displacement variables u₀ and v₀ can be eliminated from equation (18). The transverse shear stresses can be rewritten as

\begin{matrix} {\tilde{τ}}_{xz}^{k} (x, y, z) = E_{x} (z) {\bar{U}}_{x} \\ {\tilde{τ}}_{yz}^{k} (x, y, z) = E_{y} (z) {\bar{U}}_{y} \end{matrix}

(19)

where

\begin{matrix} {\bar{U}}_{x} = \frac{\partial^{2}}{\partial x^{2}} [\begin{matrix} u_{1} & u_{z} \end{matrix}]^{T} \\ {\bar{U}}_{y} = \frac{\partial^{2}}{\partial y^{2}} [\begin{matrix} v_{1} & v_{z} \end{matrix}]^{T} \end{matrix}

(20)

\begin{matrix} E_{x} (z) = \int_{- h / 2}^{z} Q_{11}^{k} [\begin{matrix} A_{1} & A_{2} \end{matrix}] dz \\ E_{y} (z) = \int_{- h / 2}^{z} Q_{22}^{k} [\begin{matrix} B_{1} & B_{2} \end{matrix}] dz \end{matrix}

(21)

in which

\begin{matrix} A_{1} = \frac{\int_{- h / 2}^{h / 2} Q_{11}^{k} z dz}{\int_{- h / 2}^{h / 2} Q_{11}^{k} dz} - z, A_{2} = \frac{\int_{- h / 2}^{h / 2} Q_{11}^{k} Φ^{k} dz}{\int_{- h / 2}^{h / 2} Q_{11}^{k} dz} - Φ^{k} \\ B_{1} = \frac{\int_{- h / 2}^{h / 2} Q_{22}^{k} z dz}{\int_{- h / 2}^{h / 2} Q_{22}^{k} dz} - z, B_{2} = \frac{\int_{- h / 2}^{h / 2} Q_{22}^{k} Ψ^{k} dz}{\int_{- h / 2}^{h / 2} Q_{22}^{k} dz} - Ψ^{k} \end{matrix}

(22)

By means of the nodal degrees of freedom and interpolation functions of the four-node quadrilateral element, the transverse shear stresses of equation (19) within an element can be expressed as

{\begin{matrix} {\tilde{τ}}_{xz}^{k} \\ {\tilde{τ}}_{yz}^{k} \end{matrix}} = {\bar{B} δ}^{e}

(23)

where

\bar{B} = [\begin{matrix} {\bar{B}}_{1} & {\bar{B}}_{2} & {\bar{B}}_{3} & {\bar{B}}_{4} \end{matrix}]

(24)

{\begin{matrix} \bar{B} \end{matrix}}_{i} = [\begin{matrix} 0 & 0 & 0 & E_{x}^{(1)} (z) \frac{\partial^{2} N_{i}}{\partial x^{2}} & E_{x}^{(2)} (z) \frac{\partial^{2} N_{i}}{\partial x^{2}} & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & E_{y}^{(1)} (z) \frac{\partial^{2} N_{i}}{\partial y^{2}} & E_{y}^{(2)} (z) \frac{\partial^{2} N_{i}}{\partial y^{2}} \end{matrix}]

(25)

where N_i is the shape function of the four-node quadrilateral element,

E_{x}^{(m)} (z)

is the mth component of the

E_{x} (z)

;

E_{y}^{(n)} (z)

is the nth component of the

E_{y} (z)

Transverse shear stresses ${\tilde{τ}}_{xz}^{k}$ and ${\tilde{τ}}_{yz}^{k}$ in equation (19) can a priori satisfy the transverse shear stress continuity conditions as well as the zero shear stress conditions on outer plate surfaces. In order to obtain the improved results, the transverse shear stresses ${\tilde{τ}}_{xz}^{k}$ and ${\tilde{τ}}_{yz}^{k}$ derived from 3D equilibrium equations are a natural choice. However, it is observed from equation (25) that matrix ${\bar{B}}_{i}$ is zero as the shape functions N_i are linear functions within one element, so that the transverse shear stresses obtained by integrating 3D equilibrium equations are also zero. For the four-node quadrilateral element, the post-processing technique based on 3D equilibrium equations is not adaptive to predict the transverse shear stresses. The motivation of the present work is to produce accurate transverse shear stresses by using the four-node quadrilateral element without any post-processing procedures, so the RMVT³⁶ is employed to eliminate the second-order derivatives of displacement variables in equation (19).

Employing the assumed displacement field given by equation (1), and the transverse shear stresses given by equation (19), the Reissner mixed variational theorem can be stated as

\begin{matrix} δ Π_{R} = \int_{S_{m}} \int_{- h / 2}^{h / 2} [σ_{x}^{k} δ ɛ_{x}^{k} + σ_{y}^{k} δ ɛ_{y}^{k} + τ_{xy}^{k} δ γ_{xy}^{k} + {\tilde{τ}}_{xz}^{k} δ γ_{xz}^{k} \\ + {\tilde{τ}}_{yz}^{k} δ γ_{yz}^{k} + δ {\tilde{τ}}_{xz}^{k} (γ_{xz}^{k} - {\tilde{γ}}_{xz}^{k}) + δ {\tilde{τ}}_{yz}^{k} (γ_{yz}^{k} - {\tilde{γ}}_{yz}^{k})] \\ d z d x d y - δ W_{e} = 0 \end{matrix}

(26)

where W_e denotes the work of external forces. The

{\tilde{τ}}_{xz}^{k}

and

{\tilde{τ}}_{yz}^{k}

are the assumed transverse shear stresses given by equation (24), and

{\tilde{γ}}_{xz}^{k}

and

{\tilde{γ}}_{yz}^{k}

are the transverse shear strains that is conjugate to

{\tilde{τ}}_{xz}^{k}

and

{\tilde{τ}}_{yz}^{k}

{\tilde{γ}}_{xz}^{k} = {\tilde{τ}}_{xz}^{k} / Q_{44}^{k}

{\tilde{γ}}_{yz}^{k} = {\tilde{τ}}_{yz}^{k} / Q_{55}^{k}

. Since the assumed transverse shear stresses

{\tilde{τ}}_{xz}^{k}

and

{\tilde{τ}}_{yz}^{k}

do not invoke the displacement variables, the variational principle given by equation (26) can be implemented in a two-step procedure

\begin{matrix} \int_{S_{m}} \int_{- h / 2}^{h / 2} (σ_{x}^{k} δ ɛ_{x}^{k} + σ_{y}^{k} δ ɛ_{y}^{k} + τ_{xy}^{k} δ γ_{xy}^{k} + {\tilde{τ}}_{xz}^{k} δ γ_{xz}^{k} + {\tilde{τ}}_{yz}^{k} δ γ_{yz}^{k}) \\ \times d z d x d y - δ W_{e} = 0 \end{matrix}

(27)

\int_{S_{m}} \int_{- h / 2}^{h / 2} (δ {\tilde{τ}}_{xz}^{k} (γ_{xz}^{k} - {\tilde{γ}}_{xz}^{k}) + δ {\tilde{τ}}_{yz}^{k} (γ_{yz}^{k} - {\tilde{γ}}_{yz}^{k})) d z d x d y = 0

(28)

Substituting equations (10) and (19) into equation (28), the following equations can be obtained

\begin{matrix} \int_{- h / 2}^{h / 2} δ {\bar{U}}_{x}^{T} E_{x}^{T} (B_{TXZ} E_{TXZ} - E_{x} {\bar{U}}_{x} / Q_{44}^{k}) dz = 0 \\ \int_{- h / 2}^{h / 2} δ {\bar{U}}_{y}^{T} E_{y}^{T} (B_{TYZ} E_{TYZ} - E_{y} {\bar{U}}_{y} / Q_{55}^{k}) dz = 0 \end{matrix}

(29)

Therefore, the higher order derivatives of displacement variables can be rewritten as

{\bar{U}}_{x} = C_{EXZ} E_{TXZ}, {\bar{U}}_{y} = C_{EYZ} E_{TYZ}

(30)

where

\begin{matrix} C_{EXZ} = [\int_{- h / 2}^{h / 2} (E_{x}^{T} E_{x} / Q_{44}^{k}) dz]^{- 1} [\int_{- h / 2}^{h / 2} E_{x}^{T} B_{TXZ} dz] \\ C_{EYZ} = [\int_{- h / 2}^{h / 2} (E_{y}^{T} E_{y} / Q_{55}^{k}) dz]^{- 1} [\int_{- h / 2}^{h / 2} E_{y}^{T} B_{TYZ} dz] \end{matrix}

(31)

Substituting equation (30) into equation (19), the improved transverse shear stresses can finally be expressed as

\begin{matrix} {\tilde{τ}}_{xz}^{k} (x, y, z) = E_{x} C_{EXZ} E_{TXZ} = Θ_{1}^{k} u_{1} + Θ_{2}^{k} u_{z} + Θ_{3}^{k} \frac{\partial w_{0}}{\partial x} \\ {\tilde{τ}}_{yz}^{k} (x, y, z) = E_{y} C_{EYZ} E_{TYZ} = ϒ_{1}^{k} v_{1} + ϒ_{2}^{k} v_{z} + ϒ_{3}^{k} \frac{\partial w_{0}}{\partial y} \end{matrix}

(32)

It can be seen from equation (32) that the second-order derivatives of in-plane displacements have been replaced by parameters u₁,u_z, $\partial w_{0} / \partial x$ , v₁, v_z and $\partial w_{0} / \partial y$ .

Finite element formulation

A C⁰-type four-node quadrilateral element

In the present work, a four-node quadrilateral isoperimetric element requiring only C⁰-continuity is presented. By means of nodal variables, the displacement parameters within one element can be expressed as follows

\begin{matrix} u_{l} = \sum_{i = 1}^{4} N_{i} u_{li}, v_{l} = \sum_{i = 1}^{4} N_{i} v_{li} \\ u_{z} = \sum_{i = 1}^{4} N_{i} u_{zi}, v_{z} = \sum_{i = 1}^{4} N_{i} v_{zi} \\ w_{0} = \sum_{i = 1}^{4} N_{i} w_{0 i} \end{matrix}

(33)

where

l = 0, 1

;

N_{i} = \frac{1}{4} (1 + ξ ξ_{i}) (1 + η η_{i}), i = 1 \sim 4

Based on equations (10) to (12), the strain vector can be expressed as

^{k} = \partial u^{k} = {B δ}^{e}

(34)

where

u^{k} = [\begin{matrix} u & v & w \end{matrix}]^{T}

B = [\begin{matrix} B_{1} & B_{2} & B_{3} & B_{4} \end{matrix}]

δ^{e} = [\begin{matrix} δ_{1}^{e} & δ_{2}^{e} & δ_{3}^{e} & δ_{4}^{e} \end{matrix}]^{T}

, and

δ_{i}^{e}

are the degrees of freedom at the ith node. The degrees of freedom at each node can be written as

δ_{i}^{e} = [\begin{matrix} u_{0 i} & v_{0 i} & w_{0 i} & u_{1 i} & u_{zi} & v_{1 i} & v_{zi} \end{matrix}]^{T}

(35)

\begin{matrix} \partial = [\begin{matrix} \partial / \partial x 00 \\ 0 \partial / \partial y 0 \\ \partial / \partial y \partial / \partial x 0 \\ \partial / \partial z 0 \partial / \partial x \\ 0 \partial / \partial z \partial / \partial y \end{matrix}] \end{matrix}

(36)

{\begin{matrix} B \end{matrix}}_{i} = [\begin{matrix} \frac{\partial N_{i}}{\partial x} & 0 & 0 & z \frac{\partial N_{i}}{\partial x} & Φ^{k} \frac{\partial N_{i}}{\partial x} & 0 & 0 \\ 0 & \frac{\partial N_{i}}{\partial y} & 0 & 0 & 0 & z \frac{\partial N_{i}}{\partial y} & Ψ^{k} \frac{\partial N_{i}}{\partial y} \\ \frac{\partial N_{i}}{\partial y} & \frac{\partial N_{i}}{\partial x} & 0 & z \frac{\partial N_{i}}{\partial y} & Φ^{k} \frac{\partial N_{i}}{\partial y} & z \frac{\partial N_{i}}{\partial x} & Ψ^{k} \frac{\partial N_{i}}{\partial x} \\ 0 & 0 & \frac{\partial N_{i}}{\partial x} & N_{i} & \frac{\partial Φ^{k}}{\partial z} N_{i} & 0 & 0 \\ 0 & 0 & \frac{\partial N_{i}}{\partial y} & 0 & 0 & N_{i} & \frac{\partial Ψ^{k}}{\partial z} N_{i} \end{matrix}]

(37)

The element stiffness matrix

The stiffness matrix of present quadrilateral element can be expressed on the basis of equation (27). The following equation within the ith element can be given by

\begin{matrix} \int_{S_{i}} \int_{- h / 2}^{h / 2} (σ_{x}^{k} δ ɛ_{x}^{k} + σ_{y}^{k} δ ɛ_{y}^{k} + τ_{xy}^{k} δ γ_{xy}^{k} + {\tilde{τ}}_{xz}^{k} δ γ_{xz}^{k} + {\tilde{τ}}_{yz}^{k} δ γ_{yz}^{k}) \\ \times d z d x d y - δ W_{i} = 0 \end{matrix}

(38)

where S_i is the area of ith element, and W_i is the work of external forces that act on ith element.

Substituting equations (10), (13), and (37) into equation (38), the element stiffness matrix ${\hat{K}}^{e}$ can be given by

{\hat{K}}^{e} = \int_{e} B^{T} Q \hat{B} d x d y d z

(39)

where Q is the transformed material constant matrix, the strain matrix B has been obtained in equation (37), and strain matrix

\hat{B}

can be written as

\hat{B} = [\begin{matrix} {\hat{B}}_{1} & {\hat{B}}_{2} & {\hat{B}}_{3} & {\hat{B}}_{4} \end{matrix}]

(40)

in which

{\begin{matrix} \hat{B} \end{matrix}}_{i} = [\begin{matrix} \frac{\partial N_{i}}{\partial x} & 0 & 0 & z \frac{\partial N_{i}}{\partial x} & Φ^{k} \frac{\partial N_{i}}{\partial x} & 0 & 0 \\ 0 & \frac{\partial N_{i}}{\partial y} & 0 & 0 & 0 & z \frac{\partial N_{i}}{\partial y} & Ψ^{k} \frac{\partial N_{i}}{\partial y} \\ \frac{\partial N_{i}}{\partial y} & \frac{\partial N_{i}}{\partial x} & 0 & z \frac{\partial N_{i}}{\partial y} & Φ^{k} \frac{\partial N_{i}}{\partial y} & z \frac{\partial N_{i}}{\partial x} & Ψ^{k} \frac{\partial N_{i}}{\partial x} \\ 0 & 0 & Θ_{3}^{k} \frac{\partial N_{i}}{\partial x} & S_{xz}^{k} Θ_{1}^{k} N_{i} & S_{xz}^{k} Θ_{2}^{k} N_{i} & 0 & 0 \\ 0 & 0 & ϒ_{3}^{k} \frac{\partial N_{i}}{\partial y} & 0 & 0 & S_{yz}^{k} ϒ_{1}^{k} N_{i} & S_{yz}^{k} ϒ_{2}^{k} N_{i} \end{matrix}]

(41)

where

S_{xz}^{k} = 1 / Q_{44}^{k}

S_{yz}^{k} = 1 / Q_{55}^{k}

Θ_{l}^{k}

and

ϒ_{l}^{k}

have been defined in equation (32), and l = 1∼3.

Based on the following equation, the vector of nodal displacement $δ^{e}$ can be given by

{\hat{K}}^{e} δ^{e} = P^{e}

(42)

where

P^{e}

is the generalized nodal load.

Geometric stiffness matrix of the element

For the buckling analysis, the geometric stiffness matrix of the element is introduced. According to equation (33), the first derivatives of transverse displacement can be written as

\begin{matrix} \frac{\partial w}{\partial x} = \sum_{i = 1}^{4} N_{i, x} w_{0 i} = G_{x} δ^{e} \\ \frac{\partial w}{\partial y} = \sum_{i = 1}^{4} N_{i, y} w_{0 i} = G_{y} δ^{e} \end{matrix}

(43)

where

[\begin{matrix} G_{x} \\ G_{y} \end{matrix}] = [\begin{matrix} G_{x 1} & G_{x 2} & G_{x 3} & G_{x 4} \\ G_{y 1} & G_{y 2} & G_{y 3} & G_{y 4} \end{matrix}]

(44)

\begin{matrix} \begin{matrix} \begin{matrix} G_{xi} = [\begin{matrix} 0 & 0 & N_{i, x} & 0 & 0 & 0 & 0 \end{matrix}] \\ G_{yi} = [\begin{matrix} 0 & 0 & N_{i, y} & 0 & 0 & 0 & 0 \end{matrix}] \end{matrix} & i = 1 \sim 4 \end{matrix} \end{matrix}

(45)

Therefore, the geometric stiffness matrix of the element is given as follows

{\hat{K}}_{G}^{e} = \int_{V_{e}} G^{T} NG d v

(46)

where

G = [\begin{matrix} G_{x} \\ G_{y} \end{matrix}]

(47)

N = [\begin{matrix} N_{xx} & N_{xy} \\ N_{xy} & N_{yy} \end{matrix}]

(48)

and

N_{xx}

N_{yy}

are the tractions in x and y directions, respectively,

N_{xy}

is the shear force. The solutions for the critical buckling loads are given by

(\sum_{e} {\hat{K}}^{e} - λ \sum_{e} {\hat{K}}_{G}^{e}) Γ = 0

(49)

where λ is the critical load and

Γ

is the vector of the mode shape of buckling.

Numerical examples and discussion

To assess the performance of the proposed zigzag model, a number of numerical examples of laminated sandwich plates with different lamination sequences, material characteristics, boundary conditions and span-to-thickness ratios will be investigated. The 3D elasticity solutions and the solutions calculated from the layerwise theory are employed as benchmark solutions. The material constants used in all examples are given by

Material (1), sandwich plate⁴⁸

Face sheets:

\begin{matrix} E_{1} = 172.4 GPa, E_{2} = E_{3} = 6.89 GPa, \\ G_{12} = G_{13} = 3.45 GPa, G_{23} = 1.378 GPa, \\ v_{12} = v_{13} = v_{23} = 0.25 \end{matrix}

Core:

\begin{matrix} E_{1}^{c} = E_{2}^{c} = 0.276 GPa, E_{3}^{c} = 3.45 GPa, \\ G_{13} = G_{23} = 0.413 GPa, G_{12}^{c} = 0.1102 GPa, \\ v_{12}^{c} = v_{13}^{c} = v_{23}^{c} = 0.25 \end{matrix}

Material (2), sandwich plate⁴⁹

Face sheets:

\begin{matrix} E_{1} = 25 GPa, E_{2} = E_{3} = 1 GPa, G_{12} = G_{13} = 0.5 GPa, \\ G_{23} = 0.2 GPa, v_{12} = v_{13} = v_{23} = 0.25 \end{matrix}

Core:

\begin{matrix} E_{1}^{c} = E_{2}^{c} = 0.04 GPa, E_{3}^{c} = 0.5 GPa, G_{12}^{c} = 0.016 GPa, \\ G_{13}^{c} = G_{23}^{c} = 0.06 GPa, v_{12}^{c} = v_{13}^{c} = v_{23}^{c} = 0.25 \end{matrix}

Material (3), laminated plate⁵⁰

\begin{matrix} E_{1} = 40 GPa, E_{2} = E_{3} = 1 GPa, G_{12} = G_{13} = 0.6 GPa, \\ G_{23} = 0.5 GPa, v_{12} = v_{13} = v_{23} = 0.25 \end{matrix}

Material (4), sandwich plate⁵¹

Face sheets:

E = 70 GPa, v = 0.3

Core:

\begin{matrix} E_{1}^{c} = 1 \times 10^{- 8} GPa, E_{2}^{c} = 0.109 GPa, \\ G_{12}^{c} = 0.0266 GPa, v_{12}^{c} = 1 \times 10^{- 5} \end{matrix}

Material (5), sandwich plate⁵²

Face sheets:

\begin{matrix} E_{1} / E_{2} = 19, E_{2} = E_{3}, G_{12} / E_{2} = G_{13} / E_{2} = 0.52, \\ G_{23} / E_{2} = 0.338, v_{12} = v_{13} = 0.32, v_{23} = 0.49 \end{matrix}

Core:

\begin{matrix} E_{1}^{c} / E_{2} = 3.2 \times 10^{- 5}, E_{2}^{c} / E_{2} = 2.9 \times 10^{- 5}, \\ E_{3}^{c} / E_{2} = 0.4, G_{12}^{c} / E_{2} = 2.4 \times 10^{- 3}, \\ G_{13}^{c} / E_{2} = 7.9 \times 10^{- 2}, G_{23}^{c} / E_{2} = 6.6 \times 10^{- 2}, \\ v_{12}^{c} = 0.99, v_{13}^{c} = v_{23}^{c} = 3 \times 10^{- 5} \end{matrix}

Example 1: Square sandwich plates under sinusoidal lateral loading q = q₀sin(πx/a). This plate is simply supported along the edges parallel to the y-axis.

Results reported in Example 1 are normalized according to the following relations

\bar{u} = \frac{E_{2}}{q_{0} h} u, \bar{w} = \frac{100 E_{2} h^{3}}{q_{0} a^{4}} w, ({\bar{σ}}_{x}, {\bar{τ}}_{xz}) = \frac{1}{q_{0}} (σ_{x}, τ_{xz})

For all the examples, a and b denote the width and length of the laminated plates, and h is the total lamination thickness. The finite element mesh is shown in Figure 2. The symmetric 11-ply double-core sandwich plate [0°/90°/0°/C/0°/90°/0°/C/0°/90°/0°] made of material (1) is firstly analyzed, where the each core thickness is 7h/20, and the thickness of each ply in the sheets is h/30. In Table 1, the displacements and stresses obtained by different mesh sizes for the 11-layer sandwich plate are presented. The values in Table 1 reveal that the converged results computed from the proposed MZZT model are close to the 3D elasticity solutions (Exact).⁵³ Moreover, the values of displacements and stresses converge at the mesh size of $28 \times 28$ , so that a mesh size of $28 \times 28$ may be considered to be adequate for the thick sandwich plates. Table 2 presents the comparison of displacements and stresses obtained from different models in detail. The errors of various models with respect to the exact solutions are also computed and presented in the table. Small errors can be observed in the results obtained from the proposed MZZT model, while large errors are there in the GLHT,¹⁸ ZZTC,⁵⁴ HSDT-Reddy⁶ and HSDT-Thai et al.⁵⁵ models. The through-the-thickness distributions of in-plane displacement and stresses obtained from different models for the 11-ply sandwich plate (a/h = 8) are depicted in Figure 3. It is found that the proposed MZZT model is in good agreement with the 3D elasticity solutions (Exact). However, the GLHT, ZZTC, HSDT-Reddy and HSDT-Thai et al. models are apparently less accurate. Figure 4 considers a symmetric 17-ply double-core sandwich plate [0°/90°/0°/90°/0°/C/0°/90°/0°/90°/0°/C/0°/90°/0°/90°/0°] made of material (1) with a/h = 8, where the each core thickness is 7 h/20, and the thickness of each ply in the face sheets is h/50. The results depicted in Figure 4 show that the proposed MZZT model matches the exact solutions quite well.

Table 1.

Convergence of displacements and stresses for the 11-ply [0°/90°/0°/C/0°/90°/0°/C/ 0°/90°/0°] sandwich plate (a/h = 5).

Models	$\bar{w} (\frac{a}{2}, \frac{b}{2}, 0)$	$\bar{u} (0, \frac{b}{2}, - \frac{h}{2})$	${\bar{σ}}_{x} (\frac{a}{2}, \frac{b}{2}, - \frac{h}{2})$	${\bar{τ}}_{xz} (0, \frac{b}{2}, 0)$
MZZT $(8 \times 8)$	7.1887	3.6986	−56.813	2.2468
MZZT $(16 \times 16)$	7.2823	3.7673	−58.990	2.0306
MZZT $(24 \times 24)$	7.2998	3.7802	−59.400	1.9893
MZZT $(28 \times 28)$	7.3035	3.7829	−59.488	1.9805
Exact⁵³	7.2069	3.6968	−58.142	1.9533

Table 2.

Comparison of displacements and stresses obtained from different models.

a/h	Models	$\bar{w} (\frac{a}{2}, \frac{b}{2}, 0)$	$\bar{u} (0, \frac{b}{2}, - \frac{h}{2})$	${\bar{σ}}_{x} (\frac{a}{2}, \frac{b}{2}, - \frac{h}{2})$	${\bar{τ}}_{xz} (0, \frac{b}{2}, 0)$
5	Exact⁵³	7.2069	3.6968	−58.142	1.9533
	MZZT	7.3035 (1.34%)	3.7829 (2.33%)	−59.488 (2.32%)	1.9805 (1.39%)
	GLHT¹⁸	5.7148 (20.7%)	3.4072 (7.83%)	−53.700 (7.64%)	6.3343 (224%)
	ZZTC⁵⁴	7.4981 (4.04%)	3.8610 (4.44%)	−61.302 (5.43%)	1.9328 (1.05%)
	HSDT-Reddy⁶	4.5181 (37.3%)	3.3540 (9.27%)	−52.863 (9.08%)	3.3103 (69.5%)
	HSDT-Thai et al.⁵⁵	4.5927 (36.3%)	3.2509 (12.1%)	−51.237 (11.9%)	3.1675 (62.2%)
10	Exact⁵³	2.9510	24.985	−196.48	3.7117
	MZZT	2.9389 (0.41%)	24.808 (0.71%)	−195.07 (0.72%)	3.7897 (2.10%)
	GLHT¹⁸	2.5354 (14.1%)	24.144 (3.49%)	−190.27 (3.16%)	12.797 (245%)
	ZZTC⁵⁴	3.1170 (5.63%)	24.795 (0.76%)	−202.77 (3.20%)	4.2556 (14.7%)
	HSDT-Reddy⁶	2.2360 (24.2%)	24.023 (3.85%)	−189.32 (3.64%)	6.6592 (79.4%)
	HSDT-Thai et al.⁵⁵	2.2529 (23.7%)	23.824 (4.25%)	−187.75 (4.44%)	6.3772 (71.8%)

Note: The numbers in bracket are the errors of various models relative to the exact solutions.

In order to further assess the accuracy of the proposed model, a thick three-ply single-core sandwich plate [0°/C/0°] made of material (1) is investigated, where the core thickness is 4 h/5, and the thickness of each face sheet is h/10. Table 3 presents the comparison of displacements and stresses of three-ply sandwich plates with different span-to-thickness ratios (a/h). The errors of present MZZT model respect to the exact solutions are also computed and presented in the table. The maximum errors of the present results with respect to the 3D elasticity solutions (Exact) are less than 4%. In addition, the influence of the stiffness ratio between the face layer and the core on the effectiveness of the present model is discussed. For the sandwich plate [0°/C/0°] made of material (1), the modulus parameter ratios (E₁ ^f / E₁ ^c ) between the skins and core is 6.2464 × 10². The core properties are divided by 10, 100 and 1000, so the modulus parameter ratios (E₁ ^f /E₁ ^c ) between the skins and core are: 6.2464 × 10³, 6.2464 × 10⁴ and 6.2464 × 10⁵, respectively. Table 4 presents the comparison of lateral deflections and transverse shear stresses of three-ply sandwich plates with different modulus parameter ratios (E₁ ^f /E₁ ^c ). The errors of present MZZT model relative to the exact solutions are also computed and presented in the table. Small errors can be observed in the results obtained from the proposed MZZT model, even for the sandwich plate with a very soft core.

Table 3.

Comparison of displacements and stresses of three-ply sandwich plate [0°/C/0°] with different span-to-thickness ratios (a/h).

a/h	Models	$\bar{w} (\frac{a}{2}, \frac{b}{2}, 0)$	$\bar{u} (0, \frac{b}{2}, - \frac{h}{2})$	${\bar{σ}}_{x} (\frac{a}{2}, \frac{b}{2}, \frac{h}{2})$	${\bar{τ}}_{xz} (0, \frac{b}{2}, 0)$
5	Exact⁵³	7.4289	3.0745	50.939	1.7258
	MZZT	7.5009 (0.97%)	3.1121 (1.22%)	48.924 (3.96%)	1.7046 (1.23%)
10	Exact⁵³	2.6612	18.211	142.84	3.5079
	MZZT	2.6405 (0.78%)	17.978 (1.28%)	141.35 (1.04%)	3.4652 (1.22%)
20	Exact⁵³	1.4228	131.32	515.53	7.0441
	MZZT	1.4001 (1.60%)	128.97 (1.79%)	507.08 (1.64%)	6.9777 (0.94%)
50	Exact⁵³	1.0726	1987.2	3126.7	17.627
	MZZT	1.0439 (2.68%)	1933.4 (2.71%)	3044.0 (2.64%)	17.762 (0.77%)
100	Exact⁵³	1.0225	15823	12453	35.266
	MZZT	0.9960 (2.59%)	15412 (2.60%)	12142 (2.50%)	35.315 (0.14%)

Note: The numbers in bracket are the errors of present MZZT model relative to the exact solutions.

Table 4.

Comparison of displacements and stresses of three-ply sandwich plate [0°/C/0°] with different modulus ratios.

E₁ ^f / E₁ ^c	Models	$\bar{w} (\frac{a}{2}, \frac{b}{2}, 0)$	${\bar{τ}}_{xz} (0, \frac{b}{2}, 0)$
6.2464 × 10²	Exact⁵³	2.6612	3.5079
	MZZT	2.6405 (0.78%)	3.4652 (1.22%)
6.2464 × 10³	Exact⁵³	16.478	3.3036
	MZZT	16.275 (1.23%)	3.2385 (1.97%)
6.2464 × 10⁴	Exact⁵³	99.369	2.1021
	MZZT	98.070 (1.31%)	2.0626 (1.88%)
6.2464 × 10⁵	Exact⁵³	213.15	0.4534
	MZZT	209.43 (1.75%)	0.4430 (2.29%)

Note: The numbers in bracket are the errors of present MZZT model relative to the exact solutions.

Example 2: Sandwich plates under double sinusoidal lateral loading $q = q_{0} sin (π x / a) sin (π y / b)$ on the top surface of the plate.

Results reported in Example 2 are normalized according to the following relations

\begin{matrix} \bar{u} = \frac{E_{2} h^{2}}{q_{0} a^{3}} u, ({\bar{σ}}_{x}, {\bar{σ}}_{y}, {\bar{τ}}_{xy}) = \frac{h^{2}}{q_{0} a^{2}} (σ_{x}, σ_{y}, τ_{xy}), ({\bar{τ}}_{xz}, {\bar{τ}}_{yz}) \\ = \frac{h}{q_{0} a} (τ_{xz}, τ_{yz}) \end{matrix}

Firstly, a simply supported (SSSS) three-ply single-core [0°/C/0°] square sandwich plate under double sinusoidal transverse loading is analyzed. The face sheets have equal thicknesses h/10, and the core thickness is 4 h/5. The material properties of the sandwich plate are the same as those used in Example 1. The comparison of displacement and stresses is shown in Figure 5. It can be seen from Figure 5 that the results obtained from the proposed MZZT model match well with the 3D elasticity solutions (Exact).⁵⁶ In order to further assess the accuracy of the proposed theory, a thick 5-ply [90°/0°/C/90°/0°] antisymmetric sandwich plate made of material (2) is studied, where the core thickness is 4 h/5, and the thickness of each ply in the face sheets is h/20. Three types of boundary conditions (BCs) are used in this study. These are SSSS (all the edges are simply supported), CCCC (all the edges are clamped) and SCSC (the edges parallel to the y-axis are simply supported and the other two edges are clamped). Comparison of transverse displacement and stresses obtained in the present analysis with the available solutions in the literature^49,56 is shown in Table 5. The results in Table 5 show that the proposed MZZT model agrees well with the exact solutions⁵⁶ and the results obtained by Khandelwal et al.⁴⁹ The through-the-thickness distributions of transverse shear stresses are plotted in Figures 6 and 7. Figures 6 and 7 show that the transverse shear stresses obtained from the MZZT model are in good agreement with the exact solutions.⁵⁶ Nevertheless, the results obtained by the refined higher-order zigzag model (RHZZT)⁵⁷ are less accurate.

Table 5.

Comparison of displacement and stresses obtained from different models (a/h = 4).

BCs	Models	$\bar{w} (\frac{a}{2}, \frac{b}{2}, 0)$	${\bar{σ}}_{x} (\frac{a}{2}, \frac{b}{2}, \frac{h}{2})$	${\bar{τ}}_{xz} (0, \frac{b}{2}, 0)$	${\bar{τ}}_{yz} (\frac{a}{2}, 0, 0)$
SSSS	Exact⁵⁶	6.1105	1.3770	0.1756	0.1751
	MZZT	6.1247	1.3652	0.1740	0.1728
	Khandelwal et al.⁴⁹	6.1342	1.4086	0.1718	0.1692
CCCC	MZZT	5.2225	0.7335	0.1182	0.1153
	Khandelwal et al.⁴⁹	5.2470	0.7442	0.1236	0.1218
SCSC	MZZT	5.6417	1.2765	0.1596	0.1211
	Khandelwal et al.⁴⁹	5.6638	1.3197	0.1576	0.1319

Example 3: Buckling analysis of a three-layer [0°/90°/0°] laminated plate with different boundary conditions.

In this example, a three-layer [0°/90°/0°] cross-ply square laminated plate made of materials (3) is chosen to calculate the biaxial critical loads. The critical buckling loads are normalized as

{\bar{N}}_{xx} = λ a^{2} / (E_{2} h^{3}), {\bar{N}}_{yy} = {\bar{N}}_{xx}

. The span-to-thickness ratio a/h is taken to be 10. Three types of boundary conditions SSSS (all the edges are simply supported), SFSF (the edges parallel to the y-axis are simply supported and the other two edges are free) and SCSC (the edges parallel to the y-axis are simply supported and the other two edges are clamped) are used in this study. Table 6 tabulates the biaxial critical buckling loads obtained by the present MZZT model together with the FSDT solutions in literature^50,58,59 and the HSDT solutions in Liu et al.⁶⁰ and the ZZT solutions in Ferreira et al.⁶¹ Numerical results show that the present results (MZZT) are in good agreement with the referenced results.

Table 6.

Normalized biaxial critical loads of three-layer [0°/90°/0°] laminated plate (a/h = 10).

Models	SSSS	SFSF	SCSC
MZZT	10.190	1.896	13.411
Liu et al.⁶⁰	10.091	1.928	12.952
Khdeir and Librescu⁵⁸	10.259	1.937	13.288
Liew and Huang⁵⁹	10.120	1.926	13.225
Nguyen-Van et al.⁵⁰	10.252	1.858	13.249
Ferreira et al.⁶¹	10.209	–	13.658

Example 4: Buckling behaviors of the soft core sandwich plates subjected to uniaxial compression are analyzed.

Firstly, a three-layer square sandwich plate made of material (4) is studied. This plate is simply supported along the edges parallel to the y-axis and the other edges are free (SFSF). The critical loads are normalized as

{\bar{N}}_{xx} = λ a^{2} / (E_{2} h^{3})

, t_c and t_f represent the thickness of core and face sheet, respectively. Table 7 tabulates the normalized uniaxial critical loads of a three-layer plate with different ratios (t_c/t_f) and span-to-thickness ratios (a/h). The errors of various theories respect to the layerwise theory are presented in Table 7. Numerical results show that the present results (MZZT) match well with the layerwise solutions (MLW).⁵¹ The maximum errors of the critical loads for the zigzag model EDZ3⁶² and the ninth-order model (NSDT) are 29.23% and 34.72%, respectively, while a very large value (1653%) of the maximum error in the critical loads is found for other global higher order models.^6,63–65 It needs to be shown that the proposed MZZT model only includes seven displacement variables and a good agreement between the present results (MZZT) and the layerwise solutions (MLW) can be observed, which demonstrate the accuracy and effectiveness of the proposed zigzag model. In order to further evaluate the accuracy of the proposed zigzag model, a square simply supported symmetric sandwich plate of material (5) is considered for buckling analysis. The plate [(0°/90°)₅/C/(90°/0°)₅] is made of equal thickness cross-ply laminated face sheets with 10 layers and a soft core. The critical buckling loads are normalized as

{\bar{N}}_{xx} = λ a^{2} / (E_{2} h^{3})

. Comparison of the normalized uniaxial critical loads is shown in Table 8. The errors of proposed zigzag model with respect to the 3D elasticity solutions⁵² are also computed and presented in the table. Small errors can be observed in the results obtained from the proposed MZZT model. It is also seen from Table 8 that the critical buckling loads increase with the ratios t_f/h and the span-to-thickness ratios a/h.

Table 7.

Comparison of normalized uniaxial critical loads obtained from various models.

		a/h
t_c/t_f	Models	2	5	10
5	MLW⁵¹	0.006222	0.01432	0.04108
	MZZT	0.006298 (1.22%)	0.01474 (2.93%)	0.04212 (2.53%)
	EDZ3⁶²	0.007854 (29.23%)	0.01597 (11.52%)	0.04285 (4.30%)
	NSDT⁶⁶	0.007216 (15.98%)	0.01864 (30.17%)	0.05535 (34.72%)
	HSDT-Reddy⁶	0.06434 (934.1%)	0.2486 (1636%)	0.4313 (949.8%)
	HSDT-Soldatos⁶³	0.06532 (949.8%)	0.2510 (1653%)	0.4332 (954.4%)
	HSDT-Ferreira et al.⁶⁴	0.05391 (766.4%)	0.2203 (1438%)	0.4083 (893.8%)
	HSDT-Akavci⁶⁵	0.05702 (816.4%)	0.2291 (1500%)	0.4158 (912.1%)
25	MLW⁵¹	0.0015299	0.0090314	0.031096
	MZZT	0.0016016 (4.69%)	0.0092038 (1.91%)	0.031394 (0.958%)
	NSDT⁶⁶	0.001777 (16.15%)	0.01017 (12.61%)	0.03477 (11.82%)
	HSDT-Reddy⁶	0.003343 (118.5%)	0.01889 (109.2%)	0.05785 (86.04%)
	HSDT-Soldatos⁶³	0.003385 (121.3%)	0.01912 (111.7%)	0.05837 (87.71%)
	HSDT-Ferreira et al.⁶⁴	0.002929 (91.45%)	0.01671 (85.02%)	0.05418 (69.06%)
	HSDT-Akavci⁶⁵	0.003052 (99.49%)	0.01736 (92.22%)	0.05418 (74.24%)

Note: The numbers in bracket are the errors of various models relative to the layerwise model.⁵¹

Table 8.

Normalized uni-axial critical loads for a 21-ply [(0°/90°)₅/core/(90°/0°)₅] sandwich plate.

		t_f/h
a/h	Models	0.025	0.05	0.075	0.01
5	3D⁵²	1.5027	2.0816	2.3973	2.6051
	MZZT	1.5420 (2.62%)	2.1167 (1.69%)	2.4234 (1.09%)	2.6245 (0.745%)
20/3	3D⁵²	1.8438	2.7911	3.3729	3.7883
	MZZT	1.8986 (2.97%)	2.8460 (1.97%)	3.4186 (1.35%)	3.8098 (0.568%)
10	3D⁵²	2.2376	3.7375	4.7637	5.6081
	MZZT	2.2803 (1.87%)	3.7827 (1.21%)	4.8178 (1.14%)	5.6478 (0.708%)

Note: The numbers in bracket are the errors of various models relative to the 3D elasticity solutions.⁵²

Conclusions

An efficient zigzag theory (MZZT) is proposed for bending and buckling analysis of thick laminated soft core sandwich plates. Based on the 3D equilibrium equations and the RMVT, the accurate and layer interface-continuous transverse shear stresses can be obtained. The second-order derivatives of in-plane displacement parameters have been eliminated from the transverse shear stresses, so a simple four-node C⁰ quadrilateral plate element is developed. By investigating the bending and buckling problems of laminated sandwich plates, the following conclusions can be drawn:

Comparison of the present results with the 3D elasticity solutions and other models reported in the literature shows that the proposed MZZT model can accurately predict the displacements and stresses of thick soft core sandwich plates without any post-processing procedures.

The proposed zigzag model can accurately predict the critical buckling loads of laminated composite and sandwich plates

The interface-continuous transverse shear stresses are included in the strain energy which can improve the accuracy of critical buckling loads.

Footnotes

Declaration of Conflicting Interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: The work described in this paper was supported by the National Key Research and Development Program of China [No. 2016YFB0200702].

ORCID iD

Qilin Jin

References

Reissner

Stavsky

. Bending and stretching of certain types of heterogeneous aelotropic elastic plates. J Appl Mech 1961; 28: 402–408.

Whitney

Pagano

. Shear deformation in heterogeneous anisotropic plates. J Appl Mech 1970; 37: 1031–1036.

Sofyev

. The vibration and buckling of sandwich cylindrical shells covered by different coatings subjected to the hydrostatic pressure. Compos Struct 2014; 117: 124–134.

Khalili

SMR

Botshekanan Dehkordi

Carrera

, et al. Non-linear dynamic analysis of a sandwich beam with pseudoelastic SMA hybrid composite faces based on higher order finite element theory. Compos Struct 2013; 96: 243–255.

Kant

Swaminathan

. Analytical solutions for the static analysis of laminated composite and sandwich plates based on a higher order refined theory. Compos Struct 2002; 56: 329–344.

Reddy

. A simple higher-order theory for laminated composite plates. J Appl Mech 1984; 51: 745–752.

Robbins

Reddy

. Modelling of thick composites using a layer wise laminate theory. J Numer Methods Eng 1993; 36: 665–677.

Shariyat

. Thermal buckling analysis of rectangular composite plates with temperature-dependent properties based on a layerwise theory. Thin Wall Struct 2007; 45: 439–452.

Dafedar

Desai

Mufti

. Stability of sandwich plates by mixed, higher-order analytical formulation. Int J Solids Struct 2003; 40: 4501–4517.

10.

Rao

Desai

. Analytical solutions for vibrations of laminated and sandwich plates using mixed theory. Compos Struct 2004; 63: 361–373.

11.

Liu

Ferreira

AJM

Xing

, et al. Analysis of composite plates using a layerwise theory and a differential quadrature finite element method. Compos Struct 2016; 156: 393–398.

12.

Yan

Pagani

Carrera

, et al. Exact solutions for the macro-, meso- and micro-scale analysis of composite laminates and sandwich structures. J Compos Mater 2018; 52: 3109–3124.

13.

Murakami

. Laminated composite plate theory with improved in-plane response. J Appl Mech 1986; 53: 661–666.

14.

Carrera

. On the use of the Murakami’s zig-zag function in the modeling of layered plates and shell. Comput Struct 2004; 82: 541–554.

15.

Di Sciuva

. Multilayered anisotropic plate models with continuous interlaminar stresses. Compos Struct 1992; 22: 149–167.

16.

Cho

Parmerter

. An efficient higher order plate theory for laminated composites. Compos Struct 1992; 20: 113–123.

17.

Cho

Parmerter

. Efficient higher-order plate theory for general lamination configuration. AIAA J 1993; 31: 1299–1306.

18.

Liu

. Generalized laminate theories based on double superposition hypothesis. Int J Numer Methods Eng 1997; 40: 1197–1212.

19.

Chen

. Refined laminated composite plate element based on global-local higher-order shear deformation theory. Compos Struct 2005; 70: 135–152.

20.

Chen

. Buckling analysis of angle-ply composite and sandwich plates by combination of geometric stiffness matrix. Comput Mech 2007; 39: 839–848.

21.

Ren

. Finite element analysis of hygrothermal effects on angle-ply composite plates. J Compos Mater 2016; 50: 2215–2233.

22.

Chen

. A refined three-node triangular element based on the HW variational theorem for multilayered composite plates. Compos Struct 2017; 161: 132–144.

23.

Chen

. Interlaminar stress analysis of multilayered composites based on the Hu-Washizu variational theorem. J Compos Mater 2018; 52: 1765–1779.

24.

Nath

Mishra

. Enhanced zigzag theory for static analysis of composite plates. Arch Appl Mech 2018; 88: 1595–1615.

25.

Chen

. A C⁰-type higher-order theory for bending analysis of laminated composite and sandwich plates. Compos Struct 2010; 92: 653–661.

26.

Chen

. A global-local higher order theory including interlaminar stress continuity and C⁰ plate bending element for cross-ply laminated composite plates. Comput Mech 2010; 45: 387–400.

27.

Sze

, et al. An improved in-plane displacement model for the stability analysis of laminated composites with general lamination configurations. Compos Struct 2011; 93: 1584–1594.

28.

Chen

Ren

. An accurate higher-order theory and C⁰ finite element for free vibration analysis of laminated composite and sandwich plates. Compos Struct 2010; 92: 1299–1307.

29.

Jin

Yao

. Hygrothermal analysis of laminated composite plates in terms of an improved C⁰ -type global-local model. Aerosp Sci Technol 2017; 63: 328–343.

30.

Shariyat

Khalili

SMR

Rajabi

. A global-local theory with stress recovery and a new post-processing technique for stress analysis of asymmetric orthotropic sandwich plates with single/dual cores. Comput Methods Appl Mech Eng 2015; 286: 192–215.

31.

Ren

Chen

. A new zig-zag theory and C⁰ plate bending element for composite and sandwich plates. Arch Appl Mech 2011; 81: 185–197.

32.

Sahoo

Singh

. Assessment of inverse trigonometric zigzag theory for stability analysis of laminated composite and sandwich plates. Int J Mech Sci 2015; 101–102: 145–154.

33.

Chalak

Chakrabarti

Iqbal

, et al. An improved C⁰ FE model for the analysis of laminated sandwich plate with soft core. Finite Elem Anal Des 2012; 56: 20–31.

34.

Chalak

Chakrabarti

Sheikh

, et al. Stability analysis of laminated soft core sandwich plates using higher order zig-zag plate theory. Mech Adv Mater Struct 2015; 22: 897–907.

35.

Murakami

. A laminated beam theory with interlayer slip. J Appl Mech 1984; 51: 551–559.

36.

Reissner

. On a certain mixed variational theorem and a proposed application. Int J Numer Methods Eng 1984; 20: 1366–1368.

37.

Carrera

. Developments, ideas, and evaluations based upon Reissner’s Mixed Variational Theorem in the modeling of multilayered plates and shells. Appl Mech Rev 2001; 54: 301–329.

38.

Messina

. Free vibrations of multilayered plates based on a mixed variational approach in conjunction with global piecewise-smooth functions. J Sound Vib 2002; 256: 103–129.

39.

Kim

Cho

. Enhanced first-order theory based on mixed formulation and transverse normal effect. Int J Solids Struct 2007; 44: 1256–1276.

40.

Lezgy-Nazargah

. A high-performance parametrized mixed finite element model for bending and vibration analyses of thick plates. Acta Mech 2016; 227: 3429–3450.

41.

Lezgy-Nazargah

Meshkani

. An efficient partial mixed finite element model for static and free vibration analyses of FGM plates rested on two-parameter elastic foundations. Struct Eng Mech 2018; 66: 665–676.

42.

Lezgy-Nazargah

Salahshuran

. A new mixed-field theory for bending and vibration analysis of multi-layered composite plate. Arch Civil Mech Eng 2018; 18: 818–832.

43.

Tessler

. Refined zigzag theory for homogeneous, laminated composite, and sandwich beams derived from Reissner’s mixed variational principle. Meccanica 2015; 50: 2621–2648.

44.

Iurlaro

Gherlone

Di Sciuva

, et al. Refined zigzag theory for laminated composite and sandwich plates derived from Reissner’s mixed variational theorem. Compos Struct 2015; 133: 809–817.

45.

Groh

RMJ

Tessler

. Computationally efficient beam elements for accurate stresses in sandwich laminates and laminated composites with delaminations. Comput Methods App Mech Eng 2017; 320: 369–395.

46.

Tessler A, Di Sciuva M and Gherlone M. Refinement of Timoshenko beam theory for composite and sandwich beams using zigzag kinematics. Technical Report NASA/TP-2007-215086.

47.

Chen

. An assessment of several displacement based theories for the vibration and stability analysis of laminated composite and sandwich beams. Compos Struct 2008; 84: 337–349.

48.

Matsunaga

. Assessment of a global higher-order deformation theory for laminated composite and sandwich plates. Compos Struct 2002; 56: 279–291.

49.

Khandelwal

Chakrabarti

Bhargava

. An efficient FE model based on combined theory for the analysis of soft core sandwich plate. Comput Mech 2013; 51: 673–697.

50.

Nguyen-Van

Mai-Duy

Karunasena

, et al. Buckling and vibration analysis of laminated composite plate/shell structures via a smoothed quadrilateral flat shell element with in-plane rotations. Comput Struct 2011; 89: 612–625.

51.

Dafedar

Desai

. Stability of composite and sandwich struts by mixed formulation. J Eng Mech 2004; 130: 762–770.

52.

Noor

Peters

Burton

. Three-dimensional solutions for initially stressed structural sandwiches. J Eng Mech ASCE 1994; 120: 284–303.

53.

Pagano

. Exact solutions for composite laminates in cylindrical bending. J Compos Mater 1969; 3: 398–411.

54.

Cho

. Higher order zig-zag plate theory under thermo-electric- mechanical loads combined. Composites 2003; 34: 67–82.

55.

Thai

Tran

, et al. Analysis of laminated composite plates using higher order shear deformation plate theory and node-based smoothed discrete shear gap method. Appl Math Model 2012; 36: 5657–5677.

56.

Pagano

. Exact solutions for rectangular bidirectional composites. J Compos Mater 1970; 4: 20–34.

57.

Khandelwal

Chakrabarti

Bhargava

. Accurate calculation of transverse shear stresses for soft-core sandwich laminates. Acta Mech 2014; 225: 2877–2891.

58.

Khdeir

Librescu

. Analysis of symmetric cross-ply elastic plates using a higher-order theory. Part II: buckling and free vibration. Compos Struct 1988; 9: 259–277.

59.

Liew

Huang

. Bending and buckling of thick symmetric rectangular laminates using the moving least-squares differential quadrature method. Int J Mech Sci 2003; 45: 95–114.

60.

Liu

Chua

Ghista

. Mesh-free radial basis function method for static, free vibration and buckling analysis of shear deformable composite laminates. Compos Struct 2007; 78: 58–69.

61.

Ferreira

AJM

Roque

CMC

Carrera

, et al. Radial basis functions collocation and a unified formulation for bending, vibration and buckling analysis of laminated plates, according to variation of Murakami’s zigzag theory. Eur J Mech A 2011; 30: 559–570.

62.

Carrera

Ciuffreda

. Closed-form solutions to assess multilayered-plate theories for various thermal stress problems. J Therm Stress 2004; 27: 1001–1031.

63.

Soldatos

. A transverse shear deformation theory for homogeneous monoclinic plates. Acta Mech 1992; 94: 195–200.

64.

Ferreira

AJM

Roque

CMC

Jorge

RMN

. Analysis of composite plates by trigonometric shear deformation theory and multiquadrics. Comput Struct 2005; 83: 2225–2237.

65.

Akavci

. Two new hyperbolic shear displacement models for orthotropic laminated composite plates. Mech Compos Mater 2010; 46: 215–226.

66.

Matsunaga

. Vibration and buckling of multilayered composite beams according to higher-order deformation theories. J Sound Vib 2001; 246: 47–62.

An accurate zigzag theory for bending and buckling analysis of thick laminated sandwich plates with soft core

Abstract

Keywords

Introduction

Kinematic assumptions

Displacement field

Improved transverse shear stresses

Finite element formulation

A C0-type four-node quadrilateral element

The element stiffness matrix

Geometric stiffness matrix of the element

Numerical examples and discussion

Conclusions

Footnotes

Declaration of Conflicting Interests

Funding

ORCID iD

References

A C⁰-type four-node quadrilateral element