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

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

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

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

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

The deformation gradient can be written as

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

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

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

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

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

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

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

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

The relationship between the Hencky and Piola-Kirchhoff stresses is

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

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.

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

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

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

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

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

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

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

we define now

which gives

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

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:

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

where the Jacobian matrix at the original configuration is

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

The covariant metric tensor can be computed from above expression as

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

In the usual FEM matrix notation

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

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

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

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:

The three distinct components of the curvature tensor are

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.

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)

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:

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

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:

This gives

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

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

where

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

leading to

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

We then have

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

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

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

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.