==Improvements in the membrane behaviour of the three node rotation-free BST shell triangle using an assumed strain approach==

​

'''Fernando G. Flores

​

National University of Córdoba 

​

Casilla de Correo 916 

​

5000 Córdoba - Argentina 

​

and 

​

Eugenio Oñate 

​

International Centre for Numerical Methods in Engineering 

​

Edificio C1, Gran Capitán s/n 

​

08034 Barcelona - Spain'''

-->

​

==Abstract==

​

In this paper an assumed strain approach is presented in order to improve the  membrane behaviour of a thin shell triangular element. The so called Basic Shell Triangle (BST) has three  nodes with only translational degrees of freedom and is based on a Total  Lagrangian Formulation. As in the original BST element the curvatures are computed  resorting to the surrounding elements (patch of four elements). Membrane  strains are now also computed from the same patch of elements which leads to a  non-conforming  membrane behaviour. Despite this non-conformity the element  passes the patch test. Large strain plasticity is considered using a logarithmic  strain-stress pair. A plane stress behaviour with an additive decomposition of  elastic and plastic strains is assumed. A hyperplastic law is considered for  the elastic part while for the plastic part an anisotropic quadratic (Hill)  yield function with non-linear isotropic hardening is adopted. The element,  termed CBST, has been implemented  in an explicit (hydro-)code adequate to simulate sheet-stamping processes and  in an implicit static/dynamic code. Several examples are given showing the good performance of the enhnaced rotation-free shell triangle.

​

16.0cm   22.0cm  0cm  -1cm

​

==1 Introduction==

​

Development of efficient and robust shell finite elements for modelling large-scale  industrial applications has been a field of constant research in the last  thirty years. The continuining efforts of many scientists has focused in the direction of deriving simple elements (preferably triangles) applicable to large scale computation, typical of practical shell problems. Some recent successful shell triangles are reported in [1,2]. For a comprehensive review of shell elements see [3]. From the theoretical point of view Kirchhoff hypothesis are  enough for the simulation of the essential aspects of most shells problems,  including structures in civil, mechanical and naval engineering, sheet metal  forming, crash-worthiness analysis and others. Due the well-known difficulties  associated to <math display="inline">C^{1}</math> continuity, “thick shell” approaches that need <math display="inline">C^{0}</math>  continuity only have been extensively used, although they are computationally  more expensive due to the number of degrees of freedom necessary to obtain thin   solutions with similar precision. The development of shell elements without  transversal shear strains then collides with the need for C<math display="inline">^{1}</math> continuity, which is not  simple to satisfy in general conditions and has led to the development of  non-conforming approaches. Among the many options of this kind, the use of shell elements with only  displacements as degrees of freedom (rotation-free) is very attractive from  the computational point of view.

​

Few shell elements exists in the literature with only translations as degrees of  freedom. They are generally based on non-conforming approaches. A state of the  art can be found in Ref.<span id='citeF-4'></span>[[#cite-4|[4]]]. Recently Cirak y  Ortiz<sup><span id='citeF-5'></span>[[#cite-5|[5]]]</sup> have developed a conforming thin shell element,  but it departs in some aspects from the standard finite element formulation. The  use of rotation-free shell elements in large scale simulations of sheet metal forming processes  is growing <sup><span id='citeF-6'></span>[[#cite-6|[6]]]-<span id='citeF-9'></span>[[#cite-9|[9]]]</sup>, specially in  explicit integration codes.

​

In Ref.<span id='citeF-10'></span>[[#cite-10|[10]]] we presented a thin shell triangular element with only  displacements as degrees of freedom, based on a total lagrangian formulation.  The element is an extension of the rotation-free Basic Shell Triangle (BST) originally developed by Oñate y  Zárate<sup><span id='citeF-4'></span>[[#cite-4|[4]]]</sup> using an updated lagrangian formulation  The element has three nodes and a linear approximation of the displacement field. The membrane part is identical to the standard  constant strain triangle. The bending part is also constant and is computed  from an integration over the element boundary using information from the  deflection gradients of the adjacent elements. This leads to a non-conforming element for  the bending part, but it converges to the right solution and is robust. The BST element  has  a very good performance in bending for structured meshes, but it is less accurate  for non-structured grids. On the other hand, for membrane dominated problems,  as is the case for most sheet metal forming simulations, it requires  fine  discretizations associated to the constant strain triangle behaviour.

​

In this paper an extension of the rotation-free BST element is presented.  The displacement gradient terms needed for the membrane and bending strains are  computed from the nodal displacements of a patch of element using a new procedure. This allows to improve the membrane  behaviour and obtain a smoother curvature field. The outline of this  paper is as follows. In sections 2 and 3 we introduce the kinematics of the  shell and the constitutive models used in the numerical simulations. Sections  4 to 6 describe the finite element approximation, including the new  proposal for the gradient evaluation from a patch of elements, the computation of  the metric tensor and the curvatures and the derivation of the necessary  element expressions for the computer implementation. Sections 7 and 8 sumarize  the numerical experiments performed in the linear and non-linear range. Numerical results   demonstrate the excellent performance of the enhanced three node rotation-free shell triangle.

​

==2 Shell kinematics==

​

A  summary of the most relevant hypothesis related to the kinematic  behaviour of the shell are presented. Further details may be found in the wide  literature dedicated to this field <span id='citeF-3'></span>[[#cite-3|[3]]].

​

Consider a shell with undeformed middle surface occupying the domain <math display="inline">\Omega ^{0}</math> in <math display="inline">R^{3}</math> with a boundary <math display="inline">\Gamma ^{0}</math>. At each point of the middle  surface a thickness <math display="inline">h^{0}</math> is defined. The positions <math display="inline">\mathbf{x}^{0}</math> and <math display="inline">\mathbf{x}</math> of any point in the undeformed and the deformed configurations  can be respectively written as a function of the position of the middle  surface <math display="inline">\mathbf{\varphi }</math> and the normal <math display="inline">\mathbf{t}_{3}</math> at the point as

​

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

|-

| 

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

|-

| style="text-align: center;" | <math>\mathbf{x}^{0}\left( \xi _{1},\xi _{2},\zeta \right)    =\mathbf{\varphi }^{0}\left( \xi _{1},\xi _{2}\right) +\lambda \mathbf{t}_{3}^{0}</math>

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

|-

| style="text-align: center;" | <math> \mathbf{x}\left( \xi _{1},\xi _{2},\zeta \right)    =\mathbf{\varphi }\left( \xi _{1},\xi _{2}\right) +\zeta \lambda \mathbf{t}_{3} </math>

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

|}

|}

​

where <math display="inline">\xi _{1},\xi _{2}</math> are curvilinear local coordinates defined over the  middle surface of the shell, and <math display="inline">\zeta </math> is the distance in the undeformed  configuration of the point to the middle surface. The product <math display="inline">\zeta \lambda </math>  is the distance of the point to the middle surface measured on the deformed  configuration. This implies a constant strain in the normal direction  associated to the parameter <math display="inline">\lambda </math> relating the thickness at the present and initial   configurations, i.e.

​

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

|-

| 

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

|-

| style="text-align: center;" | <math>  \lambda =\frac{h}{h^{0}} </math>

|}

|}

​

A convective system is defined at each point as

​

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

|-

| 

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

|-

| style="text-align: center;" | <math>\mathbf{g}_{i}\left( \mathbf{\xi }\right) =\frac{\partial \mathbf{x}}{\partial \xi _{i}} </math>

|}

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

|}

