Extended rotation-free plate and beam elements with shear deformation effects



==Abstract==

The paper describes a methodology for extending rotation-free plate and beam elements in order to accounting for transverse shear deformation effects. The ingredients for the element formulation are: a Hu-Washizu type mixed functional, a linear interpolation for the deflection and the  shear angles over standard finite elements and a finite volume approach for computing the bending moments and the curvatures over a patch of elements. As a first application of the general procedure we present an extension of the 3-noded rotation-free basic plate triangle (BPT) originally developed for thin plate analysis to accounting for shear deformation effects of relevance for thick plates and composite laminated plates. The nodal deflection degrees of freedom (DOFs) of the original BPT element are enhanced with the two shear deformation angles. This allows to computing the bending and shear deformation energies leading to a simple triangular plate element with 3 DOFs per node (termed BPT+ element). For the thin plate case the shear angles vanish and the element reproduces the good behaviour of the original thin BPT element. As a consequence the element is applicable to thick and thin plate situations without exhibiting shear locking effects. The numerical solution for the thick case can be found iteratively starting from the deflection values for the Kirchhoff theory using the original thin BPT element. A 2-noded rotation-free beam element termed CCB+ applicable to slender and thick beams is derived as a particular case of the plate formulation. The examples presented show the robustness and accuracy of the BPT+ and the CCB+ elements for thick and thin plate and beam problems.

'''keywords''' Rotation-free triangle, rotation-free beam, thick and thin plates and beams, finite elements, shear deformation


==1 INTRODUCTION==

