A family of higher-order single layer plate models meeting $C^0_z$ -- requirements for arbitrary laminates

TL;DR

This paper introduces a generic method to extend higher-order shear deformation theories to $C^0_z$ models for arbitrary laminates, improving accuracy in predicting laminated plate behavior.

Contribution

It presents an automatic approach to develop $C^0_z$ ESL plate models with enhanced cross-sectional warping description, ensuring interlaminar continuity and homogeneity conditions.

Findings

Zig-zag models outperform basis models in deflection predictions.

Models accurately satisfy interlaminar continuity and boundary conditions.

Improved results in frequency and stress response for laminated plates.

Abstract

In the framework of displacement-based equivalent single layer (ESL) plate theories for laminates, this paper presents a generic and automatic method to extend a basis higher-order shear deformation theory (polynomial, trigonometric, hyperbolic, ...) to a multilayer $C^0_z$ higher-order shear deformation theory. The key idea is to enhance the description of the cross-sectional warping: the odd high-order $C^1_z$ function of the basis model is replaced by one odd and one even high-order function and including the characteristic zig-zag behaviour by means of piecewise linear functions. In order to account for arbitrary lamination schemes, four such piecewise continuous functions are considered. The coefficients of these four warping functions are determined in such a manner that the interlaminar continuity as well as the homogeneity conditions at the plate's top and bottom surfaces are {\em a priori} exactly verified by the transverse shear stress field. These $C_z^0$ ESL models all have the same number of DOF as the original basis HSDT. Numerical assessments are presented by referring to a strong-form Navier-type solution for laminates with arbitrary stacking sequences as well for a sandwich plate. In all practically relevant configurations for which laminated plate models are usually applied, the results obtained in terms of deflection, fundamental frequency and local stress response show that the proposed zig-zag models give better results than the basis models they are issued from.

Tables

Table 1. Table 1: Original bHSDT ϕ ​ ( z ) italic-ϕ 𝑧 \phi(z) functions. The parameter k 𝑘 k of the hyperbolic function allows to generalise the original model with k = 1 𝑘 1 k=1 [ 35 ] .

Model	Type	$\varphi \left(z\right)$
Reddy[33]	Polynomial	$z-\frac{4z^3}{3h^2}$
Touratier[34]	Trigonometric	$\frac{h}{\pi }\sin \left(\pi \frac{z}{h}\right)$
Soldatos[35]	Hyperbolic	$\frac{\cosh \left(\frac{k}{2}\right)z-\frac{h}{k}\sinh \left(k\frac{z}{h}\right)}{\cosh \left(\frac{k}{2}\right)-1}$

Table 2. Table 2: Nature of the original bHSDT and corresponding couple of functions f ​ ( z ) 𝑓 𝑧 f(z) and g ​ ( z ) 𝑔 𝑧 g(z) used to build the mHSDT.

Nature	$f\left(z\right)$	$g\left(z\right)$	Name
Polynomial	$z^3$	$z^2$	ToZZ
Trigonometric	$\sin \left(\pi \frac{z}{h}\right)$	$\cos \left(\pi \frac{z}{h}\right)$	SiZZ
Hyperbolic	$\sinh \left(k\frac{z}{h}\right)$	$\cosh \left(k\frac{z}{h}\right)$	HyZZ

Table 3. Table 3: Material properties ( α = 1 , 2 𝛼 1 2 \alpha=1,2 )

	$E_1$	$E_2$	$E_3$	$G_{23}$	$G_{13}$	$G_{12}$	${\nu }_{\alpha 3}$	${\nu }_{12}$	$\rho$
Composite (c)	$25E_2^c$	$E_2^c$	$E_2^c$	$0.2E_2^c$	$0.5E_2^c$	$0.5E_2^c$	$0.25$	$0.25$	${\rho }^c$
Honeycomb (h)	$E_2^h$	$E_2^c/25$	$12.5E_2^h$	$1.5E_2^h$	$1.5E_2^h$	$0.4E_2^h$	$0.02$	$0.25$	${\rho }^c/15$

Table 4. Table 4: Comparison between the Sin and the SiZZ models for the [ 0 / 90 ] n subscript delimited-[] 0 90 𝑛 [0/90]_{n} with n = 1 , 2 , 3 , 4 , 5 , 10 𝑛 1 2 3 4 5 10 n=1,2,3,4,5,10 ( a / h = 10 𝑎 ℎ 10 a/h=10 ).

Seq.

Model

$\overline{w}\left(C\right)$

%

${\overline{\sigma }}_{13}\left(B\right)$

%

${\overline{\sigma }}_{22}\left(C^{+}\right)$

%

$\overline{\omega }$

%

Equations

$$u_{\alpha }(z)$$

$$u_3(z)$$

$$u_{\alpha }(z)$$

$$u_3(z)$$

$${\phi }_{\alpha \beta }(0)=0$$

$${\epsilon }_{\alpha \beta }(z)$$