​

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

|-

| 

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

|-

| style="text-align: center;" | <math>\mathbf{g}_{\alpha }\left( \mathbf{\xi }\right)    =\frac{\partial \left( \mathbf{\varphi }\left( \xi _{1},\xi _{2}\right) +\zeta \lambda \mathbf{t}_{3}\right) }{\partial \xi _{\alpha }}=\mathbf{\varphi }_{^{\prime }\alpha } +\zeta \left( \lambda \mathbf{t}_{3}\right) _{^{\prime }\alpha }</math>

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

|-

| style="text-align: center;" | <math> \mathbf{g}_{3}\left( \mathbf{\xi }\right)    =\frac{\partial \left( \mathbf{\varphi }\left( \xi _{1},\xi _{2}\right) +\zeta \lambda \mathbf{t}_{3}\right) }{\partial \zeta }=\lambda \mathbf{t}_{3} </math>

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

|}

|}

​

This can be particularized for the points on the middle surface as

​

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

|-

| 

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

|-

| style="text-align: center;" | <math>\mathbf{a}_{\alpha }    =\mathbf{g}_{\alpha }\left( \zeta=0\right) =\mathbf{\varphi }_{^{\prime }\alpha }</math>

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

|-

| style="text-align: center;" | <math> \mathbf{a}_{3}    =\mathbf{g}_{3}\left( \zeta=0\right) =\lambda  \mathbf{t}_{3} </math>

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

|}

|}

​

This allows to introduce the covariant (first fundamental form) and  contravariant metric tensors of the middle surface as

​

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

|-

| 

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

|-

| style="text-align: center;" | <math>a_{\alpha \beta }=\mathbf{a}_{\alpha }\cdot \mathbf{a}_{\beta } </math>

|}

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

|}

​

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

|-

| 

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

|-

| style="text-align: center;" | <math>a^{\alpha \beta }=\mathbf{a}^{\alpha }\cdot \mathbf{a}^{\beta }=\mathbf{\tilde  {\varphi }}_{^{\prime }\alpha }\cdot \mathbf{\tilde{\varphi }}_{^{\prime }\beta } </math>

|}

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

|}

​

and the curvatures (second fundamental form) of the middle surface as

​

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

|-

| 

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

|-

| style="text-align: center;" | <math>\kappa _{\alpha \beta }=\frac{1}{2}\left( \mathbf{\varphi }_{^{\prime }\alpha  }\cdot \mathbf{t}_{3^{\prime }\beta }+\mathbf{\varphi }_{^{\prime }\beta } \cdot \mathbf{t}_{3^{\prime }\alpha }\right) </math>

|}

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

|}

​

The deformation gradient can be written as

​

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

|-

| 

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

|-

| style="text-align: center;" | <math>\mathbf{F=}\left[ \begin{array}{c}[c]{ccc} \mathbf{\varphi }_{^{\prime }1}+\zeta \left( \lambda \mathbf{t}_{3}\right) _{^{\prime }1} & \mathbf{\varphi }_{^{\prime }2}+\zeta \left( \lambda  \mathbf{t}_{3}\right) _{^{\prime }2} & \lambda \mathbf{t}_{3} \end{array}  \right] </math>

|}

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

|}

​

The product <math display="inline">\mathbf{F}^{T}\mathbf{F=U}^{2}=\mathbf{C}</math> (where <math display="inline">\mathbf{U}</math> is the right stretch tensor, and <math display="inline">\mathbf{C}</math> the right Cauchy-Green  deformation tensor) is

​

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

|-

| 

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

|-

| style="text-align: center;" | <math>\mathbf{U}^{2}    =\left[ \begin{array}{c}[c]{c} \mathbf{\varphi }_{^{\prime }1}^{T}+\zeta \left( \lambda \mathbf{t}_{3}\right) _{^{\prime }1}^{T}\\

 \mathbf{\varphi }_{^{\prime }2}^{T}+\zeta \left( \lambda \mathbf{t}_{3}\right) _{^{\prime }2}^{T}\\

 \lambda \mathbf{t}_{3}^{T} \end{array}  \right] \left[ \begin{array}{c}[c]{ccc} \mathbf{\varphi }_{^{\prime }1}+\zeta \left( \lambda \mathbf{t}_{3}\right) _{^{\prime }1} & \mathbf{\varphi }_{^{\prime }2}+\zeta \left( \lambda  \mathbf{t}_{3}\right) _{^{\prime }2} & \lambda \mathbf{t}_{3} \end{array}  \right]</math>

|-

| style="text-align: center;" | <math>   =\left[ \begin{array}{c}[c]{ccc} \mathbf{\varphi }_{^{\prime }1}\cdot \mathbf{\varphi }_{^{\prime }1} &  \mathbf{\varphi }_{^{\prime }1}\cdot \mathbf{\varphi }_{^{\prime }2} & 0\\

 \mathbf{\varphi }_{^{\prime }1}\cdot \mathbf{\varphi }_{^{\prime }2} &  \mathbf{\varphi }_{^{\prime }2}\cdot \mathbf{\varphi }_{^{\prime }2} & 0\\

 0 & 0 & \lambda ^{2} \end{array}  \right]</math>

|-

| style="text-align: center;" | <math>   +\zeta \lambda \left[ \begin{array}{c}[c]{ccc} 2\mathbf{\varphi }_{^{\prime }1}\cdot \mathbf{t}_{3^{\prime }1} & \mathbf{\varphi  }_{^{\prime }1}\cdot \mathbf{t}_{3^{\prime }2}+\mathbf{\varphi }_{^{\prime }2} \cdot \mathbf{t}_{3^{\prime }1} & 0 \\

 \mathbf{\varphi }_{^{\prime }1}\cdot \mathbf{t}_{3^{\prime }2}+\mathbf{\varphi  }_{^{\prime }2}\cdot \mathbf{t}_{3^{\prime }1} & \mathbf{\varphi }_{^{\prime } 2}\cdot \mathbf{t}_{3^{\prime }2} & 0\\

 0 & 0 & 0  \end{array}  \right]</math>

|-

| style="text-align: center;" | <math>   +\zeta ^{2}\lambda ^{2}\left[ \begin{array}{c}[c]{ccc} \mathbf{t}_{3^{\prime }1}\cdot \mathbf{t}_{3^{\prime }1} & \mathbf{t}_{3^{\prime  }1}\cdot \mathbf{t}_{3^{\prime }2} & 0\\

 \mathbf{t}_{3^{\prime }1}\cdot \mathbf{t}_{3^{\prime }2} & \mathbf{t}_{3^{\prime  }2}\cdot \mathbf{t}_{3^{\prime }2} & 0\\

 0 & 0 & 0  \end{array}  \right] </math>

|}

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

|}

​

where the derivatives of the thickness ratio <math display="inline">\lambda _{^{\prime }a}</math> have been  neglected. Neglecting also the term associated to <math display="inline">\zeta ^{2}</math> and introducing  the definition of the covariant metric tensor <math display="inline">a_{\alpha \beta }</math> at the middle  surface and the curvatures <math display="inline">k_{\alpha \beta }</math> gives

​

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

|-

| 

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

|-