In past years some authors have derived a number of thin plate and shell triangular elements free of rotational degrees of freedom (DOFs) based on Kirchhoff theory <span id='citeF-1'></span>[[#cite-1|[1]]]&#8211;<span id='citeF-24'></span>[[#cite-24|[24]]]. This idea goes back to the original attempts to solve thin plate bending problems using finite difference schemes with the deflection as the only nodal variable <span id='citeF-25'></span>[[#cite-25|[25]]]&#8211;<span id='citeF-27'></span>[[#cite-27|[27]]]. In essence all methods attempt to express the curvature field over an element in terms of the displacements of a collection of nodes belonging to a patch of adjacent elements. Oñate and Cervera <span id='citeF-6'></span>[[#cite-6|[6]]] proposed a general procedure of this kind combining finite element and finite volume concepts for deriving thin plate triangles and quadrilaterals with the deflection as the only nodal variable and presented a simple and competitive rotation-free three DOFs triangular element termed BPT (for Basic Plate Triangle). These ideas were extended  in <span id='citeF-12'></span>[[#cite-12|[12]]] to derive a number of rotation-free thin plate and shell triangles. The basic ingredients of the method are a mixed Hu-Washizu formulation, a standard discretization into three-noded triangles, a linear finite element interpolation of the displacement field within each triangle and a finite volume type approach for computing constant curvature and bending moment fields within appropriate non-overlapping control domains. The so called ''cell-centered'' and ''cell-vertex'' triangular domains yield different families of rotation-free plate and shell triangles. Both the BPT plate element and its extension to shell analysis (termed BST for Basic Shell Triangle) can be derived from the ''cell-centered'' formulation. Here the control domain is an individual triangle. The constant curvature field within a triangle is computed in terms of the displacements of the six nodes belonging to the four element patch formed by the chosen triangle and the three adjacent triangles. In the ''cell-vertex'' approach the control volume is the tributary domain of a node and this leads to a different family of rotation-free plate and shell triangles. The detailed derivation of the different rotation-free triangular plate and shell element families can be found in <span id='citeF-12'></span>[[#cite-12|[12]]]. The three-noded BST element has been successfully extended to non-linear shell problems involving frictional-contact situations and dynamics <span id='citeF-15'></span><span id='citeF-19'></span><span id='citeF-20'></span>[[#cite-15|[15,19,20]]]. Practical applications of the BST element to sheet stamping analysis are reported in <span id='citeF-9'></span><span id='citeF-10'></span><span id='citeF-11'></span><span id='citeF-24'></span>[[#cite-9|[9,10,11,24]]].

The paper describes an extension of the original rotation-free thin BPT element to accounting for transverse shear deformation effects of relevance for thick plates and composite laminated plates. The nodal deflection DOFs of the original BPT element are enhanced with the two shear deformation angles. This allows to computing the bending and shear deformation energies leading to a simple triangular plate element with 3 DOFs per node (termed BPT+ element). For the thin plate case the shear angles vanish and the element reproduces the good behaviour of the original thin BPT element. As a consequence the element is applicable to thick and thin plate situations without exhibiting shear locking. It is interesting that the thick plate solution can be found iteratively starting from the deflection values obtained using the standard Kirchhoff theory and the original thin BPT element. The ingredients of the formulation are: a Hu-Washizu type mixed functional, a linear interpolation for the deflection and the shear angles over 3-noded triangles  and a finite volume approach for computing the bending moments and the curvatures over a patch of elements. Details of the element formulation are given in the paper.

The rotation-free formulation described for the BPT+ element is taken as the starting point for deriving a two-noded rotation-free beam element with shear deformation effects (termed CCB+ element). The examples presented in the last part of the paper show the robustness and accuracy of the BPT+ and the CCB+ elements for thick and thin plate and beam problems.



==2 BASIC THEORY==

===2.1 Reissner-Mindlin plate theory===

Let us consider the plate of Figure [[#img-1|1]]. We will assume Reissner-Mindlin conditions to hold, i.e.

<span id="eq-1"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>u(x,y,z)=-z \theta _x (x,y)\quad ,\quad v(x,y,z)=-z \theta _y (x,y)\quad ,\quad  w(x,y,z)=w (x,y)  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (1)
|}

<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_940614781-Fig1.png|480px|Sign convenion for the deflection and the rotations in a plate]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 1:''' Sign convenion for the deflection and the rotations in a plate
|}

with

<span id="eq-2.a"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\theta _x = {\partial w \over \partial x} +\phi _x   \quad ,\quad \theta _y = {\partial w \over \partial y} +\phi _y   </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (2.a)
|}

or

<span id="eq-2.b"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{\boldsymbol \theta } = [\theta _x,\theta _y]^T = {\boldsymbol \nabla } w + {\boldsymbol \phi }   </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (2.b)
|}

with

<span id="eq-2.c"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{\boldsymbol \nabla } = \left[{\partial  \over \partial x},{\partial  \over \partial y}\right]^T \quad \hbox{and} \quad   {\boldsymbol \phi }=[\phi _x,\phi _y]^T   </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (2.c)
|}

where <math display="inline">u,v,w</math> are the cartesian displacements <math display="inline">\theta _x</math>, <math display="inline">\theta _y</math> are the   rotations and <math display="inline">\phi _x</math>, <math display="inline">\phi _y</math> are angles coinciding with the  transverse shear deformations as shown below.

The generalized bending and shear strain vectors  are defined as <span id='citeF-29'></span><span id='citeF-30'></span>[[#cite-29|[29,30]]]

<span id="eq-3.a"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{\boldsymbol \varepsilon }_b = [\kappa _x,\kappa _y,\kappa _{xy}]^T =   \begin{bmatrix}\displaystyle                           {\partial \theta _x \over \partial x}, & \displaystyle  {\partial \theta _y \over \partial y}, &                          \displaystyle  \left({\partial \theta _x \over \partial y}+{\partial \theta _y \over \partial x} \right)                        \end{bmatrix}^T ={\boldsymbol L}_b {\boldsymbol \theta }  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (3.a)
|}

<span id="eq-3.b"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{\boldsymbol \varepsilon }_s =[\gamma _{xz},\gamma _{yz}]^T= \begin{bmatrix}\displaystyle  {\partial w \over \partial x}-\theta _x , & \displaystyle  {\partial w \over \partial y}  -\theta _y                        \end{bmatrix}^T ={\boldsymbol L}_s {\boldsymbol u}  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (3.b)
|}

where <math display="inline"> [\kappa _x,\kappa _y,\kappa _{xy}]</math> and <math display="inline"> [\gamma _{xz},\gamma _{yz}]</math> denote the bending  strains (typically called curvatures) and the transverse shear strains, respectively and

<span id="eq-4"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{\boldsymbol L}_b = \begin{bmatrix}\displaystyle {\partial  \over \partial x} & 0 \\[.3cm]                0 & \displaystyle {\partial  \over \partial y} \\[.3cm]                \displaystyle {\partial  \over \partial y} & \displaystyle {\partial  \over \partial x} \\              \end{bmatrix}  \quad ,\quad  {\boldsymbol L}_s = \begin{bmatrix}\displaystyle {\partial  \over \partial x} & -1 & 0 \\[.3cm]\displaystyle {\partial  \over \partial y} & 0 & -1 \\               \end{bmatrix} \quad ,\quad  {\boldsymbol u} =           \left\{\begin{array}{c}w \\                    \theta _x \\ \theta _y                  \end{array}             \right\}   </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (4)
|}

Substituting Eq.([[#eq-2.a|2.a]]) into ([[#eq-3.b|3.b]]) gives

<span id="eq-5"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\gamma _{xz} =-\phi _x \quad \quad ,\quad \quad \gamma _{yz}= -\phi _y  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (5)
|}

i.e.  the angles <math display="inline">\phi _x</math>, <math display="inline">\phi _y</math> coincide (with opposite sign) with  the transverse shear deformations. In the following <math display="inline">\phi _x</math> and <math display="inline">\phi _y</math> will be called“shear angles”.

The set of governing equations can be expressed in integral form starting from the following Hu-Washizu type functional <span id='citeF-12'></span><span id='citeF-28'></span>[[#cite-12|[12,28]]]

<span id="eq-6"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\begin{array}{ll}\Pi = & \displaystyle {1\over 2} \left[\iint _A {\boldsymbol \varepsilon }_b^T {\boldsymbol D}_b {\boldsymbol \varepsilon }_b \,dA+ \iint _A {\boldsymbol \varepsilon }_s^T {\boldsymbol \sigma }_s \,dA  \right]+ \iint _A [ {\boldsymbol L}_b {\boldsymbol \theta } -  {\boldsymbol \varepsilon }_b]^T {\boldsymbol  \sigma }_b\,dA -\\ 
&- \iint _A  (qw+{\boldsymbol \theta }^T {\boldsymbol m}) dA - \sum \limits _{i=1}^{n_c}  w_i P_i    \end{array} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (6)
|}

where <math display="inline">q</math> is the distributed vertical loading, <math display="inline"> {\boldsymbol m}=[m_x,m_y]^T</math> are distributed bending moments, <math display="inline">P_i</math> is the vertical force acting at point <math display="inline">i</math>, <math display="inline">n_c</math> is the number of points with external concentrated forces and <math display="inline">A</math> is the area of the plate. Concentrated bending moments have been excluded for simplicity.

In Eq.([[#eq-6|6]]) <math display="inline"> {\boldsymbol \sigma }_b = [M_x,M_y,M_{xy}]^T</math> is the bending moment vector, <math display="inline">{\boldsymbol \sigma }_s = [Q_x,Q_y]^T </math> is the shear force vector and <math display="inline">{\boldsymbol D}_b</math> is the bending constitutive matrix given by (for the isotropy case)

<span id="eq-7.a"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{\boldsymbol D}_b =\displaystyle{{E}t^3\over   12(1-\nu ^2)} \left[             \begin{array}{ccc}1 & \nu & 0 \\               \nu & 1 & 0 \\               0 & 0 & {1-\nu \over 2} \\             \end{array}           \right]  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (7.a)
|}

The shear forces and the shear strains are assumed to be related point-wise by the standard constitutive equation

<span id="eq-7.b"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{\boldsymbol \sigma }_s =  {\boldsymbol D}_s {\boldsymbol \varepsilon }_s \quad \hbox{with} \quad   {\boldsymbol D}_s = \alpha Gt \left[  \begin{array}{cc}1 & 0 \\  0 & 1 \\  \end{array} \right]  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (7.b)
|}

In Eqs.([[#2.1 Reissner-Mindlin plate theory|2.1]]) <math display="inline">t</math> is the plate thickness and <math display="inline">E,\nu </math> and <math display="inline">G</math> are the Young modulus, the Poisson's ratio and the shear modulus, respectively and <math display="inline">\alpha </math> is the shear correction factor; <math display="inline">\alpha =5/6</math> is taken for an isotropic plate.

===2.2 Using the deflection and the shear angles as main variables===

Let us express Eqs.([[#eq-3|3]])&#8211;([[#eq-6|6]]) in terms of the deflection <math display="inline">w</math> and the shear angles <math display="inline">\phi _x</math> and <math display="inline">\phi _y</math>. The resulting Reissner-Mindlin theory contains the standard expressions of Kirchhoff thin plate theory plus additional transverse shear deformation terms.

Substituting Eq.([[#eq-2|2]]) into Eqs.([[#eq-3|3]]) gives

<span id="eq-8"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{\boldsymbol \varepsilon }_b={\boldsymbol L}_b {\boldsymbol\theta }= \left[  \frac{\partial^2 w}{\partial x^2}+{\partial \phi _x \over \partial x},  \frac{\partial^2 w}{\partial y^2}  + {\partial \phi _y \over \partial y}, 2 \frac{\partial^2 w} {x\partial y} + \left({\partial \phi _x \over \partial y}  +{\partial \phi _y \over \partial x}  \right)\right]^T = {\boldsymbol L}_w w +  {\boldsymbol L}_b {\boldsymbol \phi }</math>
|-
| style="text-align: center;" | <math>  {\boldsymbol \varepsilon }_s={\boldsymbol L}_s {\boldsymbol u} =  [-\phi _x , -\phi _y]^T =  -{\boldsymbol \phi }  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (8)
|}

with

<span id="eq-9"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{\boldsymbol L}_w = \left[ \frac{\partial^2}{\partial x^2},  \frac{\partial^2}{\partial y^2},  2 \frac {\partial^2} {x\partial y}\right]^T  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (9)
|}

Substituting Eqs.([[#eq-2.b|2.b]]) and ([[#eq-8|8]]) into the functional of Eq.([[#eq-6|6]]) gives

<span id="eq-10"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\begin{array}{ll}\Pi = &\!\!\!\! \displaystyle {1\over 2} \left[\iint _A {\boldsymbol \varepsilon }_b^T {\boldsymbol D}_b {\boldsymbol \varepsilon }_b \,dA - \iint _A {\boldsymbol \phi }^T {\boldsymbol \sigma }_s \right] + \iint _A [ {\boldsymbol L}_w w + {\boldsymbol L}_b {\boldsymbol \phi } - {\boldsymbol \varepsilon }_b]^T  {\boldsymbol \sigma }_b\,dA +\\[.4cm] &\displaystyle -\!\! \iint _A  [qw+({\boldsymbol \nabla }w+ {\boldsymbol \phi })^T {\boldsymbol m}]dA - \!\sum \limits _{i=1}^{n_c} w_i P_i  \end{array}   </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (10)
|}

Note that by making the shear angles <math display="inline">{\boldsymbol \phi }</math> equal to zero we recover precisely the variational form of standard Kirchhoff thin plate theory <span id='citeF-29'></span>[[#cite-29|[29]]].

Variation of <math display="inline">\Pi </math> with respect to <math display="inline">{\boldsymbol \varepsilon }_b</math>, <math display="inline">{\boldsymbol \sigma }_b</math>, <math display="inline">w</math> and <math display="inline">{\boldsymbol \phi }</math> leads to the following equations:

===Bending constitutive equation===

<span id="eq-11"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\iint _A \delta {\boldsymbol \varepsilon }_b^T [{\boldsymbol D}_b {\boldsymbol \varepsilon }_b -{\boldsymbol \sigma }_b ]\, dA =0  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (11)
|}

===Relationship between bending strains, deflection and shear angles===

<span id="eq-12"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\iint _A \delta {\boldsymbol \sigma }_b^T [{\boldsymbol L}_w w + {\boldsymbol L}_b {\boldsymbol \phi }-{\boldsymbol \varepsilon }_b ]\, dA =0</math>
|-
| style="text-align: center;" | 
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (12)
|}

===Equilibrium equations===

<span id="eq-13.a"></span>
<span id="eq-13.b"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\iint _A [{\boldsymbol L}_w \delta w]^T {\boldsymbol \sigma }_b \, dA - \iint _A [\delta w q+ ({\boldsymbol \nabla } \delta w)^T {\boldsymbol m}] \, dA - \sum \limits _{i=1}^{n_c} \delta w_i P_i =0</math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (13.a)
|-
| style="text-align: center;" | <math> \iint _A \left[({\boldsymbol L}_b \delta {\boldsymbol \phi })^T {\boldsymbol \sigma }_b  - \delta  {\boldsymbol \phi }^T {\boldsymbol \sigma }_s\right]dA - \iint _A \delta {\boldsymbol \phi }^T {\boldsymbol m} \, dA  =0  </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (13.b)
|}
|}

These integral expressions, completed with the constitutive equation for the shear forces (Eq.([[#eq-7.b|7.b]])), are the basis for deriving the element equations.

==3 THREE-NODED BASIC ROTATION-FREE PLATE TRIANGLE WITH SHEAR DEFORMATION EFFECTS==

Let us consider an arbitrary discretization of the plate into standard three-noded triangles. We  assume a linear interpolation of the deflection <math display="inline">w</math> and the shear angles <math display="inline">{\boldsymbol \phi }</math> within each element in terms of the nodal values in the standard manner <span id='citeF-28'></span>[[#cite-28|[28]]]

<span id="eq-14"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>w =\sum \limits _{i=1}^3 N_i \bar w_i \quad ,\quad {\boldsymbol \phi }= \left\{\begin{array}{c}\phi _x \\                                     \phi _y                                   \end{array}  \right\}= \sum \limits _{i=1}^3 \mathbf{N}_i \bar {\boldsymbol \phi }_i={\boldsymbol N}_\phi \bar  {\boldsymbol \phi }^e   </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (14)
|}

where <math display="inline">N_i</math> are the  linear shape functions of the three-noded triangle <span id='citeF-28'></span><span id='citeF-30'></span>[[#cite-28|[28,30]]], <math display="inline">\bar{(\cdot )}</math> denotes nodal variables and

<span id="eq-15.a"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{N}_\phi = [\mathbf{N}_1,\mathbf{N}_2,\mathbf{N}_3]\quad ,\quad  \mathbf{N}_i=N_i {\boldsymbol I}_2 \quad ,\quad \mathbf{I}_2 = \left[\begin{array}{cc}1 & 0 \\                                                                 0 & 1 \\                                                               \end{array}  \right]  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (15.a)
|}

<span id="eq-15.b"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\bar {\boldsymbol \phi }^e= \left\{ \begin{array}{c}\bar {\boldsymbol \phi }_1 \\ \bar {\boldsymbol \phi }_2 \\ \bar {\boldsymbol \phi }_3 \\ \end{array} \right\} \quad ,\quad \bar {\boldsymbol \phi }_i= \left\{\begin{array}{c}\bar{\phi }_{x_i} \\                                     \bar{\phi }_{y_i}                                   \end{array} \right\}   </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (15.b)
|}

with

<span id="eq-15.c"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math> \begin{array}{c}   N_i = \displaystyle{1\over 2A^e} (a_i^e+b_i^e x+c_i^e y),  \\[.3cm]   a_i^e=x_j^e y_k^e -x_k^e y_j^e \, ,\, b_i^e =y_{jk}\, ,\,  c_i^e = x_{kj}\quad \hbox{with }y_{jk}= y_j^e -y_k^e\, ,\,x_{kj}=x_k^e -x_j^e \quad  i,j,k=1,2,3  \end{array}    </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (15.c)
|}

and <math display="inline">A^e</math> is the element area <span id='citeF-28'></span>[[#cite-28|[28]]]. In the above expressions and in the following super-index <math display="inline">e</math> denotes element values.

Note that the introduction of the linear approximation for <math display="inline">w</math> into the variational form ([[#eq-13.a|13a]]) will give a vanishing of the virtual bending energy term involving second derivatives of <math display="inline">w</math>. This problem is overcome by “relaxing” further the weak form using an assumed constant value for the curvatures and bending moments over the triangle and computing the curvatures from the integral of the slopes along the element sides as explained next.

===3.1 Assumed curvatures and bending moments field===

We describe the bending moments <math display="inline">{\boldsymbol \sigma }_b</math> and the curvatures <math display="inline">{\boldsymbol \varepsilon }_b </math> (and their virtual values)   by ''constant fields'' within the triangle, i.e.

<span id="eq-16"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\begin{array}{ccc}{\boldsymbol \sigma }_b =  {\boldsymbol \sigma }_b^e & , &   {\boldsymbol \varepsilon }_b =  {\boldsymbol \varepsilon }_b^e\\   \delta{\boldsymbol \sigma }_b =  \delta{\boldsymbol \sigma }_b^e & , &  \delta {\boldsymbol \varepsilon }_b =  \delta{\boldsymbol \varepsilon }_b^e \end{array}  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (16)
|}

where <math display="inline">(\cdot )^e</math> denotes constant values within the triangle.

We define the term “patch of triangular elements”. This is typically formed by four elements: a central triangle and the three adjacent triangles. Central triangles adjacent to a boundary lack  the element adjacent to the boundary side (Figure [[#img-2|2]]).

Eqs.([[#eq-16|16]]) are introduced into the governing integral equations ([[#eq-11|11]])&#8211;([[#eq-13|13]]) which are modified as follows.

===Constitutive equation for the bending moments===

<span id="eq-17"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\sum \limits _e \iint _{A^e} (\delta {\boldsymbol \varepsilon }_b^e)^T [{\boldsymbol D}_b {\boldsymbol \varepsilon }_b^e -{\boldsymbol \sigma }_b^e ]\, dA =0</math>
|-
| style="text-align: center;" | 
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (17)
|}

Recalling that the virtual curvatures are arbitrary we obtain the following constitutive equations for the assumed (constant) bending moments for each element

<span id="eq-18"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{\boldsymbol \sigma }_b^e ={\boldsymbol D}_b^e {\boldsymbol \varepsilon }_b^e \quad ,\quad {\boldsymbol  D}_b^e = {1\over A^e}\iint _A {\boldsymbol D}_b\, dA  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (18)
|}

where <math display="inline">{\boldsymbol D}_b^e</math> is the average bending constitutive matrix for the triangular element.

===Relationship between the bending strains, the deflection and the shear angles===

Substituting Eqs.([[#eq-16|16]]) into Eq.([[#eq-12|12]]) gives

<span id="eq-19"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\sum \limits _e \iint _{A^e} (\delta {\boldsymbol \sigma }_b^e)^T [{\boldsymbol L}_w w + {\boldsymbol L}_b {\boldsymbol  \phi }-{\boldsymbol \varepsilon }_b^e ]\, dA =0  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (19)
|}

<div id='img-2a'></div>
<div id='img-2b'></div>
<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_940614781-fig2a.png|348px|]]
|[[Image:Draft_Samper_940614781-fig2b.png|450px|]]
|- style="text-align: center; font-size: 75%;"
| (a) 
| (b) 
|- style="text-align: center; font-size: 75%;"
| colspan="2" | '''Figure 2:''' Triangular element patches. (a) Standard four element patch; (b) Three element boundary patch. Numbers 1,2,3 in brackets denote local node numbers for the element
|}

As the virtual bending moments are arbitrary and <math display="inline">{\boldsymbol \varepsilon }_b^e</math> is constant within each element we get

<span id="eq-20"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{\boldsymbol \varepsilon }_b^e =   {1\over A^e} \iint _{A^e}[{\boldsymbol L}_w w + {\boldsymbol L}_b {\boldsymbol \phi }]\,dA  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (20)
|}

The term <math display="inline"> {\boldsymbol L}_w w</math> in Eq.([[#eq-20|20]]) is ''integrated by parts ''to give

<span id="eq-21"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\displaystyle {\boldsymbol \varepsilon }_b^e = {1\over A^e} \left[\int _{\Gamma ^e} {\boldsymbol T} {\boldsymbol \nabla } w \, d\Gamma +  \iint _{A^e} {\boldsymbol L}_b  {\boldsymbol \phi }\, dA   \right]  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (21)
|}

where

<span id="eq-22"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{\boldsymbol T}= \begin{bmatrix}n_x & 0\\             0 & n_y\\ n_y & n_x           \end{bmatrix}            </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (22)
|}

and <math display="inline">n_x,n_y</math> are the components to the outward unit normal <math display="inline">\boldsymbol n</math> to the boundary of the element <math display="inline">\Gamma ^e</math> (Figure [[#img-2|2]]). The transformation of the integral of the curvature field over the element domain into the integral of the deflection gradient along the element boundary is a distinct feature of the BPT formulation <span id='citeF-6'></span><span id='citeF-12'></span><span id='citeF-15'></span><span id='citeF-19'></span><span id='citeF-20'></span>[[#cite-6|[6,12,15,19,20]]].

Eq.([[#eq-21|21]]) defines the bending strains as the sum of the integral of the deflection gradient along the boundary of the element and the integral over the element of the term <math display="inline">{\boldsymbol L}_b {\boldsymbol \phi }</math> including the gradients of the shear deformation angles. This term is constant within each element for a linear interpolation of <math display="inline">\boldsymbol \phi </math>.

===Shear forces-shear angles relationship===

Making <math display="inline">{\boldsymbol \varepsilon }_s= -{\boldsymbol \phi }</math> (Eq.(8)) and substituting the approximation for <math display="inline">{\boldsymbol \phi }</math> of Eq.([[#eq-14|14]]) into Eq.([[#eq-7.b|7.b]]) gives

<span id="eq-23"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{\boldsymbol \sigma }_s = {\boldsymbol D}_s {\boldsymbol \varepsilon }_s =- {\boldsymbol D}_s{\boldsymbol \phi } = - {\boldsymbol D}_s {\boldsymbol N}_\phi \bar {\boldsymbol \phi }^e  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (23)
|}

===Equilibrium equations===

Integrating by parts the first integral in Eq.([[#eq-13.a|13.a]]) and recalling that the bending moments are constant within each element and the additive property of the element integrals in the FEM, allows us to write Eqs.([[#eq-13|13]]) as

<span id="eq-24"></span>
<span id="eq-25"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\sum \limits _e  \left(\int _{\Gamma ^e} {\boldsymbol T} ({\boldsymbol \nabla }\delta w)^T\,d\Gamma \right){\boldsymbol \sigma }_b^e  - \iint _A \left[\delta w q+ \left({\boldsymbol \nabla } \delta w\right)^T {\boldsymbol m}\right]\, dA + \sum \limits _{i=1}^{n_c} \delta w_i P_i =0</math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (24)
|-
| style="text-align: center;" | <math> \sum \limits _e \left\{\left(\iint _{A^e} ({\boldsymbol L}_b {\boldsymbol \phi })^T dA\right){\boldsymbol \sigma }_b^e - \iint _{A^e} \delta  {\boldsymbol \phi }^T {\boldsymbol \sigma }_s dA \right\} - \iint _A \delta {\boldsymbol \phi }^T {\boldsymbol m}  \, dA  =0  </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (25)
|}
|}

The sum in above expressions extends over all the elements in the mesh.

===3.2 Element matrices and vectors===

===Bending strain matrices===

Substituting the linear interpolation for <math display="inline">w</math> and <math display="inline">\boldsymbol \phi </math> of Eqs.([[#eq-14|14]]) into the integral expressions for the constant bending strain field of Eq.([[#eq-21|21]]) gives

<span id="eq-26.a"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{1\over A^e}\int _{\Gamma ^e} {\boldsymbol T}{\boldsymbol \nabla }w d\Gamma ={\boldsymbol B}_w \bar {\boldsymbol w}^e  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (26.a)
|}

<span id="eq-26.b"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{1\over A^e}\iint _{A^e}{\boldsymbol L}_b {\boldsymbol \phi } dA ={\boldsymbol B}_{\phi }\bar{\boldsymbol \phi }^e  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (26.b)
|}

and, therefore

<span id="eq-27"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{\boldsymbol \varepsilon }_b^e = {\boldsymbol B}_w\bar {\boldsymbol w}^e + {\boldsymbol  B}_{\phi } \bar{\boldsymbol \phi }^e  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (27)
|}

where <math display="inline"> {\boldsymbol  B}_w</math> and <math display="inline"> {\boldsymbol B}_{\phi }</math>  are generalized bending strain matrices  and <math display="inline"> \bar {\boldsymbol w}^e</math> and <math display="inline">\bar{\boldsymbol \phi }^e </math> are nodal deflection and nodal shear angles vectors for the element. The expression for <math display="inline">\bar{\boldsymbol \phi }^e</math> is given in Eq.([[#eq-15.b|15.b]]). The expression for <math display="inline">\bar {\boldsymbol w}^e</math> is given below (Eq.([[#eq-31|31]])).

The computation of <math display="inline">{\boldsymbol B}_\phi </math> is straightforward from the second integral of Eq.([[#eq-21|21]]) as

<span id="eq-28"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{\boldsymbol B}_{\phi } = [ {\boldsymbol B}_{\phi _1},{\boldsymbol B}_{\phi _2},{\boldsymbol   B}_{\phi _3}]\qquad \hbox{with } \qquad {\boldsymbol B}_{\phi _i} = {1\over A^e}\iint _{A^e} {\boldsymbol L}_b  {N}_i\, d\Omega    </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (28)
|}

Recalling that <math display="inline">N_i</math> are linear shape functions gives (using Eq.(15c))

<span id="eq-29"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{\boldsymbol B}_{\phi _i} = \begin{bmatrix}\displaystyle {\partial N_i \over \partial x} & 0 \\[.25cm]                        0 & \displaystyle {\partial N_i \over \partial x} \\[.45cm]                       \displaystyle {\partial N_i \over \partial y} & \displaystyle {\partial N_i \over \partial x} \\                      \end{bmatrix}   = {1\over 2A^e} \begin{bmatrix}b_i^e & 0 \\                                                        0 & c_i^e \\                                                        c_i^e & b_i^e \\                                                      \end{bmatrix}   </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (29)
|}

The computation of <math display="inline">{\boldsymbol B}_w^e</math> in Eq.([[#eq-26.a|26.a]]) requires integrating the product  of the components of the normal vector and the deflection slopes along the element sides.  For the 3-noded triangle this is not so straightforward as <math display="inline">{\boldsymbol \nabla }w</math> is discontinuous at the element sides for a linear approximation of <math display="inline">w</math>. A simple way to overcoming this problem is to compute the deflection gradients at the 

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{\boldsymbol B}_w={1\over A^e} \left[\begin{matrix} y_{ij}\bar b_i^{b}+ y_{ki}\bar b_i^{d}& y_{ij}\bar  b_j^{b}+ y_{jk}\bar b_j^{c}& y_{jk}\bar b_k^{c}+ y_{ki}\bar b_k^{d}\\ -x_{ij}\bar c_i^{b}-x_{ki}\bar c_i^{d}& -x_{ij}\bar c_j^{b}-x_{jk}\bar c_j^{c}& -x_{jk}\bar c_k^{c}-x_{ki}\bar c_k^{d}\\[] [y_{ij}\bar c_i^{b}- x_{ij}\bar  b_i^{b}& [y_{ij}\bar c_j^{b}- x_{jk}\bar b_j^{b}& [y_{jk}\bar  c_k^{c}-x_{jk}\bar b_k^{c}\\ +y_{ki}\bar c_i^{d}-x_{ki}\bar b_i^{d}]& +y_{jk}\bar c_j^{c}-x_{jk}\bar b_j^{c}]& +y_{ki}\bar c_k^{d}-x_{ki}\bar  b_k^{d}]\end{matrix}\right.</math>
|}
|}

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math> \left.\begin{matrix} y_{ij}\bar b_l^{b}& y_{jk}\bar b_m^{c}& y_{ki}\bar b_n^{d}\\  -x_{ij}\bar c_l^{b}& -x_{jk}\bar c_m^{c}& -x_{ki}\bar c_n^{d}\\  y_{ij}\bar c_l^{b}-x_{ij}\bar b_l^{b}& y_{jk}\bar c_m^{c}-x_{jk}\bar b_m^{c}& y_{ki}\bar c_n^{d}-x_{ki}\bar c_n^{d}\\ \end{matrix}\right]</math>
|}
|}

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\bar b_i^{p}={{b_i^p}\over 2A^{p}},\quad \bar  c_i^{p}={{c_i^{p}}\over 2A^{p}},\quad  b_i^{p}=y_j^{p}-y_k^{p}\quad c_i^p=x_k^{p}-x_j^{p}, \hbox{etc.},\quad  p=a,b,c\quad i,j,k=1,2,3</math>
|}
|}


element sides as the average value of the gradients contributed by the two triangles adjacent to the side <span id='citeF-6'></span><span id='citeF-12'></span>[[#cite-6|[6,12]]]. This gives

<span id="eq-30"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{1\over A^e} \int _{\Gamma ^e} {\boldsymbol T} {\boldsymbol \nabla } w d\Gamma = {1\over A^e} \sum \limits _{j=1}^3 {l_j^e\over 2} {\boldsymbol T}_j \left[     \sum \limits _{i=1}^3   {\boldsymbol \nabla } N_i^e w_i^e + \sum \limits _{i=1}^3     {\boldsymbol \nabla } N_i^p w_i^p \right]_j =  {\boldsymbol B}_w \bar{\boldsymbol w}^e   </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (30)
|}

with

<span id="eq-31"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\bar{\boldsymbol w}^e =[\bar{w}^e_i, \bar{w}^e_j,\bar{w}^e_k,\bar{w}^e_l,\bar{w}^e_m,  \bar{w}^e_n]^T    </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (31)
|}

The sum in Eq.([[#eq-30|30]]) extends over the three sides of an element <math display="inline">e</math>, <math display="inline">{\boldsymbol T}_j</math> is the transformation matrix of Eq.([[#eq-22|22]]) for side <math display="inline">j</math>, <math display="inline">l^e_j</math> are the lengths of the element sides and superindex <math display="inline">p</math> refers to each of the three triangles <math display="inline">a,b,c</math> adjacent to the central triangle <math display="inline">e</math> with <math display="inline">p=a,b,c</math> for <math display="inline">j=1,2,3</math> (Figure [[#img-2|2]]). Sides are assigned the number of the opposite node. Thus, side <math display="inline">l^e_3</math> is opposite to node 3, etc.

Note that <math display="inline">{\boldsymbol B}_w</math> is a <math display="inline">3\times 6</math> matrix relating the three bending strains with the deflections at the six nodes of the patch of triangles linked to element <math display="inline">e</math>. The explicit form of matrix <math display="inline">{\boldsymbol B}_w</math> is given in Box I. This matrix coincides with the curvature matrix of the original BPT element, as presented in <span id='citeF-12'></span>[[#cite-12|[12]]].

Substituting Eq.([[#eq-27|27]]) into ([[#eq-18|18]]) gives the relationship between the bending moments and the nodal values of the deflection and the shear angles for the element as

<span id="eq-32"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{\boldsymbol \sigma }_b^e =  {\boldsymbol D}_b^e [ {\boldsymbol B}_w \bar {\boldsymbol w}^e+   {\boldsymbol B}_{\phi }  \bar{\boldsymbol \phi }^e ]  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (32)
|}

===Stiffness matrices and nodal force vectors===

Substituting Eqs.([[#eq-14|14]]), ([[#eq-23|23]]), ([[#eq-26|26]]) and ([[#eq-32|32]]) into the equilibrium equations ([[#eq-24|24]]) and ([[#eq-25|25]]) yields

<span id="eq-33.a"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\sum \limits _e (\delta \bar {\boldsymbol w}^e)^T {\boldsymbol B}_w^T  {\boldsymbol  D}_b^e \left[{\boldsymbol B}_w \bar {\boldsymbol w}^e +{\boldsymbol B}_\phi \bar{\boldsymbol \phi }^e\right]A^e -</math>
|-
| style="text-align: center;" | <math> - \iint _A \left[\left(\sum \limits _{i=1}^3 N_i \delta \bar w_i\right)q + \left( \sum \limits _{i=1}^3 {\boldsymbol \nabla } N_i \delta \bar{w}_i\right)^T {\boldsymbol  m}\right]\,dA + \sum \limits _{i=1}^{n_c} \delta \bar w_i P_i=0 </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (33.a)
|}

<span id="eq-33.b"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\sum \limits _e (\delta \bar {\boldsymbol \phi }^e)^T \left\{{\boldsymbol B}_\phi ^T {\boldsymbol D}_b^e \left[{\boldsymbol B}_w \bar {\boldsymbol w}^e +{\boldsymbol  B}_\phi \bar{\boldsymbol \phi }^e\right]A^e +\left(\iint _{A^e} {\boldsymbol N}_\phi ^T {\boldsymbol D}_s {\boldsymbol N}_\phi  dA \right)\bar {\boldsymbol \phi }^e \right\}-</math>
|-
| style="text-align: center;" | <math> -\iint _A \left(\sum \limits _{i=1}^3 {\boldsymbol N}_i \delta \bar{\boldsymbol \phi }_i\right)^T {\boldsymbol m}\,dA =0  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (33.b)
|}

After simplification of the virtual nodal deflections and the virtual nodal shear   angles we obtain finally the matrix system of equilibrium equations

<span id="eq-34"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\begin{array}{c} {\boldsymbol K}_w \bar{\boldsymbol w} + {\boldsymbol K}_{w\phi } \bar{\boldsymbol \phi } ={\boldsymbol f}_w\\[.4cm] {\boldsymbol K}_{w\phi }^T \bar{\boldsymbol w} + {\boldsymbol K}_{\phi } \bar{\boldsymbol \phi } = {\boldsymbol f}_\phi  \end{array}   </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (34)
|}

where vectors <math display="inline">\bar {\boldsymbol w}</math> and <math display="inline">\bar {\boldsymbol \phi }</math> contain the nodal deflections and the nodal shear angles for the whole mesh and the rest of the  matrices and vectors are assembled from the element contributions given by

<span id="eq-35"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\begin{array}{l}{\boldsymbol K}_{w_{ij}}^e = \displaystyle  {\boldsymbol B}_{w_i}^T {\boldsymbol D}_b^e {\boldsymbol B}_{w_j} A^e  \quad ,\quad    {\boldsymbol K}_{w\phi _{ij}}^e  =  {\boldsymbol B}_{w_i}^T {\boldsymbol D}_b^e {\boldsymbol    B}_{\phi _j} A^e\\[.4cm] {\boldsymbol K}_{\phi _{ij}}^e= \displaystyle {\boldsymbol B}_{\phi _i}^T  {\boldsymbol  D}_b^e {\boldsymbol    B}_{\phi _j} A^e + \iint _{A^e} {\boldsymbol N}_{i}^T {\boldsymbol D}_s {\boldsymbol N}_{j}\,dA   \end{array}     </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (35)
|}

The element stiffness matrices <math display="inline">{\boldsymbol K}_{w}^e</math> and <math display="inline">{\boldsymbol K}_{w\phi }^e</math>   can be explicitly computed from the strain and constitutive matrices for the element. The exact computation of the integral in the expression of <math display="inline">{\boldsymbol K}_\phi </math> requires a 3 Gauss point quadrature. Excellent results have been obtained in all examples analyzed using a reduced one point integration rule for <math display="inline">{\boldsymbol K}_{\phi }</math> which indicates no advantage in  using a full quadrature for <math display="inline">\mathbf{K}_\phi </math>. The simple one point quadrature allows computing all the element stiffness matrices in Eq.([[#eq-34|34]]) explicitly.

The equivalent nodal force vectors for the element are

<span id="eq-36"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{f}^e_{w_i}= \iint _{A^e} [ N_i q + ({\boldsymbol \nabla }N_i)^T {\boldsymbol m}]dA \quad ; \quad \displaystyle {\boldsymbol f}^e_{\phi _i}= \iint _{A^e} \mathbf{N}_i {\boldsymbol m}dA     </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (36)
|}

For a uniform distribution of <math display="inline">q</math> and <math display="inline">{\boldsymbol m}</math>

<span id="eq-37"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\begin{array}{l}\displaystyle {f}^e_{w_i}={qA^e\over 3} + {1\over 2} (b_i^e m_x + c^e_i m_y)\\[.4cm] \displaystyle {\boldsymbol f}^e_{\phi _i}= {A^e\over 3} [m_x,m_y]^T  \end{array}    </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (37)
|}

where <math display="inline">b_i^e</math> and <math display="inline">c_i^e</math> are given in Eq.([[#eq-15|15]]).

The  vertical load <math display="inline">P_i</math> acting at a node <math display="inline">i</math> contributes the following terms to the <math display="inline">i</math>th component of the global force vector <math display="inline">{\boldsymbol f}_w</math>

<span id="eq-38"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{f}_{w_i} =  P_i  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (38)
|}

The reactions (the vertical force and the bending moment) can be computed at the prescribed nodes from the nodal displacement and rotations field in the usual manner.

'''Remark 1''' . Each node has three DOFs: the deflection <math display="inline">\bar w_i</math> and the two shear angles <math display="inline">\bar {\phi }_{x_i}</math> and <math display="inline">\bar {\phi }_{y_i}</math> at the node. The size of all the element stiffness matrices in Eqs.([[#eq-35|35]]) is <math display="inline">6\times 6</math>. This means that the effective DOFs for each element are in fact the nodal deflections of the patch of four elements assigned to each element (typically six nodal deflections except for boundary elements) and the six nodal shear angles, i.e. two shear angles for each of the three nodes of the element.

===3.3 Iterative computation of the deflection and the shear angles at the nodes===

Despite that the solution for the <math display="inline">\bar{\boldsymbol w}</math> and <math display="inline">\bar {\boldsymbol \phi }</math> variables can be found simultaneously by solving Eqs.([[#eq-34|34]]), the following iterative algorithm is recommended for computing <math display="inline">\bar{\boldsymbol w}</math> and <math display="inline">\bar  {\boldsymbol \phi }</math>

<span id="eq-39.a"></span>
<span id="eq-39.b"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\displaystyle {\boldsymbol K}_{w} \bar{\boldsymbol w}^i \!\!\!=\!\!\!{\boldsymbol f}_w - {\boldsymbol  K}_{w\phi }\bar{\boldsymbol \phi }^{i-1} \rightarrow \bar {\boldsymbol w}^i</math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (39.a)
|-
| style="text-align: center;" | <math> \displaystyle {\boldsymbol K}_{\phi } \bar {\boldsymbol \phi }^i \!\!\!=\!\!\!{\boldsymbol  f}_\phi - {\boldsymbol K}_{w\phi }^T \bar{\boldsymbol w}^{i} \rightarrow \bar  {\boldsymbol \phi }^i      </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (39.b)
|}
|}

where superindex <math display="inline">i</math> denotes the number of iterations. The iterative scheme of Eqs.([[#eq-39|39]]) continues until convergence for the nodal deflection and the nodal shear deformation angles is achieved. Convergence is typically measured by the <math display="inline">L_2</math> norm of vectors <math display="inline">\bar{\boldsymbol w}</math> and <math display="inline">\bar  {\boldsymbol \phi }</math> (Eqs.([[#eq-58|58]]) and ([[#eq-59|59]])). An advantage of the above iterative scheme is that for <math display="inline">i=1</math> and <math display="inline">{\boldsymbol \phi }^0={\boldsymbol 0}</math> the value of <math display="inline">\bar{\boldsymbol w}^1</math> corresponds to the thin plate solution of Kirchhoff theory which is accurate enough for many practical cases. ''The effect of shear deformation is  introduced progressively with the number of iterations''. Shear deformation effects are negligible for thin plates and hence the <math display="inline">\bar{\boldsymbol \phi }</math> variables tend rapidly to zero in this case.

===3.4 Improved iterative scheme===

An enhanced iterative scheme can be devised by performing a smoothing of the shear angles field as follows.

===''Step 1 Computation of the nodal deflections'' ̄w¹===

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{\boldsymbol K}_w\bar{\boldsymbol w}^1 = {\boldsymbol f}_w \to \bar {\boldsymbol w}^1 \quad  \hbox{(Kirchhoff thin plate solution)}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (40)
|}

===''Step 2 Compute'' ̄ϕⁱ, i ≥ 1===

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\displaystyle {\boldsymbol K}_{\phi } \bar {\boldsymbol \phi }^i ={\boldsymbol f}_\phi  - {\boldsymbol K}_{w\phi }^T \bar{\boldsymbol w}^{i} \qquad \rightarrow \bar  {\boldsymbol \phi }^i </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (41)
|}

===''Step 3 Compute the element shear angles ̄ϕ<sup>e</sup>''===

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\bar{\boldsymbol \phi }^{e} = \frac{1}{3} \left(\bar{\boldsymbol \phi }_a + \bar{\boldsymbol \phi }_b+\bar{\boldsymbol \phi }_c   \right) </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (42)
|}

where <math display="inline">\bar{\boldsymbol \phi }^{e}</math> are mean shear angles for element <math display="inline">e</math> for the <math display="inline">i</math>th iteration  and <math display="inline">a,b,c</math> are the ''global numbers'' of the three nodes of element <math display="inline">e</math>.

===''Step 4 Compute the smoothed nodal shear angles ̃ϕₖ''===

The smoothed nodal values <math display="inline">\tilde{\boldsymbol \phi }_k^i</math> are computed as

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\tilde{\boldsymbol \phi }_k^i = \frac{1}{n_k} \sum \limits _{j=1}^{n_k} \bar{\boldsymbol \phi }^{j} \quad , \quad k =1,N </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (43)
|}

where <math display="inline">n_k</math> is the number of elements sharing the node with global number <math display="inline">k</math> and <math display="inline">N</math> is the total number of nodes in the mesh.

===''Step 5 Compute ̄wⁱ, i > 1''===

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{\boldsymbol K}_w\bar{\boldsymbol w}^i ={\boldsymbol f}_w -{\boldsymbol K}_{w\phi }\tilde{\boldsymbol  \phi }^{i-1} \qquad \rightarrow \bar{\boldsymbol w}^i </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (44)
|}

===''Return to step 2.''===

Convergence of the improved iterative scheme for the deflection field  is quite fast (2&#8211;4 iterations),  even for thick plates (see Section [[#6 EXAMPLES|6]]).

The accuracy and efficiency of the so called BPT+ triangle is shown in the examples presented in the paper.

==4 BOUNDARY CONDITIONS==

A BPT element with a side along a boundary edge has one of the triangles belonging to the patch missing (Figure [[#img-2|2]]). This is taken into account by ignoring the contribution of this element when performing the average of the deflection gradient in Eq.([[#eq-30|30]]) <span id='citeF-12'></span><span id='citeF-15'></span><span id='citeF-19'></span><span id='citeF-20'></span>[[#cite-12|[12,15,19,20]]].

===4.1 <span id='lb-4.1'></span>Clamped edge (w=0, θ=0)===

The condition <math display="inline">w=0</math> is directly imposed at the edge nodes at the solution level when solving the system of equations for <math display="inline">\bar {\boldsymbol w}</math> in the standard manner.

The condition of zero rotations at a clamped edge introduces additional terms in the system of Eqs.([[#eq-34|34]]).

In order to explain the process, let us consider for example a clamped edge corresponding to side 3 linking nodes 1 and 2 (with global numbers <math display="inline">i</math> and <math display="inline">j</math>) in the boundary element of Figure [[#img-2|2]]b. At the clamped edge

<span id="eq-45"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{\boldsymbol \theta }\vert _3^e = {\boldsymbol \nabla } w\vert _3^e + {\boldsymbol \phi }\vert _3^e =0 \quad \to \quad {\boldsymbol \nabla } w\vert _3^e = -{\boldsymbol \phi }\vert _3^e   </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (45)
|}

where <math display="inline">{\boldsymbol \theta }\vert _3^e</math> are the rotations at the edge mid point.

The shear  angles at the edge mid point are approximated by the average of the nodal values along the edge, i.e.

<span id="eq-46"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{\boldsymbol \phi }\vert _3^e = \frac{1}{2} \left({\boldsymbol \phi }_2 + {\boldsymbol \phi }_3 \right)   </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (46)
|}

Introducing Eqs.([[#eq-45|45]]) and ([[#eq-46|46]]) into the definition of matrix <math display="inline">{\boldsymbol B}_w</math> in Eq.([[#eq-30|30]]) gives

<span id="eq-47"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{1\over A^e} \int _{\Gamma ^e} \!\! {\boldsymbol \nabla } w d\Gamma \!\!=\!\! {1\over A^e} \sum \limits _{j=1}^2 {l_j^e\over 2} {\boldsymbol T}_j \left[     \sum \limits _{i=1}^3   {\boldsymbol \nabla } N_i^e w_i^e + \sum \limits _{i=1}^3     {\boldsymbol \nabla } N_i^p w_i^p \right]_j -{1\over A^e}{l_3^e\over 2}{\boldsymbol T}_3 [{\boldsymbol \phi }_2 + {\boldsymbol \phi }_3]</math>
|-
| style="text-align: center;" | <math> \!\!=\!\! {\boldsymbol B}_w \bar{\boldsymbol w}^e + \Delta {\boldsymbol B}_\phi \bar{\boldsymbol \phi }^e   </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (47)
|}

where <math display="inline">{\boldsymbol B}_w</math> is obtained by disregarding the contributions from the clamped side in the sum along the element sides in the expression of Box I and

<span id="eq-48"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\Delta {\boldsymbol B}_\phi ^e = -{l_3^e\over 2A^e}{\boldsymbol T}_3 [{\boldsymbol 0},{\boldsymbol  I}_2,{\boldsymbol I}_2]   </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (48)
|}

The bending strain field of Eq.([[#eq-27|27]]) is now modified as

<span id="eq-49"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{\boldsymbol \varepsilon }_b^e = {\boldsymbol B}_w \bar {\boldsymbol w}^e + ({\boldsymbol B}_\phi +\Delta {\boldsymbol  B}_\phi )\bar{\boldsymbol \phi }^e={\boldsymbol B}_w \bar {\boldsymbol w}^e + \bar{\boldsymbol B}_\phi  \bar{\boldsymbol \phi }^e </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (49)
|}

where

<span id="eq-50"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\bar{\boldsymbol B}_\phi = {\boldsymbol B}_\phi +\Delta {\boldsymbol B}_\phi  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (50)
|}

The new matrix <math display="inline">\bar{\boldsymbol B}_\phi </math> substitutes matrix <math display="inline">{\boldsymbol  B}_\phi </math> in the expression for <math display="inline">{\boldsymbol K}^e_\phi </math> of Eq.([[#eq-35|35]]).

We recall that the above modifications are only needed for elements with a clamped edge.

The  process is repeated twice if the element has two clamped edges.

===4.2 <span id='lb-4.2'></span>Simply supported edge (w=0, θₛ=0)===

The condition <math display="inline">w=0</math> at the nodes laying on a simply supported (SS) edge is prescribed when solving the global system of equations as for the clamped case. Prescribing <math display="inline">w=0</math> at an edge node also automatically implies that <math display="inline">{\partial  w\over  \partial s}</math> is zero along the SS edge direction <math display="inline">s</math> and this is the option taken for the standard rotation-free thin BPT element. Note that this is equivalent to assuming a “soft” simply support condition in general plate theory (i.e. <math display="inline">w=M_n=M_{ns}=0</math> at the SS edge) <span id='citeF-29'></span><span id='citeF-30'></span>[[#cite-29|[29,30]]].

The “hard” support condition requires prescribing the tangential rotation <math display="inline">\theta _s</math> to a zero value. This means

<span id="eq-51"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\theta _s = \frac{\partial w}{\partial s}+\phi _s =0 \quad \rightarrow \quad \phi _s = - \frac{\partial w}{\partial s} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (51)
|}

Prescribing <math display="inline">w</math> at the SS nodes gives automatically <math display="inline">\left(\frac{\partial  w}{\partial s} \right)_i=0</math>  and hence the nodal variable <math display="inline">\bar \phi _{s_i}</math> must be also prescribed to a zero value at these nodes. This condition can be imposed by transforming the cartesian shear strains to the boundary axes <math display="inline">s,n</math> and making <math display="inline">\bar \phi _{s_i}=0</math> at each simply supported boundary node.

Indeed, prescribing just <math display="inline">w_i =0</math> at the support nodes and letting <math display="inline">\bar \phi _{s_i}</math> free reproduces the soft support condition for Reissner-Mindlin theory. This is the approach chosen in the examples presented in the paper ( Figures [[#img-5|5]], [[#img-6|6]], [[#img-9|9]], [[#img-10|10]]).

It is interesting that the BPT element does not suffer from the difficulties associated to prescribing the tangential rotation at SS nodes which occur in standard thin plate elements <span id='citeF-29'></span><span id='citeF-32'></span><span id='citeF-33'></span>[[#cite-29|[29,32,33]]]. This is due to the fact that the rotations do not appear explicitly as variables in the BPT formulation and the SS condition is directly imposed by prescribing the deflection at the support nodes.

===4.3 <span id='lb-4.3'></span>Symmetry edge (θₙ=0)===

The condition of zero normal rotation (<math display="inline">\theta _n=0</math>) is imposed by neglecting the contributions from the normal rotation <math display="inline">\frac{\partial w}{\partial n}</math> at the symmetry edge when computing Eq.([[#eq-30|30]]). The condition <math display="inline">\theta _n=\frac{\partial w}{\partial n}+\phi _n= 0</math> for the thick case introduces additional terms in the stiffness equations. The procedure is identical as explained for the clamped edge.

'''Remark 2''' . The condition <math display="inline">\bar{\phi }_{n_i}=0</math> can also be imposed at free edge nodes. This however has not shown to be improve the quality of the solution or to lead to any computational advantage, other that the reduction in the number of DOFs, in those cases.

==5 TWO-NODED ROTATION-FREE BEAM ELEMENT==

A two-noded rotation-free beam element with shear deformation effects can be simply derived as a particular case of the formulation for the BPT element previously described. The resulting beam element is termed CCB+ as a reference to the ''cell-centered'' approach used to compute the curvature at the element center.

Figure [[#img-3|3]] shows the patch of three beam elements needed for computing the constant curvature at the central element using a cell-centered finite-volume type scheme. The starting point in the formulation of the CCB+ element is the standard Timoshenko beam theory <span id='citeF-29'></span><span id='citeF-30'></span>[[#cite-29|[29,30]]]. The relevant expressions are:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\begin{array}{l}\boldsymbol {Deflection}:\, w \\[.2cm] \boldsymbol {Rotation}:\, \theta = \displaystyle \frac{\partial w}{\partial x}+\phi \\[.2cm] \boldsymbol {Displacement~vector}:\, {\boldsymbol u} =[w,\theta ]^T \\[.2cm] \boldsymbol {Curvature}:\, {\boldsymbol \varepsilon }_b = [{\kappa }] = \left[\displaystyle \frac{\partial \theta }{\partial x}\right]= {\boldsymbol L}_b \theta \quad ,\quad {\boldsymbol L}_b =\left[\displaystyle \frac{\partial }{\partial x}\right]\\[.2cm] \boldsymbol {Transverse~shear~deformation}:\, {\boldsymbol \sigma }_s = [\gamma ] = \left[\displaystyle \frac{\partial w}{\partial x}-\theta \right]= {\boldsymbol L}_s {\boldsymbol u} \quad ,\quad {\boldsymbol L}_s =\left[\displaystyle \frac{\partial }{\partial  x},-1\right]\qquad \qquad \qquad \\ \end{array} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (52)
|}

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\begin{array}{lll}{\boldsymbol {Constitutive~equations}:\, \\ {\boldsymbol \sigma }_b=[M] = {\boldsymbol D}_b{\boldsymbol \varepsilon }_b & , & {\boldsymbol D}_b =[EI] \\ {\boldsymbol \sigma }_s =[Q] = {\boldsymbol D}_s{\boldsymbol \varepsilon }_s & , & {\boldsymbol D}_s =[\alpha GA]  \\ \end{array} &</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (53)
|}

where <math display="inline">I</math> and <math display="inline">A</math> are the inertia modulus and the area of the transverse cross section of the beam and the rest of the terms have been defined previously.

The Hu-Washizu functional for the beam has the same form as in Eq.([[#eq-10|10]]) with <math display="inline">{\boldsymbol m} = [m]</math>.

Similarly, the governing equations have identical expressions as for the plate problem Eqs.([[#eq-11|11]]&#8211;[[#Equilibrium equations|2.2]]).

A linear interpolation is  chosen for the deflection <math display="inline">w</math> and the shear angle <math display="inline">\phi </math> as

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>w =\sum \limits _{i=1}^2 N_i\bar w_i \quad ,\quad \phi = \sum \limits _{i=1}^2 N_i \bar \phi _i </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (54)
|}

<div id='img-3'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_940614781-Fig3.png|420px|(a) Patch of three rotation-free CCB+ elements. (b) Patch   of element adjacents  to a clamped boundary]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 3:''' (a) Patch of three rotation-free CCB+ elements. (b) Patch   of element adjacents  to a clamped boundary
|}

where <math display="inline">N_i = 1/2 (1+\xi \xi _i)</math> with <math display="inline">\xi _1=-1</math> and <math display="inline">\xi _2=1</math>  are the standard linear shape functions for the 2-noded Lagrange element <span id='citeF-28'></span><span id='citeF-30'></span>[[#cite-28|[28,30]]].

Following an identical process as for the 3-noded BPT+ element, the same system of Eqs.([[#eq-34|34]]) is found. The expressions for the stiffness matrices <math display="inline">{\boldsymbol  K}_w^e</math>, <math display="inline">{\boldsymbol K}_\phi ^e</math> and <math display="inline">{\boldsymbol K}_{w_\phi }^e</math> coincide with those  given in Eq.([[#eq-35|35]]) simply by noting that the integrals change from area to line  ones and substituting the area <math display="inline">A^e</math> by the element length <math display="inline">l^e</math>. The form of  matrices <math display="inline">{\boldsymbol  B}_w</math> and <math display="inline">{\boldsymbol B}_\phi </math> is given in Box II.

The expressions for the equivalent force vector are similar to Eqs.([[#eq-36|36]])&#8211;([[#eq-38|38]]), substituting the area integrals by integrals along the element length. For a uniform distribution of <math display="inline">q</math> and <math display="inline">m</math>

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\begin{array}{l}f_{w_i}= \displaystyle \frac{ql^e}{2}+\frac{2m}{l^e}\xi _i     \\[.2cm]     f_{\phi _i}= \displaystyle \frac{A^e}{2}m\\   \end{array}   </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (55)
|}


{|  class="floating_tableSCP" style="text-align: left; margin: 1em auto;border-top: 2px solid;border-bottom: 2px solid;min-width:50%;"
|-
| <math>  \begin{array}{l} \displaystyle  \mathbf{B}_w = \frac{1}{2 l^{a} l^{e} l^{b}} \left[ l^{b},-l^{b},-l^{a},l^{a} \right]\\[.3cm] \displaystyle  \mathbf{B}_\phi = \frac{1}{l^{e}}[1,-1] \end{array}</math>'''Box II'''. Matrices <math>{\boldsymbol B}^e_{w}</math> and <math>{\boldsymbol B}^e_{\phi }</math> for the 2-noded rotation-free CCB+ element
|}
where <math display="inline">\xi _i</math> is the value of the natural coordinate <math display="inline">\xi </math> at node <math display="inline">i</math>.

The procedure for imposing the boundary conditions follows precisely the lines described for the BPT+ element in Section 4. The <math display="inline">\Delta{\boldsymbol B}_\phi </math> matrix for an element with a clamped node is (Eq.(48) and Figure 3b)

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\hbox{Left clamped node: }\Delta{\boldsymbol B}_\phi = \frac{1}{l^{e}} [1,0]  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (56)
|}

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\hbox{Right clamped node: } \Delta{\boldsymbol B}_\phi = \frac{1}{l^{e}} [0,-1] </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (57)
|}

A similar expression is obtained for an element with a node on the symmetry axis.

The CCB+ element is an extension of the 2-noded rotation-free CCB element based on Euler-Bernouilli beam theory. A description of the CCB element can be found in <span id='citeF-11'></span><span id='citeF-30'></span>[[#cite-11|[11,30]]].

==6 EXAMPLES==

===6.1 Square and circular plates===

The efficiency and accuracy of the BPT+ element has been tested in the analysis of a number of square plates of side <math display="inline">L</math> and circular plates of diameter <math display="inline">2L</math>  under a uniformly distributed loading and a central point load for different thicknesses ranging from to <math display="inline">t/L =10^{-3}</math> (very thin plate) to <math display="inline">t/L= 0.1</math> (thick plate). Simply supported and clamped boundary conditions have been considered. ''For rectangular plates the soft SS condition has been assumed''. For circular plates <math display="inline">\phi _s</math> is automatically zero at the SS boundary for the problems studied due to symmetry and hence results correspond to the hard SS case. Figure 4 shows some of the different meshes  used for the analysis.

<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_940614781-placas_fig1.png|219px|Some of meshes of BPT+ elements used for analysis of square and circular plates]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 4:''' Some of meshes of BPT+ elements used for analysis of square and circular plates
|}

<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_940614781-placas_fig2.png|364px|Simple supported square plate (soft support) under uniform load. Convergence of  central   deflection for different thicknesses. Upper curves show convergence of the   vertical deflection and the shear angles for a thick plate with the number   of iterations. Lower curves show the distribution of Mₓ and Qₓ along the   central line and their contours ]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 5:''' Simple supported square plate (soft support) under uniform load. Convergence of  central   deflection for different thicknesses. Upper curves show convergence of the   vertical deflection and the shear angles for a thick plate with the number   of iterations. Lower curves show the distribution of <math>M_x</math> and <math>Q_x</math> along the   central line and their contours 
|}

<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_940614781-placas_fig3.png|364px|Simple supported square plate (soft support) under central point load. Convergence of  central   deflection for different thicknesses. Upper curves show convergence of the   vertical deflection and the shear angles for a thick plate with the number   of iterations. Lower curves show the distribution of Mₓ and Qₓ along the   central line and their contours]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 6:''' Simple supported square plate (soft support) under central point load. Convergence of  central   deflection for different thicknesses. Upper curves show convergence of the   vertical deflection and the shear angles for a thick plate with the number   of iterations. Lower curves show the distribution of <math>M_x</math> and <math>Q_x</math> along the   central line and their contours
|}

<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_940614781-placas_fig4.png|371px|Clampled square plate under uniform load. Convergence of  central   deflection for different thicknesses. Upper curves show convergence of the   vertical deflection and the shear angles for a thick plate with the number   of iterations. Lower curves show the distribution of Mₓ and Qₓ along the   central line and their contours]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 7:''' Clampled square plate under uniform load. Convergence of  central   deflection for different thicknesses. Upper curves show convergence of the   vertical deflection and the shear angles for a thick plate with the number   of iterations. Lower curves show the distribution of <math>M_x</math> and <math>Q_x</math> along the   central line and their contours
|}

<div id='img-8'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_940614781-placas_fig5.png|365px|Clampled square plate under central point load. Convergence of  central   deflection for different thicknesses. Upper curves show convergence of the   vertical deflection and the shear angles for a thick plate with the number   of iterations. Lower curves show the distribution of Mₓ and Qₓ along the   central line and their contours ]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 8:''' Clampled square plate under central point load. Convergence of  central   deflection for different thicknesses. Upper curves show convergence of the   vertical deflection and the shear angles for a thick plate with the number   of iterations. Lower curves show the distribution of <math>M_x</math> and <math>Q_x</math> along the   central line and their contours 
|}

<div id='img-9'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_940614781-placas_fig6.png|371px|Simple supported circular plate under uniform load. Convergence of  central   deflection for different thicknesses. Upper curves show convergence of the   vertical deflection and the shear angles for a thick plate with the number   of iterations. Lower curves show the distribution of Mₓ and Qₓ along the   central line and their contours ]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 9:''' Simple supported circular plate under uniform load. Convergence of  central   deflection for different thicknesses. Upper curves show convergence of the   vertical deflection and the shear angles for a thick plate with the number   of iterations. Lower curves show the distribution of <math>M_x</math> and <math>Q_x</math> along the   central line and their contours 
|}

<div id='img-10'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_940614781-placas_fig7.png|364px|Simple supported circular plate under central point load. Convergence of  central   deflection for different thicknesses. Upper curves show convergence of the   vertical deflection and the shear angles for a thick plate with the number   of iterations. Lower curves show the distribution of Mₓ and Qₓ along the   central line and their contours ]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 10:''' Simple supported circular plate under central point load. Convergence of  central   deflection for different thicknesses. Upper curves show convergence of the   vertical deflection and the shear angles for a thick plate with the number   of iterations. Lower curves show the distribution of <math>M_x</math> and <math>Q_x</math> along the   central line and their contours 
|}

<div id='img-11'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_940614781-placas_fig8.png|363px|Clampled circular plate under uniform load. Convergence of  central   deflection for different thicknesses. Upper curves show convergence of the   vertical deflection and the shear angles for a thick plate with the number   of iterations. Lower curves show the distribution of Mₓ and Qₓ along the   central line and their contours ]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 11:''' Clampled circular plate under uniform load. Convergence of  central   deflection for different thicknesses. Upper curves show convergence of the   vertical deflection and the shear angles for a thick plate with the number   of iterations. Lower curves show the distribution of <math>M_x</math> and <math>Q_x</math> along the   central line and their contours 
|}

