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.'' <span id='citeF-22'></span>[[#cite-22|[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 <span id='citeF-1'></span>[[#cite-1|[1]]] and the more advanced Timoshenko theory  <span id='citeF-2'></span>[[#cite-2|[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 <span id='citeF-3'></span>[[#cite-3|[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 <span id='citeF-3'></span>[[#cite-3|[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 <span id='citeF-3'></span><span id='citeF-4'></span>[[#cite-3|[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 <span id='citeF-5'></span><span id='citeF-6'></span>[[#cite-5|[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 <span id='citeF-7'></span>[[#cite-7|[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 <span id='citeF-7'></span>[[#cite-7|[7,8]]], a quadratic displacement field <span id='citeF-10'></span>[[#cite-7|[10,11]]] or a cubic displacement field <span id='citeF-12'></span>[[#cite-12|[12,13,14]]] through the thickness of the laminate.

​

Many zigzag theories  require <math display="inline">C^1</math> continuity for the deflection field, which is a drawback versus simpler <math display="inline">C^\circ </math> continuous FEM approximations. Also many zigzag theories run into theoretical difficulties to satisfy equilibrium of forces at a clamped support.

​

Averill et al. <span id='citeF-9'></span><span id='citeF-10'></span><span id='citeF-16'></span><span id='citeF-17'></span>[[#cite-9|[9,10,16,17]]] developed linear, quadratic and cubic zigzag beam theories that overcame the need for <math display="inline">C^1</math> continuity. The shear strain angle is introduced as a kinematic variable together with the deflection, the rotation and the zigzag function. A <math display="inline">C^\circ </math> 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 <span id='citeF-16'></span><span id='citeF-17'></span>[[#cite-16|[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 <span id='citeF-18'></span>[[#cite-18|[18]]]. Good results are reported for cantilever and clamped composite and sandwich beams.

​

An assessment of different zigzag theories for beam is reported in <span id='citeF-19'></span><span id='citeF-20'></span>[[#cite-19|[19,20]]]. A review of zigzag theory for plate analysis can be found in <span id='citeF-21'></span>[[#cite-21|[21]]].

​

Tessler ''et al.'' <span id='citeF-22'></span>[[#cite-22|[22,23]]] have developed a refined zigzag (RZ) theory  starting from the standard Timoshenko kinematic assumptions. This allows one using <math display="inline">C^\circ </math> 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. <span id='citeF-24'></span>[[#cite-24|[24]]] have  developed two and three-noded <math display="inline">C^\circ </math> 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 <span id='citeF-24'></span>[[#cite-24|[24]]].

​

Quite simultaneously to the above work, Oñate et al. <span id='citeF-25'></span>[[#cite-25|[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  <span id='citeF-25'></span>[[#cite-25|[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

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>u^k (x,z)  = u_0 (x) - z\theta (x) + \bar u^k (x,z)</math>

|-

| style="text-align: center;" | <math> w(x,z)  =w_0(x) </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (1.a)

|}

​

where

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>\bar u^k = \phi ^k (z)\Psi (x) \quad ,\quad k = 1, N </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (1.b)

|}

​

is the ''zigzag displacement function'' (Figure [[#img-1|1]]).

​

In Eqs.(1) <math display="inline">N</math> is the number of layers, superscript <math display="inline">k</math> indicates quantities within the <math display="inline">k</math>th layer with <math display="inline">z_k \le z\le z_{k+1}</math> and <math display="inline">z_k</math> is the vertical coordinate of the <math display="inline">k</math>th interface. In Eq.(1a) the uniform axial displacement <math display="inline">u_0(x)</math>, the rotation <math display="inline">\theta (x)</math> and the transverse deflection <math display="inline">w_0(x)</math> are the primary kinematic  variables of the underlying equivalent single-layer Timoshenko beam theory. In Eq.(1b) function <math display="inline">\phi ^k (z)</math> denotes ''a piecewise linear zigzag function'', yet to be established, and <math display="inline">\Psi (x)</math> 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 [[#img-1|1]]c.

​

<div id='img-1'></div>

{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"

|-

|[[Image:Draft_Samper_624436055-test-Fig1.png|600px|Thickness distribution of the the zigzag function ϕ<sup>k</sup> (a), the zigzag displacement function ̄u<sup>k</sup> (b),  and the axial displacement (c) in refined zigzag beam theory]]

|- style="text-align: center; font-size: 75%;"

| colspan="1" | '''Figure 1:''' Thickness distribution of the the zigzag function <math>\phi ^k</math> (a), the zigzag displacement function <math>\bar u^k</math> (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.

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>\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 </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (2.a)

|}

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>\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 </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (2.b)

|}

​

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

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>\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} </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (2.c)

|}

​

where <math display="inline">\hat {\boldsymbol \varepsilon }_p</math> and <math display="inline">\hat{\boldsymbol \varepsilon }_t</math> are the generalized in-plane and transverse shear strain vectors, respectively. Vector <math display="inline">\hat {\boldsymbol \varepsilon }_p</math> contains the axial elongation <math display="inline">\left(\frac{\partial u_0}{\partial x}\right)</math>, the pseudo-curvature <math display="inline">\left(\frac{\partial \theta }{\partial x}  \right)</math> and the derivatives of the amplitude of the zigzag function <math display="inline">\left(\frac{\partial \Psi }{\partial x} \right)</math>. In <math display="inline">\hat{\boldsymbol \varepsilon }_t</math>, <math display="inline">\gamma = \frac{\partial w_0}{\partial x} -\theta </math> is the average transverse shear strain of Timoshenko beam theory and <math display="inline">\beta ^k = \frac{\partial \phi ^k}{\partial z}</math>. Note that since <math display="inline">\phi ^k (z)</math> is piecewise linear, <math display="inline">\beta ^k </math> is constant across each layer.

​

For major principal material axes that are coincident with the beam <math display="inline">x</math> axis, Hooke stress-strain relations for the <math display="inline">k</math>th orthotropic layer have the standard form

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>\sigma _x^k = E^k \varepsilon _x^k = E^k  {\boldsymbol S}_p\hat {\boldsymbol \varepsilon }_p</math>

| style="width: 5px;text-align: right;white-space: nowrap;" | (3.a)

|-

| style="text-align: center;" | <math> \tau _{xz}^k  = G^k \gamma _{xz}^k = G^k {\boldsymbol S}_t \hat{\boldsymbol \varepsilon }_t </math>

| style="width: 5px;text-align: right;white-space: nowrap;" | (3.b)

|}

|}

​

where <math display="inline">E^k</math> and <math display="inline">G^k</math> are the axial and shear moduli for the <math display="inline">k</math>th layer, respectively.

​

In the above equations we have distinguished all variables within a layer with superscript <math display="inline">k</math>.

​

==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, <math display="inline">C^\circ </math> continuity is only required for the FEM approximation of the kinematic variables.

​

Within each layer the zigzag function is expressed as

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>\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 </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (4)

|}

​

where <math display="inline">\bar \phi ^k</math> and <math display="inline">\bar \phi ^{k-1}</math> are the zigzag function values of the <math display="inline">k</math> and <math display="inline">k-1</math> interface, respectively with <math display="inline">\bar \phi ^0 = \bar \phi ^N =0</math> and <math display="inline">\zeta ^k = \frac{2(z - z^{k-1})}{h^k}-1</math>.

​

Collectively, function <math display="inline">\phi ^k</math> has the zigzag distribution shown in Figure [[#img-1|1]]a. Due to the dependence between the zigzag displacement function <math display="inline">\bar u^k</math> and <math display="inline">\bar \phi ^k</math> (see Eq.(1b)), <math display="inline">\bar u^k</math> also vanishes at the top and bottom layers. The axial displacement field is plotted in Figure [[#img-1|1]]c.

​

The above form of <math display="inline">\phi ^k</math> gives the constant value of <math display="inline">\beta ^k</math> for each layer as

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>\beta ^k = \frac{\partial \phi ^k}{\partial z} = \frac{\bar \phi ^k - \bar \phi ^{k-1}}{h^k} </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (5.a)

|}

​

and

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>\iint _A \beta ^k dA =0 </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (5.b)

|}

​

The <math display="inline">\beta ^k</math> 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 <math display="inline">\Psi </math> is independent of <math display="inline">z</math> yields

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>\gamma = \frac{1}{A} \iint _A \gamma _{xz}^k dA </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (6)

|}

​

i.e. <math display="inline">\gamma </math> represents the average shear strain of the cross section, as expected.

​

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

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>\tau _{xz}^k =G^k \eta + G^k (1+\beta ^k) \Psi  </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (7)

|}

​

where <math display="inline">\eta =\gamma - \Psi </math> is a difference function.

​

'''Remark 1'''. Function <math display="inline">\Psi </math>  can be interpreted as a weighted-average shear strain angle <span id='citeF-22'></span>[[#cite-22|[22]]]. The value of <math display="inline">\Psi </math> 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 <math display="inline">\tau _{xz}^k</math> within each layer is constant, as <math display="inline">\eta </math> is independent of the zigzag function and <math display="inline">\beta ^k </math> is constant.

​

The distribution of <math display="inline">\tau _{xz}^k</math> is now ''enforced to be independent of the zigzag function''. This can be achieved by constraining the term multiplying <math display="inline">\Psi </math> in Eq.(7) to be constant, i.e.

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>G^k (1+\beta ^k) = G^{k+1} (1+\beta ^{k+1})=G, \quad \hbox{constant} </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (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

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>\beta ^k = \frac{G}{G^k} -1 </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (9)

|}

​

Substituting <math display="inline">\beta ^k</math> in the integral of Eq.(5b) gives

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>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} </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (10)

|}

​

where <math display="inline">h</math> 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

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>\bar u_k = \sum \limits _{i=1}^k h^i \beta ^i \quad \hbox{with}\quad u^0 =u^N=0 </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (11)

|}

​

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

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>\phi ^k = \frac{h^k\beta ^k}{2} (\zeta ^k -1)+  \sum \limits _{i=1}^k  h^i \beta ^i </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (12)

|}

​

Recall that superindex <math display="inline">k</math> denotes the number of each material layer.

​

'''Remark 3.''' For homogeneous material <math display="inline">G^k= G</math> and <math display="inline">\beta ^k =0</math>. Hence, the zigzag function <math display="inline">\phi ^k</math> 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

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>\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</math>

| style="width: 5px;text-align: right;white-space: nowrap;" | (13)

|-

| style="text-align: center;" | <math> \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  </math>

| style="width: 5px;text-align: right;white-space: nowrap;" | (14)

|}

|}

​

In vectors <math display="inline">\hat{\boldsymbol \sigma }_p</math> and <math display="inline">\hat{\boldsymbol \sigma }_t</math>, <math display="inline">N,M</math> and <math display="inline">Q</math> are respectively the axial force, the bending moment and the transverse shear force of standard beam theory, whereas <math display="inline">M_\phi </math> and <math display="inline">Q_\phi </math> are an additional bending moment and an additional shear force which are conjugate to the new generalized strains <math display="inline">{\partial \Psi  \over \partial x}</math> and <math display="inline">\Psi </math>, respectively.

​

The generalized constitutive matrices <math display="inline">\hat {\boldsymbol D}_p</math> and <math display="inline">\hat {\boldsymbol D}_t</math> are

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>\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] </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (15.a)

|}