| style="text-align: center;" | <math>\mathbf{U}^{2}=\left[ \begin{array}{c}[c]{ccc} a_{11}+2k_{11}\zeta \lambda & a_{12}+2k_{12}\zeta \lambda & 0\\

 a_{12}+2k_{12}\zeta \lambda & a_{22}+2k_{22}\zeta \lambda & 0\\

 0 & 0 & \lambda ^{2} \end{array}  \right] </math>

|}

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

|}

​

Above expression shows that <math display="inline">\mathbf{U}^{2}</math>   is non-zero at the original  configuration for curved surfaces (<math display="inline">\kappa _{ij}^{0}\neq{0}</math>). The changes of curvature of the middle surface are computed by

​

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

|-

| 

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

|-

| style="text-align: center;" | <math>\chi _{ij}=\kappa _{ij}-\kappa _{ij}^{0} </math>

|}

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

|}

​

For computational convenience the following approximate expression (which is  exact for initially flat surfaces) will be adopted

​

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

|-

| 

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

|-

| style="text-align: center;" | <math>\mathbf{U}^{2}=\left[ \begin{array}{c}[c]{ccc} a_{11}+2\chi _{11}\zeta \lambda & a_{12}+2\chi _{12}\zeta \lambda & 0\\

 a_{12}+2\chi _{12}\zeta \lambda & a_{22}+2\chi _{22}\zeta \lambda & 0\\

 0 & 0 & \lambda ^{2} \end{array}  \right] </math>

|}

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

|}

​

Above expression of <math display="inline">\mathbf{U}^{2}</math> is useful  to compute different Lagrangian strain  measures. An advantage of these measures is that they are associated to  material fibres, what makes it easy to take into account  material anisotropy. It is also useful to compute  the eigen decomposition of <math display="inline">\mathbf{U}</math> as

​

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

|-

| 

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

|-

| style="text-align: center;" | <math>\mathbf{U=}\sum _{\alpha=1}^{3}\lambda _{\alpha } \mathbf{r}_{\alpha } \otimes \mathbf{r}_{\alpha } </math>

|}

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

|}

​

where <math display="inline">\lambda _{\alpha }</math> and <math display="inline">\mathbf{r}_{\alpha }</math> are the eigenvalues and  eigenvectors of the right stretch tensor <math display="inline">\mathbf{U}</math>.

​

In order to treat plasticity at finite strains an adequate stress-strain pair must  be used. The Hencky measures will be adopted here. The (logarithmic)  strains are defined as

​

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

|-

| 

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

|-

| style="text-align: center;" | <math>\mathbf{E}_{\ln }\mathbf{=}\left[ \begin{array}{c}[c]{ccc} \varepsilon _{11} & \varepsilon _{21} & 0\\

 \varepsilon _{12} & \varepsilon _{22} & 0\\

 0 & 0 & \varepsilon _{33} \end{array}  \right] =\sum _{\alpha=1}^{3}\ln \left( \lambda _{\alpha }\right) \mathbf{r}_{\alpha }\otimes \mathbf{r}_{\alpha } </math>

|}

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

|}

​

Consistently, the Hencky stress tensor <math display="inline">\mathbf{T}</math>  will be  used as the stress measure. Using a total lagrangian formulation it is convenient to work with  the the second Piola-Kirchhoff stresses (<math display="inline">\mathbf{S}</math>) for the evaluation of  the residual forces. We define the rotated tensors as

​

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

|-

| 

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

|-

| style="text-align: center;" | <math>\mathbf{T}_{L}    =\mathbf{R}_{L}^{T}\;\mathbf{T\;R}_{L}</math>

|-

| style="text-align: center;" | <math> \mathbf{S}_{L}    =\mathbf{R}_{L}^{T}\;\mathbf{S\;R}_{L} </math>

|}

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

|}

​

where <math display="inline">\mathbf{R}_{L}</math> is the rotation tensor obtained from the eigenvectors  of <math display="inline">\mathbf{U}</math> given by

​

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

|-

| 

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

|-

| style="text-align: center;" | <math>\mathbf{R}_{L}=\left[ \begin{array}{c}[c]{ccc} \mathbf{r}_{1} & \mathbf{r}_{2} & \mathbf{r}_{3} \end{array}  \right] </math>

|}

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

|}

​

The relationship between the Hencky and Piola-Kirchhoff stresses is

​

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

|-

| 

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

|-

| style="text-align: center;" | <math>\left[ S_{L}\right] _{\alpha \alpha }    =\frac{1}{\lambda _{\alpha }^{2} }\left[ T_{L}\right] _{\alpha \alpha }</math>

|-

| style="text-align: center;" | <math> \left[ S_{L}\right] _{\alpha \beta }    =\frac{\ln \left( \lambda _{\alpha  }/\lambda _{\beta }\right) }{\frac{1}{2}\left( \lambda _{\alpha }^{2} -\lambda _{\beta }^{2}\right) }\left[ T_{L}\right] _{\alpha \beta } </math>

|}

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

|}

​

The second Piola-Kirchhoff stress tensor can be finally computed by

​

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

|-

| 

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

|-

| style="text-align: center;" | <math>\mathbf{S}=\mathbf{R}_{L}\;\mathbf{S}_{L}\mathbf{\;R}_{L}^{T} </math>

|}

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

|}

​

The generalized stresses (forces and moments) are obtained by integrating across the  original thickness the second Piola-Kirchhoff stress tensor using the actual  distance to the middle surface for the evaluation of the bending moments, i.e.

​

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

|-

| 

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

|-

| style="text-align: center;" | <math>\mathbf{N}=\int _{h^{0}}\mathbf{S} d\zeta </math>

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

|-

| style="text-align: center;" | <math> \mathbf{M}=\int _{h^{0}}\mathbf{S }\lambda \zeta  d\zeta  </math>

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

|}

|}

​

With these values the weak form of the equilibrium equations can be written as

​

<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>\delta \Pi =\int _{\Omega ^{0}}\left[ \delta \mathbf{E}:\mathbf{N} +\delta {\boldsymbol K}:\mathbf{M}\right] d\Omega ^{0}+\delta \Pi _{ext}=0   </math>

|}

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

|}

​

where <math display="inline">{\boldsymbol K}</math> is the virtual curvature vector and <math display="inline">\mathbf{E}</math> is the virtual  Green-Lagrange strain vector of the middle surface, with

​

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

|-

| 

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

|-

| style="text-align: center;" | <math>E_{ij}={1\over 2} (a_{ij}-\delta _{ij})  </math>

|}

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

|}

​

==3 Constitutive models==

​

In this section, a brief description  of the constitutive models used  in the numerical examples presented in Sections 7 and 8 is given. Two types of  materials are described: an elastic-plastic material associated to thin rolled  metal sheets and a hyper-elastic material for rubbers.

​

In the case of metals, where the elastic strains are small, the use of a  logarithmic strain measure reasonably allows to adopt an additive  decomposition of elastic and plastic components as

​

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

|-

| 

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

|-

| style="text-align: center;" | <math>\mathbf{E}_{\ln }\mathbf{=E}_{\ln }^{e}+\mathbf{E}_{\ln }^{p} </math>

|}

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

|}

​

A constant linear relationship between (plane) stresses and elastic strains is  also adopted giving

​

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

|-

| 

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

|-