<div id='img-12'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_940614781-placas_fig9.png|373px|Clampled circular plate under central point load. Convergence of  central   deflection for different thicknesses. Upper curves show convergence of the   vertical deflection and the shear angles for a thick plate with the number   of iterations. Lower curves show the distribution of Mₓ and Qₓ along the   central line and their contours ]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 12:''' Clampled circular plate under central point load. Convergence of  central   deflection for different thicknesses. Upper curves show convergence of the   vertical deflection and the shear angles for a thick plate with the number   of iterations. Lower curves show the distribution of <math>M_x</math> and <math>Q_x</math> along the   central line and their contours 
|}

Results for the  cases studied using the enhanced iterative scheme of Section 3.4 are presented in Figures 5&#8211;12. Each figure shows:

* The convergence of the vertical deflection values and the shear     angles with the number of iterations measured as

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>

L_2^w = \left[\sum \limits _{j=1}^N \frac{(\bar w_j^i - \bar w_j^{i-1})^2}{     (\bar w_j^i)^2}\right]^{1/2}     </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (58)
|}

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>

L_2^{\boldsymbol \phi } = \left[\sum \limits _{j=1}^N     \frac{[\bar {\boldsymbol \phi }_j^i - \bar {\boldsymbol \phi }_j^{i-1}]^T(\bar {\boldsymbol \phi }_j^i -     \bar {\boldsymbol \phi }_j^{i-1})}{     [{\boldsymbol \phi }_j^i]^T {\boldsymbol \phi }_j^i}\right]^{1/2}     </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (59)
|}

* where <math display="inline">N</math> is the number of nodes in the mesh and an upper index denotes     the iteration number. A value of <math display="inline">\bar w_j^0=0</math> and <math display="inline">\bar {\boldsymbol \phi }_j^0=\mathbf{0}</math> has been taken. For the examples     considered, the iterative scheme stops when <math display="inline">L_2^w< 10^{-3}</math>.

* The value of the normalized central deflection     (<math display="inline">w_c = \frac{D\bar w_c\cdot 10^5}{qL^4}</math> for uniform load and <math display="inline">w_c = \frac{D\bar w_c\cdot 10^4}{PL^2}</math> for point load with <math display="inline">D=\frac{Et^3}{12(1-\nu ^2)}</math>)     for the range of     thickness ratios <math display="inline">10^{-3}-10^{-1}</math> for each of the  meshes considered. Results are     compared with analytical and series values for the thin and thick cases     <span id='citeF-26'></span><span id='citeF-29'></span><span id='citeF-30'></span><span id='citeF-31'></span>[[#cite-26|[26,29,30,31]]] when available, or, alternatively,  with     FEM results obtained for the deflection at the center of the midle plane using a mesh of <math display="inline">40\times 40 \times 6</math> eight-noded hexahedra in a     quarter of plate.<p>     For the point load case the analytical value for the deflection     under the load     given by thick plate theory is infinity. Hence, results for the deflection     for thick plates are compared     at the mid-point along a central line in this case.

</p>
* The distribution of the bending moment <math display="inline">M_x</math> and the shear force <math display="inline">Q_x</math> along the central line for the thick case (<math display="inline">t/L=0.10</math>) for each of the     five meshes considered. The isovalues of <math display="inline">M_x</math> and <math display="inline">Q_x</math> over a     quarter of the plate are also shown for the finer mesh.

The following conclusions are drawn from the examples:

<ol>

<li>The BPT+ element reproduces accurately the expected results for the     deflection field for thin and thick plates.  </li>

<li>A converged solution of <math display="inline">L_2^w<10^{-2}</math> for  the deflection      field using the iterative algorithm of Section 3.4 was obtained     in a maximum of four iterations for the thick  case for all the meshes     considered. The number of iterations reduced to two for     thin plates, due to the less importance of shear effects in these cases.  </li>

<li>The convergence of the shear angles is slightly slower than for the deflection field. However an error norm of <math display="inline">L_2^\phi < 10^{-1}</math> obtained in 3-4 iterations was found to give accurate results for the shear strains and the shear forces distribution for the thick case.  </li>

<li>The distribution of the bending moments and the shear forces (obtained     directly from the shear angles by Eq.(23)) was good and in accordance with     the expected results for the thick case.  </li>

<li>  For thin plates the distribution of bending moments is also very good. However the     distribution of the shear forces deteriorates slightly if computed      via Eq.(23). This is due to the fact that the shear angles  tend to zero as the plate is thinner. It     is more appropriate in theses cases to compute the shear forces from the     bending moment distribution as in standard thin plate     theory <span id='citeF-29'></span><span id='citeF-30'></span>[[#cite-29|[29,30]]].     </li>
<li>Similar good results were obtained for SS square plate problems     solved with the “hard” SS condition obtained by prescribing <math display="inline">\bar \phi _{s_i} =0</math>     at the support nodes (see Table [[#table-1|1]]). </li>

</ol>

Indeed for problems with a discontinuous shear force field (such as the case of internal point forces), the nodal continuity of the shear strains which is intrinsic to the element formulation is a drawback to accurately capturing shear force jumps. This can be overcome by computing the shear forces at the center of the elements adjacent to the point load and extrapolating the solution  within each element so as to reproduce the shear force jump at the common node.


{|  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 Simply supported square thick plate (hard support) under uniform load. Normalized central deflection values for <math>t/L =0.10</math> for the meshes of Figure [[#img-5|5]]
|- style="border-top: 2px solid;"
| colspan='7' | SS (hard) square thick plate. Uniform load, <math>t/L =0.10</math>

|- style="border-top: 2px solid;"
| style="text-align: left;" | 
| Mesh 1
| Mesh 2 
| Mesh 3 
| Mesh 4 
| Mesh 5 
| Mesh 6 

|- style="border-top: 2px solid;"
| style="text-align: left;" | <math>w_c</math>
| -381,73 
| -404,59 
| -416,31
| -421,98 
| -424,65 
| -425,59 
|-
| style="text-align: left;" | <math>w_c/w_c^a</math>
| 0,893
| 0,947 
| 0,974 
| 0,988 
| 0,994
| 0,996   

|- style="border-top: 2px solid;"
| colspan='7' style="text-align: left;" | <math>w_c^a</math> (Series solution): -427.28 <span id='citeF-29'></span><span id='citeF-30'></span>[[#cite-29|[29,30]]]
|-
| colspan='7' style="text-align: left;" | Convergence achieved in a maximum of 4 iterations for each mesh
|}

===6.2 Simple supported and cantilever beams===

The accuracy of the CCB+ beam element of Section 5 was tested in the analysis of simple supported and cantilever thick beams under distributed and point loads. Results of the study plotted in Figures [[#img-13|13]]&#8211;[[#img-16|16]] show the normalized distribution of the  deflection along the beams and the normalized value of the central deflection (for the clamped beam) and the end deflection (for the cantilever beam) with the number of elements. Good results (error less than 10%) are obtained with a relatively coarse mesh (8 elements).

The distribution of the bending moment and the shear force along the beam is also plotted for the 40 element mesh. Results are practically coincident with the analytical values. Similar good behaviour was obtained for other thick and thin beam problems studied with the CCB+ element.

==7 CONCLUDING REMARKS==

A methodology for extending the rotation-free plate and beam elements initially designed for thin/slender situations so as to account with shear deformation effects has been presented. The method allows one to introduce the effect of shear deformation in a progressive (iterative) manner starting from the initial thin solution.

The formulation is useful for analysis of plates and beams of a variety of thickness and materials. A particular interesting application of the new elements is the analysis of composite laminated plates and beams for which shear deformation effects are relevant.

The formulation  is also suitable for implementing an adaptive solution scheme where the shear angles are introduced in structures (or zones of a structure) where the effect of shear deformation is relevant.

<div id='img-13'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_940614781-vigas_fig1.png|600px|Simple supported thick beam under uniform load. Convergence of central   deflection and distribution of the deflection   for different meshes of CCB+ elements. Bending moment and shear   force diagrams for 40 element mesh]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 13:''' Simple supported thick beam under uniform load. Convergence of central   deflection and distribution of the deflection   for different meshes of CCB+ elements. Bending moment and shear   force diagrams for 40 element mesh
|}