​

with

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>D_s =\iint _A G^k dA \quad ,\quad g = D_s - GA </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (15.b)

|}

​

In the derivation of the expression for <math display="inline">\hat {\boldsymbol D}_s</math> we have used the definition of <math display="inline">\beta ^k</math> of Eq.(9).

​

The generalized constitutive equation can be written as

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>\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] </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (16)

|}

​

===3.4 Virtual work expression===

​

The virtual work expression for a distributed load <math display="inline">q</math> is

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>\iiint _V (\delta \varepsilon _x^k \sigma _x^k + \delta \gamma _{xz}^k \tau _{xz}^k)dV - \int _l \delta w q dA =0 </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (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. <math display="inline">V</math> and <math display="inline">l</math> 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

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>\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 =</math>

|-

| style="text-align: center;" | <math> =\!\! \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 </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (18)

|}

​

The virtual work is therefore written as

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>\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 w q dx=0 </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (19)

|}

​

==4 TWO-NODED LRZ BEAM ELEMENT==

​

The  four kinematic variables are <math display="inline">u_0,w_0, \theta </math> and <math display="inline">\Psi </math>. They can be discretized using  2-noded linear <math display="inline">C^\circ </math> beam elements of length <math display="inline">l^e</math> in the standard form as

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>{\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  </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (20)

|}