| style="text-align: center;" | <math>\mathbf{T}=\mathbf{CE}_{\ln }^{e} </math>

|}

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

|}

​

These constitutive equations are integrated using a standard return algorithm.  The following Mises-Hill<sup><span id='citeF-11'></span>[[#cite-11|[11]]]</sup> yield function with non-linear isotropic  hardening is adopted

​

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

|-

| 

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

|-

| style="text-align: center;" | <math>\left( G+H\right) \;T_{11}^{2}+\left( F+H\right) \;T_{22}^{2} -2H\;T_{11}T_{22}+2N\;T_{12}^{2}=\kappa \left( \varepsilon _{0}+e^{p}\right) ^{n} </math>

|}

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

|}

​

where <math display="inline">F, G, H</math> and <math display="inline">N</math> define the non-isotropic shape of the yield  surface and the parameters <math display="inline">\kappa ,</math> <math display="inline">\varepsilon _{0}</math> and <math display="inline">n</math> define its size  as a function of the effective plastic strain <math display="inline">e^{p}</math>.

​

The Mises-Hill yield function has the advantage of simplicity and it allows, as a  first approximation, to treat rolled thin metal sheets with  planar and transversal anisotropy.

​

For the case of rubbers, the Ogden<sup><span id='citeF-12'></span>[[#cite-12|[12]]]</sup> model extended to  the compressible range is considered. The material behaviour is characterized  by the strain energy density per unit undeformed volume of the form

​

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

|-

| 

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

|-

| style="text-align: center;" | <math>\psi =\frac{K}{2}\left( \ln J\right) ^{2}+\sum _{p=1}^{N}\frac{\mu _{p}}{\alpha _{p}}\left[ J^{-\dfrac{\alpha _{p}}{3}}\left( \sum _{i=1}^{3} \lambda _{i}^{\alpha _{p}-1}\right) -3\right] </math>

|}

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

|}

​

where <math display="inline">K</math> is the bulk modulus of the material, <math display="inline">J</math> is the determinant of <math display="inline">\mathbf{U}</math>, <math display="inline">N</math>, <math display="inline">\mu _{i}</math> and <math display="inline">\alpha _{i}</math> are material parameters, <math display="inline">\mu _{i}\,,\,\alpha _{i}</math> are real numbers such that <math display="inline">\mu _{i}\alpha _{i}>0</math> <math display="inline"> (\forall i=1,N)</math> and <math display="inline">N</math> is a positive integer.

​

The stress measures associated to the principal logarithmic strains are  denoted by <math display="inline">\beta _{i}</math>. They can be computed noting that

​

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

|-

| 

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

|-

| style="text-align: center;" | <math>\beta _{i}=\frac{\partial \psi \left( \Lambda \right) }{\partial \left( \ln \lambda _{i}\right) }=K\left( \ln J\right) +\lambda _{i}\sum _{p=1}^{N} \mu _{p}J^{-\dfrac{\alpha _{p}}{3}}\left( \lambda _{i}^{\alpha _{p}-1}-\frac{1}{3}\frac{1}{\lambda _{i}}\sum _{j=1}^{3}\lambda _{j}^{\alpha _{p}}\right) </math>

|}

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

|}

​

we define now

​

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

|-

| 

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

|-

| style="text-align: center;" | <math>a_{p}=\sum _{j=1}^{3}\lambda _{j}^{\alpha _{p}} </math>

|}

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

|}

​

which gives

​

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

|-

| 

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

|-

| style="text-align: center;" | <math>\beta _{i}=K\left( \ln J\right) +\sum _{p=1}^{N}\mu _{p}J^{-\dfrac{\alpha _{p} }{3}}\left( \lambda _{i}^{\alpha _{p}}-\frac{1}{3}a_{p}\right) </math>

|}

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

|}

​

The values of <math display="inline">\beta _i</math>, expressed in the principal strains directions, allow to evaluate  the stresses  in the convective coordinate system as

​

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

|-

| 

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

|-

| style="text-align: center;" | <math>\mathbf{T}=\sum _{i=1}^{3}\beta _{i}\;\mathbf{r}_{i}\otimes \mathbf{r}_{i} </math>

|}

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

|}

​

We note that the Hencky stress tensor <math display="inline">\mathbf{T}</math> can be easily  particularized for the plane stress case.

​

==4 Interpolation functions and gradient evaluation==

​

The starting point of the enhanced rotation-free of BST element is  to discretize the shell surface with a standard 3-node  triangular mesh. The difference with a standard finite element method is that for the  computation of strains within an element, the configuration of the three adjacent triangular  elements is  used. Then at each triangle a patch of four elements, formed by the central triangle and the three adjacent ones is considered (see Figure [[#img-1|1]].a). This allows to define a  quadratic interpolation of the geometry from the position  of the 6 nodes. In the isoparametric space, we keep the vertices of the  3-node master element (standard linear triangle) with the positions

​

* node 1: <math display="inline">\left( \xi ,\eta \right) =\left( 0,0\right) </math>

​

* node 2: <math display="inline">\left( \xi ,\eta \right) =\left( 1,0\right) </math>

​

* node 3: <math display="inline">\left( \xi ,\eta \right) =\left( 0,1\right) </math>

​

and the three additional nodes that complete the patch with the positions

​

* node 4: <math display="inline">\left( \xi ,\eta \right) =\left( 1,1\right) </math>

​

* node 5: <math display="inline">\left( \xi ,\eta \right) =\left( -1,1\right) </math>

​

* node 6: <math display="inline">\left( \xi ,\eta \right) =\left( 1,-1\right) </math>

​

<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_340530075-test-fig1.png|600px|Patch of four elements  (a)in spacial coordinates (b)in natural coordinates]]

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

| colspan="1" | '''Figure 1:''' Patch of four elements  (a)in spacial coordinates (b)in natural coordinates

|}

​

It is possible to define now a set of (non-standard) quadratic shape functions over  over the six node element as:

​

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

|-

| 

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

|-

| style="text-align: center;" | <math>\begin{array}{c}[c]{ccc} N^{1}=\zeta{+\xi}\eta &  & N^{4}=\frac{\zeta }{2}\left( \zeta{-1}\right)\\

 N^{2}=\xi{+\eta}\zeta &  & N^{5}=\frac{\xi }{2}\left( \xi{-1}\right)\\

 N^{3}=\eta{+\zeta}\xi &  & N^{6}=\frac{\eta }{2}\left( \eta{-1}\right) \end{array}  </math>

|}

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

|}

​

This interpolation allows to compute the displacement gradients at selected  points in order to use an assumed strain approach. The computation of the  gradients is performed at three points over the boundary of the central  element of the patch. These points are located at the mid-point of each side,  denoted by G1, G2 and G3 in Figure [[#img-1|1]].b. This election has the  following characteristics.

​

* Gradients at these points depend only on the nodes pertaining to  the two elements adjacent to the side. This can be easily verified by sampling the derivatives of  the shape functions  at each mid-side point.

​

* When gradients are computed at the common mid-side point from two  adjacent elements, the same values are  obtained, as the coordinates of  the same four points are used. That is, the gradients at the mid-side points are  independent of the element where they are computed. A side-oriented  implementation of the finite element could lead to a unique evaluation of the  gradients per side.

​

Next, some details of the implementation chosen are presented. We denote by <math display="inline">\mathbf{t}_{1}</math> and <math display="inline">\mathbf{t}_{2}</math> the orthogonal unit vectors over the  tangent plane (undeformed configuration) associated to a conveniently selected  local cartesian system. The cartesian derivatives of the shape functions are computed at the original configuration by the standard expression

​

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

|-

| 

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

|-

| style="text-align: center;" | <math>\left[ \begin{array}{c}[c]{c} N_{^{\prime }1}^{I}\\

 N_{^{\prime }2}^{I} \end{array}  \right] =\mathbf{J}^{-1}\left[ \begin{array}{c}[c]{c} N_{^{\prime }\xi }^{I}\\

 N_{^{\prime }\eta }^{I} \end{array}  \right] </math>

|}

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

|}

