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

Latest revision as of 16:01, 29 October 2018

Published in Comput. Methods Appl. Mech. Engrg., Vol. 194, pp. 907-932, 2005;
doi: 10.1016/j.cma.2003.08.012

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.

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 ${\textstyle C^{1}}$ continuity, “thick shell” approaches that need ${\textstyle C^{0}}$ 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 ${\textstyle ^{1}}$ 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. [4]. Recently Cirak y Ortiz [5] 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 [6]-[9], specially in explicit integration codes.

In Ref. [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 [4] 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 [3].

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

$\mathbf {x} ^{0}\left(\xi _{1},\xi _{2},\zeta \right)=\mathbf {\varphi } ^{0}\left(\xi _{1},\xi _{2}\right)+\lambda \mathbf {t} _{3}^{0}$	(1)
$\mathbf {x} \left(\xi _{1},\xi _{2},\zeta \right)=\mathbf {\varphi } \left(\xi _{1},\xi _{2}\right)+\zeta \lambda \mathbf {t} _{3}$	(2)

where ${\textstyle \xi _{1},\xi _{2}}$ are curvilinear local coordinates defined over the middle surface of the shell, and ${\textstyle \zeta }$ is the distance in the undeformed configuration of the point to the middle surface. The product ${\textstyle \zeta \lambda }$ 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 ${\textstyle \lambda }$ relating the thickness at the present and initial configurations, i.e.

A convective system is defined at each point as

(3)

$\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 }$	(4)
$\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}$	(5)

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

$\mathbf {a} _{\alpha }=\mathbf {g} _{\alpha }\left(\zeta =0\right)=\mathbf {\varphi } _{^{\prime }\alpha }$	(6)
$\mathbf {a} _{3}=\mathbf {g} _{3}\left(\zeta =0\right)=\lambda \mathbf {t} _{3}$	(7)

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

(8)

(9)

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

(10)

The deformation gradient can be written as