​

with

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>{\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\} </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (21)

|}

​

where <math display="inline">N_i = \frac{1}{2} (1+\xi \xi _i)</math> with <math display="inline">\xi = 1 - \frac{2x}{l^e}</math> are the standard one-dimensional linear shape functions, <math display="inline">{\boldsymbol a}_i</math> is the vector of nodal kinematic variables and <math display="inline">{\boldsymbol I}_4</math> is the <math display="inline">4\times 4</math> unit matrix.

​

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

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>\hat {\boldsymbol \varepsilon }_p = {\boldsymbol B}_p {\boldsymbol a}^e \quad ,\quad \hat {\boldsymbol \varepsilon }_t = {\boldsymbol B}_t {\boldsymbol a}^e </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (22)

|}

​

The  generalized strain matrices <math display="inline">{\boldsymbol B}_p</math> and <math display="inline">{\boldsymbol B}_t </math> are

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>{\boldsymbol B}_p =[{\boldsymbol B}_{p_1},{\boldsymbol B}_{p_2} ]\quad ,\quad {\boldsymbol B}_t =[{\boldsymbol B}_{t_1},{\boldsymbol B}_{t_2} ] </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (23.a)

|}

​

with

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>{\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} </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (23.b)

|}

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>\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} </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (23.c)