​

where the Jacobian matrix at the original configuration is

​

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

|-

| 

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

|-

| style="text-align: center;" | <math>\mathbf{J=}\left[ \begin{array}{c}[c]{cc} \mathbf{\varphi }_{^{\prime }\xi }^{\left( 0\right) }\cdot \mathbf{t}_{1} &  \mathbf{\varphi }_{^{\prime }\eta }^{\left( 0\right) }\cdot \mathbf{t}_{1}\\

 \mathbf{\varphi }_{^{\prime }\xi }^{\left( 0\right) }\cdot \mathbf{t}_{2} &  \mathbf{\varphi }_{^{\prime }\eta }^{\left( 0\right) }\cdot \mathbf{t}_{2} \end{array}  \right] </math>

|}

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

|}

​

With the previous definitions the deformation gradient on the middle surface,  associated to an arbitrary spatial cartesian system and to the material  cartesian system defined on the middle surface, is

​

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

|-

| 

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

|-

| style="text-align: center;" | <math>\left[ \mathbf{\varphi }_{^{\prime }1},\mathbf{\varphi }_{^{\prime }2}\right] =\left[ \mathbf{\varphi }_{^{\prime }\xi },\mathbf{\varphi }_{^{\prime }\eta  }\right]  \mathbf{J}^{-1} </math>

|}

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

|}

​

The covariant metric tensor can be computed from above expression as

​

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

|-

| 

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

|-

| style="text-align: center;" | <math>\mathbf{g=}\left[ \begin{array}{c}[c]{cc} a_{11} & a_{12}\\

 a_{21} & a_{22} \end{array}  \right] =\left[ \begin{array}{c}[c]{cc} \mathbf{\varphi }_{^{\prime }1}\cdot \mathbf{\varphi }_{^{\prime }1} &  \mathbf{\varphi }_{^{\prime }1}\cdot \mathbf{\varphi }_{^{\prime }2}\\

 \mathbf{\varphi }_{^{\prime }2}\cdot \mathbf{\varphi }_{^{\prime }1} &  \mathbf{\varphi }_{^{\prime }2}\cdot \mathbf{\varphi }_{^{\prime }2} \end{array}  \right] </math>

|}

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

|}

​

The component of <math display="inline">\boldsymbol g</math> can be used to compute any convenient membrane strain measures. For example the  Green-Lagrange strain tensor is obtained by

​

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

|-

| 

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

|-

| style="text-align: center;" | <math>\mathbf{E}_{GL}\mathbf{=}\frac{1}{2}\left[ \begin{array}{c}[c]{cc} \mathbf{\varphi }_{^{\prime }1}\cdot \mathbf{\varphi }_{^{\prime }1}-1 &  \mathbf{\varphi }_{^{\prime }1}\cdot \mathbf{\varphi }_{^{\prime }2}\\

 \mathbf{\varphi }_{^{\prime }2}\cdot \mathbf{\varphi }_{^{\prime }1} &  \mathbf{\varphi }_{^{\prime }2}\cdot \mathbf{\varphi }_{^{\prime }2}-1  \end{array}  \right] =\frac{1}{2}\left[ \begin{array}{c}[c]{cc} a_{11}-1 & a_{12}\\

 a_{21} & a_{22}-1  \end{array}  \right] </math>

|}

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

|}

​

In the usual FEM matrix notation

​

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

|-

| 

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

|-

| style="text-align: center;" | <math>\left[ \begin{array}{c}[c]{c} E_{11}\\

 E_{22}\\

 2E_{12} \end{array}  \right] =\frac{1}{2}\left[ \begin{array}{c}[c]{c} \mathbf{\varphi }_{^{\prime }1}\cdot \mathbf{\varphi }_{^{\prime }1}-1\\

 \mathbf{\varphi }_{^{\prime }2}\cdot \mathbf{\varphi }_{^{\prime }2}-1\\

 2\mathbf{\varphi }_{^{\prime }1}\cdot \mathbf{\varphi }_{^{\prime }2} \end{array}  \right] </math>

|}

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

|}

​

The virtual strains are obtained by the variation of above expression as

​

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

|-

| 

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

|-

| style="text-align: center;" | <math>\delta \left[ \begin{array}{c}[c]{c} E_{11}\\

 E_{22}\\

 2E_{12} \end{array}  \right] =\left[ \begin{array}{c}[c]{c} \mathbf{\varphi }_{^{\prime }1}\cdot \delta \mathbf{\varphi }_{^{\prime }1}\\

 \mathbf{\varphi }_{^{\prime }2}\cdot \delta \mathbf{\varphi }_{^{\prime }2}\\

 \mathbf{\varphi }_{^{\prime }1}\cdot \delta \mathbf{\varphi }_{^{\prime }2} +\delta \mathbf{\varphi }_{^{\prime }1}\cdot \mathbf{\varphi }_{^{\prime }2} \end{array}  \right] </math>

|}

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

|}

​

In sides of elements where the adjacent element does not exist  (i.e. at the  boundaries of the shell), the gradient on that side is made equal to the  gradient computed over the 3 nodes of the central element of the patch. From  the definition of the gradients at the three mid-side points of the triangle,  the element formulation can be interpreted as an ''assumed strain'' approach, where  the metric tensor is interpolated over the element using the values computed  on the sides by

​

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

|-

| 

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

|-

| style="text-align: center;" | <math>\mathbf{g}\left( \xi ,\eta \right) =\left( 1-2\zeta \right) \mathbf{g}^{1}  +\left( 1-2\xi \right) \mathbf{g}^{2} +\left( 1-2\eta \right) \mathbf{g}^{3} </math>

|}

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

|}

​

==5 Computation of curvatures==

​

Curvatures are assumed to be constant over the element. They are computed by performing  an average over the central element as

​

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

|-

| 

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

|-

| style="text-align: center;" | <math>\kappa _{\alpha \beta }\mathbf{ =}\frac{-1}{A^{M}}\int _{A^{M}}\mathbf{t}_{3}\cdot \mathbf{\varphi }_{^{\prime }\beta \alpha }\;dA^{M} </math>

|}

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

|}

​

where (<math display="inline">A^{M}</math> is the original area of the element).

​

Integrating by parts Eq.(43), the following integral over the element boundary <math display="inline">\Gamma ^{M}</math> is  obtained:

​

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

|-

| 

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

|-

| style="text-align: center;" | <math>\kappa _{\alpha \beta }\mathbf{ =}\frac{-1}{A^{M}}{\displaystyle \oint  _{\Gamma ^{M}}}\mathbf{t}_{3}\cdot \mathbf{\varphi }_{^{\prime }\alpha }\;n_{\beta  }\;d\Gamma ^{M} </math>