\mathbf {F=} \left[{\begin{array}{cc}\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]

(11)

The product ${\textstyle \mathbf {F} ^{T}\mathbf {F=U} ^{2}=\mathbf {C} }$ (where ${\textstyle \mathbf {U} }$ is the right stretch tensor, and ${\textstyle \mathbf {C} }$ the right Cauchy-Green deformation tensor) is

\mathbf {U} ^{2}=\left[{\begin{array}{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}{cc}\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]

=\left[{\begin{array}{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]

+\zeta \lambda \left[{\begin{array}{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]

+\zeta ^{2}\lambda ^{2}\left[{\begin{array}{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]

(12)

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

(13)

Above expression shows that ${\textstyle \mathbf {U} ^{2}}$ is non-zero at the original configuration for curved surfaces ( ${\textstyle \kappa _{ij}^{0}\neq {0}}$ ). The changes of curvature of the middle surface are computed by

(14)

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

(15)

Above expression of ${\textstyle \mathbf {U} ^{2}}$ 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 ${\textstyle \mathbf {U} }$ as

(16)

where ${\textstyle \lambda _{\alpha }}$ and ${\textstyle \mathbf {r} _{\alpha }}$ are the eigenvalues and eigenvectors of the right stretch tensor ${\textstyle \mathbf {U} }$ .

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

(17)

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

(18)

where ${\textstyle \mathbf {R} _{L}}$ is the rotation tensor obtained from the eigenvectors of ${\textstyle \mathbf {U} }$ given by

(19)

The relationship between the Hencky and Piola-Kirchhoff stresses is

(20)

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

(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.

$\mathbf {N} =\int _{h^{0}}\mathbf {S} d\zeta$	(22.a)
$\mathbf {M} =\int _{h^{0}}\mathbf {S} \lambda \zeta d\zeta$	(22.b)

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

(23)

where ${\textstyle {\boldsymbol {K}}}$ is the virtual curvature vector and ${\textstyle \mathbf {E} }$ is the virtual Green-Lagrange strain vector of the middle surface, with

(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

(25)

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

(26)

These constitutive equations are integrated using a standard return algorithm. The following Mises-Hill [11] yield function with non-linear isotropic hardening is adopted

(27)

where ${\textstyle F,G,H}$ and ${\textstyle N}$ define the non-isotropic shape of the yield surface and the parameters ${\textstyle \kappa ,}$ ${\textstyle \varepsilon _{0}}$ and ${\textstyle n}$ define its size as a function of the effective plastic strain ${\textstyle e^{p}}$ .

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 [12] 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

(28)

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

The stress measures associated to the principal logarithmic strains are denoted by ${\textstyle \beta _{i}}$ . They can be computed noting that

(29)

we define now

(30)

which gives

(31)

The values of ${\textstyle \beta _{i}}$ , expressed in the principal strains directions, allow to evaluate the stresses in the convective coordinate system as

(32)

We note that the Hencky stress tensor ${\textstyle \mathbf {T} }$ 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 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: ${\textstyle \left(\xi ,\eta \right)=\left(0,0\right)}$

node 2: ${\textstyle \left(\xi ,\eta \right)=\left(1,0\right)}$

node 3: ${\textstyle \left(\xi ,\eta \right)=\left(0,1\right)}$

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

node 4: ${\textstyle \left(\xi ,\eta \right)=\left(1,1\right)}$

node 5: ${\textstyle \left(\xi ,\eta \right)=\left(-1,1\right)}$

node 6: ${\textstyle \left(\xi ,\eta \right)=\left(1,-1\right)}$

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:

{\begin{array}{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}}

(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 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 ${\textstyle \mathbf {t} _{1}}$ and ${\textstyle \mathbf {t} _{2}}$ 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

(34)

where the Jacobian matrix at the original configuration is

\mathbf {J=} \left[{\begin{array}{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]

(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

(36)

The covariant metric tensor can be computed from above expression as

\mathbf {g=} \left[{\begin{array}{cc}a_{11}&a_{12}\\a_{21}&a_{22}\end{array}}\right]=\left[{\begin{array}{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]

(37)

The component of ${\textstyle {\boldsymbol {g}}}$ can be used to compute any convenient membrane strain measures. For example the Green-Lagrange strain tensor is obtained by

\mathbf {E} _{GL}\mathbf {=} {\frac {1}{2}}\left[{\begin{array}{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}{cc}a_{11}-1&a_{12}\\a_{21}&a_{22}-1\end{array}}\right]

(38)

In the usual FEM matrix notation

\left[{\begin{array}{c}E_{11}\\E_{22}\\2E_{12}\end{array}}\right]={\frac {1}{2}}\left[{\begin{array}{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]

(39)

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

\delta \left[{\begin{array}{c}E_{11}\\E_{22}\\2E_{12}\end{array}}\right]=\left[{\begin{array}{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]

(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

(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

(42)

where ( ${\textstyle A^{M}}$ is the original area of the element).

Integrating by parts Eq.(43), the following integral over the element boundary ${\textstyle \Gamma ^{M}}$ is obtained:

(43)

The three distinct components of the curvature tensor are

\left[{\begin{array}{c}\kappa _{11}\\\kappa _{22}\\2\kappa _{12}\end{array}}\right]={\frac {-1}{A^{M}}}{\displaystyle \oint _{\Gamma ^{M}}}\left[{\begin{array}{cc}n_{1}&0\\0&n_{2}\\n_{2}&n_{1}\end{array}}\right]\left[{\begin{array}{c}\mathbf {t} _{3}\cdot \mathbf {\varphi } _{^{\prime }1}\\\mathbf {t} _{3}\cdot \mathbf {\varphi } _{^{\prime }2}\end{array}}\right]d\Gamma ^{M}

(44)

The numerical evaluation of the boundary integral in Eq.(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.

\left[{\begin{array}{c}\kappa _{11}\\\kappa _{22}\\2\kappa _{12}\end{array}}\right]={\frac {-1}{A^{M}}}\sum _{G=1}^{3}l_{G}\left[{\begin{array}{cc}n_{1}&0\\0&n_{2}\\n_{2}&n_{1}\end{array}}\right]^{G}\left[{\begin{array}{c}\mathbf {\varphi } _{1}^{G}\cdot \mathbf {t} _{3}\\\mathbf {\varphi } _{2}^{G}\cdot \mathbf {t} _{3}\end{array}}\right]

(45)

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

\left[{\begin{array}{c}\mathbf {\varphi } _{1}\\\mathbf {\varphi } _{2}\end{array}}\right]^{M}={\frac {1}{2A^{M}}}\sum _{i=1}^{3}\left[{\begin{array}{c}-b_{i}\\a_{i}\end{array}}\right]\mathbf {\varphi } ^{i}=\left[{\begin{array}{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}\mathbf {\varphi } ^{1}\\\mathbf {\varphi } ^{2}\\\mathbf {\varphi } ^{3}\end{array}}\right]

(46)

(47)

From these expressions it is also possible to compute in the original configuration the element area ${\textstyle A^{M}}$ , the outer normals ${\textstyle \left(n_{1},n_{2}\right)^{i}}$ at each side and its lengths ${\textstyle l_{i}}$ . Eq.(47) also allows to evaluate the thickness ratio ${\textstyle \lambda }$ in the deformed configuration and the actual normal ${\textstyle \mathbf {t} _{3}}$ . In Eq.(46), ${\textstyle {\bar {N}}^{I}}$ are the linear shape functions of the three node linear triangle [3].

As one integration point is used over each side, it is not necessary to distinguish between sides ( ${\textstyle I}$ ) and integrations points ( ${\textstyle G}$ ) any more. In this way the curvatures can be computed by:

$\left[{\begin{array}{c}\kappa _{11}\\\kappa _{22}\\2\kappa _{12}\end{array}}\right]={\frac {-1}{A^{M}}}\sum _{I=1}^{3}l_{I}\left[{\begin{array}{cc}[n_{1}&0\\0&n_{2}\\n_{2}&n_{1}\end{array}}\right]^{G}\left[{\begin{array}{c}\mathbf {\varphi } _{1}^{I}\cdot \mathbf {t} _{3}\\\mathbf {\varphi } _{2}^{I}\cdot \mathbf {t} _{3}\end{array}}\right]$	(48)
$=2\sum _{I=1}^{3}\left[{\begin{array}{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}\mathbf {\varphi } _{1}^{I}\cdot \mathbf {t} _{3}\\\mathbf {\varphi } _{2}^{I}\cdot \mathbf {t} _{3}\end{array}}\right]$	(49)

In the original BST element [4,10] the gradient ${\textstyle \mathbf {\varphi } _{\alpha }^{I}}$ 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 ${\textstyle I}$ from the quadratic interpolation

\left[{\begin{array}{c}\mathbf {\varphi } _{1}\\\mathbf {\varphi } _{2}\end{array}}\right]^{I}=\left[{\begin{array}{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}\mathbf {\varphi } ^{1}\\\mathbf {\varphi } ^{2}\\\mathbf {\varphi } ^{3}\\\mathbf {\varphi } ^{I+3}\end{array}}\right]

(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 ( ${\textstyle I+3}$ ), associated precisely to the side ( ${\textstyle I}$ ) where the gradient is computed.

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

(51)

This gives

(52)

This last expression allows to interpret the curvatures as the projections of the vectors ${\textstyle \mathbf {h} _{ij}}$ over the normal of the central element. Direction t ${\textstyle _{3}}$ can be seen as a reference direction. If a different direction is chosen, at an angle ${\textstyle \theta }$ with the former, this has an influence of order ${\textstyle \theta ^{2}}$ in the projection. This justifies Eq.(47) for the definition of t ${\textstyle _{3}}$ 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.(52) the variation of the curvatures is obtained as

(53)

The first term is crucial in the evaluation of the curvatures variation. The variations of ${\textstyle \mathbf {h} _{ij}}$ are:

(54)

where

\left[{\begin{array}{c}\delta \mathbf {\varphi } _{1}\\\delta \mathbf {\varphi } _{2}\end{array}}\right]^{I}=\left[{\begin{array}{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}\delta \mathbf {u} ^{1}\\\delta \mathbf {u} ^{2}\\\delta \mathbf {u} ^{3}\\\delta \mathbf {u} ^{I+3}\end{array}}\right]

(55)

The variations of t ${\textstyle _{3}}$ can be computed as shown in the original BST element [10].

(56)

leading to

(57)

where ${\textstyle \mathbf {\tilde {\varphi }} _{^{\prime }\alpha }}$ are the contravariant base vectors in the central triangle

{\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}}

(58)

We then have

(59)

Substituting the last expression in Eq.(53), results

\left[{\begin{array}{c}\delta \chi _{11}\\\delta \chi _{22}\\2\delta \chi _{12}\end{array}}\right]=2\sum \limits _{I=1}^{3}\left[{\begin{array}{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\{\sum \limits _{J=1}^{3}\left[{\begin{array}{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}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\}

-2\sum \limits _{I=1}^{3}\left[{\begin{array}{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)

(60)

where the projections of the vectors ${\textstyle \mathbf {h} _{ij}}$ over the contravariant base vectors ${\textstyle \mathbf {\tilde {\varphi }} _{^{\prime }\alpha }}$ have been included

(61)

These projections together with Eq.(52) allows to write

(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 ${\textstyle \mathbf {t} _{3}}$ 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 ${\textstyle \mathbf {\varphi } _{^{\prime }n}^{\left(0\right)}}$ that does not change during the process.

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 2) and defined as

(63)

where vector ${\textstyle \mathbf {\varphi } _{^{\prime }n}^{\left(0\right)}}$ is fixed while direction ${\textstyle \mathbf {\bar {\varphi }} _{^{\prime }s}}$ emerges as the intersection of the symmetry plane with the plane defined by the central element ( ${\textstyle M}$ ). The plane (gradient) defined by the central element is

(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 ( ${\textstyle J}$ and ${\textstyle K}$ ) and the original length of the side ${\textstyle L_{0}}$ , i.e.

(65)

That together with the outer normal to the side ${\textstyle \mathbf {n=} \left[\nu _{1},\nu _{2}\right]^{T}}$ (resolved in the selected original convective cartesian system) leads to

(66)

where the normal component of the gradient ${\textstyle \mathbf {\varphi } _{^{\prime }n}}$ is

(67)

In this way the contribution of the gradient at side ${\textstyle I}$ to vectors ${\textstyle \mathbf {h} _{ij}}$ results in

\left[{\begin{array}{c}\mathbf {h} _{11}^{T}\\\mathbf {h} _{22}^{T}\\2\mathbf {h} _{12}^{T}\end{array}}\right]^{I}=2\left[{\begin{array}{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}{cc}\nu _{1}&-\nu _{2}\\\nu _{2}&\nu _{1}\end{array}}\right]\left[{\begin{array}{c}\mathbf {\varphi } _{^{\prime }n}^{T}\\\mathbf {\varphi } _{^{\prime }s}^{T}\end{array}}\right]

(68)

The only difference with Eq.(54) necessary for the computation of the curvature variations, is that now the contribution from the gradient at side ${\textstyle I}$ is

\delta \left[{\begin{array}{c}\mathbf {h} _{11}^{T}\\\mathbf {h} _{22}^{T}\\2\mathbf {h} _{12}^{T}\end{array}}\right]^{I}=2\left[{\begin{array}{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}{cc}\nu _{1}&-\nu _{2}\\\nu _{2}&\nu _{1}\end{array}}\right]\left[{\begin{array}{c}\mathbf {0} \\{\frac {1}{L_{o}}}\left[\delta \mathbf {u} ^{K}-\delta \mathbf {u} ^{J}\right]^{T}\end{array}}\right]

(69)

={\frac {2}{L_{0}}}\left[{\begin{array}{c}-{\bar {N}}_{^{\prime }1}^{I}\nu _{2}\\{\bar {N}}_{^{\prime }2}^{I}\nu _{1}\\{\bar {N}}_{^{\prime }1}^{I}\nu _{1}-{\bar {N}}_{^{\prime }2}^{I}\nu _{2}\end{array}}\right]\left[\delta \mathbf {u} ^{K}-\delta \mathbf {u} ^{J}\right]^{T}

(70)

where the influence of variations in the length of vector ${\textstyle \mathbf {\varphi } _{^{\prime }n}}$ has been neglected.

In the case of natural boundary conditions (either free or simple supported) the gradient at the sides is supposed to be equal to the gradient in the central element, similarly as in the membrane approach, i.e.

(71)

Using Eq.(71), vectors ${\textstyle \mathbf {h} _{ij}}$ and their variations can be easily computed. A second possibility is to make use of the natural boundary condition constraining the normal curvature to a zero value. In simple supported sides where the curvature along the side is zero, this leads to zero values for both bending moments.

6 Stiffness matrix evaluation

When a predictor-corrector scheme is used to trace the movement of the shell, the derivative (stiffness matrix) of the weak form of the momentum equations (23) is needed. As usual, material and geometric components are computed separately. The material part does not offer difficulties and it is computed from the integral

(72)

where matrix ${\textstyle \mathbf {B} ={\boldsymbol {B}}_{m}+{\boldsymbol {B}}_{\phi }}$ includes:

a) a membrane part ${\textstyle {\boldsymbol {B}}_{m}}$ computed at each mid side point ${\textstyle I}$ from

\delta \left[{\begin{array}{c}E_{11}\\E_{22}\\2E_{12}\end{array}}\right]^{I}=\left[{\begin{array}{cc}\mathbf {\varphi } _{^{\prime }1}^{T}&\mathbf {0} _{3\times {1}}^{T}\\\mathbf {0} _{3\times {1}}^{T}&\mathbf {\varphi } _{^{\prime }2}^{T}\\\mathbf {\varphi } _{^{\prime }2}^{T}&\mathbf {\varphi } _{^{\prime }1}^{T}\end{array}}\right]^{I}\left[{\begin{array}{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]\left[{\begin{array}{c}\delta \mathbf {u} ^{1}\\\delta \mathbf {u} ^{2}\\\delta \mathbf {u} ^{3}\\\delta \mathbf {u} ^{I+3}\end{array}}\right]

=\left[{\begin{array}{cc}\mathbf {\varphi } _{^{\prime }1}^{T}&\mathbf {0} _{3\times {1}}^{T}\\\mathbf {0} _{3\times {1}}^{T}&\mathbf {\varphi } _{^{\prime }2}^{T}\\\mathbf {\varphi } _{^{\prime }2}^{T}&\mathbf {\varphi } _{^{\prime }1}^{T}\end{array}}\right]^{I}\left[{\begin{array}{cccc}N_{^{\prime }1}^{1\left(I\right)}&N_{^{\prime }1}^{2\left(I\right)}&N_{^{\prime }1}^{3\left(I\right)}&N_{^{\prime }1}^{4\left(I\right)}\\N_{^{\prime }2}^{1\left(I\right)}&N_{^{\prime }2}^{2\left(I\right)}&N_{^{\prime }2}^{3\left(I\right)}&N_{^{\prime }2}^{4\left(I\right)}\end{array}}\right]\left[{\begin{array}{c}\delta \mathbf {u} ^{1}\\\delta \mathbf {u} ^{2}\\\delta \mathbf {u} ^{3}\\\delta \mathbf {u} ^{4}={\boldsymbol {B}}\delta {\boldsymbol {u}}\end{array}}\right]^{\left(I\right)}

(73)

Typically one integration point is used for computing the terms in ${\textstyle {\boldsymbol {B}}_{m}}$ .

b) a bending part ${\textstyle {\boldsymbol {B}}_{\phi }}$ which results from Eq.(60) that is constant over the element.

In Eq.(72), C is the elasticity matrix integrated in the thickness. If material non-linearity is considered, a layer wise numerical integration across the thickness is required This leads to the tangent or algorithmic elastic-plastic constitutive matrix ${\textstyle \mathbf {C} _{ep}}$ .

6.1 Geometric stiffness (membrane part)

The geometric stiffness due to membrane forces results from computing

(74)

This can be written as the sum of the contributions of the three sides, i.e.

={\frac {A^{M}}{3}}\sum _{K=1}^{3}\sum _{I=1}^{4}\sum _{J=1}^{4}\left[N_{^{\prime }1}^{I\left(K\right)}N_{^{\prime }1}^{J\left(K\right)}\;N_{11}^{\left(K\right)}+N_{^{\prime }2}^{I\left(K\right)}N_{^{\prime }2}^{J\left(K\right)}\;N_{22}^{\left(K\right)}\right.

={\frac {A^{M}}{3}}\sum _{K=1}^{3}\sum _{I=1}^{4}\sum _{J=1}^{4}\left\{\delta \mathbf {u} ^{I}\;\left[{\begin{array}{cc}N_{^{\prime }1}^{J}&N_{^{\prime }2}^{J}\end{array}}\right]\left[{\begin{array}{c}[c]{cc}N_{11}&N_{12}\\N_{21}&N_{22}\end{array}}\right]\left[{\begin{array}{c}N_{^{\prime }1}^{J}\\N_{^{\prime }2}^{J}\end{array}}\right]\Delta \mathbf {u} ^{J}\right\}^{\left(K\right)}

(75)

6.2 Geometric stiffness (bending part)

The geometric stiffness associated to bending moments is more involved. Its expression stems from

(76)

Recalling expressions (51) and (52) it can be written

(77)

where

(78)

The first two terms lead to symmetric components. The second one ${\textstyle \left(\delta \mathbf {t} _{3}\cdot \Delta \mathbf {h} _{ij}\right)}$ can be expressed as

=\left\{-\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)\right\}\cdot \left\{\sum _{I=1}^{3}\sum _{K=1}^{4}\left({\bar {N}}_{^{\prime }i}^{I}\;N_{^{\prime }j}^{K}+{\bar {N}}_{^{\prime }j}^{I}\;N_{i}^{K}\right)\Delta \mathbf {u} ^{K\left(I\right)}\right\}

=-\sum _{J=1}^{3}\sum _{I=1}^{3}\sum _{K=1}^{4}\left(\delta \mathbf {u} ^{J}\right)^{T}\left[{\bar {N}}_{^{\prime }1}^{J}\left(\mathbf {t} _{3}\otimes \mathbf {\tilde {\varphi }} _{^{\prime }1}\right)+{\bar {N}}_{^{\prime }2}^{J}\left(\mathbf {t} _{3}\otimes \mathbf {\tilde {\varphi }} _{^{\prime }2}\right)\right]\left({\bar {N}}_{^{\prime }i}^{I}\;N_{^{\prime }j}^{K}+{\bar {N}}_{^{\prime }j}^{I}\;N_{i}^{K}\right)\Delta \mathbf {u} ^{K\left(I\right)}

(79)

Then, a first component of the geometric bending stiffness can be written as

\delta \mathbf {u} ^{T}\mathbf {K} _{b1}^{G}\Delta \mathbf {u} =A^{M}\left\{-\sum _{J=1}^{3}\left(\delta \mathbf {u} ^{J}\right)^{T}\left[{\begin{array}{c}[c]{cc}{\bar {N}}_{^{\prime }1}^{J}\mathbf {1} &{\bar {N}}_{^{\prime }2}^{J}\mathbf {1} \end{array}}\right]\left[{\begin{array}{c}\mathbf {t} _{3}\otimes \mathbf {\tilde {\varphi }} _{^{\prime }1}\\\mathbf {t} _{3}\otimes \mathbf {\tilde {\varphi }} _{^{\prime }2}\end{array}}\right]\right\}

\sum _{I=1}^{3}\left[{\begin{array}{cc}{\bar {N}}_{^{\prime }1}^{I}&{\bar {N}}_{^{\prime }2}^{I}\end{array}}\right]\left[{\begin{array}{cc}M_{11}&M_{12}\\M_{12}&M_{22}\end{array}}\right]\sum _{K=1}^{4}\left[{\begin{array}{c}N_{^{\prime }1}^{K\left(I\right)}\\N_{^{\prime }2}^{K\left(I\right)}\end{array}}\right]\Delta \mathbf {u} ^{K\left(I\right)}

(80)

The last term in (78), can be obtained observing that (sum in ${\textstyle \alpha }$ and ${\textstyle \beta }$ )

$\Delta \left(\delta \mathbf {t} _{3}\right)=\left[\Delta \left(\delta \mathbf {t} _{3}\right)\right]_{1}\mathbf {\tilde {\varphi }} _{^{\prime }1}+\left[\Delta \left(\delta \mathbf {t} _{3}\right)\right]_{2}\mathbf {\tilde {\varphi }} _{^{\prime }2}+\left[\Delta \left(\delta \mathbf {t} _{3}\right)\right]_{3}\mathbf {t} _{3}$
$=-\left(\mathbf {\tilde {\varphi }} _{^{\prime }\alpha }\cdot \mathbf {\tilde {\varphi }} _{^{\prime }\beta }\right)\left[\delta \mathbf {\varphi } _{^{\prime }\alpha }^{T}\left(\mathbf {t} _{3}\otimes \mathbf {t} _{3}\right)\Delta \mathbf {\varphi } _{^{\prime }\beta }\right]\mathbf {t} _{3}+$	(81)
$+\left[\delta \mathbf {\varphi } _{^{\prime }\alpha }^{T}\left(\mathbf {t} _{3}\otimes \mathbf {\tilde {\varphi }} _{^{\prime }\alpha }\right)\Delta \mathbf {\varphi } _{^{\prime }\beta }\right]\mathbf {\tilde {\varphi }} _{^{\prime }\beta }+\left[\delta \mathbf {\varphi } _{^{\prime }\alpha }^{T}\left(\mathbf {\tilde {\varphi }} _{^{\prime }\beta }\otimes \mathbf {t} _{3}\right)\Delta \mathbf {\varphi } _{^{\prime }\beta }\right]\mathbf {\tilde {\varphi }} _{^{\prime }\alpha }$

then

\Delta \left(\delta \mathbf {t} _{3}\right)\cdot \mathbf {h} _{ij}=-\left(\mathbf {\tilde {\varphi }} _{^{\prime }\alpha }\cdot \mathbf {\tilde {\varphi }} _{^{\prime }\beta }\right)\left[\delta \mathbf {\varphi } _{^{\prime }\alpha }^{T}\left(\mathbf {t} _{3}\otimes \mathbf {t} _{3}\right)\Delta \mathbf {\varphi } _{^{\prime }\beta }\right]\left(\mathbf {t} _{3}\cdot \mathbf {h} _{ij}\right)+

=\delta \mathbf {\varphi } _{^{\prime }\alpha }^{T}\left\{\rho _{ij}^{\beta }\left(\mathbf {t} _{3}\otimes \mathbf {\tilde {\varphi }} _{^{\prime }\alpha }\right)+\rho _{ij}^{\alpha }\left(\mathbf {\tilde {\varphi }} _{^{\prime }\beta }\otimes \mathbf {t} _{3}\right)-g^{\alpha \beta }\kappa _{ij}\left(\mathbf {t} _{3}\otimes \mathbf {t} _{3}\right)\right\}\Delta \mathbf {\varphi } _{^{\prime }\beta }

(82)

(83)

with

(84)

where the covariant metric tensor at the middle surface ${\textstyle a^{\alpha \beta }}$ has been used.

Denoting the sum

(85)

the term

(86)

results in

(87)

The last expression (87) has components only in the nodes of the central element, which stems from the definition used for t ${\textstyle _{3}}$ .

Numerical experiments have shown that the bending part of the geometric stiffness is not so important and can be disregarded in the iterative process.

7 Numerical examples in the linear range

In this section, a summary of the results obtained in the assessment of the proposed element in the linear range are presented. In this first part of the numerical experiments, results have been obtained with a static/dynamic code that uses a implicit solution of the discretized equilibrium equations. In the examples, the original BST element [10] is compared with the enhanced version here proposed, denoted by CBST when a numerical integration is performed with three points, and CBST1 when only one integration point is used. Note that one integration quadrature is equivalent to averaging the metric tensors computed at each side.

7.1 Membrane patch test

The main objective of the present formulation is to obtain a (non conforming) membrane element with a similar performance than the linear strain triangle, that satisfies the patch test. To assess this, a square domain subjected to nodal forces associated to a constant stress state (in both directions and in shear) is considered. Figure 3 shows the patch of elements and the necessary nodal forces to obtain a uniform tensile stress in direction ${\textstyle x_{1}}$ . Note that the nodal forces used correspond to the constant strain triangle, not to the linear one. For the distorted mesh shown, correct results are obtained using either 1 or 3 integration points.

Figure 3: Membrane patch test

7.2 Bending patch test (torsion)

The element bending formulation does not allow to apply external bending moments (there are not rotational DOFs). Hence it is not possible to analyse a patch of elements under loads leading to a uniform bending moment. A uniform torsion can be considered if a point load is applied at the corner of a rectangular plate with two consecutive free sides and two simple supported sides. Figure 4 shows three patches leading to correct results both in displacements and stresses. All three patches are structured meshes. When the central node in the third patch is displaced from the center, the results obtained with the CBST element are not correct. The original element BST gives correct results in all cases, if natural boundary conditions are imposed in the formulation. If this is not the case, incorrent results are obtained even with structured meshes.

Figure 4: Patch test for uniform torsion

7.3 Cook's membrane problem

This example is used to assess the membrane performance of the CBST element and to compare it with the linear triangle (constant strain) and the quadratic triangle (linear strain). This example involves important shear energy and was proposed to assess the distortion capability of elements. Figure 5.a shows the geometry and the applied load. Figure 5.b plots the vertical displacement of the upper vertex as a function of the number of nodes in the mesh. In this figure results obtained with other isoparametric elements have been also plotted for comparison. They include the constain strain triangle (CST), the bilinear quadrilateral (QUAD4) and the linear strain triangle (LST). Note that in this membrane problem both the BST and the CST elements give identical results.

Figure 5: Cook membrane problem (a) Geometry (b) Results

From the plot shown it can be seen that the new element with three integration points (CBST) gives values slightly better that the constant strain triangle for the most coarse mesh (only two elements). However, when the mesh is refined, a performance similar to the linear strain triangle is obtained, that is dramatically superior that the former. On the other hand, if a one point quadrature is used (CBST1) the convergence in the reported displacement is notably better that for the rest of the elements.

7.4 Cylindrical roof

In this example an effective membrane interpolation is of primary importance. Hence this is good test to assess the new element. The geometry is a cylindrical roof supported by a rigid diaphragm at both ends and it is loaded by a uniform dead weight (see 6.a.). Only one quarter of the structure is modelled due to symmetry conditions. Unstructured and structured meshes are considered. In the latter case two orientations are possible (Figure 6 shows orientation B).

Tables 1, 2 and 3 present the normalized vertical displacements at the crown (point A) and in the midpoint of the free side (point B) for the two orientations of structured meshes and non-structured meshes. Values used for normalization are ${\textstyle u_{A}=0.5407}$ y ${\textstyle u_{B}=-3.610}$ that are quoted in reference [13].

Figure 6: Cylindrical roof under dead weight. E=3 E

,

, Thickness =3.0, shell weight =0.625 per unit area.

Table. 1 Cylindrical roof under dead weigth. Normalized displacements for mesh orientation A
	Point-A			Point-B
NDOFs	CBST	CBST1	BST	CBST	CBST1	BST
16	0.65724	0.91855	0.74161	0.40950	0.70100	1.35230
56	0.53790	1.03331	0.74006	0.54859	1.00759	0.75590
208	0.89588	1.04374	0.88491	0.91612	1.02155	0.88269
800	0.99658	1.01391	0.96521	0.99263	1.00607	0.96393
3136	1.00142	1.00385	0.99105	0.99881	1.00102	0.98992

Table. 2 Cylindrical roof under dead weigth. Normalized displacements for mesh orientation A
	Point-A			Point-B
NDOFs	CBST	CBST1	BST	CBST	CBST1	BST
16	0.26029	0.83917	0.40416	0.52601	0.86133	0.89778
56	0.81274	1.10368	0.61642	0.67898	0.93931	0.68238
208	0.97651	1.04256	0.85010	0.93704	1.00255	0.86366
800	1.00085	1.01195	0.95626	0.99194	1.00211	0.95864
3136	1.00129	1.00337	0.98879	0.99828	1.00017	0.98848

Table. 3 Cylindrical roof under dead weigth. Normalized displacements for non-structured mesh
	Point-A			Point-B
NDOFs	CBST	CBST1	BST	CBST	CBST1	BST
851	0.97546	0.8581	0.97598	0.97662	1.0027	0.97194
3311	0.98729	0.9682	0.98968	0.98476	1.0083	0.98598
13536	0.99582	0.9992	1.00057	0.99316	0.9973	0.99596

Figure 6.b shows the normalized displacement of point-B over the structured meshes as a function of the number of degrees of freedom. An excelent convergence the CBST element can be seen. The version with only one integration point (CBST1) presents a behavior a little more flexible and converges from above. Table 3 shows that for non-structured meshes the element converges to the correct value but more slowly.

7.5 Open semi-espherical dome with point loads

The main problem of finite elements with initially curved geometry is the so called membrane locking. The CBST element proposed has a quadratic interpolation of the geometry, then it may suffer from this problem. To assess this we resort to an example of inextensional bending. This is an hemispherical shell of radius ${\textstyle r=10}$ and thickness ${\textstyle h=0.04}$ with an 18 ${\textstyle ^{o}}$ hole in the pole and free at all boundaries, subjected to two inward and two outward forces 90 ${\textstyle ^{o}}$ apart. Material properties are ${\textstyle E=6.825\times {10}^{7}}$ and ${\textstyle \nu =0.3}$ . Figure 7.a shows the discretized geometry (only one quarter of the geometry is considered due to symmetry).

Figure 7: Pinched hemispherical shell with a hole, (a)geometry, (b)normalized displacement

In Figure 7.b the displacements of the points under the loads have been plotted versus the number of nodes used in the discretization. Due to the orientation of the meshes chosen, the displacement of the point under the inward load is not the same as the displacement under the outward load, so in the figure an average (the absolute values) has been used. Results obtained with other elements have been included for comparison: two membrane locking free elements namely the original linear BST and a transverse shear-deformable quadrilateral (QUAD); a transverse shear deformable quadratic triangle (TRIA) [1] that is vulnerable to locking and an assumed strain version of this triangle (TRIA-1)^[1] that alleviates membrane locking but does not eliminate it.

From the plotted results it can be seen that the CBST element presents membrane locking in bending dominated problems with initially curved geometries. This locking is much less severe than in a standard quadratic triangle and even than for the assumed strain version used for comparison. When only one integration point is used (CBST1 element) membrane locking disappears.

8 Non linear numerical examples

Results for examples with geometric and material non-linearities are presented next using the CBST1 element. Due to the features of the modelized problems, with strong nonlinearities associated to instabilities and frictional contact conditions, a code with explicit integration of the dynamic equilibrium equations has been used [14]. This code allows to obtain pseudo-static solutions throught dynamic relaxation.

8.1 Inflation of a sphere

The example is the inflation of a spherical shell under internal pressure. An incompressible Mooney-Rivlin constitutive material have been considered. The Ogden parameters are ${\textstyle N=2}$ , ${\textstyle \alpha _{1}=2}$ , ${\textstyle \mu _{1}=40}$ , ${\textstyle \alpha _{2}=-2}$ , ${\textstyle \mu _{2}=-20}$ . Due to the simple geometry an analytical solution exists [15] (with ${\textstyle \gamma =R/R^{\left(0\right)}}$ ):

In this numerical simulation the same geometric and material parameters used in Ref. [16] have been adopted: ${\textstyle R^{\left(0\right)}=1}$ and ${\textstyle h^{\left(0\right)}=0.02}$ . The three meshes considered to evaluate convergence are shown in Figure 8.a. The value of the actual radius as a function of the internal pressure is plotted in Figure 8.b for the different meshes and is also compared with the analytical solution. It can be seen that with a few degrees of freedom it is possible to obtain an excellent agreement for the range of strains considered. The final value corresponds to a thickness radius ratio of ${\textstyle h/R=0.00024}$ .


Figure 8: Inflation of sphere of Mooney-Rivlin material. (a) Meshes used in the analysis (b) Change of radius as a function of the internal pressure.

8.2 Inflation of an air-bag

This example has also been taken from Ref.[16] where it is shown that the final configuration is mesh dependent due to the strong instabilities leading to a non-uniqueness of the solution. In [16] it is also discussed the important regularizing properties of the bending energy, that when disregarded leads to massive wrinkling in the compressed zones.

The air bag geometry is initially square with an undeformed diagonal of ${\textstyle 1.20}$ . The constitutive material is a linear isotropic elastic one with modulus of elasticity ${\textstyle E=5.88\times {10}^{8}}$ and Poisson's ratio ${\textstyle \nu =0.4}$ . Only one quarter of the geometry has been modelled due to symmetry. Only the normal to the original plane is constrained along the boundaries. The thickness considered is ${\textstyle h=0.001}$ and the inflation pressure is ${\textstyle 5000}$ . Using a density ${\textstyle \delta =1000}$ , pressure is linearly incremented from 0 to the final value in ${\textstyle t=0.1}$ .

With comparative purposes and also to backup the comments in Ref. [16] two analyses have been performed, a purely membranal one and one including bending effects. Figure 9 shows the final deformed configurations for three meshes with 289, 1089 and 4225 nodes. The top row corresponds to a full analysis including bending and the central row is a pure membrane analysis. The bottom row is also an analysis including bending where the mesh orientation has been changed.

Figure 9: Inflation of a square air-bag . Deformed configurations for three different meshes with 800, 3136 and 12416 degrees of freedom.

The top and bottom lines show the final shapes change according to the degree of discretization and mesh orientation due to instabilities and non uniqueness of the solution. The central row shows the noteworthy increment of wrinkles as the mesh is refined. Finally in Figure 10 presents the convergence of the maximum displacement of the central point as the mesh is refined.

Figure 10: Inflation of a square air-bag. Convergence of the maximum displacement versus the number of degrees of freedom in the model.

8.3 Deep drawing of a rolled sheet

The element performance in problems involving large strains and anisotropic plastic behaviour is assesed in a benchmark proposed in a recent NUMISHEET [17] meeting. This is the deep drawing of a circular mild steel sheet with a spherical punch of radius ${\textstyle 50mm}$ ( ${\textstyle R_{p}}$ ). The original radius of the sheet ( ${\textstyle R_{s}}$ ) is ${\textstyle 100mm}$ ( ${\textstyle {\frac {R_{L}}{R_{P}}}=\beta =2}$ ). The drawing depth is ( ${\textstyle DD}$ ) is 85 mm ( ${\textstyle {\frac {DD}{R_{P}}}=1.7}$ ) and the force over the blank holder is ${\textstyle 80kN}$ . The constitutive material of the rolled sheet is characterized by three tensile tests along three different directions respect to rolling direction ( ${\textstyle RD}$ ) (Table 4). Note that the modellization requires to treat the material as transversally anisotropic, as this improves the deep drawability and if not considered leads to an early localized necking of the material. The yield function originally proposed by Hill^[11] has been chosen for both associative and non-associative potential functions. For elastic-plastic problems a numerical integration along the thickness direction is necessary. This is accomplished here using four integration points.

Table. 4 Plastic characterization of the constitutive material
Thickness	Orientation	Yield	Tensile	$\varepsilon _{u}$	$\varepsilon _{t}$	n	r
	respect to RD	stres	strength	(uniform)	(total)
[mm]	[ ${\textstyle ^{o}}$ ]	[N/mm ${\textstyle ^{2}}$ ]	[N/mm ${\textstyle ^{2}}$ ]	[%]	[%]	-	-
	0	176	322	24	40	0.214	1.73
0.98	45	185	333	22	39	0.203	1.23
	90	180	319	23	44	0.206	2.02

Only one quarter of the geometry has been considered due to symmetry. Tools are modelized as rigid surfaces, the punch has been defined by 1439 points and 2730 triangles. For the die 744 points and 490 quadrilaterals have been used. Finally the blank holder has been discretized with 155 points and 120 quadrilaterals. Figure 11.a shows a perspective of the tools and Figure 11.a the final deformed sheet. The sheet has been modelled with 6370 3-node triangles and 3284 nodes (9274 DOFs).

Figure 11: Deep drawing of a circular sheet. (a)geometry of the tools, (b)final deformed shape of the sheet.

The reported results are associated to three different meridiands: A at 90 ${\textstyle ^{0}}$ of the rolling direction, B: at 45 ${\textstyle ^{0}}$ of the rolling direction and C: in the rolling direction. The numerical results are compared with a set of experimental values reported by Thyssen Krupp Stahl AG (who proposed the benchmark and supplied the sheet samples) [17]. It must be noted that there was a great dispersion in the experimental data send to NUMISHEET, so it is difficult to extract conclusions from a unique comparison. In any case the experimental data supplied seems to be consistent enought to have an idea if the numerical model gives reasonable results.

The first element of comparison associated to the planar anisotropy, are the different displacements (draw in) of the points in the boundary of the sheet. Table 5 shows the values obtained for the three reference meridians chosen.

Table. 5 Flange draw in at three reference meridians
	Draw in [mm]
Model	Section A	Section B	Section C
CBST1-Non Associative	29.06	31.91	30.77
CBST1-Associative	27.08	31.36	28.95
Experimental	30.75	32.30	30.00

The values for the non-associative model are much like the experimental ones, those obtained with the associative model are a little lower. The largest draw-in is at the meridian at 45 ${\textstyle ^{o}}$ of rolling direction ( B). In the numerial simulation the draw-in at meridian C is larger than in meridian A, while in the experimental results the latter are larger than the former. It should be reminded that other numerical results presented at NUMISHEET and most of the numerical simulations showed the same tendency of present values.

The experimental data supplied are the three logarithmic principal strains associated to the three meridians defined above. Here one strain at each meridian has been choosen for comparison. Figure 12.a shows the meridian strain ( ${\textstyle E_{1}}$ ) along meridian A (direction transversal to rolling); Figure 12.b shows the circunferencial strain ( ${\textstyle E_{2}}$ ) along meridian B and in Figure 12.c the thickness strain ( ${\textstyle E_{3}}$ ) along meridian C (rolling direction) has been plotted.

The experimental data show a saw-tooth profile (specially for the thickness strain) that is difficult to accept. These values are therefore not reliable point-wise but as a whole. Based on this fact and on the above mentioned differences between experimental values, it can be said that the present simulations agree quite well with the experimental values, specially when the non-associative plasticity model was used.

In the central part of the sheet the experimental data gives strains lower than the numerical simulation. It can be said that experimentally the sheet has better drawability. In the external zone there are less differences for in plane strains but for thickness strains they are larger. Experimentally thickness strains are more or less uniform between ${\textstyle -0.1}$ and ${\textstyle -0.2}$ what differs from present numerical simulations and also from other numerical simulations and experimental data presented at NUMISHEET [17].

Finally Figure 12 shows resistance to punch as a funtion of punch travel. In the final part simulation forces are smaller than experimental values. This may be due to an incorrect definition of friction with the tools.

Figure 12: Results from NUMISHEET benchmark. (a) Meridional strain

along transversal direction to rolling (meridian A); (b) hoop strain

along direction at 45

of rolling (meridian B); (c) Thickness strain

along rolling direction (meridian C); (d) Punch force versus Punch travel;

From the present results it is difficult to draw definite conclusions about what should be improved in the numerical model due to all the factors involved, including the constitutive model, the thin shell theory, the contact treatment, the reliability of the experimental data of both the material characterization and the deep drawing, etc..

9 Conclusions

A membrane and bending non-conforming rotation-free shell finite element has been presented. The element pass the membrane patch test and the numerical experiments performed do not show problems. From the bending point of view the element is a little stiffer than the original BST element, but presents a smoother continuity of displacement gradients. Convergence rate in membrane dominated problems is similar to the linear strain triangle. Using three integration points the elements is vulnerable to membrane locking for initially curved surfaces. This locking effect dissapears when only one integration point is used, which the usual case for elastic-plastic problems and for any large scale simulation.

Acknowledgments

The first author is a member of the scientific staff of the Science Research Council of Argentina (CONICET). The support provided by grants of CONICET and Agencia Córdoba Ciencia S.E. is gratefully acknowledged.

Bibliography

[1] F.G. Flores, E. Oñate and F. Zárate “New assumed strain triangles for non-linear shell analysis”, Computational Mechanics, 17, 107-114, 1995.

[2] J.H. Argyris, M. Papadrakakis, C. Aportolopoulun and S. Koutsourelakis. The TRIC element. Theoretical and numerical investigation. Comput. Meth. Appl. Mech. Engrg., 182, 217–245, 2000.

[3] O.C. Zienkiewicz and R.L. Taylor. The finite element method. Solid Mechanics. Vol II, Butterworth-Heinemann, 2000.

[4] E. Oñate and F. Zárate, “Rotation-free plate and shell triangles”, Int. J. Num. Meth. Engng, 47, pp. 557–603, 2000.

[5] F. Cirak and M. Ortiz, Subdivision surfaces: A new paradigm for thin-shell finite element analysis, Int. J. Num. Meths in Engng, vol. 47, 2000, págs. 2039-2072.

[6] D.Y. Yang, D.W. Jung, L.S. Song, D.J. Yoo and J.H. Lee, “Comparative investigation into implicit, explicit and iterative implic/explicit schemes for simulation of sheet metal forming processes” NUMISHEET'93, A. Makinouchi, E. Nakamachi, E. Oñate and R.H. Wagoner (Eds.), RIKEN, pp. 35–42, Tokyo, 1993.

[7] M. Brunet and F. Sabourin, “Prediction of necking and wrinkles with a simplified shell element in sheet forming”, Int. Conf. of Metal Forming Simulation in Industry, Vol. II, pp. 27–48, B. Kröplin (Ed.), 1994.

[8] G. Rio, B. Tathi and H. Laurent, “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 Technoloy, Ed. M.J.M. Barata Marques, 18th IDDRG Biennial Congress, Lisbon, 1994.

[9] Rojek J., Oñate E. and Postek E ., Application of explicit FE codes to simulation of sheet and bulk metal forming processes Journal of the Materials Processing Thecnology, 41, págs. 620-627, 1998.

[10] F.G. Flores F.G. and E. Oñate A basic thin shell triangle with only translational DOFs for large strain plasticity, Int. J. Num. Meths in Engng, vol. 51, pp 57-83, 2001.

[11] R. Hill, “A Theory of the Yielding and Plastic Flow of Anisotropic Metals”, Proc. Royal Society London, vol. A193, pp. 281, 1948.

[12] R.W. Ogden “Large deformation isotropic elasticity: on the correlation of theory and experiments for incompressible rubberlike solids”, Proceedings of the Royal Society of London, vol. A326, pp. 565-584, 1972.

[13] H.C. Huang, Static and Dynamic Analysis of Plates and Shells, page 40, Springer-Verlag, Berlin, 1989.

[14] STAMPACK, A General Finite Element System for Sheet Stamping and Forming Problems, Quantech ATZ, Barcelona, Spain, 2003 (www.quantech.es).

[15] A. Needleman, Inflation of spherical rubber ballons, Int. J. of Solids and Structures, 13, 409–421, 1977.

[16] F. Cirak and M. Ortiz, Fully ${\textstyle C^{1}}$ -conforming subdivision elements for finite deformations thin-shell analysis , Int. J. Num. Meths in Engng, vol. 51, 2001, 813-833.

[17] NUMISHEET'99, Fourth International Conference and Workshop on Numerical Simulation of 3D Sheet Forming Processes, NUMISHEET'99, J.C. Gelin and P. Picart (Eds.), BURS Edition, Besancon, France, 1999.