<div id='img-14'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_940614781-vigas_fig2.png|600px|Simple supported thick beam under central point load. Convergence of central deflection  and distribution of the deflection for different meshes of CCB+ elements. Bending moment and shear   force diagrams for  40 element mesh]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 14:''' Simple supported thick beam under central point load. Convergence of central deflection  and distribution of the deflection for different meshes of CCB+ elements. Bending moment and shear   force diagrams for  40 element mesh
|}

<div id='img-15'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_940614781-vigas_fig3.png|600px|Cantilever thick beam under uniform load. Convergence of end deflection and   distribution of the deflection   for different meshes of CCB+ elements. Bending moment and shear   force diagrams for 40 element mesh]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 15:''' Cantilever thick beam under uniform load. Convergence of end deflection and   distribution of the deflection   for different meshes of CCB+ elements. Bending moment and shear   force diagrams for 40 element mesh
|}

<div id='img-16'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_940614781-vigas_fig4.png|600px|Cantilever thick beam under central point load. Convergence of end deflection and distribution of deflection  for different meshes of CCB+ elements. Bending moment and shear   force diagrams for  40 element mesh]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 16:''' Cantilever thick beam under central point load. Convergence of end deflection and distribution of deflection  for different meshes of CCB+ elements. Bending moment and shear   force diagrams for  40 element mesh
|}

==ACKNOWLEDGEMENTS==

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

===References===

<div id="cite-1"></div>
'''[[#citeF-1|[1]]]'''  Nay RA, Utku S. An alternative to the finite element method. '' Variational Methods Eng.'' 1972; Vol. 1:3/63&#8211;3/74.

<div id="cite-2"></div>
'''[2]'''  Chan ASL, Davies GAO Davies. A simplified finite element model for the impact of thin shell. In Bulson PS (ed.), ''Structures under Shock and Impact'', ''Proc. 1st. Int. Conf.'', Cambridge, MA, USA, 365&#8211;380, 1989.

<div id="cite-3"></div>
'''[3]''' Hampshire JK,  Topping BHV,  Chan HC. Three node triangular elements with one degree of freedom per node. ''Engng. Comput.'' 1992; '''9''':49&#8211;62.

<div id="cite-4"></div>
'''[4]''' Phaal R, Calladine CR. A simple class of finite elements for plate and shell problems. I: Elements for beams and thin plates. ''Int. J. Num. Meth. Engng.'' 1992; '''35''':955&#8211;977.

<div id="cite-5"></div>
'''[5]''' Phaal R, Calladine CR. A simple class of finite elements for plate and shell problems. II: An element for thin shells with only translational degrees of freedom. ''Int. J. Num. Meth. Engng.'' 1992;  '''35''':979&#8211;996.

<div id="cite-6"></div>
'''[[#citeF-6|[6]]]''' Oñate E, Cervera M. Derivation of thin plate bending elements with one degree of freedom per node. ''Engineering Computations'' 1993; '''10''':553&#8211;561.

<div id="cite-7"></div>
'''[7]''' Rio G, Tathi B, Laurent H. A new efficient finite element model of shell with only three degrees of freedom per node. Applications to industrial deep drawing test. in ''Recent Developments in Sheet Metal Forming Technology'', MJM Barata Marques (ed), 18th IDDRG Biennial Congress, Lisbon, 1994.

<div id="cite-8"></div>
'''[8]'''  Sabourini F, Bounet M. Analysis of plates and shells with a simplified three noded triangular element. '' Thin Walled Structures'' 1995; '''21''':209&#8211;223.

<div id="cite-9"></div>
'''[[#citeF-9|[9]]]'''  Rojek J, Oñate E. Sheet springback analysis using a simple shell triangle with translational degrees of freedom only. ''Int. J. of Forming Processes'' 1998; '''1'''(3):275&#8211;296.

<div id="cite-10"></div>
'''[[#citeF-10|[10]]]'''  Rojek J, Oñate E, Postek E. Application of explicit FE codes to simulation of sheet and bulk forming processes. ''J. of Materials Processing Technology'' 1998; '''80-81''':620&#8211;627.

<div id="cite-11"></div>
'''[[#citeF-11|[11]]]'''  Jovicevic J, Oñate E. ''Analysis of beams and shells using a rotation-free finite element-finite volume formulation'', Monograph 43, CIMNE, Barcelona, 1999.

<div id="cite-12"></div>
'''[[#citeF-12|[12]]]''' Oñate E, Zárate F. Rotation-free plate and shell triangles. ''Int. J. Num. Meth. Engng.'' 2000; '''47''':557&#8211;603.

<div id="cite-13"></div>
'''[13]''' Cirak F, Ortiz M. Subdivision surfaces: A new paradigm for thin-shell finite element analysis. ''Int. J. Num. Meths in Engng'' 2000; '''47''':2039&#8211;2072.

<div id="cite-14"></div>
'''[14]''' Cirak F, Ortiz M.  Fully <math display="inline">C^{1}</math>-conforming subdivision elements for finite deformations thin-shell analysis. ''Int. J. Num. Meths in Engng'' 2001; '''51''':813&#8211;833.

<div id="cite-15"></div>
'''[[#citeF-15|[15]]]''' Flores FG, Oñate E. A basic thin shell triangle with only translational DOFs for large strain plasticity. '' Int. J. Num. Meths in Engng'' 2001; '''51''':57&#8211;83.

<div id="cite-16"></div>
'''[16]''' Guo YQ, Gatis W, Naceur H, Batoz JL. An efficient dkt rotation free shell element for spring simulations in sheet metal forming. ''Computers & Structures'' 2002; '''80''':2299&#8211;2312.

<div id="cite-17"></div>
'''[17]''' Engel G, Garikipati K,  Hughes TJR,  Larson MG, Mazzei L, Taylor RL. Continuous/discontinuous finite element approximation of fourth-order elliptic problems in structural and continuum mechanics with applications to thin beams and plates, and strain gradient elasticity. ''Comput. Methods Appl. Mech. Engrg.'' 2002; '''191''':3669&#8211;3750.

<div id="cite-18"></div>
'''[18]'''  Oñate E, Cendoya P,   Miquel J. Non linear explicit dynamic analysis of shells using the BST rotation-free triangle. ''Engineering Computations'' 2002; '''19'''(6):662&#8211;706.

<div id="cite-19"></div>
'''[[#citeF-19|[19]]]'''  Oñate E, Flores FG. Advances in the formulation of the rotation-free shell triangle. ''Comput. Meth. in Appl. Mech. Engng.'' 2005; '''194''':2406&#8211;2443.

<div id="cite-20"></div>
'''[[#citeF-20|[20]]]''' Flores FG, Oñate E. Improvements in the membrane behaviour of the three node rotation-free BST shell triangle using an assumed strain approach. ''Computer Methods in Applied Mechanics and Engineering'' 2005; '''194''':907&#8211;932.

<div id="cite-21"></div>
'''[21]'''  Brunet M, Sabourin F. Analysis of a rotation-free 4-node shell element. ''Int. J. Num. Meth. Engng.'' 2006; '''66''':1483&#8211;1510.

<div id="cite-22"></div>
'''[22]'''  Sabourini F, Brunet M. Detailed formulation of the rotation-free triangular element “S3” for general purpose shell analysis. ''Engrg. Comput.'' 2006; '''23''':469&#8211;502.

<div id="cite-23"></div>
'''[23]''' Gärdsback M, Tibert G. A comparison of rotation-free triangular shell elements for unstructured meshes. ''Comput. Meth. Appl. Mech. Engng.'' 2007; ''' 196''':5001&#8211;5015.

<div id="cite-24"></div>
'''[[#citeF-24|[24]]]'''  Oñate E, Flores FG, Neamtu L. Enhanced rotation-free basic shell triangle. Applications to sheet metal forming. E. Oñate and R. Owen (eds.), ''Computational Plasticity'', 239-265, 2007 Springer, The Netherlands.

<div id="cite-25"></div>
'''[[#citeF-25|[25]]]'''  Bushnell D, Almroth BO. Finite difference energy method for non linear shell analysis. ''J. Computers and Structures'' 1971; '''1''':361.

<div id="cite-26"></div>
'''[[#citeF-26|[26]]]'''  Timoshenko SP. ''Theory of Plates and Shells'', McGraw Hill, New York, 1971.

<div id="cite-27"></div>
'''[[#citeF-27|[27]]]'''  Ugural AC. ''Stresses in  Plates and Shells'', McGraw Hill, New York, 1981.

<div id="cite-28"></div>
'''[[#citeF-28|[28]]]'''  Zienkiewicz OC, Taylor RL, Zhu JZ. ''The Finite Element Method. Its basis and fundamentals''. 6th Edition, Elsevier, 2005.

<div id="cite-29"></div>
'''[[#citeF-29|[29]]]'''  Zienkiewicz OC, Taylor RL. ''The Finite Element Method for Solid and Structural Mechanics''. 6th Edition, Elsevier, 2005.

<div id="cite-30"></div>
'''[[#citeF-30|[30]]]'''  Oñate E, ''Structural Analysis with the Finite Element Method. Linear Statics. Vol.2: Beams, Plates and Shells''. CIMNE-Springer, 2009.

<div id="cite-31"></div>
'''[[#citeF-31|[31]]]'''  Batoz, JL, Dhatt G. ''Modelisation des structures par élements finis''.  ''Vol. 2: Poutres et Plaques''. Hermes, Paris, 1990.

<div id="cite-32"></div>
'''[[#citeF-32|[32]]]'''  Hugues TJR. ''The finite element method. Linear static and dynamic finite element analysis''. (p. 324), Prentice-Hall, 1987.

<div id="cite-33"></div>
'''[[#citeF-33|[33]]]'''  Babuŝka I, Scapolla T. Benchmark computation and performance evaluation for a rhombic plate bending problem. ''International Journal for Numerical Methods in Engineering'' 1989; '''28''':155&#8211;180.

<div id="cite-34"></div>
'''[34]'''  Salerno VL, Golberg MA. Effect of shear deformation on the bending of rectangular plates. ''Journal Applied Mechanics'' 1960; '''27''': 54&#8211;58.

<div id="cite-35"></div>
'''[35]'''  Shrinivas S, Rao AK. Flexure of thick rectangular plates. ''Journal Applied Mechanics'' 1973; 298&#8211;299.