|}

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

|}

​

The three distinct components of the curvature tensor are

​

<span id="eq-44"></span>

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

|-

| 

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

|-

| style="text-align: center;" | <math>\left[ \begin{array}{c}[c]{c} \kappa _{11}\\

 \kappa _{22}\\

 2\kappa _{12} \end{array}  \right] =\frac{-1}{A^{M}}{\displaystyle \oint _{\Gamma ^{M}}}\left[ \begin{array}{c}[c]{cc} n_{1} & 0\\

 0 & n_{2}\\

 n_{2} & n_{1} \end{array}  \right] \left[ \begin{array}{c}[c]{c} \mathbf{t}_{3}\cdot \mathbf{\varphi }_{^{\prime }1}\\

 \mathbf{t}_{3}\cdot \mathbf{\varphi }_{^{\prime }2} \end{array}  \right] d\Gamma ^{M}  </math>

|}

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

|}

​

The numerical evaluation of the boundary integral in Eq.([[#eq-44|44]]) results in a sum over the integration  points on the element boundary which are, in fact, the same points used for the computation of  the gradients in the membrane formulations, i.e.

​

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

|-

| 

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

|-

| style="text-align: center;" | <math>\left[ \begin{array}{c}[c]{c} \kappa _{11}\\

 \kappa _{22}\\

 2\kappa _{12} \end{array}  \right] =\frac{-1}{A^{M}}\sum _{G=1}^{3}l_{G}\left[ \begin{array}{c}[c]{rr} n_{1} & 0\\

 0 & n_{2}\\

 n_{2} & n_{1} \end{array}  \right] ^{G}\left[ \begin{array}{c}[c]{c} \mathbf{\varphi }_{1}^{G}\cdot \mathbf{t}_{3}\\

 \mathbf{\varphi }_{2}^{G}\cdot \mathbf{t}_{3} \end{array}  \right] </math>

|}

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

|}

​

For the definition of <math display="inline">\mathbf{t}_{3}</math>, the linear interpolation over the  central element is used. In this case the tangent plane components are (with <math display="inline">a_{i}</math>  and <math display="inline">b_{i}</math> being the usual cartesian projections of the sides when area coordinates  are used)

​

<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>\left[ \begin{array}{c}[c]{c} \mathbf{\varphi }_{1}\\

 \mathbf{\varphi }_{2} \end{array}  \right] ^{M}    =\frac{1}{2A^{M}}\sum _{i=1}^{3}\left[ \begin{array}{c}[c]{c} -b_{i}\\

 a_{i} \end{array}  \right] \mathbf{\varphi }^{i}=\left[ \begin{array}{c}[c]{ccc} \bar{N}_{^{\prime }1}^{1} & \bar{N}_{^{\prime }1}^{2} & \bar{N}_{^{\prime }1} ^{3}\\

 \bar{N}_{^{\prime }2}^{1} & \bar{N}_{^{\prime }2}^{2} & \bar{N}_{^{\prime }2}^{3} \end{array}  \right] \left[ \begin{array}{c}[c]{c} \mathbf{\varphi }^{1}\\

 \mathbf{\varphi }^{2}\\

 \mathbf{\varphi }^{3} \end{array}  \right]  </math>

|}

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

|}

​

<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>\mathbf{t}_{3}    =\frac{\mathbf{\varphi }_{1}^{M}\times \mathbf{\varphi }_{2}^{M}}{\left\vert \mathbf{\varphi }_{1}^{M}\times \mathbf{\varphi }_{2} ^{M}\right\vert }=\lambda \;\mathbf{\varphi }_{1}^{M}\times \mathbf{\varphi }_{2}^{M}  </math>

|}

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

|}

​

From these expressions it is also possible to compute in the original  configuration the element area <math display="inline">A^{M}</math>, the outer normals <math display="inline">\left( n_{1} ,n_{2}\right) ^{i}</math> at each side and its lengths <math display="inline">l_{i}</math>. Eq.([[#eq-47|47]]) also allows to evaluate the thickness ratio <math display="inline">\lambda </math> in the deformed configuration and the actual normal <math display="inline">\mathbf{t}_{3}</math>. In Eq.([[#eq-46|46]]), <math display="inline">\bar{N}^{I}</math> are the linear shape functions of the three node linear triangle <span id='citeF-3'></span>[[#cite-3|[3]]].

​

As one integration point is used over each side, it is not necessary to  distinguish  between sides (<math display="inline">I</math>) and integrations points (<math display="inline">G</math>) any more. In  this way the curvatures can be computed by:

​

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

|-

| 

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

|-

| style="text-align: center;" | <math>\left[ \begin{array}{c}[c]{c} \kappa _{11}\\

 \kappa _{22}\\

 2\kappa _{12} \end{array}  \right]    =\frac{-1}{A^{M}}\sum _{I=1}^{3}l_{I}\left[ \begin{array}{c}[c]{rr} n_{1} & 0\\

 0 & n_{2}\\

 n_{2} & n_{1} \end{array}  \right] ^{G}\left[ \begin{array}{c}[c]{c} \mathbf{\varphi }_{1}^{I}\cdot \mathbf{t}_{3}\\

 \mathbf{\varphi }_{2}^{I}\cdot \mathbf{t}_{3} \end{array}  \right]</math>

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

|-

| style="text-align: center;" | <math>   =2\sum _{I=1}^{3}\left[ \begin{array}{c}[c]{rr} \bar{N}_{^{\prime }1}^{I} & 0\\

 0 & \bar{N}_{^{\prime }2}^{I}\\

 \bar{N}_{^{\prime }2}^{I} & \bar{N}_{^{\prime }1}^{I} \end{array}  \right] \left[ \begin{array}{c}[c]{c} \mathbf{\varphi }_{1}^{I}\cdot \mathbf{t}_{3}\\

 \mathbf{\varphi }_{2}^{I}\cdot \mathbf{t}_{3} \end{array}  \right] </math>

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

|}

|}

​

In the original BST element<sup><span id='citeF-4'></span><span id='citeF-10'></span>[[#cite-4|[4,10]]]</sup> the gradient <math display="inline">\mathbf{\varphi  }_{\alpha }^{I}</math> was computed as the average of the linear approximations over  the two adjacent elements. In the new version, the gradient is  computed at each side <math display="inline">I</math> from the quadratic interpolation

​

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

|-

| 

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

|-

| style="text-align: center;" | <math>\left[ \begin{array}{c}[c]{c} \mathbf{\varphi }_{1}\\

 \mathbf{\varphi }_{2} \end{array}  \right] ^{I}=\left[ \begin{array}{c}[c]{cccc} N_{^{\prime }1}^{1} & N_{^{\prime }1}^{2} & N_{^{\prime }1}^{3} & N_{^{\prime } 1}^{I+3}\\

 N_{^{\prime }2}^{1} & N_{^{\prime }2}^{2} & N_{^{\prime }2}^{3} & N_{^{\prime } 2}^{I+3} \end{array}  \right] ^{I}\left[ \begin{array}{c}[c]{c} \mathbf{\varphi }^{1}\\

 \mathbf{\varphi }^{2}\\

 \mathbf{\varphi }^{3}\\

 \mathbf{\varphi }^{I+3} \end{array}  \right] </math>

|}

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

|}

​

Note again than at each side the gradients  depend only on the positions of the three nodes of the central triangle and of an  extra node (<math display="inline">I+3</math>), associated precisely to the side (<math display="inline">I</math>) where the gradient is computed.

​

An alternative form to express the curvatures, which is useful when their  variations are needed, is to define:

​

<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>\mathbf{h}_{ij}=\sum _{I=1}^{3}\left( \bar{N}_{^{\prime }i}^{I} \;\mathbf{\varphi }_{j}^{I}+\bar{N}_{^{\prime }j}^{I}\;\mathbf{\varphi }_{i} ^{I}\right)  </math>

|}

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

