Simple and accurate two-noded beam element for composite laminated beams using a refined zigzag theory

Published in Comput. Methods Appl. Mech. Engrg. Vol. 213–216, pp. 362–382, 2012
doi: 10.1016/j.cma.2011.11.023

Abstract

In this work we present a new simple linear two-noded beam element adequate for the analysis of composite laminated and sandwich beams based on the combination of classical Timoshenko beam theory and the refined zigzag kinematics proposed by Tessler et al. [22]. The element has just four kinematic variables per node. Shear locking is eliminated by reduced integration. The accuracy of the new beam element is tested in a number of applications to the analysis of composite laminated beams with simple supported and clamped ends under point loads and uniformly distributed loads. An example showing the capability of the new element for accurately reproducing delamination effects is also presented.

Keywords: Two-noded beam element, zigzag kinematics, Timoshenko theory, composite, sandwich beams

1 INTRODUCTION

It is well known that both the classical Euler-Bernouilli beam theory [1] and the more advanced Timoshenko theory [2] produce inadequate predictions when applied to relatively thick composite laminated beams with material layers that have highly dissimilar stiffness characteristics. Even with a judiciously chosen shear correction factor, Timoshenko theory tends to underestimate the axial stress at the top and bottom outer fibers of a beam. Also, along the layer interfaces of a laminated beam the transverse shear stresses predicted exhibit erroneous discontinuities. These difficulties are due to the higher complexity of the “true” variation of the axial displacement field across a highly heterogeneous beam cross-section.

Indeed to achieve accurate computational results, 3D finite element analyses are often preferred over beam models. For composite laminates with hundred of layers, however, 3D modelling becomes prohibitely expensive, specially for non linear and progressive failure analyses.

Improvements to the classical beam theories have been obtained by the so called equivalent single layer (ESL) theories that assume a priori the behaviour of the displacement and/or the stress through the laminate thickness [3,4]. Despite of being computational efficient, ESL theories often produce innacurate distributions for the stresses and strains (in particular the transverse shear stress) across the thickness.

The need for composite laminated beam theories with better predictive capabilities has led to the development of the so-called higher order theories. In these theories higher-order kinematic terms with respect to the beam depth are added to the expression for the axial displacement and, in some cases, to the expressions for the deflection. A review of these theories can be found in [3,4].

Accurate predictions of the correct shear and axial stresses for thick and highly heterogenous composite laminated and sandwich beams can be obtained by using layer-wise theory. In this theory the thickness coordinate is split into a number of analysis layers that may or not coincide with the number of laminate plies. The kinematics are independently described within each layer and certain physical continuity requirements are enforced [3,4].

A drawback of layer-wise theory is that the number of kinematic variables depends on the number of analysis layers. The layer displacements can be condensed at each section in terms of the axial displacement for the top layer during the equation solution process [5,6]. The displacement condensation processes can be however expensive for problems involving many analysis layers.

Discrete layer theories in which the number of unknowns in the model does not depend on the number of layers in the laminate are described in [7,8,9]. In this class of discrete layerwise theories (called zigzag theories) a piewise linear in-plane displacement function (the zigzag function) is superimposed over a linear displacement field [7,8], a quadratic displacement field [10,11] or a cubic displacement field [12,13,14] through the thickness of the laminate.

Many zigzag theories require ${\textstyle C^{1}}$ continuity for the deflection field, which is a drawback versus simpler ${\textstyle C^{\circ }}$ continuous FEM approximations. Also many zigzag theories run into theoretical difficulties to satisfy equilibrium of forces at a clamped support.

Averill et al. [9,10,16,17] developed linear, quadratic and cubic zigzag beam theories that overcame the need for ${\textstyle C^{1}}$ continuity. The shear strain angle is introduced as a kinematic variable together with the deflection, the rotation and the zigzag function. A ${\textstyle C^{\circ }}$ interpolation can be used for all these variables. The relationship between the shear angle, the deflection and the rotation of each layer is introduced as a constraint via a penalty method. This also ensures the continuity of the transverse shear stress across the laminate depth and the satisfaction of the shear traction boundary conditions. However, Averill theories have difficulties to model correctly clamped boundary conditions. For this reason, analytical and numerical (FEM) studies based on Averill theory have mainly focused on simple supported beams [16,17].

A 2-noded beam element based on Euler-Bernouilli beam theory and an extension of Averill's zigzag theory including a cubic in-plane displacement field within each layer has been recently proposed by Alam and Upadhyay [18]. Good results are reported for cantilever and clamped composite and sandwich beams.

An assessment of different zigzag theories for beam is reported in [19,20]. A review of zigzag theory for plate analysis can be found in [21].

Tessler et al. [22,23] have developed a refined zigzag (RZ) theory starting from the standard Timoshenko kinematic assumptions. This allows one using ${\textstyle C^{\circ }}$ continuous interpolation for all the kinematic variables. Timoshenko beam theory also introduces naturally shear deformation effect for the homogeneous material case, which is advantageous for many problems. The zigzag functions chosen have the property of vanishing on the top and bottom surfaces of a laminate. A particular feature of this zigzag theory is that the transverse shear stresses are not required to be continuous at the layer interfaces. This results in simple piewice-constant functions that approximate the true shear stress distribution. An accurate continuous thickness distribution of the transverse shear stress can be obtained “a posteriori” in terms of the axial stress by integrating the equilibrium equations. This theory also provides good results for clamped supports.

Gherlone et al. [24] have developed two and three-noded ${\textstyle C^{\circ }}$ beam elements based on the RZ theory for analysis of multilayered composite and sandwich beams. Locking-free elements are obtained by using anisoparametric interpolations that are adapted to approximate the four independent kinematic variables that model the beam deformation. A family of beam elements is achieved by imposing different constraints on the original displacement approximation. The constraint conditions requiring a constant variation of the transverse shear force provide an accurate 2-noded beam element [24].

Quite simultaneously to the above work, Oñate et al. [25] proposed a simple 2-noded beam element for composite laminated beams based on the RZ theory. A standard linear displacement field is used to model the four variables of the so called LRZ element. Shear locking is avoided by using reduced integration on selected terms of the shear stiffness matrix.

In this paper we present in detail the formulation of the LRZ beam element originally reported in [25] and explore the capabilities of the new element for multilayered beams and delamination analysis. A study of the locking-free behaviour of the LRZ element for slender beams is presented. The good performance of the element is demonstrated for simply supported and clamped composite laminated beams with different layers under point load and uniformly distributed loads.

Finally, an example showing the capability of the LRZ element to model delamination effects is presented.

2 GENERAL CONCEPTS OF ZIGZAG BEAM THEORY

The kinematic field in zigzag beam theory is generally written as

$u^{k}(x,z)=u_{0}(x)-z\theta (x)+{\bar {u}}^{k}(x,z)$ $w(x,z)=w_{0}(x)$

(1.a)

where

${\bar {u}}^{k}=\phi ^{k}(z)\Psi (x)\quad ,\quad k=1,N$

(1.b)

is the zigzag displacement function (Figure 1).

In Eqs.(1) ${\textstyle N}$ is the number of layers, superscript ${\textstyle k}$ indicates quantities within the ${\textstyle k}$th layer with ${\textstyle z_{k}\leq z\leq z_{k+1}}$ and ${\textstyle z_{k}}$ is the vertical coordinate of the ${\textstyle k}$th interface. In Eq.(1a) the uniform axial displacement ${\textstyle u_{0}(x)}$, the rotation ${\textstyle \theta (x)}$ and the transverse deflection ${\textstyle w_{0}(x)}$ are the primary kinematic variables of the underlying equivalent single-layer Timoshenko beam theory. In Eq.(1b) function ${\textstyle \phi ^{k}(z)}$ denotes a piecewise linear zigzag function, yet to be established, and ${\textstyle \Psi (x)}$ is a primary kinematic variable that defines the amplitude of the zigzag function along the beam. Collectively, the interfacial axial displacement field has a zigzag distribution, as shown in Figure 1c.

Figure 1: Thickness distribution of the the zigzag function $\phi ^{k}$ (a), the zigzag displacement function ${\bar {u}}^{k}$ (b), and the axial displacement (c) in refined zigzag beam theory

The strain-displacement relations are derived by substituting Eq.(1a) into the expressions of classical Timoshenko beam theory, i.e.

$\varepsilon _{x}^{k}={\frac {\partial u^{k}}{\partial x}}={\frac {\partial u_{0}}{\partial x}}-z{\frac {\partial \theta }{\partial x}}+\phi ^{k}{\frac {\partial \Psi }{\partial x}}=[1,-z,\phi ^{k}]\left\{{\begin{matrix}\displaystyle {\frac {\partial u_{0}}{\partial x}}\\[.3cm]\displaystyle {\frac {\partial \theta }{\partial x}}\\[.3cm]\displaystyle {\frac {\partial \Psi }{\partial x}}\end{matrix}}\right\}={\boldsymbol {S}}_{p}{\hat {\boldsymbol {\varepsilon }}}_{p}$