|}

​

where <math display="inline">{\boldsymbol B}_{p_i}</math> and <math display="inline">{\boldsymbol B}_{t_i}</math> are the in-plane and transverse shear strain matrices for node <math display="inline">i</math>.

​

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

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>\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 </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (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

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>{\boldsymbol K}{\boldsymbol a}-{\boldsymbol f}=0 </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (25)

|}

​

where <math display="inline">{\boldsymbol a}</math> is the vector of nodal kinematic variables for the  whole mesh.

​

The stiffness matrix <math display="inline">{\boldsymbol K}</math> and the equivalent nodal force vector <math display="inline">{\boldsymbol f}</math> are obtained by assembling the element contributions <math display="inline">{\boldsymbol K}^e</math> and <math display="inline">{\boldsymbol f}^e</math> given by

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>{\boldsymbol K}^e = {\boldsymbol K}^e_p + {\boldsymbol K}^e_t </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (26)

|}

​

with

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>{\boldsymbol K}^e_{p_{ij}} = \int _{l^e} {\boldsymbol B}^T_{p_i} \hat {\boldsymbol D}_p {\boldsymbol B}_{p_j}dx ~~,~~ {\boldsymbol K}^e_{t_{ij}} = \int _{l^e} {\boldsymbol B}^T_{t_i} \hat {\boldsymbol D}_t {\boldsymbol B}_{t_j}dx </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (27)

|}

​

and

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>{\boldsymbol f}^e_i = \int _{l^e} N_i q [1,0,0,0]^T dx </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (28)