$${\epsilon }_{\alpha 3}(z)$$

$${\epsilon }_{33}(z)$$

$${\phi }_{\alpha \beta }^{\prime }(0^{-})={\delta }_{\alpha \beta }$$

$${\sigma }_{\alpha 3}(z)=C_{\alpha 3\beta 3}(z){\phi }_{\beta \gamma }^{\prime }(z){\gamma }_{\gamma 3}^{0^{-}}$$

$$C_{\alpha 3\gamma 3}({\zeta }_i^{-}){\phi }_{\gamma \beta }^{\prime }({\zeta }_i^{-})=C_{\alpha 3\gamma 3}({\zeta }_i^{+}){\phi }_{\gamma \beta }^{\prime }({\zeta }_i^{+})(i\in \{1,N-1\})$$

$$C_{\alpha 3\gamma 3}(-h/2){\phi }_{\gamma \beta }^{\prime }(-h/2)=0$$

$$C_{\alpha 3\gamma 3}(+h/2){\phi }_{\gamma \beta }^{\prime }(+h/2)=0$$

$${\phi }_{\delta \beta }^{\prime }(-h/2)=0$$

$${\phi }_{\delta \beta }^{\prime }(+h/2)=0$$

$${\phi }_{\alpha \beta }(0)=0$$

$${\phi }_{\alpha \beta }({\zeta }_i^{-})={\phi }_{\alpha \beta }({\zeta }_i^{+}),(i\in \{1,N-1\})$$

$${\phi }_{\alpha \beta }^{\prime }(0^{-})={\delta }_{\alpha \beta }$$

$$C_{\alpha 3\gamma 3}({\zeta }_i^{-}){\phi }_{\gamma \beta }^{\prime }({\zeta }_i^{-})=C_{\alpha 3\gamma 3}({\zeta }_i^{+}){\phi }_{\gamma \beta }^{\prime }({\zeta }_i^{+})(i\in \{1,N-1\})$$

$${\phi }_{\alpha \beta }^{\prime }(-h/2)=0$$

$${\phi }_{\alpha \beta }^{\prime }(+h/2)=0$$

$$(\varphi (z),\gamma (z),{\text{Z}}_i(z),{\text{1}}_i(z))$$

$$f(z)=-f(-z)$$

$$g(z)=g(-z)$$

$$g^{\prime }(\pm h/2)\ne 0$$

$${\phi }_{\alpha \beta }(z)=a_{\alpha \beta }f(z)+b_{\alpha \beta }g(z)+c_{\alpha \beta }^i{\text{Z}}_i(z)+d_{\alpha \beta }^i{\text{1}}_i(z)$$

$${\phi }_{\alpha \beta }^{\prime }(z)=a_{\alpha \beta }{f}^{\prime }(z)+b_{\alpha \beta }{g}^{\prime }(z)+c_{\alpha \beta }^i{\text{1}}_i(z)$$

$$d_{\alpha \beta }^{i_0}=0$$

$$c_{\alpha \beta }^i{\zeta }_i+d_{\alpha \beta }^i=c_{\alpha \beta }^{i+1}{\zeta }_i+d_{\alpha \beta }^{i+1}(i\in \{1,N-1\})$$

$$a_{\alpha \beta }{f}^{\prime }(0^{-})+b_{\alpha \beta }{g}^{\prime }(0^{-})+c_{\alpha \beta }^{i_0}={\delta }_{\alpha \beta }$$

$$C_{\alpha 3\gamma 3}^i\left(a_{\gamma \beta }{f}^{\prime }({\zeta }_i)+b_{\gamma \beta }{g}^{\prime }({\zeta }_i)+c_{\gamma \beta }^i\right)$$

$$=C_{\alpha 3\gamma 3}^{i+1}\left(a_{\gamma \beta }{f}^{\prime }({\zeta }_i)+b_{\gamma \beta }{g}^{\prime }({\zeta }_i)+c_{\gamma \beta }^{i+1}\right)(i\in \{1,N-1\})$$

$$a_{\alpha \beta }{f}^{\prime }(-h/2)+b_{\alpha \beta }{g}^{\prime }(-h/2)+c_{\alpha \beta }^1=0$$

$$a_{\alpha \beta }{f}^{\prime }(+h/2)+b_{\alpha \beta }{g}^{\prime }(+h/2)+c_{\alpha \beta }^N=0$$

$${\sigma }_{\alpha 3}(z)=C_{\alpha 3\beta 3}(z){\phi }_{\beta \gamma }^{\prime }(z){\gamma }_{\gamma 3}^{0^{-}}=4C_{\alpha 3\beta 3}(z){\phi }_{\beta \gamma }^{\prime }(z)S_{\gamma 3\delta 3}(0^{-}){\sigma }_{\delta 3}^0$$

$${\psi }_{\alpha \beta }(z)=4C_{\alpha 3\gamma 3}(z){\phi }_{\gamma \delta }^{\prime }(z)S_{\delta 3\beta 3}(0^{-})$$