(2.a)

$\gamma _{xz}^{k}={\frac {\partial u^{k}}{\partial z}}+{\frac {\partial w}{\partial x}}={\frac {\partial w_{0}}{\partial x}}-\theta +{\frac {\partial \phi ^{k}}{\partial z}}\Psi =\gamma +\beta ^{k}\Psi =[1,\beta ^{k}]\left\{{\begin{matrix}\gamma \\\Psi \end{matrix}}\right\}={\boldsymbol {S}}_{t}{\hat {\boldsymbol {\varepsilon }}}_{t}$

(2.b)

In Eq.s (2a) and (2b)

${\begin{array}{ll}{\boldsymbol {S}}_{p}=[1,-z,\phi ^{k}]\quad ,&{\hat {\boldsymbol {\varepsilon }}}_{p}=\displaystyle \left[{\frac {\partial u_{0}}{\partial x}},{\frac {\partial \theta }{\partial x}},{\frac {\partial \Psi }{\partial x}}\right]^{T}\\{\boldsymbol {S}}_{t}=[1,\beta ^{k}]\quad ,&{\hat {\boldsymbol {\varepsilon }}}_{t}=\left[\gamma ,\Psi \right]^{T}\end{array}}$

(2.c)

where ${\textstyle {\hat {\boldsymbol {\varepsilon }}}_{p}}$ and ${\textstyle {\hat {\boldsymbol {\varepsilon }}}_{t}}$ are the generalized in-plane and transverse shear strain vectors, respectively. Vector ${\textstyle {\hat {\boldsymbol {\varepsilon }}}_{p}}$ contains the axial elongation ${\textstyle \left({\frac {\partial u_{0}}{\partial x}}\right)}$, the pseudo-curvature ${\textstyle \left({\frac {\partial \theta }{\partial x}}\right)}$ and the derivatives of the amplitude of the zigzag function ${\textstyle \left({\frac {\partial \Psi }{\partial x}}\right)}$. In ${\textstyle {\hat {\boldsymbol {\varepsilon }}}_{t}}$, ${\textstyle \gamma ={\frac {\partial w_{0}}{\partial x}}-\theta }$ is the average transverse shear strain of Timoshenko beam theory and ${\textstyle \beta ^{k}={\frac {\partial \phi ^{k}}{\partial z}}}$. Note that since ${\textstyle \phi ^{k}(z)}$ is piecewise linear, ${\textstyle \beta ^{k}}$ is constant across each layer.

For major principal material axes that are coincident with the beam ${\textstyle x}$ axis, Hooke stress-strain relations for the ${\textstyle k}$th orthotropic layer have the standard form

$\sigma _{x}^{k}=E^{k}\varepsilon _{x}^{k}=E^{k}{\boldsymbol {S}}_{p}{\hat {\boldsymbol {\varepsilon }}}_{p}$ (3.a) $\tau _{xz}^{k}=G^{k}\gamma _{xz}^{k}=G^{k}{\boldsymbol {S}}_{t}{\hat {\boldsymbol {\varepsilon }}}_{t}$ (3.b)

where ${\textstyle E^{k}}$ and ${\textstyle G^{k}}$ are the axial and shear moduli for the ${\textstyle k}$th layer, respectively.

In the above equations we have distinguished all variables within a layer with superscript ${\textstyle k}$.

3 REFINED ZIGZAG THEORY

3.1 Zigzag kinematics

The key attributes of the refined zigzag (RZ) theory proposed by Tessler et al. [22] are, first, the zigzag function vanishes at the top and bottom surfaces of the beam section and does not require full shear-stress continuity across the laminated-beam depth. Second, all boundary conditions can be modelled adequately. And third, ${\textstyle C^{\circ }}$ continuity is only required for the FEM approximation of the kinematic variables.

Within each layer the zigzag function is expressed as

$\phi ^{k}={\frac {1}{2}}(1-\zeta ^{k}){\bar {\phi }}^{k-1}+{\frac {1}{2}}(1+\zeta ^{k}){\bar {\phi }}^{k}={\frac {{\bar {\phi }}^{k}+{\bar {\phi }}^{k-1}}{2}}+{\frac {{\bar {\phi }}^{k}-{\bar {\phi }}^{k-1}}{2}}\zeta ^{k}$

(4)

where ${\textstyle {\bar {\phi }}^{k}}$ and ${\textstyle {\bar {\phi }}^{k-1}}$ are the zigzag function values of the ${\textstyle k}$ and ${\textstyle k-1}$ interface, respectively with ${\textstyle {\bar {\phi }}^{0}={\bar {\phi }}^{N}=0}$ and ${\textstyle \zeta ^{k}={\frac {2(z-z^{k-1})}{h^{k}}}-1}$.

Collectively, function ${\textstyle \phi ^{k}}$ has the zigzag distribution shown in Figure 1a. Due to the dependence between the zigzag displacement function ${\textstyle {\bar {u}}^{k}}$ and ${\textstyle {\bar {\phi }}^{k}}$ (see Eq.(1b)), ${\textstyle {\bar {u}}^{k}}$ also vanishes at the top and bottom layers. The axial displacement field is plotted in Figure 1c.

The above form of ${\textstyle \phi ^{k}}$ gives the constant value of ${\textstyle \beta ^{k}}$ for each layer as

$\beta ^{k}={\frac {\partial \phi ^{k}}{\partial z}}={\frac {{\bar {\phi }}^{k}-{\bar {\phi }}^{k-1}}{h^{k}}}$

(5.a)

and

$\iint _{A}\beta ^{k}dA=0$

(5.b)

The ${\textstyle \beta ^{k}}$ parameter is useful for computing the zigzag function as explained in the next section.

3.2 Computation of the zigzag function

Integrating Eq.(2b) over the cross section and using Eq.(5b) and the fact that ${\textstyle \Psi }$ is independent of ${\textstyle z}$ yields

$\gamma ={\frac {1}{A}}\iint _{A}\gamma _{xz}^{k}dA$

(6)

i.e. ${\textstyle \gamma }$ represents the average shear strain of the cross section, as expected.

The shear strain-shear stress relationship of Eq.(3b) is written as

$\tau _{xz}^{k}=G^{k}\eta +G^{k}(1+\beta ^{k})\Psi$

(7)

where ${\textstyle \eta =\gamma -\Psi }$ is a difference function.

Remark 1. Function ${\textstyle \Psi }$ can be interpreted as a weighted-average shear strain angle [22]. The value of ${\textstyle \Psi }$ should be prescribed to zero at a clamped edge and left unprescribed at free and simply supported edges.

Eq.(7) shows that the distribution of ${\textstyle \tau _{xz}^{k}}$ within each layer is constant, as ${\textstyle \eta }$ is independent of the zigzag function and ${\textstyle \beta ^{k}}$ is constant.

The distribution of ${\textstyle \tau _{xz}^{k}}$ is now enforced to be independent of the zigzag function. This can be achieved by constraining the term multiplying ${\textstyle \Psi }$ in Eq.(7) to be constant, i.e.

$G^{k}(1+\beta ^{k})=G^{k+1}(1+\beta ^{k+1})=G,\quad {\hbox{constant}}$

(8)

This is equivalent to enforcing the interfacial continuity of the second term in the r.h.s. of Eq.(7).

Remark 2. We emphasize that this zigzag theory does not enforce the continuity of the transverse shear stresses across the section. This is consistent with the kinematic freedom inherent in the lower order kinematic approximation of the underlying beam theory. An accurate continuous distribution of the transverse shear stress across the thickness of the laminate can be obtained “a posteriori” in terms of axial stresses by integrating the equilibrium equations as explained in Section 7.3.

From Eq.(8) we deduce

$\beta ^{k}={\frac {G}{G^{k}}}-1$

(9)

Substituting ${\textstyle \beta ^{k}}$ in the integral of Eq.(5b) gives

$G=\left[{\frac {1}{A}}\iint _{A}{\frac {dA}{G^{k}}}\right]^{-1}=\left[h\sum \limits _{k=1}^{N}{\frac {h^{k}}{G^{k}}}\right]^{-1}$

(10)

where ${\textstyle h}$ is the section depth. Substituting Eq.(9) into Eq.(5a) gives the following recursion relation for the zigzag displacement function values at the layer interfaces

${\bar {u}}_{k}=\sum \limits _{i=1}^{k}h^{i}\beta ^{i}\quad {\hbox{with}}\quad u^{0}=u^{N}=0$

(11)

Introducing Eq.(11) into (4) gives the expression for the zigzag function as

$\phi ^{k}={\frac {h^{k}\beta ^{k}}{2}}(\zeta ^{k}-1)+\sum \limits _{i=1}^{k}h^{i}\beta ^{i}$

(12)

Recall that superindex ${\textstyle k}$ denotes the number of each material layer.