|}

​

Matrix <math display="inline">{\boldsymbol K}^e_p</math> is integrated with a one-point numerical quadrature which is exact in this case. Full integration of matrix <math display="inline">{\boldsymbol K}^e_t</math> 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 <math display="inline">{\boldsymbol K}^e_t</math> for overcoming the  shear locking problem we split <math display="inline">{\boldsymbol K}^e_t</math> as follows

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>{\boldsymbol K}^e_t = {\boldsymbol K}^e_s +{\boldsymbol K}^e_\psi + {\boldsymbol K}^e_{s\psi } + [{\boldsymbol K}^e_{s\psi }]^T </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (29.a)

|}

​

with

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>{\boldsymbol K}^e_{s_{ij}} = \int _{l^e}  D_s {\boldsymbol B}^T_{s_i} {\boldsymbol B}_{s_j}dx ~~,~~ {\boldsymbol K}^e_{\psi _{ij}} = \int _{l^e} g {\boldsymbol B}_{\psi _i}^T {\boldsymbol B}_{\psi _j}dx </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (29.b)

|}

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math> {\boldsymbol K}^e_{s\psi _{ij}}= \int _{l^e} (-g ) {\boldsymbol B}^T_{s_i}{\boldsymbol B}_{\psi _j}dx </math>

|}

|}

​