|}

​

This gives

​

<span id="eq-52"></span>

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

|-

| 

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

|-

| style="text-align: center;" | <math>\kappa _{ij}=\mathbf{h}_{ij}\cdot \mathbf{t}_{3}  </math>

|}

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

|}

​

This last expression allows to interpret the curvatures as the projections of the  vectors <math display="inline">\mathbf{h}_{ij}</math> over the normal of the central element. Direction  '''t'''<math display="inline">_{3}</math> can be seen as a reference direction.  If a different  direction  is chosen, at an angle <math display="inline">\theta </math> with the former, this has an  influence of order <math display="inline">\theta ^{2}</math> in the projection. This justifies Eq.([[#eq-47|47]]) for the  definition of '''t'''<math display="inline">_{3}</math> as a function exclusively of the three nodes of  the central triangle,  instead of the 6-node  isoparametric interpolation.

​

===5.1 Variation of the curvatures===

​

From Eq.([[#eq-52|52]]) the variation of the curvatures is obtained as

​

<span id="eq-53"></span>

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

|-

| 

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

|-

| style="text-align: center;" | <math>\delta \kappa _{ij}=\delta \mathbf{h}_{ij}\cdot \mathbf{t}_{3}+\mathbf{h}_{ij}\cdot \delta \mathbf{t}_{3}  </math>

|}

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

|}

​

The first term is crucial in the evaluation of the curvatures variation. The variations of <math display="inline">\mathbf{h}_{ij}</math> are:

​

<span id="eq-54"></span>

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

|-

| 

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

|-

| style="text-align: center;" | <math>\delta \mathbf{h}_{ij}=\sum _{I=1}^{3}\left( \bar{N}_{^{\prime }i}^{I} \;\delta \mathbf{\varphi }_{j}^{I}+\bar{N}_{^{\prime }j}^{I}\;\delta  \mathbf{\varphi }_{i}^{I}\right)  </math>

|}

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

|}

​

where

​

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

|-

| 

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

|-

| style="text-align: center;" | <math>\left[ \begin{array}{c}[c]{c} \delta \mathbf{\varphi }_{1}\\

 \delta \mathbf{\varphi }_{2} \end{array}  \right] ^{I}=\left[ \begin{array}{c}[c]{cccc} N_{^{\prime }1}^{1} & N_{^{\prime }1}^{2} & N_{^{\prime }1}^{3} & N_{^{\prime } 1}^{I+3}\\

 N_{^{\prime }2}^{1} & N_{^{\prime }2}^{2} & N_{^{\prime }2}^{3} & N_{^{\prime } 2}^{I+3} \end{array}  \right] ^{I}\left[ \begin{array}{c}[c]{c} \delta \mathbf{u}^{1}\\

 \delta \mathbf{u}^{2}\\

 \delta \mathbf{u}^{3}\\

 \delta \mathbf{u}^{I+3} \end{array}  \right] </math>

|}

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

|}

​

The variations of '''t'''<math display="inline">_{3}</math> can be computed as shown in the  original BST element <sup><span id='citeF-10'></span>[[#cite-10|[10]]]</sup>.

​

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

|-

| 

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

|-

| style="text-align: center;" | <math>\begin{array}{l}\delta t_{31}    =-\mathbf{t}_{3}\cdot \delta \mathbf{\varphi }_{^{\prime }1} ^{M}\\

 \\

 \delta t_{32}    =-\mathbf{t}_{3}\cdot \delta \mathbf{\varphi }_{^{\prime }2}^{M} \end{array}  </math>

|}

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

|}

​

leading to

​

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

|-

| 

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

|-

| style="text-align: center;" | <math>\delta \mathbf{t}_{3}=\delta t_{31}\mathbf{\tilde{\varphi }}_{^{\prime }1}+\delta  t_{32}\mathbf{\tilde{\varphi }}_{^{\prime }2} </math>

|}

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

|}

​

where <math display="inline">\mathbf{\tilde{\varphi }}_{^{\prime }\alpha }</math> are the  contravariant base vectors in the central triangle

​

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

|-

| 

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

|-

| style="text-align: center;" | <math>\begin{array}{l}\mathbf{\tilde{\varphi }}_{^{\prime }1}    =\lambda \;\mathbf{\varphi }_{^{\prime }2}^{M}\times \mathbf{t}_{3}\\

 \\

 \mathbf{\tilde{\varphi }}_{^{\prime }2}    =-\lambda \;\mathbf{\varphi  }_{^{\prime }1}^{M}\times \mathbf{t}_{3} \end{array}  </math>

|}

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

|}

​

We then have

​

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

|-

| 

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

|-

| style="text-align: center;" | <math>\delta \mathbf{t}_{3}    =\left( -\mathbf{t}_{3}\cdot \delta \mathbf{\varphi  }_{^{\prime }1}^{M}\right) \mathbf{\tilde{\varphi }}_{^{\prime }1}+\left( -\mathbf{t}_{3}\cdot \delta \mathbf{\varphi }_{^{\prime }2}^{M}\right) \mathbf{\tilde{\varphi }}_{^{\prime }2}</math>

|-

| style="text-align: center;" | <math>   =-\sum _{J=1}^{3}\left[ \bar{N}_{^{\prime }1}^{J}\mathbf{\tilde{\varphi }}_{^{\prime }1}+\bar{N}_{^{\prime }2}^{J}\mathbf{\tilde{\varphi }}_{^{\prime } 2}\right] \left( \mathbf{t}_{3}\cdot \delta \mathbf{u}^{J}\right) </math>

|}

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

|}

​