Remark 3. For homogeneous material ${\textstyle G^{k}=G}$ and ${\textstyle \beta ^{k}=0}$. Hence, the zigzag function ${\textstyle \phi ^{k}}$ vanishes and we recover the kinematics and constitutive expressions of the standard Timoshenko composite laminated beam theory. This is a particular feature of this zigzag theory.

Remark 4. Note that differently from standard Timoshenko beam theory, a shear correction parameter is not needed in the RZ theory.

3.3 Constitutive relationship

The in-plane bending and transverse shear resultant stresses are defined as

${\hat {\boldsymbol {\sigma }}}_{p}=\left\{{\begin{array}{c}N\\M\\M_{\phi }\end{array}}\right\}=\iint _{A}{\boldsymbol {S}}_{p}^{T}\sigma _{x}^{k}dA=\left(\iint _{A}{\boldsymbol {S}}_{p}^{T}{\boldsymbol {S}}_{p}E^{k}dA\right){\hat {\boldsymbol {\varepsilon }}}_{p}={\hat {\boldsymbol {D}}}_{p}{\hat {\boldsymbol {\varepsilon }}}_{p}$ (13) ${\hat {\boldsymbol {\sigma }}}_{t}=\left\{{\begin{array}{c}Q\\Q_{\phi }\end{array}}\right\}=\iint _{A}{\boldsymbol {S}}_{t}^{T}\tau _{xz}^{k}dA=\left(\iint _{A}{\boldsymbol {S}}_{t}^{T}{\boldsymbol {S}}_{t}G^{k}dA\right){\hat {\boldsymbol {\varepsilon }}}_{t}={\hat {\boldsymbol {D}}}_{t}{\hat {\boldsymbol {\varepsilon }}}_{t}$ (14)

In vectors ${\textstyle {\hat {\boldsymbol {\sigma }}}_{p}}$ and ${\textstyle {\hat {\boldsymbol {\sigma }}}_{t}}$, ${\textstyle N,M}$ and ${\textstyle Q}$ are respectively the axial force, the bending moment and the transverse shear force of standard beam theory, whereas ${\textstyle M_{\phi }}$ and ${\textstyle Q_{\phi }}$ are an additional bending moment and an additional shear force which are conjugate to the new generalized strains ${\textstyle {\partial \Psi \over \partial x}}$ and ${\textstyle \Psi }$, respectively.

The generalized constitutive matrices ${\textstyle {\hat {\boldsymbol {D}}}_{p}}$ and ${\textstyle {\hat {\boldsymbol {D}}}_{t}}$ are

${\hat {\boldsymbol {D}}}_{p}=\iint _{A}E^{k}\left[{\begin{matrix}1&-z&\phi ^{k}\\-z&z^{2}&-z\phi ^{k}\\\phi ^{k}&-z\phi ^{k}&(\phi ^{k})^{2}\end{matrix}}\right]dA\quad ,\quad {\hat {\boldsymbol {D}}}_{t}=\left[{\begin{matrix}D_{s}&-g\\-g&g\end{matrix}}\right]$

(15.a)

with

$D_{s}=\iint _{A}G^{k}dA\quad ,\quad g=D_{s}-GA$

(15.b)

In the derivation of the expression for ${\textstyle {\hat {\boldsymbol {D}}}_{s}}$ we have used the definition of ${\textstyle \beta ^{k}}$ of Eq.(9).

The generalized constitutive equation can be written as

${\hat {\boldsymbol {\sigma }}}=\left\{{\begin{array}{c}{\hat {\boldsymbol {\sigma }}}_{p}\\{\hat {\boldsymbol {\sigma }}}_{t}\end{array}}\right\}={\hat {\boldsymbol {D}}}{\hat {\boldsymbol {\varepsilon }}}={\hat {\boldsymbol {D}}}\left\{{\begin{array}{c}{\hat {\boldsymbol {\varepsilon }}}_{p}\\{\hat {\boldsymbol {\varepsilon }}}_{t}\end{array}}\right\}\quad \quad {\hbox{with }}\quad {\hat {\boldsymbol {D}}}=\left[{\begin{matrix}{\hat {\boldsymbol {D}}}_{p}&{\boldsymbol {0}}\\{\boldsymbol {0}}&{\hat {\boldsymbol {D}}}_{t}\end{matrix}}\right]$

(16)

3.4 Virtual work expression

The virtual work expression for a distributed load ${\textstyle q}$ is

$\iiint _{V}(\delta \varepsilon _{x}^{k}\sigma _{x}^{k}+\delta \gamma _{xz}^{k}\tau _{xz}^{k})dV-\int _{l}\delta wqdA=0$

(17)

The l.h.s. of Eq.(17) contains the internal virtual work performed by the axial and tangential stresses and the r.h.s. is the external virtual work carried out by the distributed load. ${\textstyle V}$ and ${\textstyle l}$ are the volume and length of the beam, respectively.

Substituting Eqs.(3) into the expression for the virtual internal work and using Eqs.(13) and (14) gives

$\iiint _{V}\left(\delta \varepsilon _{x}^{k}\sigma _{x}^{k}+\delta \gamma _{xz}^{k}\tau _{xz}^{k}\right)dV=\!\!\iiint _{V}\left(\delta {\hat {\boldsymbol {\varepsilon }}}_{p}^{T}{\boldsymbol {S}}_{p}^{T}\sigma _{x}^{k}+\delta {\hat {\boldsymbol {\varepsilon }}}_{t}^{T}{\boldsymbol {S}}_{t}^{T}\tau _{xz}^{k}\right)dV=$ $=\!\!\int _{l}\left(\delta {\hat {\boldsymbol {\varepsilon }}}_{p}^{T}{\hat {\boldsymbol {\sigma }}}_{p}+\delta {\hat {\boldsymbol {\varepsilon }}}_{t}^{T}{\hat {\boldsymbol {\sigma }}}_{t}\right)dx$

(18)

The virtual work is therefore written as

$\int _{l}\left(\delta {\hat {\boldsymbol {\varepsilon }}}_{p}^{T}{\hat {\boldsymbol {\sigma }}}_{p}+\delta {\hat {\boldsymbol {\varepsilon }}}_{t}^{T}{\hat {\boldsymbol {\sigma }}}_{t}\right)dx-\int _{l}\delta wqdx=0$

(19)

4 TWO-NODED LRZ BEAM ELEMENT

The four kinematic variables are ${\textstyle u_{0},w_{0},\theta }$ and ${\textstyle \Psi }$. They can be discretized using 2-noded linear ${\textstyle C^{\circ }}$ beam elements of length ${\textstyle l^{e}}$ in the standard form as

${\boldsymbol {u}}=\left\{{\begin{array}{c}u_{0}\\w_{0}\\\theta \\\Psi \\\end{array}}\right\}=\sum \limits _{i=1}^{2}N_{i}{\boldsymbol {a}}_{i}={\boldsymbol {N}}{\boldsymbol {a}}^{e}$

(20)

with

${\boldsymbol {N}}=[N_{1}{\boldsymbol {I}}_{4},N_{2}{\boldsymbol {I}}_{4}]~~,~~{\boldsymbol {a}}^{e}=\left\{{\begin{array}{c}{\boldsymbol {a}}_{1}\\{\boldsymbol {a}}_{2}\end{array}}\right\}~~,~~{\boldsymbol {a}}_{i}=\left\{{\begin{array}{c}u_{0_{i}}\\w_{0_{i}}\\\theta _{i}\\\Psi _{i}\\\end{array}}\right\}$

(21)

where ${\textstyle N_{i}={\frac {1}{2}}(1+\xi \xi _{i})}$ with ${\textstyle \xi =1-{\frac {2x}{l^{e}}}}$ are the standard one-dimensional linear shape functions, ${\textstyle {\boldsymbol {a}}_{i}}$ is the vector of nodal kinematic variables and ${\textstyle {\boldsymbol {I}}_{4}}$ is the ${\textstyle 4\times 4}$ unit matrix.

Substituting Eq.(20) into the generalized strain vectors in Eq.(2c) gives

${\hat {\boldsymbol {\varepsilon }}}_{p}={\boldsymbol {B}}_{p}{\boldsymbol {a}}^{e}\quad ,\quad {\hat {\boldsymbol {\varepsilon }}}_{t}={\boldsymbol {B}}_{t}{\boldsymbol {a}}^{e}$

(22)

The generalized strain matrices ${\textstyle {\boldsymbol {B}}_{p}}$ and ${\textstyle {\boldsymbol {B}}_{t}}$ are

${\boldsymbol {B}}_{p}=[{\boldsymbol {B}}_{p_{1}},{\boldsymbol {B}}_{p_{2}}]\quad ,\quad {\boldsymbol {B}}_{t}=[{\boldsymbol {B}}_{t_{1}},{\boldsymbol {B}}_{t_{2}}]$

(23.a)

with