where <math display="inline">{\boldsymbol B}_{s_i}</math> and <math display="inline">{\boldsymbol B}_{\psi _i}</math> are defined in Eq.(23c) and <math display="inline">D_s</math> and <math display="inline">g</math> 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 <math display="inline">{\boldsymbol K}^e_s</math>, <math display="inline">{\boldsymbol K}^e_\psi </math> and <math display="inline">{\boldsymbol K}^e_{s\psi }</math>  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 <math display="inline">L</math> under an end point load of value <math display="inline">F=1</math> (Figure [[#img-2|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: <math display="inline">\lambda = 5,10,50,100</math> (<math display="inline">\lambda = L/h</math>) 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 [[#img-3|3]]).  Results for the plane stress analysis are labeled “PS” in the figures.

​

<div id='img-2'></div>

{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"

|-

|[[Image:Draft_Samper_624436055-test-Fig3.png|400px|]]

|- style="text-align: center; font-size: 75%;"

| colspan="1" | '''Figure 2:''' Cantilever beam under point  load 

|}

​

<div id='img-3'></div>

{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"

|-

|[[File:Draft_Samper_624436055_3572_Fig3-con299.png|600px|]]

|- style="text-align: center; font-size: 75%;"

| colspan="1" | '''Figure 3:''' Mesh of 27000 4-noded plane stress rectangular elements for analysis of cantilever and simple supported beams

|}

​

​

{|  class="floating_tableSCP wikitable" style="text-align: center; margin: 1em auto;min-width:50%;"

|+ style="font-size: 75%;" |<span id='table-1'></span>Table. 1 Symmetric 3-layered cantilever beam. Material properties for shear locking study

|- style="border-top: 2px solid;"

| style="text-align: left;" |   

| colspan='3' style="text-align: left;" | '''Composite material properties'''

|- style="border-top: 2px solid;"

| style="text-align: left;" |   

| Layer 1 

| Layer 2  

| Layer 3 

|-

| style="text-align: left;" | 

| (bottom)

| (core) 

| (top)

|- style="border-top: 2px solid;"

| style="text-align: left;" |  h [mm] 

| 6.6667 

| 6.6667 

| 6.6667 

|-

| style="text-align: left;" | E [MPa] 

| 2.19E5 

| 2.19E3 

| 2.19E5 

|- style="border-bottom: 2px solid;"

| style="text-align: left;" | G [MPa] 

| 0.876E5 

| 8.80E2 

| 0.876E5 

​

|}

​

Figure [[#img-4|4]] shows the ratio <math display="inline">r</math> between the end node deflection  obtained with the LRZ element (<math display="inline">w_{zz}</math>) and with the plane stress quadrilateral (<math display="inline">w_{ps}</math>) (i.e. <math display="inline">r=\frac{w_{zz}}{w_{ps}}</math>) versus the beam span-to-thickness ratio <math display="inline">d = \frac{L}{h}</math>. Results for the LRZ element have been obtained using ''exact'' two-point integration for all terms of  matrix <math display="inline">{\boldsymbol K}^e_t</math> (Eq.(27)) and  a one-point ''reduced'' integration for the following three groups of matrices: <math display="inline">{\boldsymbol K}^e_s</math>; <math display="inline">{\boldsymbol K}^e_s</math> and <math display="inline">{\boldsymbol K}^e_{s\psi }</math>; and <math display="inline">{\boldsymbol K}_s</math>, <math display="inline">{\boldsymbol K}^e_{s\psi }</math> and <math display="inline">{\boldsymbol K}^e_{\psi }</math> (Eqs.(29b)).

​

Labels ``all'<nowiki/>''', ``S'<nowiki/>''', ``SPsi'<nowiki/>''' and `''`''Psi'''' in Figures [[#img-4|4]]&#8211;[[#img-7|7]] refer to matrices <math display="inline">{\boldsymbol K}^e_t</math>, <math display="inline">{\boldsymbol K}^e_s</math>, <math display="inline">{\boldsymbol K}^e_{s\psi }</math> and <math display="inline">{\boldsymbol K}^e_{\psi }</math> of Eq.(29a), respectively.

​

Results in Figure [[#img-4|4]] clearly show that the exact integration of <math display="inline">{\boldsymbol K}^e_t</math> 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 [[#img-5|5]]&#8211;[[#img-7|7]] show the thickness distribution of <math display="inline">\tau _{xz}</math> in sections located at distances <math display="inline">\frac{L}{20},\frac{L}{4},\frac{L}{2}</math> and <math display="inline">\frac{3}{4}L</math> from the clamped end for span-to-thickness ratios of <math display="inline">\lambda =5 , 10</math> 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 <math display="inline">\frac{5}{6}</math>. Indeed a more accurate value of the shear correction factor in TBT can be used for laminated sections <span id='citeF-28'></span>[[#cite-28|[28]]].

​

<div id='img-4'></div>

{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"

|-

|[[Image:Draft_Samper_624436055-test-Fig4a.png|600px|r ratio \left(r = \fracw<sub>zz</sub>wₚₛ\right) versus L/h for cantilever under point load analyzed with the LRZ element. Labels ``''all'''', S, SPsi and Psi refer to matrices K<sup>e</sup>ₜ, Kₛ<sup>e</sup>,K<sub>sψ</sub><sup>e</sup> and K<sub>ψ</sub><sup>e</sup>, respectively]]

|- style="text-align: center; font-size: 75%;"

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

|}

​

<div id='img-5a'></div>

<div id='img-5b'></div>

<div id='img-5c'></div>

<div id='img-5d'></div>

<div id='img-5'></div>

{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"

|-

|[[Image:Draft_Samper_624436055-test-Grafico2.png|360px|τ<sub>xz</sub>,λ= 5 , L/20]]

|[[Image:Draft_Samper_624436055-test-Grafico3.png|360px|τ<sub>xz</sub> ,λ= 5 , L/4]]

|- style="text-align: center; font-size: 75%;"

| (a) <math>\tau _{xz},\lambda = 5 , L/20</math>

| (b) <math>\tau _{xz} ,\lambda = 5 , L/4</math>

|-

|[[Image:Draft_Samper_624436055-test-Grafico4.png|360px|τ<sub>xz</sub> ,λ= 5 , L/2]]

|[[Image:Draft_Samper_624436055-test-Grafico5.png|360px|τ<sub>xz</sub> ,λ= 5 , 3L/4]]

|- style="text-align: center; font-size: 75%;"

| (c) <math>\tau _{xz} ,\lambda = 5 , L/2</math>

| (d) <math>\tau _{xz} ,\lambda = 5 , 3L/4</math>

|- style="text-align: center; font-size: 75%;"

| colspan="2" | '''Figure 5:''' Symmetric 3-layered cantilever thick beam under end point load. Thickness distribution of shear stress for <math>\lambda =5</math> at different sections

|}

​

<div id='img-6a'></div>

<div id='img-6b'></div>

<div id='img-6c'></div>

<div id='img-6d'></div>

<div id='img-6'></div>

{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"

|-

|[[Image:Draft_Samper_624436055-test-Grafico6.png|360px|τ<sub>xz</sub> ,λ= 10 , L/20]]

|[[Image:Draft_Samper_624436055-test-Grafico7.png|360px|τ<sub>xz</sub> ,λ= 10 , L/4]]

|- style="text-align: center; font-size: 75%;"

| (a) <math>\tau _{xz} ,\lambda = 10 , L/20</math>

| (b) <math>\tau _{xz} ,\lambda = 10 , L/4</math>

|-

|[[Image:Draft_Samper_624436055-test-Grafico8.png|360px|τ<sub>xz</sub> ,λ= 10 , L/2]]

|[[Image:Draft_Samper_624436055-test-Grafico9.png|360px|τ<sub>xz</sub> ,λ= 10 , 3L/4]]

|- style="text-align: center; font-size: 75%;"

| (c) <math>\tau _{xz} ,\lambda = 10 , L/2</math>

| (d) <math>\tau _{xz} ,\lambda = 10 , 3L/4</math>

|- style="text-align: center; font-size: 75%;"

| colspan="2" | '''Figure 6:''' Symmetric 3-layered cantilever thick beam under end point load.  Thickness distribution of shear stress for <math>\lambda =10</math> at different sections

|}

​

<div id='img-7a'></div>

<div id='img-7b'></div>

<div id='img-7c'></div>

<div id='img-7d'></div>

<div id='img-7'></div>

{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"

|-

|[[Image:Draft_Samper_624436055-test-Grafico10.png|360px|τ<sub>xz</sub> ,λ= 100 , L/20]]

|[[Image:Draft_Samper_624436055-test-Grafico11.png|360px|τ<sub>xz</sub> ,λ= 100 , L/4]]

|- style="text-align: center; font-size: 75%;"

| (a) <math>\tau _{xz} ,\lambda = 100 , L/20</math>

| (b) <math>\tau _{xz} ,\lambda = 100 , L/4</math>

|-

|[[Image:Draft_Samper_624436055-test-Grafico12.png|360px|τ<sub>xz</sub> ,λ= 100 , L/2]]

|[[Image:Draft_Samper_624436055-test-Grafico13.png|360px|τ<sub>xz</sub> ,λ= 100 , 3L/4]]

|- style="text-align: center; font-size: 75%;"

| (c) <math>\tau _{xz} ,\lambda = 100 , L/2</math>

| (d) <math>\tau _{xz} ,\lambda = 100 , 3L/4</math>

|- style="text-align: center; font-size: 75%;"

| colspan="2" | '''Figure 7:''' Symmetric 3-layered cantilever thick beam under end point load. Thickness distribution of transverse shear stress for <math>\lambda =100</math> at different sections

|}

​

The conclusion is that for small values of <math display="inline">\lambda </math> the reduced or exact reduced integration of matrix <math display="inline">{\boldsymbol K}^e_t</math> leads to similar results. For slender beams, however, results obtained using reduced integration for <math display="inline">{\boldsymbol K}^e_s</math>; <math display="inline">{\boldsymbol K}^e_s</math> and <math display="inline">{\boldsymbol K}^e_{s\psi }</math>; and <math display="inline">{\boldsymbol K}^e_s</math>, <math display="inline">{\boldsymbol K}^e_{s\psi }</math> and <math display="inline">{\boldsymbol K}^e_{\psi }</math> are different. Slightly more accurate results are obtained with the second choice for the section at <math display="inline">x=L/4</math> and <math display="inline">\lambda =100</math> (Figure [[#img-7|7]]b).

​

In conclusion, we recommend using a reduced one-point integration for matrices <math display="inline">{\boldsymbol K}^e_s</math> and <math display="inline">{\boldsymbol K}^e_{s\psi }</math>, while matrix <math display="inline">{\boldsymbol K}^e_\psi </math> should be integrated with a 2-point quadrature.

​

==6 CONVERGENCE STUDY==

​

The same three-layered cantilever beam of Figure [[#img-2|2]] was studied next for three different set of thickness and material properties for the three layers as listed in Table [[#table-2|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 [[#table-3|3]]&#8211;[[#table-5|5]] show the convergence with the number of elements for the deflection and function <math display="inline">\Psi </math> at the beam end,  the maximum axial stress <math display="inline">\sigma _x</math> at the end section and the maximum shear stress <math display="inline">\tau _{xz}</math> at the mid section.

​

​

{|  class="floating_tableSCP wikitable" style="text-align: center; margin: 1em auto;min-width:50%;"

|+ style="font-size: 75%;" |<span id='table-2'></span>Table. 2 Non symmetric 3-layered cantilever beams. Material properties for convergence analysis

​

|-

| 

| 

| colspan='3' | '''Material properties'''

​

|-

| 

| 

| '''Layer 1'''(bottom) 

| '''Layer 2''' (core) 

| '''Layer 3''' (top) 

​

|- style="border-top: 2px solid;"

| '''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 

​

|- style="border-top: 2px solid;"

| '''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 

​

|- style="border-top: 2px solid;"

| '''Composite C''' 

| h [mm] 

| 2 

| 16 

| 2 

|-

| 

| E [MPa] 

| 7.30E5 

| 7.30E2 

| 2.19E5 

|- style="border-bottom: 2px solid;"

| 

| G [MPa] 

| 2.92E5 

| 2.20E2 

| 8.76E4  

​

|}

​

​

{|  class="floating_tableSCP wikitable" style="text-align: center; margin: 1em auto;min-width:50%;"

|+ style="font-size: 75%;" |<span id='table-3'></span>Table. 3 Non symmetric 3-layered cantilever thick beams under end point load (<math>\lambda =5</math>). Relative error for <math>w</math> at <math>x=L</math>

​

|- style="border-top: 2px solid;"

| colspan='4' | <math>e_r %  - w</math> at <math>x=L</math>

​

|- style="border-top: 2px solid;"

| '''Number of''' 

| colspan='3' | '''Composites'''

​

|-

| '''elements''' 

| '''A''' 

| '''B''' 

| '''C''' 

​

|- style="border-top: 2px solid;"

| 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 

|- style="border-bottom: 2px solid;"

| 300 

| 0.0000 

| 0.0000 

| 0.0000 

​

|}

​

​

{|  class="floating_tableSCP wikitable" style="text-align: center; margin: 1em auto;min-width:50%;"

|+ style="font-size: 75%;" |<span id='table-4'></span>Table. 4 Non symmetric 3-layered cantilever thick beams under end point load (<math>\lambda =5</math>). Convergence study. Relative error for <math>\Psi </math> at <math>x=L</math>

​

|- style="border-top: 2px solid;"

| colspan='4' | <math>e_r %  - \Psi </math> at <math>x=L</math>

​

|- style="border-top: 2px solid;"

| '''Number of''' 

| colspan='3' | '''Composites'''

​

|-

| '''elements''' 

| '''A''' 

| '''B''' 

| '''C''' 

​

|- style="border-top: 2px solid;"

| 5 

| 0.040 

| 8.563 

| 36.113 

|-

| 10 

| 0.003

| 1.814

| 8.042 

|-

| 25 

| 0.000

| 0.259 

| 0.328 

|-