$${\sigma }_{\alpha 3}(z)={\psi }_{\alpha \beta }(z){\sigma }_{\beta 3}^0$$

$$\textbf{P}_{\theta}=\left[\begin{array}[]{cc}c&s\\ -s&c\end{array}\right]\quad\text{with}\leavevmode\nobreak\ c=\cos(\theta)\leavevmode\nobreak\ \text{and}\leavevmode\nobreak\ s=\sin(\theta)$$

$$\left\{\begin{array}[]{l}\bar{\varphi}_{11}=\varphi_{11}c^{2}+\varphi_{22}s^{2}+(\varphi_{12}+\varphi_{21})sc\\ \bar{\varphi}_{12}=-(\varphi_{11}-\varphi_{22})sc+\varphi_{12}c^{2}-\varphi_{21}s^{2}\\ \bar{\varphi}_{21}=-(\varphi_{11}-\varphi_{22})sc-\varphi_{12}s^{2}+\varphi_{21}c^{2}\\ \bar{\varphi}_{22}=\varphi_{11}s^{2}+\varphi_{22}c^{2}-(\varphi_{12}+\varphi_{21})sc\end{array}\right.$$

$$\begin{array}{rl}A=(a/2,0,0);& A^{\prime }=(a/2,0,h/4)\\ B=(0,a/2,0);& B^{\prime }=(0,a/2,h/4)\\ C=(a/2,a/2,0);& C^{+}=(a/2,a/2,h/2)\end{array}$$

Full text

††journal: Composite Structures, accepted

1 Introduction

Among the numerous theories that have been developed for multilayered plates, those belonging to the Equivalent-Single Layer (ESL) family are of practical interest due to their relatively small number of unknowns, that is independent of the number of layers. Within this ESL class, the classical lamination theory (CLT) which has been proposed first, does not take into account the transverse shear behaviour and is, therefore, accurate only for thin plates, for which the transverse shear deformation can be neglected. First order shear-deformation theories (FSDT) have then been proposed to overcome this problem upon retaining a transverse shear deformation that is constant throughout the plate’s thickness. Its accuracy with regard to the plate’s gross response (transverse deflection, low vibration frequencies) results nevertheless dependent on shear correction factors. Higher-order shear deformation theories (HSDT) have been subsequently proposed in order to avoid the need of these problem-dependent shear correction factors. This is accomplished upon describing the through-thickness behaviour of the in-plane displacement field by means of functions of order greater than one, which thus introduces an enhanced description of the transverse shear deformation. The most well-known HSDT is the Vlasov–Levinson–Reddy’s third order theory (ToSDT). Polynomial functions are not the unique way to enrich the kinematic field, a wide variety of functions have been used, in particular trigonometric, hyperbolic, and exponential functions, as summarized in recent review papers [1, 2].

In all the theories cited above, the transverse shear deformation is included into the kinematic field by means of functions of class $C^{1}$ along the thickness direction $z$. This leads to continuous transverse shear strains and hence to a discontinuous transverse shear stress field, which violates equilibrium conditions within multilayered structures. Authors have thus proposed to use piecewise continuous, differentiable functions of $z$, often referred to as zig-zag functions. These functions are commonly constructed in such a manner that appropriate jumps of their derivatives at the interfaces restore the transverse stress continuity. Such Zig-Zag (ZZ) theories have appeared around the half of the twentieth century with Ambartsumyan [3, 4], and Osternik and Barg [5], and have been continuously receiving attention until now. An earlier paper that belongs to this category is due to Lekhnitskii [6], but it was limited to the study of beams. The approach by Lekhnitskii has been extended to plates 50 years later by Ren [7]. While the approaches by Ambartsumyan and Lekhnitskii rely on the exact verification of the constitutive equation connecting the transverse shear stress and the kinematic fields, Murakami [8] formulated a ZZ theory by postulating these two fields in an independent manner thanks to Reissner’s mixed variational theorem (RMVT) dedicated to multilayered plates [9]. For more recent developments of RMVT-based ZZ theories, the interested reader may refer to papers by Carrera [10], Demasi [11] and Tessler [12, 13]. Murakami’s zig-zag function (MZZF) has been also extensively applied to classical displacement-based variable kinematics approaches, see, e.g., [14, 15]. Within a comprehensive discussion about ZZ theories, Groh and Weaver have recently proposed a mixed ZZ theory based on Hellinger-Reissner’s principle [16].

Among these several approaches, more details of which can be found in the review papers [17, 18, 19], only a subset of these ZZ-theories are able to satisfy the appropriate interlaminar continuity (IC) of both, the displacement and the transverse shear stress fields. These two requirements (ZZ and IC) have been summarized by the acronym $C_{z}^{0}$–requirements [20]. In the following, we shall limit our attention to those ZZ theories that satisfy exactly the $C_{z}^{0}-$requirements. A more detailed examination is next proposed of some pioneering works, with the aim of establishing the background and highlighting the differences with the family of models proposed in this paper.

Ambartsumyan’s approach is based on the early paper [3], in which it is assumed that the transverse shear stress vary along $z$ according to a quadratic parabola with nil values at the outer surfaces. On page 20 of the later book [4], this assumption is formally expressed by the expressions $\tau_{xz}=f_{1}(z)\varphi(x,y)$ and $\tau_{yz}=f_{2}(z)\psi(x,y)$. It leads to 4 kinematic functions, but only two of them are independent, based on the primitives of $f_{1}(z)$ and $f_{2}(z)$. Conceived for orthotropic shells, the theory presented in the 1958 paper is not yet a ZZ theory, but in the 1970 book [4], page 75, an extension to symmetric multilayered orthotropic plates is given, which exhibits 4 zig-zag kinematic functions issued from two modified functions $f_{1}(z)$ and $f_{2}(z)$. This extension is perhaps due to Osternik and Barg [5], which is cited in the 1970 Ambartsumyan’s book (see also Carrera’s review paper [17]).

• 2.

In 1969, Whitney extends Ambrartsumyan’s theory to anisotropic plates, more specifically to general symmetric laminates and orthotropic non-symmetric laminates [21]. Following Ambartsumyan’s approach, Whitney starts from an assumed transverse shear stress field and ends up with 4 kinematic functions that are expressed as the superimposition of a polynomial third order ToSDT function and a zig-zag linear functions, see equation 6 of [21].

• 3.

Sun and Whitney propose in 1973 a layerwise model, which is the starting point for deriving an ESL model upon eliminating the parameters of the $N-1$ upper layers [22]. The resulting ESL model is equivalent to a first order zig-zag model with 4 kinematic functions. The link between such models has been discussed in detail in [23]. It has constant transverse shear stresses, which is the drawback of first-order models.

• 4.

The 1986 paper by Ren [7] proposes 4 kinematic functions from a priori given transverse shear stress functions and applies the resulting zig-zag model to cross-ply laminates. Four displacement unknowns are introduced to take into account the transverse shear behaviour, which yields a 7-parameter model which is difficult to compare with our 5-parameter models.

• 5.

Cho and Parmerter formulate a ZZ theory for symmetric [24] and general [25] orthotropic composite plates. Starting point are 2 kinematic functions which are the superimposition of a cubic polynomial and a linear zig-zag function expressed in terms of the Heaviside function. The coefficients are determined by enforcing the transverse stress continuity, which leads to the coupling of the $x-$ and $y-$directions in the kinematic field. It therefore appears that the theory is in fact based on 4 kinematic functions (see equation 4 of [25]). In his historical review, Carrera demonstrated the equivalence of Cho and Parmerter’s model with Ambartsumyan’s model.

• 6.

It is finally worth mentioning the ZZ theory fulfilling the $C^{0}_{z}$–requirements that is based on trigonometric functions and developed by Ossadzow and coworkers [26, 27]. The construction of the kinematic zig-zag functions follows a similar path as proposed by Cho and Parmerter, but the trigonometric $\cos$ and $\sin$ functions replace the quadratic and cubic terms of the polynomial expansion, respectively.

A conforming finite element based on a trigonometric ZZ theory enhanced through a transverse normal strain has been developed for laminated plates [28]. References [29, 30] extend the procedure of Cho and Parmerter to a wider family including polynomial, trigonometric, exponential and hyperbolic functions. However, the authors only consider two kinematic functions, which reduces the applicability of their models to cross-ply laminates.

In [31, 32], corresponding polynomial and trigonometric $C^{0}_{z}$ zig-zag models have been constructed from basis polynomial and trigonometric HSDT by using four functions, which shall be hereafter referred to as warping functions. While in these works the warping functions were obtained from transverse stress fields obtained from 3D solutions or from equilibrium equations, the present paper presents a general procedure for extending a basis higher-order shear deformation theory (bHSDT) to a multilayer higher-order shear deformation theory (mHSDT) that meets the $C_{z}^{0}-$requirements. This extension consists in the construction of four $C^{0}_{z}$ warping functions starting from the native functions that characterize the basis theory: it can be applied to any couple of odd and even functions and has no limitation concerning the lamination scheme.

The paper is organized as follows. Section 2 introduces the notation and points out the properties that the four warping functions are required to fulfil. The extension of a bHSDT up to an mHSDT fulfilling all $C^{0}_{z}$–requirements is described in Section 3. Three different basis models are exemplarily considered, which pertain respectively to the polynomial, trigonometric and hyperbolic type. It is also shown that the warping functions are components of a second-order tensor, hence being covariant with rotations about the $z-$axis. Section 4 reports the numerical evaluations: the $C^{0}_{z}$ warping functions effectively increase the accuracy of the basis (non zig-zag) model and this enhancement is quite insensitive with respect to the type of functions used for the model. A discussion is proposed in Section 5 in order to substantiate the limitations of conventional ZZ models with respect to particularly “constrained” configurations with very low number of layers and length-to-thickness ratios: in these cases, accuracy may only be assured by resorting to warping functions that contain more layer-specific information, just as LayerWise models do. The main conclusions are finally summarised in Section 6.

2 Definitions and general properties

This paper deals with a generic method to extend a basis higher-order shear deformation theory (bHSDT) to a multilayer higher-order shear deformation theory (mHSDT). This Section introduces the notation employed for identifying the various plate theories along with the fundamental properties that the underlying approximating functions are required to satisfy.

2.1 The basis theories

We consider a basis high-order shear deformation theory (bHSDT) for which the kinematic field can be written in the following general form:

[TABLE]

where $u^{0}_{\alpha}$ are the membrane displacements at $z=0$, $w^{0}$ is the deflection at $z=0$, $\gamma^{0^{-}}_{\alpha 3}$ are the transverse engineering strains at $z=0^{-}$, and $\phi(z)$ is a $C^{1}$ odd function. The choice of the $0^{-}$ coordinate is a convention useful to avoid undetermined shear strains if an interface lies at $z=0$. The $z=0$ plane is assumed to be the middle plane of the plate, the lower and upper faces are respectively located at $z=-h/2$ and $z=h/2$.

Written in the form reported in equation (1), the function $\phi(z)$ must verify $\phi(0)=0$ and $\phi^{\prime}(0^{-})=\phi^{\prime}(0)=1$ to give sense to the notations, and $\phi^{\prime}(\pm h/2)=0$ to enforce null transverse shear stresses at the top and bottom of the plate. Due to the $C^{1}$ property of $\phi(z)$, such theories do not have particular abilities to deal with multilayered plates. Indeed, the continuity of $\phi^{\prime}(z)$ induces discontinuities of the transverse shear stresses $\sigma_{\alpha 3}$ at the interfaces. Table 1 summarizes the functions $\phi(z)$, along with the reference author and the type of the approximation, that will be extended to a ZZ model in Section 3.

2.2 The multilayer HSDT

A multilayer theory is a plate theory which is dedicated to composite plates upon fulfilling the $C^{0}_{z}$–requirements. The generic expression for the kinematics of a multilayer HSDT (mHSDT) is of the form:

[TABLE]

where $\varphi_{\alpha\beta}(z)$ are four piecewise $C^{1}$ functions, sometimes called warping functions, that are requested to fulfil specific properties, as it will be discussed below. Among these properties, specific jump values need to be prescribed to their derivatives for enforcing continuity of transverse stresses at the interfaces. In order to construct an mHSDT that is applicable to arbitrary laminates, it is important to consider four functions $\varphi_{\alpha\beta}(z)$, and hence to retain the coupling between $\gamma^{0^{-}}_{13}$ (resp. $\gamma^{0^{-}}_{23}$) and $u_{2}$ (resp. $u_{1}$). In fact, models written with only two functions, i.e., with $\varphi_{12}(z)=\varphi_{21}(z)=0$, are only applicable to cross-ply laminates.

2.3 Required properties for the functions $\varphi_{\alpha\beta}(z)$

Formula (2a) shows that the piecewise functions $\varphi_{\alpha\beta}(z)$ must be continuous at each interface to respect the continuity of the in-plane displacements. Since $u_{\alpha}^{0}$ denotes the membrane displacements at $z=0$, the functions $\varphi_{\alpha\beta}$ must fulfil the following homogeneity condition

[TABLE]

Figure 1 illustrates the continuity and homogeneity conditions for a practical example (a $[-10/0/40]$ laminate, trigonometric mHSDT).

The compatible strain field defined by generic mHSDT kinematic field of Eq. (2) reads

[TABLE]

The transverse shear strains $\varepsilon_{\alpha 3}(z)$ must be defined in each layer, but they also should be discontinuous at the interfaces for allowing the transverse shear stresses $\sigma_{\alpha 3}(z)$ to be continuous in order to fulfil the equilibrium condition. Indicating by $\zeta_{i}$ the $z$-coordinate of the $i$th interface, with $i=1,2\ldots N-1$, $\zeta_{0}=-h/2$ and $\zeta_{N}=h/2$, the functions $\varphi_{\alpha\beta}(z)$ are thus required to be piecewise $C^{1}$ over the intervals $]\zeta_{i-1},\zeta_{i}[$. Furthermore, since $\gamma^{0^{-}}_{\beta 3}$ represents the engineering shear strain at $z=0^{-}$, the derivatives of the functions $\varphi_{\alpha\beta}(z)$ are required to fulfil the following homogeneity conditions:

[TABLE]

Figure 2 illustrates these conditions with the same practical example as before.

Due to the continuity of the $\varphi_{\alpha\beta}(z)$, in-plane strains $\varepsilon_{\alpha\beta}(z)$ are continuous. Following the classical plate approach, the normal stress $\sigma_{33}(z)$ is set to [math], which leads to the use of reduced (in-plane) stiffnesses $Q_{\alpha\beta\gamma\delta}(z)$ in place of the stiffnesses $C_{\alpha\beta\gamma\delta}(z)$. It is further recalled that there is no physical reason for the in-plane stresses $\sigma_{\alpha\beta}(z)$ to be continuous at interfaces between adjacent layers with dissimilar stiffness coefficients.

The constitutive equation defines the transverse shear stresses in terms of strains as follows:

[TABLE]

For equilibrium reasons, these transverse stresses need to be continuous at the interfaces, and also to be null at $z=\pm h/2$ if the applied load remains normal to the plate. These conditions are expressed as

[TABLE]

We can pre-multiply equations (7b) and (7c) with the compliance tensors $S_{\delta 3\alpha 3}(-h/2)$ and $S_{\delta 3\alpha 3}(+h/2)$, respectively, to obtain:

[TABLE]

Figure 3 illustrates the top and bottom nullity conditions and the jump conditions prescribed on the $\varphi^{\prime}_{\alpha\beta}(z)$ functions. The following relations summarize all the properties that the functions $\varphi_{\alpha\beta}(z)$ are required to verify:

[TABLE]

According to the free indexes in the above formulas, we can see that there are $4+4(N-1)+4+4(N-1)+4+4=8(N+1)$ conditions for the four functions $\varphi_{\alpha\beta}$.

3 The extension process

This Section describes the procedure for extending a generic bHSDT to a corresponding mHSDT. The construction of the $C_{z}^{0}$ warping functions is described and their tensorial character highlighted.

3.1 Construction of the four $\varphi_{\alpha\beta}(z)$ functions

Given a composite plate consisting of $N$ layers, the goal is to find four $\varphi_{\alpha\beta}(z)$ functions that obey to all the properties summarized in Eq. (9). Observing that the $\phi(z)$ function of the bHSDT is an odd function, one possibility would be to use it directly and to merely find an even $\gamma(z)$ function with suitable properties, in order to build the four warping functions from the basis spanned by the following $2N+2$ elements

[TABLE]

$\textbf{Z}_{i}(z)$ and $\textbf{1}_{i}(z)$ represent the restrictions on the interval $]\zeta_{i-1},\zeta_{i}[$ of the linear and the constant (unitary) functions, respectively. Instead of $\phi(z)$ and $\gamma(z)$, the method proposed here uses more general and less constrained $f(z)$ and $g(z)$ functions, and hence is easier to use. The link that remains between the bHSDT and the corresponding mHSDT is the nature of the functions that will be used to form the basis (polynomial, trigonometric, hyperbolic…). In any case, these high-order functions are responsible for tailoring the transverse shear deformation, while the constant and linear elements introduce the characteristic zig-zag distribution of the in-plane displacements.

We need two functions of class $C^{1}$: an odd function $f(z)$, and an even function $g(z)$ verifying $g^{\prime}(\pm h/2)\neq 0$, viz.:

[TABLE]

Now consider the four piecewise functions:

[TABLE]

where summation is implied over the dummy index $i=1,2,\ldots N$. These four functions are defined with $4(2+2N)=8(N+1)$ constants. The expression for the derivatives of the four functions is

[TABLE]

Just as the transverse shear strains, these four derivatives are not defined at the $N-1$ interfaces. Although Dirac’s delta function might be used to formally write these derivatives, this is not useful because the relation in Eq. (9d) only involves their values at the layers’ limits. Substituting Eqs. (12) and (13) into Eqs. (9) yields the following system of $8(N+1)$ equations:

[TABLE]

The index $i_{0}$ corresponds to the number of the layer which contains the $z=0^{-}$ coordinate. Since it seems difficult to formulate a recursive process to determine all the coefficients, the linear system (14) is solved for the $8(N+1)$ unknown coefficients $a_{\alpha\beta},b_{\alpha\beta},c_{\alpha\beta}^{i},d_{\alpha\beta}^{i}$ $(i=1,2\ldots N)$.

Table 2 reports some functions $f(z)$ and $g(z)$ that can be chosen to build an mHSDT model. While these functions allow to accommodate the transverse shear behaviour inside each layer, the linear and constant contributions are responsible for the ZZ behaviour, that is the respect of displacement and transverse stress continuities at the layers’ interfaces. It should be noted that a “mixed” model can be constructed by using functions of different nature, for example the hyperbolic odd function $\sinh(kz/h)$ can be considered in conjunction with the even trigonometric function $\cos(\pi z/h)$. Analytical expressions for the warping functions for a single-layer plate are reported explicitly in Appendix.

3.2 Computational aspects

The construction of the linear system (14) can be automated because its structure does not depend on the choice of the functions $f(z)$ and $g(z)$. Indeed, only few values of these functions and of their derivatives, taken at specific $z$ coordinates, have to be sent to the routine. The solution of the $(8N+8)\times(8N+8)$ linear algebraic system can be carried out with a classical algorithm and provides the coefficients defining the four warping functions $\varphi_{\alpha\beta}(z)$. It may be noted that either numerical or semi-analytical versions of the warping functions can be used. Numerical versions, which consist on a sufficiently dense table of values, are more suitable for computing the numerous generalized stiffness and mass terms of the plate model within a numerical quadrature scheme.

3.3 Stress functions

Once the $\varphi_{\alpha\beta}(z)$ functions are built, one can compute the corresponding stress functions $\psi_{\alpha\beta}(z)$. They do not bring new information to the models, as these stress functions are a direct consequence of the warping functions, but they are useful to illustrate and understand the static response of the mHSDT. Let us replace in equation (6), the middle-plane transverse strains by the corresponding middle-plane stresses:

[TABLE]

Note that the $0^{-}$ exponent of $\gamma^{0^{-}}_{\gamma 3}$ is not required to appear in the interlaminar continuous stress $\sigma^{0}_{\delta 3}$, but it is found in the $S_{\gamma 3\delta 3}(0^{-})$ term. Eq. (15) permits to define the 4 stress functions of the model

[TABLE]

through which the transverse shear stresses are expressed as

[TABLE]

3.4 Tensorial character of the $\varphi_{\alpha\beta}(z)$

The tensorial character of the $\varphi_{\alpha\beta}(z)$ functions follows directly from their definition, see Eq. (2). This tensorial character concerns only the 2D $(x,y)$ space. Also, the equations of the system (9) are tensor equations, i.e., their form is invariant with respect to rotations about the $z$ axis. It implies that all the coefficients $a_{\alpha\beta}$, $b_{\alpha\beta}$, $c^{i}_{\alpha\beta}$, $d^{i}_{\alpha\beta}$ are second order tensors and must, therefore, obey to the formulas of coordinate transformation for second order tensors.

The tensorial character of the warping functions implies that the four functions $\varphi_{\alpha\beta}(z)$ of a laminate whose lamination sequence is $s=[\theta_{1}/\theta_{2}/\dots/\theta_{N}]$ must be linked to the four functions $\bar{\varphi}_{\alpha\beta}(z)$ of the laminate whose stacking sequence is $\bar{s}=s+\theta=[(\theta_{1}+\theta)/(\theta_{2}+\theta)/\dots/(\theta_{N}+\theta)]$.

In order to identify this relation, let us consider the $\bar{s}$–laminate, in a Cartesian frame $(x,y,z)$, and suppose it undergoes a pure shear deformation of its middle plane $\bar{\boldsymbol{\upgamma}}$. In this case, the kinematic field of Eq. (2a) can be written $\bar{u}_{\alpha}(z)=\bar{\varphi}_{\alpha\beta}(z)\bar{\gamma}^{0}_{\beta 3}$ or, in matrix notation, $\bar{\textbf{u}}=\bar{\boldsymbol{\upvarphi}}\bar{\boldsymbol{\upgamma}}$. Consider the matrix of change of coordinates $\textbf{P}_{\theta}$ from the Cartesian frame $(x,y,z)$ to a Cartesian frame $(x^{\prime},y^{\prime},z)$, rotated from the previous one by an angle $\theta$ about the $z$-axis. In the rotated frame, this shear strain is $\boldsymbol{\upgamma}=\textbf{P}_{\theta}\bar{\boldsymbol{\upgamma}}$ and it “acts” on the $s$–laminate producing the in-plane kinematic field u. In the original frame, the $\bar{s}$–laminate is then subjected to the kinematic field $\bar{\textbf{u}}=\textbf{P}_{\theta}^{-1}\textbf{u}$. Therefore, the following relation is established: $\bar{\boldsymbol{\upvarphi}}=\textbf{P}_{\theta}^{-1}\boldsymbol{\upvarphi}\textbf{P}_{\theta}$. Since the transformation matrix is

[TABLE]

one can compute the warping functions for the $\bar{s}$–laminate directly from those for the $s$–laminate according to

[TABLE]

Two examples of such transformations are given next for illustration purposes. Figure 4 compares the native warping functions of a $[45/{-}45]$ laminate against those obtained from a $[0/90]$ laminate after rotating them by an angle of $45^{\circ}$. The same comparison is proposed in figure 5 for the two laminates $[-25/{-}15/25]$ and $[-10/0/40]$ and with a rotation of $-15^{\circ}$.

4 Numerical results

A numerical evaluation is proposed in order to assess the accuracy of the basis models and their corresponding enhancement through ZZ warping functions with respect to the plate’s length-to-thickness ratio, number of layers, and stacking sequence. All bHSDT listed in Table 1 are compared with their corresponding enhancements defined by the functions listed in Table 2. The factor $k$ in the hyperbolic functions has been set equal to 2 in the subsequent numerical investigations. Note finally that all considered models have the same number of DOF as the bHSDT, i.e., 5 DOF.

In order to encompass a quite broad range of stiffness mismatch between adjacent layers, the study will investigate laminated as well as sandwich plates with composite skins and a honeycomb core. The material properties used for the composite and the honeycomb layers are reported in Table 3.

The numerical assessment of the different models is carried out by referring to an exact solution of the 2D differential equations governing the plate bending problem. Square, simply-supported plates are considered, for which we compute the fundamental eigenfrequency as well as the static response under bi-sinusoidal transverse pressure loads of amplitude $q/2$ acting at the top and the bottom surfaces of the plate. A Navier-type strong-form solution is found for all considered laminates, where the simply-support conditions of arbitrary, non cross-ply laminates are opportunely relaxed as discussed in [36]. The warping functions and stress functions of the bHSDT and mHSDT are compared against those that have been extracted from the 3D solutions following the procedure detailed out in [31].

The points at which quantities are output are defined as follows:

[TABLE]

Deflections $w$, first natural frequencies $\omega$ and stresses are given according to following adimensionalisation

[TABLE]

where for sandwich plates $E_{2}^{\text{ref}}$ and $\rho^{\text{ref}}$ are the values of the core material. It is important to specify that the transverse shear stress values reported in the tables and their distributions across the plate thickness plotted in the figures are obtained from the equilibrium equations upon integrating the in-plane stresses.

4.1 The $[0/90]_{n}$ laminates

The models are tested for the antisymmetric cross-ply laminates $[0/90]_{n}$, where different numbers of layers are considered with $n=1,2,3,4,5,10$. In table 4.1, non-dimensional deflection, transverse shear stresses and fundamental frequency are given for the Sin and the SiZZ models, and compared to the exact solution. The length-to-thickness ratio is set to $a/h=10$. Very similar results are obtained with polynomial and hyperbolic bHSDT/mHSDT models and are omitted from Table 4.1 for the sake of clarity. The results clearly shows the accuracy improvement introduced by the $C^{0}_{z}$ warping functions, in particular for the deflection and the fundamental frequency: the enhancement on these two quantities appears to decrease as the number of layers increases, although for $n=10$ it is still larger than 5% and 3%, respectively.

Only two warping functions are required for a cross-ply laminate because the cross-coupling functions are identically nil, $\varphi_{12}(z)=\varphi_{21}(z)=0$. The functions $\varphi_{11}(z)$ and $\varphi_{22}(z)$ of the polynomial bHSDT (ToSDT) and of the polynomial, trigonometric and hyperbolic mHSDT are compared in Figures 6 and 9 for the $n=2$ and $n=4$ configurations, respectively. The differences between the 3 mHSDT are seen to be negligible, and the curves for the trigonometric and hyperbolic bHSDT have been omitted for the sake of clarity because they are practically coincident with those of the ToSDT.

As far as the impact of warping functions on the local stress response is concerned, the values in Table 4.1 for the transverse shear stresses at the selected points do not allow to well appreciate it, but their through-the-thickness distributions obtained with the extended mHSDT model are closer to the exact solution in comparison to the bHSDT models. This can be seen by comparing the two stress functions $\psi_{11}(z)$ and $\psi_{22}(z)$ depicted in Figures 7 and 10 for the cases $n=2$ and $n=4$, respectively. On the other hand, Figures 8 and 11 report the transverse shear stress distributions computed for the two configurations $n=2$ and $n=4$, respectively, upon integrating the equilibrium equations starting from the in-plane stresses. This post-processing procedure is seen to annihilate all differences between the bHSDT and the mHSDT, thus providing distributions that very accurately recover the exact 3D solution.

All considered models for the cross-ply laminate $[0/90]_{2}$ are assessed in Table 4.1 with respect to the length-to-thickness ratio $a/h$. It is obvious that the improvement of the mHSDT over the bHSDT is more important for thick plates than for thin plates, it decreases from more than 10% for $a/h=4$ to about 0.1% for $a/h=100$. This is not surprising as it is well known that the effect of the transverse shear increases as $a/h$ diminishes.

It is worthwhile to make some comments about an expected symmetry for the warping functions of the considered antisymmetric cross-ply laminates. Indeed, one should expect the $\varphi_{11}(z)$ functions to be equal to the corresponding $\varphi_{22}(-z)$ functions, but figures 6 and 9 show that this is not the case. This is due to the fact that the $z=0^{-}$ coordinate has been chosen for prescribing the $\varphi^{\prime}_{11}(0^{-})=\varphi^{\prime}_{22}(0^{-})=1$ conditions. This choice “hides” the expected property, which can nevertheless be easily restored: dividing $\varphi_{22}(-z)$ by $\varphi^{\prime}_{22}(0^{+})$ yields in fact exactly the function $\varphi_{11}(z)$. Note that, since there are no such constraints on the stress functions, the symmetry $\psi_{11}(z)=\psi_{22}(-z)$ is immediately apparent in Figures 7 and 10.