${\boldsymbol {B}}_{p_{i}}={\begin{bmatrix}\displaystyle {\frac {\partial N_{i}}{\partial x}}&0&0&0\\0&0&\displaystyle {\frac {\partial N_{i}}{\partial x}}&0\\0&0&0&\displaystyle {\frac {\partial N_{i}}{\partial x}}\end{bmatrix}}$

(23.b)

$\displaystyle {\boldsymbol {B}}_{t_{i}}={\begin{bmatrix}0&\displaystyle {\frac {\partial N_{i}}{\partial x}}&-N_{i}&0\\--&--&--&--\\0&0&0&N_{i}\end{bmatrix}}={\begin{bmatrix}{\boldsymbol {B}}_{s_{i}}\\--\\{\boldsymbol {B}}_{\psi _{i}}\end{bmatrix}}$

(23.c)

where ${\textstyle {\boldsymbol {B}}_{p_{i}}}$ and ${\textstyle {\boldsymbol {B}}_{t_{i}}}$ are the in-plane and transverse shear strain matrices for node ${\textstyle i}$.

The virtual displacement and generalized strain fields are expressed in terms of the virtual nodal kinematic variables as

$\delta {\boldsymbol {u}}={\boldsymbol {N}}\delta {\boldsymbol {a}}^{e}~~,~~\delta {\hat {\boldsymbol {\varepsilon }}}_{p}={\boldsymbol {B}}_{p}\delta {\boldsymbol {a}}^{e}~~,~~\delta {\hat {\boldsymbol {\varepsilon }}}_{t}={\boldsymbol {B}}_{t}\delta {\boldsymbol {a}}^{e}$

(24)

The discretized equilibrium equations are obtained by substituting Eqs.(13), (14), (20), (22) and (24) into the virtual work expression (19). After simplification of the virtual nodal kinematic variables, the following standard matrix equation is obtained

${\boldsymbol {K}}{\boldsymbol {a}}-{\boldsymbol {f}}=0$

(25)

where ${\textstyle {\boldsymbol {a}}}$ is the vector of nodal kinematic variables for the whole mesh.

The stiffness matrix ${\textstyle {\boldsymbol {K}}}$ and the equivalent nodal force vector ${\textstyle {\boldsymbol {f}}}$ are obtained by assembling the element contributions ${\textstyle {\boldsymbol {K}}^{e}}$ and ${\textstyle {\boldsymbol {f}}^{e}}$ given by

${\boldsymbol {K}}^{e}={\boldsymbol {K}}_{p}^{e}+{\boldsymbol {K}}_{t}^{e}$

(26)

with

${\boldsymbol {K}}_{p_{ij}}^{e}=\int _{l^{e}}{\boldsymbol {B}}_{p_{i}}^{T}{\hat {\boldsymbol {D}}}_{p}{\boldsymbol {B}}_{p_{j}}dx~~,~~{\boldsymbol {K}}_{t_{ij}}^{e}=\int _{l^{e}}{\boldsymbol {B}}_{t_{i}}^{T}{\hat {\boldsymbol {D}}}_{t}{\boldsymbol {B}}_{t_{j}}dx$

(27)

and

${\boldsymbol {f}}_{i}^{e}=\int _{l^{e}}N_{i}q[1,0,0,0]^{T}dx$

(28)

Matrix ${\textstyle {\boldsymbol {K}}_{p}^{e}}$ is integrated with a one-point numerical quadrature which is exact in this case. Full integration of matrix ${\textstyle {\boldsymbol {K}}_{t}^{e}}$ requires a two-point Gauss quadrature. This however leads to shear locking for slender composite laminated beams (Section 5).

In order to asses the influence of the reduced integration of matrix ${\textstyle {\boldsymbol {K}}_{t}^{e}}$ for overcoming the shear locking problem we split ${\textstyle {\boldsymbol {K}}_{t}^{e}}$ as follows

${\boldsymbol {K}}_{t}^{e}={\boldsymbol {K}}_{s}^{e}+{\boldsymbol {K}}_{\psi }^{e}+{\boldsymbol {K}}_{s\psi }^{e}+[{\boldsymbol {K}}_{s\psi }^{e}]^{T}$

(29.a)

with

${\boldsymbol {K}}_{s_{ij}}^{e}=\int _{l^{e}}D_{s}{\boldsymbol {B}}_{s_{i}}^{T}{\boldsymbol {B}}_{s_{j}}dx~~,~~{\boldsymbol {K}}_{\psi _{ij}}^{e}=\int _{l^{e}}g{\boldsymbol {B}}_{\psi _{i}}^{T}{\boldsymbol {B}}_{\psi _{j}}dx$

(29.b)

${\boldsymbol {K}}_{s\psi _{ij}}^{e}=\int _{l^{e}}(-g){\boldsymbol {B}}_{s_{i}}^{T}{\boldsymbol {B}}_{\psi _{j}}dx$

where ${\textstyle {\boldsymbol {B}}_{s_{i}}}$ and ${\textstyle {\boldsymbol {B}}_{\psi _{i}}}$ are defined in Eq.(23c) and ${\textstyle D_{s}}$ and ${\textstyle g}$ are given in Eq.(15b).

The new linear beam element based on the RZ theory is termed LRZ.

A study of the accuracy of the LRZ beam element for analysis of slender laminated beams using one and two-point quadratures for integrating matrices ${\textstyle {\boldsymbol {K}}_{s}^{e}}$, ${\textstyle {\boldsymbol {K}}_{\psi }^{e}}$ and ${\textstyle {\boldsymbol {K}}_{s\psi }^{e}}$ is presented in the next section.

5 STUDY OF SHEAR LOCKING FOR THE LRZ BEAM ELEMENT

We study the performance of the LRZ beam element for the analysis of a cantilever beam of length ${\textstyle L}$ under an end point load of value ${\textstyle F=1}$ (Figure 2). The beam is formed by a symmetric three-layered material whose properties are listed on Table 1. The analysis is performed for four span-to-thickness ratios: ${\textstyle \lambda =5,10,50,100}$ (${\textstyle \lambda =L/h}$) using a mesh of 100 LRZ beam elements. Results for the LRZ element are labelled “ZZ” in the figures.

The same beam was analized using a mesh of 27000 four-noded plane stress rectangles for comparison purposes (Figure 3). Results for the plane stress analysis are labeled “PS” in the figures.

Figure 2: Cantilever beam under point load

Figure 3: Mesh of 27000 4-noded plane stress rectangular elements for analysis of cantilever and simple supported beams

Table. 1 Symmetric 3-layered cantilever beam. Material properties for shear locking study

	Composite material properties
	Layer 1	Layer 2	Layer 3
	(bottom)	(core)	(top)
h [mm]	6.6667	6.6667	6.6667
E [MPa]	2.19E5	2.19E3	2.19E5
G [MPa]	0.876E5	8.80E2	0.876E5

Figure 4 shows the ratio ${\textstyle r}$ between the end node deflection obtained with the LRZ element (${\textstyle w_{zz}}$) and with the plane stress quadrilateral (${\textstyle w_{ps}}$) (i.e. ${\textstyle r={\frac {w_{zz}}{w_{ps}}}}$) versus the beam span-to-thickness ratio ${\textstyle d={\frac {L}{h}}}$. Results for the LRZ element have been obtained using exact two-point integration for all terms of matrix ${\textstyle {\boldsymbol {K}}_{t}^{e}}$ (Eq.(27)) and a one-point reduced integration for the following three groups of matrices: ${\textstyle {\boldsymbol {K}}_{s}^{e}}$; ${\textstyle {\boldsymbol {K}}_{s}^{e}}$ and ${\textstyle {\boldsymbol {K}}_{s\psi }^{e}}$; and ${\textstyle {\boldsymbol {K}}_{s}}$, ${\textstyle {\boldsymbol {K}}_{s\psi }^{e}}$ and ${\textstyle {\boldsymbol {K}}_{\psi }^{e}}$ (Eqs.(29b)).

Labels ``all', ``S', ``SPsi' and ``Psi' in Figures 4–7 refer to matrices ${\textstyle {\boldsymbol {K}}_{t}^{e}}$, ${\textstyle {\boldsymbol {K}}_{s}^{e}}$, ${\textstyle {\boldsymbol {K}}_{s\psi }^{e}}$ and ${\textstyle {\boldsymbol {K}}_{\psi }^{e}}$ of Eq.(29a), respectively.

Results in Figure 4 clearly show that the exact integration of ${\textstyle {\boldsymbol {K}}_{t}^{e}}$ leads to shear locking as expected. Good (locking-free) results are obtained by one-point reduced integration of the three groups of matrices.