Substituting the last expression in Eq.([[#eq-53|53]]), results

​

<span id="eq-60"></span>

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

|-

| 

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

|-

| style="text-align: center;" | <math>\left[ \begin{array}{c}[c]{c} \delta \chi _{11}\\

 \delta \chi _{22}\\

 2\delta \chi _{12} \end{array}  \right]    =2\sum \limits _{I=1}^{3}\left[ \begin{array}{c}[c]{rr} \bar{N}_{^{\prime }1}^{I} & 0\\

 0 & \bar{N}_{^{\prime }2}^{I}\\

 \bar{N}_{^{\prime }2}^{I} & \bar{N}_{^{\prime }1}^{I} \end{array}  \right] \left\{ \sum \limits _{J=1}^{3}\left[ \begin{array}{c}[c]{c} N_{^{\prime }1}^{J\left( I\right) }\left( \mathbf{t}_{3}\cdot \delta  \mathbf{u}^{J}\right)\\

 N_{^{\prime }2}^{J\left( I\right) }\left( \mathbf{t}_{3}\cdot \delta  \mathbf{u}^{J}\right) \end{array}  \right] +\left[ \begin{array}{c}[c]{c} N_{^{\prime }1}^{I+3\left( I\right) }\left( \mathbf{t}_{3}\cdot  \delta \mathbf{u}^{I+3}\right)\\

 N_{^{\prime }2}^{I+3\left( I\right) }\left( \mathbf{t}_{3}\cdot  \delta \mathbf{u}^{I+3}\right) \end{array}  \right] \right\}</math>

|-

| style="text-align: center;" | <math>   -2\sum \limits _{I=1}^{3}\left[ \begin{array}{c}[c]{c} \left( \bar{N}_{^{\prime }1}^{I}\rho _{11}^{1}+\bar{N}_{^{\prime }2}^{I} \rho _{11}^{2}\right)\\

 \left( \bar{N}_{^{\prime }1}^{I}\rho _{22}^{1}+\bar{N}_{^{\prime }2}^{I} \rho _{22}^{2}\right)\\

 \left( \bar{N}_{^{\prime }1}^{I}\rho _{12}^{1}+\bar{N}_{^{\prime }2}^{I} \rho _{12}^{2}\right) \end{array}  \right] \left( \mathbf{t}_{3}\cdot \delta \mathbf{u}^{I}\right)  </math>

|}

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

|}

​

where the projections of the vectors <math display="inline">\mathbf{h}_{ij}</math> over the  contravariant base vectors <math display="inline">\mathbf{\tilde{\varphi }}_{^{\prime }\alpha }</math> have  been included

​

<span id="eq-61"></span>

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

|-

| 

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

|-

| style="text-align: center;" | <math>\rho _{ij}^{\alpha }=\mathbf{h}_{ij}\cdot \mathbf{\tilde{\varphi }}_{^{\prime  }\alpha }  </math>

|}

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

|}

​

These projections together with Eq.([[#eq-52|52]]) allows to write

​

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

|-

| 

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

|-

| style="text-align: center;" | <math>\mathbf{h}_{ij}=\sum _{\alpha=1}^{2}\rho _{ij}^{\alpha }\;\mathbf{\varphi  }_{^{\prime }\alpha }+\kappa _{ij}\mathbf{t}_{3} </math>

|}

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

|}

​

===5.2 Boundary conditions===

​

Elements at the domain boundary, where an adjacent element does not  exist, deserve a special attention. The treatment of essential boundary  conditions associated to translational constraints is straightforward, as they  are the natural degrees of freedom of the element. The conditions associated to the  normal vector are crucial in this formulation. For clamped sides or  symmetry planes, the normal vector <math display="inline">\mathbf{t}_{3}</math> must be kept fixed (clamped case), or  constrained to move in the plane of symmetry (symmetry case). The former case can be seen as a  special case of the latter, so we will consider symmetry planes only. This  restriction can be imposed through the definition of the tangent plane at the  boundary, including the normal to the plane of symmetry <math display="inline">\mathbf{\varphi  }_{^{\prime }n}^{\left( 0\right) }</math> that does not change during the process.

​

<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_340530075-test-fig2.png|300px|Local cartesian system for the treatment of boundary conditions]]

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

| colspan="1" | '''Figure 2:''' Local cartesian system for the treatment of boundary conditions

|}

​

The tangent plane at the boundary is expressed in terms of two orthogonal  unit vectors referred to as local-to-the-boundary cartesian system (see Figure  [[#img-2|2]]) and defined as

​

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

|-

| 

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

|-

| style="text-align: center;" | <math>\left[ \mathbf{\varphi }_{^{\prime }n}^{\left( 0\right) },\;\mathbf{\bar  {\varphi }}_{^{\prime }s}\right] </math>

|}

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

|}

​

where vector <math display="inline">\mathbf{\varphi }_{^{\prime }n}^{\left( 0\right) }</math> is fixed while  direction <math display="inline">\mathbf{\bar{\varphi }}_{^{\prime }s}</math> emerges as the intersection of  the symmetry plane with the plane defined by the central element (<math display="inline">M</math>). The  plane (gradient) defined by the central element is

​

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

|-

| 

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

|-

| style="text-align: center;" | <math>\left[ \mathbf{\varphi }_{^{\prime }1}^{M},\;\mathbf{\varphi }_{^{\prime }2} ^{M}\right] </math>

|}

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

|}

​

the intersection of this plane with the plane of symmetry can be written in  terms of the position of the nodes that define the side (<math display="inline">J</math> and <math display="inline">K</math>) and the  original length of the side <math display="inline">L_{0}</math>, i.e.

​

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

|-

| 

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

|-

| style="text-align: center;" | <math>\mathbf{\varphi }_{^{\prime }s}=\frac{1}{L_{0}}\left( \mathbf{\varphi }^{K}-\mathbf{\varphi }^{J}\right) </math>

|}

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

|}

​

That together with the outer normal to the side <math display="inline">\mathbf{n =}\left[ \nu  _{1},\nu _{2}\right]^T  </math>(resolved in the selected original convective cartesian  system) leads to

​

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

|-

| 

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

|-

| style="text-align: center;" | <math>\left[ \mathbf{\varphi }_{^{\prime }1}^{I},\;\mathbf{\varphi }_{^{\prime }2} ^{I}\right] =\left[ \mathbf{\varphi }_{^{\prime }n},\;\mathbf{\varphi  }_{^{\prime }s}\right] \left[ \begin{array}{c}[c]{cc} \nu _{1} & \nu _{2}\\

 -\nu _{2} & \nu _{1} \end{array}  \right] </math>

|}

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

|}

​

where the normal component of the gradient <math display="inline">\mathbf{\varphi }_{^{\prime }n}</math> is

​

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

|-

| 

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

|-

| style="text-align: center;" | <math>\mathbf{\varphi }_{^{\prime }n}=\frac{\mathbf{\varphi }_{^{\prime }n}^{\left( 0\right) }}{\lambda |\mathbf{\varphi }_{^{\prime }s}|} </math>

|}

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

|}

​

In this way the contribution of the gradient at side <math display="inline">I</math> to vectors <math display="inline">\mathbf{h}_{ij}</math> results in

​

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

|-

| 

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

|-

| style="text-align: center;" | <math>\left[ \begin{array}{c}[c]{c} \mathbf{h}_{11}^{T}\\

 \mathbf{h}_{22}^{T}\\

 2\mathbf{h}_{12}^{T} \end{array}  \right] ^{I}=2\left[ \begin{array}{c}[c]{cc} \bar{N}_{^{\prime }1}^{I} & 0\\

 0 & \bar{N}_{^{\prime }2}^{I}\\

 \bar{N}_{^{\prime }2}^{I} & \bar{N}_{^{\prime }1}^{I} \end{array}  \right] \left[ \begin{array}{c}[c]{cc} \nu _{1} & -\nu _{2}\\

 \nu _{2} & \nu _{1} \end{array}  \right] \left[ \begin{array}{c}[c]{c} \mathbf{\varphi }_{^{\prime }n}^{T}\\

 \mathbf{\varphi }_{^{\prime }s}^{T} \end{array}  \right] </math>

|}

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

|}

​

The only difference with Eq.([[#eq-54|54]]) necessary for the computation  of the curvature variations, is that now the contribution from the gradient at  side <math display="inline">I</math> is

​

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

|-

| 

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

|-

| style="text-align: center;" | <math>\delta \left[ \begin{array}{c}[c]{c} \mathbf{h}_{11}^{T}\\