The influence of reduced integration in the distribution of the transverse shear stress was studied next for the three groups of matrices. Figures 5–7 show the thickness distribution of ${\textstyle \tau _{xz}}$ in sections located at distances ${\textstyle {\frac {L}{20}},{\frac {L}{4}},{\frac {L}{2}}}$ and ${\textstyle {\frac {3}{4}}L}$ from the clamped end for span-to-thickness ratios of ${\textstyle \lambda =5,10}$ and 100. Results are compared with the plane stress solution and also with results obtained with a mesh of 300 standard 2-noded elements based on laminated Timoshenko beam theory (labelled TBT in the figures). All TBT results presented in the paper have been used with a simple shear correction factor of ${\textstyle {\frac {5}{6}}}$. Indeed a more accurate value of the shear correction factor in TBT can be used for laminated sections [28].

Figure 4:' $r$ ratio $\left(r={\frac {w_{zz}}{w_{ps}}}\right)$ versus $L/h$ for cantilever under point load analyzed with the LRZ element. Labels ``all', $S$, $SPsi$ and $Psi$ refer to matrices ${\boldsymbol {K}}_{t}^{e}$, ${\boldsymbol {K}}_{s}^{e},{\boldsymbol {K}}_{s\psi }^{e}$ and ${\boldsymbol {K}}_{\psi }^{e}$, respectively

(a) $\tau _{xz},\lambda =5,L/20$	(b) $\tau _{xz},\lambda =5,L/4$
(c) $\tau _{xz},\lambda =5,L/2$	(d) $\tau _{xz},\lambda =5,3L/4$
Figure 5: Symmetric 3-layered cantilever thick beam under end point load. Thickness distribution of shear stress for $\lambda =5$ at different sections

(a) $\tau _{xz},\lambda =10,L/20$	(b) $\tau _{xz},\lambda =10,L/4$
(c) $\tau _{xz},\lambda =10,L/2$	(d) $\tau _{xz},\lambda =10,3L/4$
Figure 6: Symmetric 3-layered cantilever thick beam under end point load. Thickness distribution of shear stress for $\lambda =10$ at different sections

(a) $\tau _{xz},\lambda =100,L/20$	(b) $\tau _{xz},\lambda =100,L/4$
(c) $\tau _{xz},\lambda =100,L/2$	(d) $\tau _{xz},\lambda =100,3L/4$
Figure 7: Symmetric 3-layered cantilever thick beam under end point load. Thickness distribution of transverse shear stress for $\lambda =100$ at different sections

The conclusion is that for small values of ${\textstyle \lambda }$ the reduced or exact reduced integration of matrix ${\textstyle {\boldsymbol {K}}_{t}^{e}}$ leads to similar results. For slender beams, however, results obtained using reduced integration for ${\textstyle {\boldsymbol {K}}_{s}^{e}}$; ${\textstyle {\boldsymbol {K}}_{s}^{e}}$ and ${\textstyle {\boldsymbol {K}}_{s\psi }^{e}}$; and ${\textstyle {\boldsymbol {K}}_{s}^{e}}$, ${\textstyle {\boldsymbol {K}}_{s\psi }^{e}}$ and ${\textstyle {\boldsymbol {K}}_{\psi }^{e}}$ are different. Slightly more accurate results are obtained with the second choice for the section at ${\textstyle x=L/4}$ and ${\textstyle \lambda =100}$ (Figure 7b).

In conclusion, we recommend using a reduced one-point integration for matrices ${\textstyle {\boldsymbol {K}}_{s}^{e}}$ and ${\textstyle {\boldsymbol {K}}_{s\psi }^{e}}$, while matrix ${\textstyle {\boldsymbol {K}}_{\psi }^{e}}$ should be integrated with a 2-point quadrature.

6 CONVERGENCE STUDY

The same three-layered cantilever beam of Figure 2 was studied next for three different set of thickness and material properties for the three layers as listed in Table 2. Material A is the more homogeneous one, while material C is clearly the more heterogeneous.

The problem was studied with six meshes of LRZ elements ranging from 5 to 300 elements. Tables 3–5 show the convergence with the number of elements for the deflection and function ${\textstyle \Psi }$ at the beam end, the maximum axial stress ${\textstyle \sigma _{x}}$ at the end section and the maximum shear stress ${\textstyle \tau _{xz}}$ at the mid section.

Table. 2 Non symmetric 3-layered cantilever beams. Material properties for convergence analysis

		Material properties
		Layer 1(bottom)	Layer 2 (core)	Layer 3 (top)
Composite A	h [mm]	6.66	6.66	6.66
	E [MPa]	4.40E5	2.19E4	2.19E5
	G [MPa]	2.00E5	8.80E3	8.76E4
Composite B	h [mm]	6.66	6.66	6.66
	E [MPa]	2.19E5	2.19E3	2.19E5
	G [MPa]	8.76E4	8.80E2	8.76E4
Composite C	h [mm]	2	16	2
	E [MPa]	7.30E5	7.30E2	2.19E5
	G [MPa]	2.92E5	2.20E2	8.76E4

Table. 3 Non symmetric 3-layered cantilever thick beams under end point load ($\lambda =5$). Relative error for $w$ at $x=L$

$e_{r}\%-w$ at $x=L$
Number of	Composites
elements	A	B	C
5	1.800	9.588	42.289
10	0.506	2.901	19.277
25	0.0860	0.499	4.913
50	0.0191	0.123	1.406
100	0.0048	0.031	0.339
300	0.0000	0.0000	0.0000

Table. 4 Non symmetric 3-layered cantilever thick beams under end point load ($\lambda =5$). Convergence study. Relative error for $\Psi$ at $x=L$

$e_{r}\%-\Psi$ at $x=L$
Number of	Composites
elements	A	B	C
5	0.040	8.563	36.113
10	0.003	1.814	8.042
25	0.000	0.259	0.328
50	0.000	0.063	0.033
100	0.000	0.016	0.007
300	0.000	0.000	0.000

Table. 5 Non symmetric 3-layered cantilever thick beams under end point load ($\lambda =5$). Convergence study. (a) Relative error for the maximum value of $\sigma _{x}$ at $x=L$ and (b) idem for $\tau _{xz}$ at $x=L/2$

(a)
$e_{r}\%-(\sigma _{x})_{\max }$ at $x=L$
Number of	Composites
elements	A	B	C
5	0.568	6.923	18.239
10	0.076	2.704	12.437
25	0.013	0.568	4.266
50	0.003	0.131	1.095
100	0.001	0.029	0.250
300	0.000	0.000	0.000

(b)
$e_{r}\%-(\tau _{xz})_{\max }$ at ${\frac {L}{2}}$
Number of	Composites
elements	A	B	C
5	7.020	19.283	50.938
10	0.352	5.176	20.602
25	0.052	0.888	3.408
50	0.010	0.210	0.707
100	0.003	0.049	0.147
300	0.000	0.000	0.000

Convergence is measured by the relative error defined (in absolute value) as

$e_{r}=\left|{\frac {v_{6}-v_{i}}{v_{6}}}\right|$

(30)

where ${\textstyle v_{6}}$ and ${\textstyle v_{i}}$ are the values of the magnitude of interest obtained using the finest grid (300 elements) and the ${\textstyle i}$th mesh (${\textstyle i=1,2,\cdots 5}$), respectively.

Results clearly show that convergence is always slower for the heterogeneous material case, as expected.

For a mesh of 25 elements the errors for all the magnitudes considered are less than 1% for materials A and B. For material C the maximum error does not exceed 5% (Table 5). For the 50 element mesh errors of the order of 1% or less were obtained in all cases.

Results for a 10 element mesh are good for material A (errors less than 0.4%), relatively good for material B (errors less than around ${\textstyle 5\%}$) and unacceptable for material C (errors ranging from around 8% to 20%

7 EXAMPLES OF APPLICATION

7.1 Three-layered thick cantilever beam with non symmetric material properties

We present results for a laminated thick cantilever beam under an end point load. The material properties are those of Composite C in Table 2. The span-to-thickness ratio is ${\textstyle \lambda =5}$.

For the laminated sandwich considered the core is eight times thicker than the face sheets. In addition, the core is three orders of magnitude more compliant than the bottom face sheet. Moreover, the top face sheet has the same thickness as the bottom face sheet, but is about three times stiffer. Note that this laminate does not possess material symmetry with respect to the mid-depth reference axis. The high heterogeneity of this stacking sequence is very challenging for the beam theories considered herein to model adequately.

As in previous section, the legend caption PS denotes the reference solution obtained with the structured mesh of 27000 four-noded plane stress quadrilaterals shown in Figure 3. TBT denotes the solution obtained with a mesh of 300 2-noded beam elements based on standard laminated Timoshenko beam theory. LRZ-300, LRZ-50, LRZ-25, LRZ-10 refer to the solution obtained with the LRZ beam element using meshes of 300, 50, 25 and 10 elements, respectively.

Figure 8 shows the deflection values along the beam length. Very good agreement with the plane stress solution is obtained already for the LRZ-50 mesh as expected from the conclusions of the previous section.

TBT results are considerable stiffer. The difference with the reference solution is about six times stiffer for the end deflection value.

Figure 9 shows the distribution of the axial displacements at the upper and lower surfaces of layer 3 (top layer) along the beam length. Excellent results are again obtained with the 50 element mesh. The TBT results are far from the correct ones.

Figure 8: Non symmetric 3-layered cantilever thick beam under end point load ($\lambda =5$). Distribution of the vertical deflection $w$ for different theories and meshes

(a)	(b)
Figure 9: Non symmetric 3-layered cantilever thick beam under end point load ($\lambda =5$). Axial displacement $u$ at the upper and lower surfaces of the top layer (layer 3)

(a)	(b)
(c)
Figure 10: Non symmetric 3-layered cantilever thick beam under end point load ($\lambda =5$). Thickness distribution of the axial displacement $u$ at $x=L/4$ (a), $x=L/2$ (b) and $x=L$ (c)

Figure 10 shows the thickness distribution for the axial displacement at sections located at distances ${\textstyle {\frac {L}{4}},{\frac {L}{2}}}$ and ${\textstyle {\frac {3L}{4}}}$ from the clamped end. Results for the LRZ element (LRZ-25, LRZ-50 and LRZ-300) are in good agreement with the reference solution. The TBT results have the standard linear distribution which is far from the correct zigzag results.

Figure 11 shows the distribution along the beam length of the axial stress ${\textstyle \sigma _{x}}$ at the top and bottom surfaces of the beam cross section. Very good agreement between the reference PS solution and the LRZ-50 and LRZ-300 results is obtained. Results for the LRZ-25 mesh compare reasonably well with the PS solution except in the vicinity of the clamped edge. This error is corrected for the LRZ-50 and LRZ-300 meshes. The TBT results yield a linear distribution of the axial stress along the beam, as expected. This introduces large errors in the axial stress values in the vicinity of the clamped edge, as clearly shown in Figure 11.

(a)	(b)
Figure 11: Non symmetric 3-layered cantilever thick beam under end point load ($\lambda =5$). Axial stress $\sigma _{x}$ at upper (a) and lower (b) surfaces of the cross section along the beam length

Figures 12 and 13 show the thickness distribution for the axial stress ${\textstyle \sigma _{x}}$ at the clamped section and at the center of the beam. LRZ results agree quite well with those of the reference solution. TBT results have an erroneous stress distribution for the top and bottom layers at the clamped end. These differences are less important for the central section.

Figure 14 shows the thickness distribution for the transverse shear stress ${\textstyle \tau _{xz}}$ at different sections (${\textstyle {\frac {L}{20}},{\frac {L}{4}},{\frac {L}{2}}}$ and ${\textstyle {\frac {3L}{4}}}$). LRZ results provide an accurate estimate of the average transverse shear stress value for each layer. The distribution of ${\textstyle \tau _{xz}}$ across the thickness can be substantially improved by using the equilibrium equations for computing ${\textstyle \tau _{xz}}$ “a posteriori” as explained in Section 7.3.

TBT results are acceptable for the central layer and clearly overestimate the transverse shear stress in sections far from the clamped end.

LRZ and TBT results for the distribution of the (constant) tangential shear stress ${\textstyle \tau _{xz}}$ for each of the three layers along the beam length are shown in Figure 15. TBT results are clearly inaccurate (except for the values at the clamped edge).

Figure 12: Non symmetric 3-layered cantilever thick beam under end point load ($\lambda =5$). Thickness distribution of the axial stress $\sigma _{x}$ at $x=0$

Figure 13: Non symmetric 3-layered cantilever thick beam under end point load ($\lambda =5$). Thickness distribution of the axial stress $\sigma _{x}$ at $x=L/2$

(a)	(b)
(c)	(d)
Figure 14: Non symmetric 3-layered cantilever thick beam under end point load ($\lambda =5$). Thickness distribution of transverse shear stress $\tau _{xz}$ at $L/20$ (a), $L/4$ (b), $L/2$ (c) and $3L/4$ (d)

(a)	(b)
(c)
Figure 15: Non symmetric 3-layered cantilever thick beam under end point load ($\lambda =5$). LRZ and TBT results for the transverse shear stress $\tau _{xz}$ along the beam. Layer 1 (a), layer 2 (b) and layer 3(c)

7.2 Three-layered simple supported (SS) thick beams under uniform load

The next example is the analysis of a three-layered simple supported thick beam under a uniformly distributed load of unit value (${\textstyle q=1}$). The material properties and the thickness for the three layers are shown in Table 6. The material has a non symmetric distribution with respect to the beam axis. An unusually low value for the shear modulus of the core layer has been taken, thus reproducing the effect of a damaged material in this zone. The span-to-thickness ratio is ${\textstyle \lambda =5}$. Results obtained with the LRZ element are once more compared with those obtained with a mesh of 300 2-noded TBT elements and with the mesh of 27000 4-noded plane stress (PS) rectangles shown in Figure 3. The PS solution has been obtained by fixing the vertical displacement of all nodes at the end sections and the horizontal displacement of the mid-line node at ${\textstyle x=0}$ and ${\textstyle x=L}$ to a zero value. This way of approximating a simple support condition leads to some discrepancies between the PS results and those obtained with beam theory.

No advantage of the symmetry of the problem for the discretization has been taken.

Figure 16 shows the distribution of the vertical deflection for the different methods. The error in the “best” maximum central deflection value versus the “exact” PS solution is ${\textstyle \simeq }$ 12% The discrepancy is due to the difference in the way the simple support condition is modelled in beam and PS theories, as well as to the limitations of beam theory to model accurately very thick beams. TBT results are inaccurate, as expected.

Figure 17 shows the distribution of the axial stress ${\textstyle \sigma _{x}}$ along the beam for the top surface of the second and third layer.

The accuracy of the LRZ results is remarkable with a maximum error of 10% despite of the modeling limitations mentioned above. TBT results are incorrect.

Figure 18 shows the thickness distribution of the axial displacement at the left end section and at ${\textstyle x={\frac {L}{4}}}$. The LRZ element captures very well the zigzag shape of the axial displacement field even for a coarse mesh of 10 elements. The TBT element yields an unrealistic linear distribution.

Figures 19 and 20 show the thickness distribution of the axial stress and the transverse shear stress at the left end and mid sections. The accuracy of the LRZ results is again noticeable (even for the coarse 10 element mesh). The TBT element fails to capture the zigzag distribution of the axial stress (Figure 19) and gives a wrong value of almost zero shear stress at the core layer for the two sections chosen (Figure 20).

Figure 21 shows the distribution of the shear stress ${\textstyle \tau _{xz}}$ along the beam for each of the three layers obtained with the LRZ and TBT elements. TBT results are accurate for the first and third layer but are wrong for the core layer.

Figure 22 shows a similar set of results for a moderately thick SS beam (${\textstyle \lambda =10}$) and the same material properties. Results shown are the distribution along the beam of the deflection and the axial stress at the top surface of layer 2. The accuracy of the LRZ element is again noticeable.

Table. 6 Thickness and material properties for 3-layered non-symmetric simple supported (SS) beam

	Thickness and material properties
	Layer 1 (bottom)	Layer 2 (core)	Layer 3 (top)
h [mm]	6.6666	6.6666	6.6666
E [MPa]	2.19E${\textstyle ^{5}}$	5.30E${\textstyle ^{5}}$	7.30E${\textstyle ^{5}}$
G [MPa]	8.76E${\textstyle ^{4}}$	2.90E${\textstyle ^{2}}$	2.92E${\textstyle ^{5}}$

Figure 16: Non symmetric 3-layered SS thick beam under uniformly distributed load $\left(\lambda =5\right)$. Distribution of vertical deflection $w$ along the beam length

(a)	(b)
Figure 17: Non symmetric 3-layered SS thick beam under uniformly distributed load $(\lambda =5)$. Distribution of axial stress $\sigma _{x}$ at upper surface of layer 2 (a) and layer 3 (b)

(a) $x=0$	(b) $x={\frac {L}{2}}$
Figure 18: Non symmetric 3-layered SS thick beam under uniformly distributed load $(\lambda =5)$. Thickness distribution of axial displacement at $x=0$ (a) and at $x=L/2$ (b)

(a)	(b)
Figure 19: Non symmetric 3-layered SS thick beam under uniformly distributed load $(\lambda =5)$. Thickness distribution of axial stress $\sigma _{x}$ at $x=0$ (a) and at $x=L/2$ (b)

(a)	(b)
Figure 20: Non symmetric 3-layered SS thick beam under uniformly distributed load $(\lambda =5)$. Thickness distribution of the shear stress at $x=L/20$ (a) and at $x=L/4$ (b)

(a)	(b)
(c)
Figure 21: Non symmetric 3-layered SS thick beam under uniformly distributed load $(\lambda =5)$. LRZ and TBT results for the distribution of $\tau _{xz}$ along the beam for layer 1 (a), layer 2 (b) and layer 3 (c)

(a)	(b)
Figure 22: Non symmetric 3-layered SS moderately thick beam under uniformly distributed load $(\lambda =10)$. Distribution along the beam length of the vertical deflection $w$ (a) and the axial stress $\sigma _{x}$ at the upper of layer 2 (b)

7.3 Non-symmetric ten-layered clamped slender beam under uniformly distributed loading

We present results for a ten-layered clamped slender rectangular beam (${\textstyle L=100}$ mm, ${\textstyle h=5}$ mm, ${\textstyle b=1}$ mm, ${\textstyle \lambda =20}$) under uniformly distributed loading (${\textstyle q=1}$ KN/mm). The composite material has the non-symmetric distribution across the thickness shown in Table 7.

Table. 7 10-layered clamped slender rectangular beam under uniformly distributed loading. (a) Thickness and material number for each of the 10 layers. (b) Properties of each material

(a)
Layer	$h_{i}$	Material
1	0.5	IV
2	0.6	I
3	0.5	V
4	0.4	III
5	0.7	IV
6	0.1	III
7	0.4	II
8	0.5	V
9	0.3	I
10	1	II

(b)
Material	E [MPa]	G [MPa]
I	2.19e5	0.876e5
II	7.3e5	2.92e5
III	0.0073e5	0.0029e5
IV	5.3e5	2.12e5
V	0.82e5	0.328e5

Figure 23 shows results for the deflection along the beam for LRZ meshes with 10 and 300 elements (LRZ-10 and LRZ-300). Results obtained with a mesh of 27.000 4-noded plane stress quadrilaterals and with a mesh of 300 TBT elements are also shown for comparison. Note the accuracy of the coarse LRZ-10 mesh and the erroneous results of the TBT solution.

Figure 24 shows the thickness distribution of the axial displacement and the axial stress (${\textstyle \sigma _{x}}$) for the section at ${\textstyle x={\frac {L}{4}}}$. The accuracy of the LRZ results is once more remarkable.

Figure 23: 10-layered clamped slender beam under uniform loading. Distribution of the deflection along the beam

(a)	(b)
Figure 24: 10-layered clamped slender beam under uniform loading. Thickness distribution of axial displacement (a) and axial stress $\sigma _{x}$ (b) for $x={\frac {L}{4}}$

Figure 25 shows the thickness distribution of the transverse shear stress at ${\textstyle x={\frac {L}{4}}}$. Results in Figure 25a show the values directly obtained with the LRZ-10 and LRZ-300 meshes. These results are clearly better than those obtained with the TBT element but only coincide in an average sense with the plane stress FEM solution.

(a)	(b)
Figure 25: 10-layered clamped slender beam under uniform loading. Thickness distribution of $\tau _{xz}$ at $x={\frac {L}{4}}$. (a) Comparison of LRZ-10 and LRZ-300 results with plane stress (PS) and TBT solutions. (b) PS solution and LRZ-10-$N_{z}$ and LRZ-300-$N_{z}$ results for $\tau _{xz}$ obtained by thickness integration of the equilibrium equation using the LRZ-10 and LRZ-300 results (Eq.(32))

LRZ results for the thickness distribution of ${\textstyle \tau _{xz}}$ can be much improved by computing ${\textstyle \tau _{xz}}$ “a posteriori” from the axial stress field using the equilibrium equation

${\frac {\partial \sigma _{x}}{\partial x}}+{\frac {\partial \tau _{xz}}{\partial z}}=0$

(31)

The transverse shear stress at a point across the thickness with coordinate ${\textstyle z}$ is computed by integrating Eq.(31) as

$\tau _{xz}(z)=-\int _{-{\frac {h}{2}}}^{z}{\frac {\partial \sigma _{x}}{\partial x}}dz=-{\frac {\partial N_{z}}{\partial x}}$

(32)

where

$N_{z}=\int _{-{\frac {h}{2}}}^{z}\sigma _{x}dz$

(33)

is the axial force (per unit width) resulting from the thickness integration of ${\textstyle \sigma _{x}}$ between the coordinates ${\textstyle -{\frac {h}{2}}}$ and ${\textstyle z}$.

The space derivative of ${\textstyle N_{z}}$ in Eq.(32) is computed at a node ${\textstyle i}$ as

${\frac {\partial N_{z}}{\partial x}}={\frac {2}{l^{e}+l^{e-1}}}(N_{z}^{e}-N_{z}^{e-1})$

(34)

where (${\textstyle l^{e},N_{z}^{e}}$) and (${\textstyle l^{e-1},N_{z}^{e-1}}$) are the element length and the value of ${\textstyle N_{z}}$ at elements ${\textstyle e}$ and ${\textstyle e-1}$ adjacent to node ${\textstyle i}$, respectively. A value of ${\textstyle \tau _{xz}(-{\frac {h}{2}})=0}$ is taken. It is remarkable that the method yields automatically ${\textstyle \tau _{xz}({\frac {h}{2}})\simeq 0}$.

Results for ${\textstyle \tau _{xz}}$ obtained with this procedure are termed LRZ-10-${\textstyle N_{z}}$ and LRZ-300-${\textstyle N_{z}}$ in Figure 25. We note the accuracy of the “recovered” thickness distribution for ${\textstyle \tau _{xz}}$, even for the coarse mesh of 10 LRZ elements.

7.4 Modeling of delamination with the LRZ element

Prediction of delamination in composite laminated beams is a challenge for all beam models. A method for predicting delamination in beams using a Hermitian zigzag theory was presented in [26,27]. A sub-laminate approach is used for which the number of kinematic unknowns depends of the number of physical layers. This increases the number of variables but it yields the correct an accurate transverse shear stress distribution without integrating the equilibrium equations.

Delamination effects in composite laminated beams can be effectively reproduced with the LRZ element without introducing additional kinematic variables. The delamination model simply implies introducing a very thin “interface layer” between adjacent material layers in the actual composite laminated section. Delamination is produced when the material properties of the interface layer are drastically reduced to almost a zero value in comparison with those of the adjacent layers due to interlamina failure. This simple delamination model allows the LRZ element to take into account the reduction of the overall beam stiffness due to the failure of the interface layer leading to an increase in the deflection and rotation field. Moreover, the LRZ element can also accurately represent the jump in the axial displacement field across the interface layer and the change in the axial and tangential stress distributions over the beam sections as delamination progresses.

Figures 26–30 show an example of the capabilities of the LRZ beam element to model delamination. The problem represents the analysis of a cantilever thick rectangular beam (${\textstyle \lambda =5}$) under an end point load. The beam section has three layers of composite material with properties shown in Table 8. Delamination between the upper and core layers has been modelled by introducing a very thin interface layer (${\textstyle h=0.01}$ mm) between these two layers (Figure 26). The initial properties of the interface layer coincide with those of the upper layer. Next, the shear modulus value for the interface layer has been progressively reduced up to 11 orders of magnitude from ${\textstyle G2=8.76\times 10^{4}}$ MPa (Model 1) to ${\textstyle G2=8.76\times 10^{-7}}$ MPa (Model 12) (Table 9).

Figure 26: Modeling of interface layer for delamination study in 3-layered thick cantilever beam ($\lambda =5$) under end point load

Table. 8 Thickness and layer properties for delamination study in a 3-layered cantilever beam under end point load. Layer 2 is the interface layer. G2 values are given in Table 9

	Composite material
	Layer 1	Layer 2	Layer 3	Layer 4
h [mm]	2	0.01	16	2
E [MPa]	2.19E5	2.19E5	0.0073E5	7.30E5
G [MPa]	0.876E5	G2	0.0029E5	2.92E5

Table. 9 Shear modulus values for the interface layer for delamination study in a 3-layered cantilever beam. Values of G2 in MPa

Model	G2	Model	G2	Model	G2
1	8.76E+004	5	8.76E+000	9	8.76E-004
2	8.76E+003	6	8.76E-001	10	8.76E-005
3	8.76E+002	7	8.76E-002	11	8.76E-006
4	8.76E+001	8	8.76E-003	12	8.76E-007

We note that the reduction of the shear modulus has been applied over the whole beam length in this case. However it can applied in selected beam regions as appropriate.

Figure 27 shows results for the end deflection in terms of the shear modulus value of the interface layer for the LRZ-100 mesh. Note that the deflection increases one order of magnitude versus the non-delaminated case. It is also interesting that the end deflection does not change after the shear modulus of the interface layer is reduced beyond eight orders of magnitude (results for Model 9 in Figure 27). Results agree reasonably well (error ${\textstyle \simeq 10\%}$) with those obtained with the plane stress model of Figure 3 introducing a similar reduction in the shear modulus of an ad hoc interface layer.

Figure 27: Delamination study in 3-layered cantilever beam under end point load. Evolution of end deflection with the shear modulus value for the interface layer LRZ-100 results and PS solution

Figure 28 shows the thickness distribution for the axial displacement at the mid section for four decreasing values of the shear modulus at the interface layer: ${\textstyle G2=8.76,8.76\times 10^{-1},8.76\times 10^{-3}}$ and ${\textstyle 8.76\times 10^{-6}}$ MPa. The jump of the axial displacement across the thickness at the interface layer during delamination is well captured. We again note that the displacement jump at the interface layer remains stationary after a reduction of the material properties in that layer of six orders of magnitude. Results agree well with the plane stress solution also shown in the figure.

Figure 28: Delamination study in 3-layered cantilever beam under end point load. Thickness distribution of axial displacement at $x={\frac {L}{2}}$ for four decreasing values of the shear modulus at the interface layer (Models 5, 6, 8 and 11, Table 9)

Figure 29 shows the thickness distribution of the axial stress (${\textstyle \sigma _{x}}$) for the same four decreasing values of the shear modulus at the interface layer. The effect of delamination in the stress distribution is clearly visible. Once again the LRZ-100 results agree well with the plane stress solution.

Figure 29: Delamination study in 3-layered cantilever beam under end point load. Thickness distribution of $\sigma _{x}$ at $x={\frac {L}{2}}$ for four decreasing values of the shear modulus at the interface layer (Models 5, 6, 8 and 11, Table 9)

Figure 30 finally shows the thickness distribution for the transverse shear stress at ${\textstyle x={\frac {L}{2}}}$ for the same four values of the shear modulus at the interface layer. The three graphs show the PS results, the LRZ-100 results and the solution obtained by integrating the equilibrium equation (via Eqs.(31)–(34)) using the LRZ-100 results. Note the accuracy of the later solution versus the standard LRZ-100 results as delamination develops and the transverse shear stress progressively vanishes at the interface layer.

Similar good results for predicting the delamination and the thickness distribution of the axial and transverse shear stresses are obtained over the entire beam length.

The example shows clearly the capability of the LRZ element to model a complex phenomenon such as delamination in composite laminated beams without introducing additional kinematic variables. More evidence of the good behaviour of the LRZ beam element for predicting delamination in beams are reported in [29].

Figure 30: Delamination study in 3-layered cantilever beam under end point load. Thickness distribution of $\tau _{xz}$ at $x={\frac {L}{2}}$ for four values of G at the interface layer (Models 5, 6, 8 and 11, Table 9). LRZ-100 results, plane stress (PS) solution and LRZ-100-Sx results obtained by integrating the equilibrium equation (Eq.(32)) using the LRZ-100 results

8 CONCLUSIONS

We have presented a simple and accurate 2-noded beam element based on the combination of the refined zigzag beam theory proposed by Tessler et al. [22]. The element has four degrees of freedom per node (the axial displacement, the deflection, the rotation and the amplitude of the zigzag function). A standard ${\textstyle C^{\circ }}$ linear interpolation is used for all variables. The resulting LRZ beam element is shear locking-free and has shown an excellent behaviour for analysis of thick and thin composite laminated beams with clamped and simple supported conditions. Numerical results agree in practically all cases with those obtained with a two-dimensional plane-stress FEM using a far larger number of degrees of freedom. It is remarkable that the zigzag distribution of the axial displacement and the axial stress across the thickness, typical of composite laminated beams, is very accurately captured with the basic approximation chosen. The possibilities of the new LRZ beam element for predicting delamination effects has been demonstrated in a simple but representative example of application.

ACKNOWLEDGEMENTS

This research was partially supported by project SEDUREC of the Consolider Programme of the Ministerio de Educación y Ciencia of Spain.

REFERENCES

[1] Timoshenko S.P. and Woinowsky-Krieger S., Theory of Plates and Shells. McGraw-Hill, New York, 3rd Edition, 1959.

[2] Timoshenko S. On the correction for shear of differential equations for transverse vibrations of prismatic bars. Philosophical Magazine Series, Vol. 41, pp. 744-746, 1921.

[3] Liu D. and Li X. An overall view of laminate theories based on displacement hypothesis. Journal of Composite Materials, Vol. 30(14), pp. 1539-1561, 1996.

[4] Reddy, J.N., Mechanics of laminated composite plates and shells. Theory and analysis. 2nd Edition, CRC Press, Boca Raton, 2004.

[5] Owen D.R.J. and Li Z.H., A refined analysis of laminated plates by finite element displacement methods. Part I. Fundamentals and static analysis. II Vibration and stability. Comp. Struct, Vol. 26, pp. 907–923, 1987.

[6] Botello, S. Oñate, E. and Canet, J.M., A layer-wise triangle for analysis of laminated composite plates and shells. Computers and Structures, Vol. 70, 635–646, 1999.

[7] Di Sciuva M. Bending, vibration and buckling of simply supported thick multilayered orthotropic plates: an evaluation of a new displacement model. Journal of Sound and Vibration, Vol. 105, pp. 425-442, 1986.

[8] Murakami H. Laminated composite plate theory with improved in-plane responses. ASME Journal of Applied Mechanics, Vol. 53, pp. 661-666, 1986.

[9] Aitharaju V.R. and Averill R.C., An assessment of zig-zag kinematic displacement models for the anlaysis of laminated composites. Mech. Composite Mat. and Struct., 6, 1–26, 1999a.

[10] Aitharaju V.R. and Averill R.C., C${\textstyle ^{\circ }}$ zig-zag finite element for analysis of laminated composite beams. J. Engineering Mechanics, ASCE, 323–330, 1999b.

[11] Li X. and Liu D. An interlaminar shear stress continuity theory for both thin and thick composite laminates. J. Appl. Mech., 59, 502–509, 1992.

[12] Di Sciuva M. An improved shear-deformation theory for moderately thick multilayered anisotropic shells and plates. J. Appl. Mech. 54, 589–594, 1987.

[13] Toledano A and Murakami H. A higher-order laminate plate theory with improved in-plane response. Int. J. Solids and Struct., 23, 111–131, 1987b.

[14] Cho M. and Parmerter R.R. Efficient higher order composite plate theory for general laminations configuration. AIAA J., 31, 1299–1306, 1993.

[15] Kapuria, S., Dumir, P.C., Ahmed, A. and Alam, N. Finite element model of efficient zigzag theory for static analysis of hybrid piezoelectric beams, Computational Mechanics, Vol. 34 (6), pp. 475–483, 2004.

[16] Averill R.C. Static and dynamic response of moderately thick laminated beams with damage. Composites Engineering, Vol. 4(4), pp. 381-395, 1994.

[17] Averill R.C. and Yuen Cheong Yip. Development of simple, robust finite elements based on refined theories for thick laminated beams. Computers & Structures, Vol. 59(3), pp. 529-546, 1996.

[18] Alam N.M. and Upadhyay N.Kr. Finite element analysis of laminated composite beams for zigzag theory using MATLAB. Int. J. of Mechanics and Solids, Vol. 5(1), pp. 1–14, 2010.

[19] Savoia M. On the accuracy of one-dimensional models for multilayered composite beams. Int. J. Solids Struct., Vol. 33, pp. 521–-44, 1996.

[20] Kapuria S, Dumir P.C. and Jain N.K. Assessment of zigzag theory for static loading, buckling, free and forced response of composite and sandwich beams. Composite Structures, Vol. 64, pp. 317–-27, 2004.

[21] Carrera, E. Historical review of zigzag theories for multilayered plate and shell. Applied Mechanics Review, 56(3), 287–308, 2003.

[22] Tessler, A., Di Sciuva M. and Gherlone, M. A refined zigzag beam theory for composite and sandwich beams. J. of Composite Materials, Vol. 43, 1051–1081, 2009.

[23] Tessler, A., Di Sciuva M. and Gherlone, M. A consistent refinement of first-order shear-deformation theory for laminated composite and sandwich plates using improved zigzag kinematics. J. of Mechanics of Materials and Structures, 5(2), 341–367, 2010.

[24] Gherlone, M., Tessler, A. and Di Sciuva M. C${\textstyle ^{\circ }}$ beam element based on the refined zigzag theory for multilayered composite and sandwich laminates. Composite Structures. Accepted for publication, May 2011.

[25] Oñate, E., Eijo, A. and Oller, S. Two-noded beam element for composite and sandwich beams using Timoshenko theory and refined zigzag kinematics. Publication CIMNE PI346, CIMNE, Barcelona, October 2010.

[26] Di Sciuva M. and Gherlone, M. A global/local third-order Hermitian displacement field with damaged interfaces and transverse extensibility: FEM formulation. Composite Structures, 59(4), 433-444, 2003.

[27] Di Sciuva M. and Gherlone, M. Quasi-3D static and dynamic analysis of undamaged and damaged sandwich beams. Journal of Sandwich Structures & Materials, 7(1), 31–52, 2005.

[28] Oñate, E. Structural Analysis with the Finite Element Method. Vol. 2 Beams, Plates and Shells. CIMNE-Springer, 2011.

[29] Oñate, E., Eijo, A. and Oller, S. Modeling of delamination in composite laminated beams using a two-noded beam element based in refined zigzag theory. Publication CIMNE, PI367, November 2011.