Rotation-free finite element for the non-linear analysis of beams and axisymmetric shells

Latest revision as of 10:35, 19 October 2018

Published in Comput. Methods Appl. Mech. Engrg., Vol. 195, pp. 5297–5315, 2006
doi:10.1016/j.cma.2005.08.021

Abstract

In this paper a finite element for the non-linear analysis of two dimensional beams and axisymmetric shells is presented. The element uses classical thin shell assumptions (no transverse shear strains). The main feature of the element is that it has no rotational degrees of freedom. Curvatures are computed using geometrical information from the patch of three elements formed by the main element and the two neighbour (adjacent) elements. Special attention is devoted to non-smooth geometries and branching shells. An elastic-plastic material law is considered. Large strains are treated using a logarithmic strain measure and a through-the-thickness numerical integration of the constitutive equations. Several examples are presented including linear problems to study convergence properties, and non-linear problems for both elastic and elastic-plastic materials and large strains.

Keywords: finite elements, beams, axisymmetric shells, rotation-free element, large strains

1 Introduction

The development of numerical techniques for shell analysis without rotational degrees of freedom has been mainly associated to the finite differences method (see for example Refs. [1,2,3]. Nevertheless, the idea of developing finite elements without rotational DOFs is not new [4,5] and different attempts have been reported in the last twenty years [6,7,8,9,10]. Despite these attempts it is just in the last few years that reliable rotation-free elements for industrial applications have been developed [11,12,13,14]. All the approximations share in common the definition of a patch (neighborhood) of elements to interpolate the geometry and the displacements. The main difference between the different approaches is the interpolation of the geometry and the theoretical basis used. One of the main aspects that remains unsolved satisfactorily is the treatment of non-smooth surfaces and specially branching shells. An adequate handling of multiple surfaces, i.e. when more than two elements share a side or edge, is mandatory if the element is to be used for the analysis of aeronautic and naval structures or frame structures typical in buildings, among others.

In this paper two-dimensional shell problems are tackled (previous work in this line can be found in Ref. [15,16]) with special focus on non-smooth shells and branching lines. This work may be seen as a first step towards the development of three dimensional rotation-free shell finite elements capable of handling kinked and branching shells.

The outline of the paper is as follows. First a summary of the most relevant equations governing the kinematic behavior of two dimensional Kirchhoff-Love shells is presented. In Section 3 a two-node rotation-free finite element for the analysis of straight but non-homogeneous (different thickness or materials) beam/axisymmetric shells is developed. In Section 4 the element is extended to deal with curved and non-smooth shells. The formulation is extended further to branching shells in Section 5. Several examples are presented in Section 6 showing convergence properties in linear and non-linear problems. At the end of this section two examples including elasto-plasticity with large strains are shown. In Section 7 some conclusions are drawn.

2 Governing equations of beams and thin two-dimensional shells

In this section a summary of the most relevant equations governing the kinematic behavior of thin shells is presented. More details can be found in the wide literature dedicated to the subject [1]. In order to introduce the problem, the kinematics of two dimensional beams are considered. These equations are latter extended to axisymmetric shells.

Figure 1: Beam geometry and notation.

2.1 Euler-Bernoulli beams

Within a classical Euler-Bernoulli beam theory (transverse shear strains disregarded) the geometry of a curved beam (see Figure 1) is completely defined by the position of central axis (line of centroids) ${\textstyle \mathbf {\varphi } }$ as a function of the arc length ${\textstyle s}$ measured from an arbitrary point on the reference (undeformed) configuration

(1)

where ${\textstyle \varphi _{i}}$ are the Cartesian components referred to the canonical base ${\textstyle \left(\mathbf {e} _{1},\mathbf {e} _{2}\right)}$ . At each material point of the line of centroids associated to the coordinate ${\textstyle s}$ it is possible to compute a convective system defined by the tangent ${\textstyle \mathbf {t} }$ to the line of centroids ${\textstyle \mathbf {\varphi } }$ and the normal (in the plane) ${\textstyle \mathbf {n} }$

(2)

where ${\textstyle \mathbf {e} _{3}}$ is the normal to the beam plane.

The position ${\textstyle \mathbf {x} }$ of an arbitrary point located on the beam cross section can be written in terms of the position of the central line ${\textstyle \mathbf {\varphi } }$ and the normal direction ${\textstyle \mathbf {n} }$ as a function of the distance ${\textstyle \zeta }$ of the point to the central axis as:

(3)

The undeformed stress-free (natural) configuration will be used as the reference configuration where the material points are defined by their coordinates ( ${\textstyle s}$ and ${\textstyle \zeta }$ ). Values measured on the reference configuration will be denoted with a left supra-index ${\textstyle o}$ .

Dropping in the notation the dependence with the arc length ${\textstyle s}$ , and denoting by ${\textstyle {\frac {\partial \left(\right)}{\partial s}}=\left(\right)_{^{\prime }s}}$ , the in-plane deformation gradient is defined by

(4)

While the right Cauchy-Green deformation tensor, restricted to the in-plane components is given by

\mathbf {C} =\mathbf {F} ^{T}\mathbf {F} =\left[{\begin{array}{c}[c]{cc}\mathbf {\varphi } _{^{\prime }s}\cdot \mathbf {\varphi } _{^{\prime }s}&0\\0&1\end{array}}\right]+\zeta \left[{\begin{array}{c}[c]{cc}2\mathbf {\varphi } _{^{\prime }s}\cdot \mathbf {n} _{^{\prime }s}&0\\0&0\end{array}}\right]+\zeta ^{2}\left[{\begin{array}{c}[c]{cc}\mathbf {n} _{^{\prime }s}\cdot \mathbf {n} _{^{\prime }s}&0\\0&0\end{array}}\right]

(5)

Disregarding terms associated with ${\textstyle \zeta ^{2}}$ the only non trivial component of ${\textstyle \mathbf {C} }$ is

(6)

We note that as the symmetric second order tensors considered are diagonal, only one index is used to alleviate notation when individual components are described.

The geometric measures of the central axis that allow to compute the strains are:

The stretch of the line of centroids (squared)

(7)

The in-plane curvature of the line of centroids

(8)

The later can also be expressed using the angle ${\textstyle \alpha }$ between the tangent ${\textstyle \mathbf {t} }$ and axis ${\textstyle \mathbf {e} _{1}}$ . Noting that

and that

then

(9)

Above expressions (starting with (3)) assume a non-deformable cross section, which is adequate for negligible elastic small strains. For large or moderately large strains associated with an elastic-plastic behavior, it is necessary to modify the previous expressions to account for the section stretching. Here it will be assumed that the section does not modify its shape but its size, in such a way that the stretching of the line of centroids does not change the volume. Then, if the same stretching is supposed in both directions of the cross section ${\textstyle \mathbf {n} }$ and ${\textstyle \mathbf {e} _{3}}$ , the normal stretching ${\textstyle \lambda _{n}}$ is

(10)

Note that this hypothesis does not imply any additional constraint on the beam theory neither an additional parameter. It is only introduced to update the distance of a point to the line of centroids (originally at a distant ${\textstyle \zeta }$ ) as ${\textstyle \zeta \lambda _{n}}$ . Consequently the present position of a point originally at ${\textstyle ^{o}\mathbf {x} \left(s,\zeta \right)}$ will be written as:

(11)

while the deformation gradient and the right Cauchy-Green deformation tensor are now

(12)

\mathbf {C} =\left[{\begin{array}{c}[c]{cc}\mathbf {\varphi } _{^{\prime }s}\cdot \mathbf {\varphi } _{^{\prime }s}&0\\0&\lambda _{n}^{2}\end{array}}\right]+2\lambda _{n}\zeta \left[{\begin{array}{c}[c]{cc}\mathbf {\varphi } _{^{\prime }s}\cdot \mathbf {n} _{^{\prime }s}&0\\0&0\end{array}}\right]

(13)

where the derivative of the normal stretching ${\textstyle \lambda _{n}}$ for the computation of the deformation gradient ${\textstyle \mathbf {F\,} }$ , and terms associated with ${\textstyle \zeta ^{2}}$ in the definition of ${\textstyle \mathbf {C} }$ , have been disregarded.

2.2 Extension to axisymmetric shells

Assuming the axis of revolution coincident with axis ${\textstyle \mathbf {e} _{2}}$ , the extension of above expressions to axisymmetric shells requires the definition of the stretch in the circumferential direction ${\textstyle \lambda _{\theta }}$ ( ${\textstyle \mathbf {e} _{\theta }=\mathbf {e} _{3}}$ ). This stretch depends only on the coordinate components over the axis ${\textstyle \mathbf {e} _{1}}$ in both present ( ${\textstyle x_{1}}$ ) and original ( ${\textstyle ^{o}x_{1}}$ ) configurations. The additional non-trivial components of the deformation gradient and the deformation tensor are:

(14)

(15)

Disregarding the influence of distance ${\textstyle \zeta }$ in the denominator and consistently the influence of terms in ${\textstyle \zeta ^{2}}$ in the numerator of Eq.(15) gives

$C_{\theta }={\frac {\left(\mathbf {\varphi } \cdot \mathbf {e} _{1}\right)^{2}+2\lambda _{n}\zeta \left(\mathbf {\varphi } \cdot \mathbf {e} _{1}\right)\left(\mathbf {n} \cdot \mathbf {e} _{1}\right)}{\left(^{o}\mathbf {\varphi } \cdot \mathbf {e} _{1}\right)^{2}}}$
$={\frac {\left(\mathbf {\varphi } \cdot \mathbf {e} _{1}\right)^{2}}{\left(^{o}\mathbf {\varphi } \cdot \mathbf {e} _{1}\right)^{2}}}+{\frac {2\lambda _{n}\zeta \left(\mathbf {\varphi } \cdot \mathbf {e} _{1}\right)\left(\mathbf {n} \cdot \mathbf {e} _{1}\right)}{\left(^{o}\mathbf {\varphi } \cdot \mathbf {e} _{1}\right)^{2}}}$	(16)
$=\left({\frac {\varphi _{1}}{^{o}\varphi _{1}}}\right)^{2}+{\frac {2\lambda _{n}\zeta \varphi _{1}\left(\mathbf {-} \sin \alpha \right)}{\left(^{o}\varphi _{1}\right)^{2}}}$	(17)

For axisymmetric shells the stretching in the normal direction ${\textstyle \mathbf {n} }$ , using the isochoric assumption mentioned above, can be computed by

(18)

The additional geometrical measures that allow to compute the strains are:

The middle surface stretching in the circumferential direction directly related with the Green strain

(19)

The second principal (hoop) curvature

(20)

In this way, the right Cauchy-Green deformation tensor can be written as (for the large strain case):

\mathbf {C} =\left[{\begin{array}{c}[c]{ccc}\mathbf {\varphi } _{^{\prime }s}\cdot \mathbf {\varphi } _{^{\prime }s}&0&0\\0&\lambda _{n}^{2}&0\\0&0&\left({\dfrac {\varphi _{1}}{^{o}\varphi _{1}}}\right)^{2}\end{array}}\right]+2\lambda _{n}\zeta \left[{\begin{array}{c}[c]{ccc}\mathbf {\varphi } _{^{\prime }s}\cdot \mathbf {n} _{^{\prime }s}&0&0\\0&0&0\\0&0&\lambda _{\theta }{\dfrac {\left(-\sin \alpha \right)}{^{o}\varphi _{1}}}\end{array}}\right]

(21)

2.3 Strain and stress measures

Note that for initially curved beams/shells, tensor ${\textstyle \mathbf {C} }$ (Eq. 21) is not the unit tensor, except for the points laying in the middle surface. Disregarding the influence of the initial curvatures in the definition of ${\textstyle ^{o}\mathbf {C} }$ and defining the curvature changes as:

(22)

it is possible to approximate (with the aim of an efficient numerical implementation) the expression of ${\textstyle \mathbf {C} }$ of Eq.(21) by

\mathbf {C} =\left[{\begin{array}{c}[c]{ccc}\mathbf {\varphi } _{^{\prime }s}\cdot \mathbf {\varphi } _{^{\prime }s}&0&0\\0&\lambda _{n}^{2}&0\\0&0&\left({\dfrac {\varphi _{1}}{^{o}\varphi _{1}}}\right)^{2}\end{array}}\right]+2\lambda _{n}\zeta \left[{\begin{array}{c}[c]{ccc}\mathbf {\varphi } _{^{\prime }s}\cdot \mathbf {n} _{^{\prime }s}-^{o}\mathbf {\varphi } _{^{\prime }s}\cdot ^{o}\mathbf {n} _{^{\prime }s}&0&0\\0&0&0\\0&0&{\dfrac {\sin ^{o}\alpha {-\lambda }_{\theta }\sin \alpha }{^{o}\varphi _{1}}}\end{array}}\right]

(23)

The Green-Lagrange strain measures can be directly obtained from ${\textstyle \mathbf {C} }$ by:

$E_{s}\left(\zeta \right)={\frac {1}{2}}\left(C_{s}-1\right)={\frac {1}{2}}\left(\lambda _{s}^{2}-1\right)+\zeta \lambda _{n}\chi _{s}={\hat {E}}_{s}+\zeta \lambda _{n}\chi _{s}$	(24.a)
$E_{\theta }\left(\zeta \right)={\frac {1}{2}}\left(C_{\theta }-1\right)={\frac {1}{2}}\left(\lambda _{\theta }^{2}-1\right)+\zeta \lambda _{n}\chi _{\theta }={\hat {E}}_{\theta }+\zeta \lambda _{n}\chi _{\theta }$	(24.b)

where the strains at points out of the middle surface can be obtained by adding the strains at the middle surface ${\textstyle {\hat {E}}_{s}}$ and ${\textstyle {\hat {E}}_{\theta }}$ and the strains due to the curvatures.

To deal with moderately large or large strains, the most convenient strain measure is the logarithmic strain. Note that tensor ${\textstyle \mathbf {C} }$ is diagonal so that the principal strain directions are fixed (in material axis). Hence, the Hencky strain tensor components for an arbitrary point with coordinate ${\textstyle \zeta }$ along the normal ${\textstyle \mathbf {n} }$ are:

$\varepsilon _{s}={\frac {1}{2}}\ln C_{s}={\frac {1}{2}}\ln \left(\lambda _{s}^{2}+2\zeta \lambda _{n}\chi _{s}\right)$	(25.a)
$\varepsilon _{\theta }={\frac {1}{2}}\ln C_{\theta }={\frac {1}{2}}\ln \left(\lambda _{\theta }^{2}+2\zeta \lambda _{n}\chi _{\theta }\right)$	(25.b)

These strain measures can be associated with the (convective) Kirchhoff stress tensor ${\textstyle \mathbf {T} }$ through a constitutive equation adequate for large strains. The usual thin shell plane stress assumptions are used. In general, it is more convenient to write the weak form of the momentum equations in terms of the second Piola-Kirchhoff stress tensor ${\textstyle \mathbf {S} }$ related to the reference configuration. The relation between both stress tensors is:

$S_{s}={\frac {T_{s}}{C_{s}}}$	(26.a)
$S_{\theta }={\frac {T_{\theta }}{C_{\theta }}}$	(26.b)

The through-the-thickness integrated forces (generalized or resultant stresses) are:

\left[{\begin{array}{c}[c]{c}N_{s}\\N_{\theta }\\M_{s}\\M_{\theta }\end{array}}\right]=b\int _{-{\dfrac {h}{2}}}^{+{\dfrac {h}{2}}}\left[{\begin{array}{c}[c]{c}S_{s}\\S_{\theta }\\S_{s}\lambda _{n}\zeta \\S_{\theta }\lambda _{n}\zeta \end{array}}\right]d\zeta

(27)

where ${\textstyle h}$ is the original thickness of the beam/shell, ${\textstyle b}$ is the original width ( ${\textstyle 2\pi ^{o}\varphi _{1}}$ for axisymmetric shells) and ${\textstyle \lambda _{n}\zeta }$ is the present distance (deformed) of the point to the middle line/surface. The weak form of the equilibrium conditions can be written as a virtual work equation. The internal virtual work (IVW) can be reduced to

(28)

and

{\hbox{IVW}}=\int _{^{o}V}{\frac {1}{2}}\delta \mathbf {C} :\mathbf {S} d^{o}V=2\pi \int _{^{o}L}\left[\lambda _{s}\delta \lambda _{s},\lambda _{\theta }\delta \lambda _{\theta },\delta \kappa _{s},\delta \kappa _{\theta }\right]\left[{\begin{array}{c}[c]{c}N_{s}\\N_{\theta }\\M_{s}\\M_{\theta }\end{array}}\right]\varphi _{1}ds

(29)

for beams and axisymmetric shells, respectively.

3 Rotation-free beam/axisymmetric shell element

A two-node element with only translational degrees of freedom (rotation-free) is proposed. The geometrical configuration is defined by the positions of the two end nodes ${\textstyle \mathbf {\varphi } ^{1}}$ and ${\textstyle \mathbf {\varphi } ^{2}}$ (see Figure 2). For each element a linear interpolation for the line of centroids is proposed

(30)

where ${\textstyle N^{I}\left(\xi \right)}$ are the standard linear Lagrangian polynomials (shape functions) defined over the interval ${\textstyle \left[0,1\right]}$ . This allows the evaluation of middle surface stretches and the hoop curvature but it is not enough for the computation of the in-plane curvature ${\textstyle \kappa _{s}}$ . The later is computed resorting to the geometry of the surrounding elements leading to an assumed strain approach as shown below.

Figure 2: Two-node element.

3.1 Axial strain

The constant axial strain is easily computed from the element configuration. The tangent vector in the original configuration is:

(31)

where ${\textstyle ^{o}\mathbf {\varphi } ^{i}}$ , ${\textstyle i=1,2}$ are the original position vectors of the element nodes and ${\textstyle ^{o}L}$ the initial element length given by

(32)

Performing the same operations for the deformed configuration, we obtain

(33)

(34)

The axial stretch is

(35)

and the Green-Lagrange strain of the line of centroids is computed as:

(36)

while its variation is:

$\delta {\hat {E}}_{s}=\lambda _{s}\delta \lambda _{s}={\frac {L}{^{o}L}}{\frac {\delta L}{^{o}L}}$	(37.a)
$={\frac {\left(\mathbf {\varphi } ^{2}-\mathbf {\varphi } ^{1}\right)\cdot \left(\delta \mathbf {\varphi } ^{2}-\delta \mathbf {\varphi } ^{1}\right)}{\left(^{o}L\right)^{2}}}$	(37.b)
$=\lambda _{s}{\frac {\mathbf {t} \cdot \left(\delta \mathbf {\varphi } ^{2}-\delta \mathbf {\varphi } ^{1}\right)}{^{o}L}}$	(37.c)

Matrix ${\textstyle \mathbf {B} _{ms}}$ relating axial strain increment with nodal incremental displacements can be expressed as:

\delta {\hat {E}}_{s}=\lambda _{s}\delta \lambda _{s}={\frac {\lambda _{s}}{^{o}L}}\left[-\mathbf {t} ^{T},\mathbf {t} ^{T}\right]\left[{\begin{array}{c}[c]{c}\delta \mathbf {u} ^{1}\\\delta \mathbf {u} ^{2}\end{array}}\right]={\underset {1\times {4}}{\mathbf {B} _{ms}}}{\underset {4\times {1}}{\delta \mathbf {u} ^{e}}}

(38)

where the incremental displacements at the nodes have been grouped in vector u ${\textstyle ^{e}}$ .

3.2 Hoop strains in axisymmetric shells

The computation of the Green hoop strain ${\textstyle {\hat {E}}_{\theta }}$ at the middle surface (2.3.b) and the hoop change of curvature ${\textstyle \chi _{\theta }}$ (22) associated to axisymmetric shells does not offer difficulties. Denoting with an upper index the node and with a lower index the component of the Cartesian system, in both the configuration ${\textstyle \mathbf {\varphi } }$ and the displacement ${\textstyle \mathbf {u} }$ :

{\hat {E}}_{\theta }={\frac {1}{2}}\left[\left(\lambda _{\theta }\right)^{2}-1\right]={\frac {1}{2}}\left\{{\frac {\left(N^{1}\varphi _{1}^{1}+N^{2}\varphi _{1}^{2}\right)^{2}}{\left(N^{1o}\varphi _{1}^{1}+N^{2o}\varphi _{1}^{2}\right)^{2}}}-1\right\}

(39)

(40)

The variations of these strain measures, to be used in the weak form of the momentum equations, are:

$\delta {\hat {E}}_{\theta }=\lambda _{\theta }\delta \lambda _{\theta }=\lambda _{\theta }{\frac {N^{1}\delta u_{1}^{1}+N^{2}\delta u_{1}^{2}}{N^{1o}\varphi _{1}^{1}+N^{2o}\varphi _{1}^{2}}}$	(41)
$={\frac {\lambda _{\theta }}{N^{1o}\varphi _{1}^{1}+N^{2o}\varphi _{1}^{2}}}\left[N^{1},0,N^{2},0\right]\left[{\begin{array}{c}\delta u_{1}^{1}\\\delta u_{2}^{1}\\\delta u_{1}^{2}\\\delta u_{2}^{2}\end{array}}\right]={\underset {1\times {4}}{\mathbf {B} _{m\theta }}}{\underset {4\times 1}{\delta \mathbf {u} ^{e}}}$	(42)

\delta \chi _{\theta }={\frac {\left(N^{1}\delta u_{1}^{1}+N^{2}\delta u_{1}^{2}\right)\left(-\sin \alpha \right)-\left(N^{1}\varphi _{1}^{1}+N^{2}\varphi _{1}^{2}\right)\cos \alpha \left(\mathbf {n} \cdot \delta \mathbf {t} \right)}{\left(N^{1o}\varphi _{1}^{1}+N^{2o}\varphi _{1}^{2}\right)^{2}}}

=-{\frac {\sin \alpha \left(N^{1}\delta u_{1}^{1}+N^{2}\delta u_{1}^{2}\right)+\left(N^{1}\varphi _{1}^{1}+N^{2}\varphi _{1}^{2}\right)\cos \alpha \left(\mathbf {n} \cdot {\dfrac {\delta \mathbf {\varphi } ^{2}-\delta \mathbf {\varphi } ^{1}}{L}}\right)}{\left(N^{1o}\varphi _{1}^{1}+N^{2o}\varphi _{1}^{2}\right)^{2}}}={\underset {1\times {4}}{\mathbf {B} _{b\theta }}}{\underset {4\times {1}}{\delta \mathbf {u} ^{e}}}

(43)

3.3 Evaluation of the in-plane curvature

For the in-plane curvature an assumed strain approach is adopted. The curvatures (9) are computed at the two nodes and a linear interpolation of these values is considered within the element

${\bar {\kappa }}_{s}\left(s\right)=\left(1-{\frac {s}{^{o}L}}\right)\kappa _{s}^{\left(1\right)}+{\frac {s}{^{o}L}}\kappa _{s}^{\left(2\right)}$	(44)
$=-N^{1}\lambda _{s}^{\left(1\right)}\alpha _{^{\prime }s}^{\left(1\right)}-N^{2}\lambda _{s}^{\left(2\right)}\alpha _{^{\prime }s}^{\left(2\right)}$	(45)

where a supra-index between brackets indicates the node of evaluation.

As in other rotation-free element, for the computation of the curvatures, the position of the nodes in the adjacent element is considered. In this case a patch of three elements (see Figure 3) is involved. There are different options for the computation of the curvatures at each element node. Initially an approximation adequate for smooth and homogeneous beam is defined. This approximation is latter modified to deal with more general geometries, including non-smooth and branching lines.

Figure 3: Patch of elements used for the evaluation of the curvature.

For the initial approximation of (9) the following (finite differences) computations can be used at each node:

a) the axial stretch ${\textstyle \lambda _{s}}$ at the node is computed as the average of the axial stretches at the adjacent elements

b) the derivative of the angle ${\textstyle \alpha }$ is computed as the difference of the angles between the tangent vectors ${\textstyle \mathbf {t} }$ divided by the average of the initial element lengths ${\textstyle ^{o}L}$

then

$\kappa _{s}^{\left(1\right)}={\frac {\lambda _{sl}+\lambda _{s}}{2}}{\frac {2\phi ^{1}}{^{o}L_{l}+^{o}L}}={\frac {\lambda _{sl}+\lambda _{s}}{^{o}L_{l}+^{o}L}}\phi ^{1}$	(46.a)
$\kappa _{s}^{\left(2\right)}={\frac {\lambda _{s}+\lambda _{sr}}{^{o}L+^{o}L_{r}}}\phi ^{2}$	(46.b)

In Eqs.(46) indexes ${\textstyle l}$ and ${\textstyle r}$ denote the left and right end nodes in the three elements patch (see Figure 3). Also ${\textstyle ^{o}L_{r}}$ and ${\textstyle ^{o}L_{l}}$ are the initial lengths of the elements respectively located at the right and the left of the central element in the patch.

For future reference the following angles must be identified

${\textstyle \alpha }$ is, as defined above, the angle between the element (vector ${\textstyle \mathbf {t} }$ of Equation (33)) and the global axe ${\textstyle \mathbf {e} _{1}}$

${\textstyle \beta }$ is the rotation of the element with respect to the original configuration (angle between vector ${\textstyle \mathbf {t} }$ and vector ${\textstyle ^{o}\mathbf {t} }$ of Equation (31))

${\textstyle \phi }$ is the angle between two adjacent elements (relative angle). For example for the first node of the element (see Figure 3)

For the initially horizontal elements of Figure 3 the angle ${\textstyle \alpha }$ is also the rotation ${\textstyle \beta }$ of the element. While the relative angle ${\textstyle \phi ^{1}}$ between the tangent vector of the central element ( ${\textstyle \mathbf {t} }$ ) and the tangent vector of the left element ( ${\textstyle \mathbf {t} _{l}}$ ) measured from ${\textstyle \mathbf {t} }$ to ${\textstyle \mathbf {t} _{l}}$ counter-clockwise (see Figure 4) is also the relative rotation ${\textstyle \beta {-\beta }_{l}}$ . Similarly ${\textstyle \phi ^{2}}$ is the angle between the tangent vector of the central element ( ${\textstyle \mathbf {t} }$ ) and the tangent vector of the right element ( ${\textstyle \mathbf {t} _{r}}$ ) measured from ${\textstyle \mathbf {t} _{r}}$ to ${\textstyle \mathbf {t} }$ counter-clockwise.

Figure 4: Change of angle between elements.

The approximation above leading to expressions (46) define a unique curvature value at each node, independently of the element from where it is computed. This represents a drawback for nonhomogeneous beams/shells, i.e. when there are changes in the material or in the transverse section properties of the beam/shell. Then at the intersection node ${\textstyle n}$ of the two elements ${\textstyle i}$ and ${\textstyle j}$ of Figure 4 with different bending stiffness ( ${\textstyle EI}$ ), equilibrium can not be fulfilled at the node as

(47)

because the nodal curvatures ${\textstyle \kappa _{si}^{\left(n\right)}}$ and ${\textstyle \kappa _{sj}^{\left(n\right)}}$ are equal.

To avoid this problem, at each node ${\textstyle n}$ shared by two elements ( ${\textstyle i}$ and ${\textstyle j}$ ), the nodal curvatures are re-defined independently for each element as:

(48)

(49)

(Note that lower indexes refer to element variables and upper indexes to nodal variables) where, with reference to Figure 4, ${\textstyle \gamma _{k}^{\left(n\right)}}$ is twice the difference between the rotation of the element ${\textstyle \beta _{k}}$ and the rotation of the node ${\textstyle {\bar {\beta }}^{n}}$ (undefined yet)

$\gamma _{i}^{\left(n\right)}=2\left(\beta _{i}-{\bar {\beta }}^{n}\right)$	(50.a)
$\gamma _{j}^{\left(n\right)}=2\left({\bar {\beta }}^{n}-\beta _{j}\right)$	(50.b)

The rotation of the node ${\textstyle {\bar {\beta }}^{n}}$ can be seen as the new position of a corrotational system from which curvatures are computed using the relative rotation ${\textstyle \gamma _{k}^{\left(n\right)}}$ . The sum of the relative rotations is of course the angle between both elements

(51)

where the relative angle at node ${\textstyle \phi ^{n}}$ was defined above (see Figure 3).

Replacing the new definitions Eqs.(51) and (52) in the nodal equilibrium condition (50), for small strains ( ${\textstyle \lambda _{si}\simeq \lambda _{sj}\simeq {1}}$ ), it results:

$\left({\frac {EI}{^{o}L}}\right)_{i}\gamma _{i}^{\left(n\right)}=\left({\frac {EI}{^{o}L}}\right)_{j}\gamma _{j}^{\left(n\right)}$	(52)
$R_{i}\gamma _{i}^{\left(n\right)}=R_{j}\gamma _{j}^{\left(n\right)}{\hbox{ with}}R_{k}=\left({\frac {EI}{^{o}L}}\right)_{k}\quad k=i,j$	(53)

Replacing Eq.(51) in the equilibrium relation (53) the values of ${\textstyle \gamma _{k}^{\left(n\right)}}$ can be expressed in terms of the relative angle ${\textstyle \phi ^{n}}$ and the stiffness measures ( ${\textstyle R_{k}}$ )

$\gamma _{i}^{\left(n\right)}={\frac {2R_{j}}{R_{i}+R_{j}}}\left(\beta _{i}-\beta _{j}\right)=c_{j}\phi ^{n}$	(54.a)
$\gamma _{j}^{\left(n\right)}={\frac {2R_{i}}{R_{i}+R_{j}}}\left(\beta _{i}-\beta _{j}\right)=c_{i}\phi ^{n}$	(54.b)

Finally, noting that the difference between the relative rotations is

(55)

replacing Eq.(54) into (55) it allows to define the rotation of the node ${\textstyle {\bar {\beta }}^{n}}$ as:

(56)

Then, the rotation of the node ${\textstyle {\bar {\beta }}^{n}}$ can be computed as the average of the rotation of the two adjacent elements weighted by their respective stiffness measure ( ${\textstyle R_{k}}$ ). Above expressions (54) allow to compute the curvature at any point of the element by Eq.(45).

The variation of the curvature at an arbitrary point can be expressed as:

(57)

From Eq. (49) the variation of the nodal curvature is composed of two parts

(58)

The variations ${\textstyle \delta \lambda _{s}}$ can be obtained as in (37c)

(59)

While the variation ${\textstyle \delta \gamma ^{\left(n\right)}}$ can be obtained from (54) as (see Appendix):

$\delta \gamma ^{\left(1\right)}={\frac {2R_{l}}{R+R_{l}}}\delta \left(\beta _{l}-\beta \right)=c_{l}\delta \phi ^{1}=c_{l}\left[{\frac {\delta \mathbf {\varphi } ^{1}-\delta \mathbf {\varphi } ^{2}}{L}}\cdot \mathbf {n} +{\frac {\delta \mathbf {\varphi } ^{1}-\delta \mathbf {\varphi } ^{l}}{L}}\cdot \mathbf {n} _{l}\right]$	(60.a)
$\delta \gamma ^{\left(2\right)}={\frac {2R_{r}}{R+R_{r}}}\delta \left(\beta {-\beta }_{r}\right)=c_{r}\delta \phi ^{2}=c_{r}\left[{\frac {\delta \mathbf {\varphi } ^{2}-\delta \mathbf {\varphi } ^{1}}{L}}\cdot \mathbf {n} +{\frac {\delta \mathbf {\varphi } ^{r}-\delta \mathbf {\varphi } ^{2}}{L_{r}}}\cdot \mathbf {n} _{r}\right]$	(60.b)

Substituting Eqs. (59) and (60) into (58) gives

$\delta \kappa _{s}^{\left(1\right)}={\frac {1}{^{o}L}}\left[\gamma ^{\left(1\right)}{\frac {\delta \mathbf {\varphi } ^{2}-\delta \mathbf {\varphi } ^{1}}{^{o}L}}\cdot \mathbf {t} +\lambda _{s}c_{l}\left({\frac {\delta \mathbf {\varphi } ^{1}-\delta \mathbf {\varphi } ^{l}}{L_{l}}}\cdot \mathbf {n} _{l}+{\frac {\delta \mathbf {\varphi } ^{1}-\delta \mathbf {\varphi } ^{2}}{L}}\cdot \mathbf {n} \right)\right]$	(61)
$={\frac {1}{^{o}L}}\left\{-{\frac {\lambda _{s}c_{l}}{L_{l}}}\mathbf {n} _{l},\left[\lambda _{s}c_{l}\left({\frac {\mathbf {n} _{l}}{L_{l}}}+{\frac {\mathbf {n} }{L}}\right)-{\frac {\gamma ^{\left(1\right)}}{^{o}L}}\mathbf {t} \right],\left[{\frac {\gamma ^{\left(1\right)}}{^{o}L}}\mathbf {t} -{\frac {\lambda _{s}c_{l}}{L}}\mathbf {n} \right],0\right\}\left[{\begin{array}{c}[c]{c}\delta \mathbf {u} ^{l}\\\delta \mathbf {u} ^{1}\\\delta \mathbf {u} ^{2}\\\delta \mathbf {u} ^{r}\end{array}}\right]$	(62)
$=\mathbf {B} _{bs}^{\left(1\right)}\delta \mathbf {u} ^{p}$	(63)

where ${\textstyle \delta \mathbf {u} ^{p}}$ gathers the displacements of all the nodes in the patch (although node ${\textstyle r}$ have no contributions in the curvature at node ${\textstyle 1}$ ).

The variation of (48) is similarly obtained as:

$\delta \kappa _{s}^{\left(2\right)}={\frac {1}{^{o}L}}\left[\gamma ^{\left(2\right)}{\frac {\delta \mathbf {\varphi } ^{2}-\delta \mathbf {\varphi } ^{1}}{^{o}L}}\cdot \mathbf {t} +\lambda _{s}c_{r}\left({\frac {\delta \mathbf {\varphi } ^{2}-\delta \mathbf {\varphi } ^{1}}{L}}\cdot \mathbf {n} +{\frac {\delta \mathbf {\varphi } ^{2}-\delta \mathbf {\varphi } ^{r}}{L_{r}}}\cdot \mathbf {n} _{r}\right)\right]$	(64)
$={\frac {1}{^{o}L}}\left\{0,\left[{\frac {\gamma ^{2}}{^{o}L}}\mathbf {t} -{\frac {\lambda _{s}c_{r}}{L}}\mathbf {n} \right],\left[\lambda _{s}c_{r}\left({\frac {\mathbf {n} _{r}}{L_{r}}}+{\frac {\mathbf {n} }{L}}\right)-{\frac {\gamma ^{\left(2\right)}}{^{o}L}}\mathbf {t} \right],-{\frac {\lambda _{s}c_{r}}{L_{r}}}\mathbf {n} _{r}\right\}\left[{\begin{array}{c}[c]{c}\delta \mathbf {u} ^{l}\\\delta \mathbf {u} ^{1}\\\delta \mathbf {u} ^{2}\\\delta \mathbf {u} ^{r}\end{array}}\right]$	(65)
$=\mathbf {B} _{bs}^{\left(2\right)}\delta \mathbf {u} ^{p}$	(66)

Finally (57) can be written as:

$\delta {\bar {\kappa }}_{s}=N^{1}\delta \kappa _{s}^{\left(1\right)}+N^{2}\delta \kappa _{s}^{\left(2\right)}=\left[N^{1}\mathbf {B} _{bs}^{\left(1\right)}+N^{2}\mathbf {B} _{bs}^{\left(2\right)}\right]\delta \mathbf {u} ^{p}$	(67)
$={\underset {1\times {8}}{\mathbf {B} _{bs}\left(s\right)}}{\underset {8\times {1}}{\delta \mathbf {u} ^{p}}}$	(68)

Note that matrix ${\textstyle \mathbf {B} _{b}}$ is linear in the coordinate ${\textstyle s}$ , then two integrations points are necessary to correctly integrate the tangent stiffness matrix or the residual force vector. If one integration point is desired some stabilization procedure is needed.

3.4 Boundary conditions

The imposition of translation boundary conditions associated with the line of centroids is straightforward as these are the element degrees of freedom. For boundary conditions affecting the rotation of the element normal, it must be noted that in such points the neighbor element does not exist. Two cases may be considered:

Clamp or symmetry condition: in this case the rotation is null, i.e. the tangent to the beam/shell at that point is constant along the deformation, and equal to the normal to the symmetry (clamping) plane ${\textstyle ^{o}\mathbf {t} }$ . The relative angle ${\textstyle \phi }$ (in this case when the symmetric point is missing) can be expressed as twice the angle change between ${\textstyle ^{o}\mathbf {t} }$ and the element tangent. If, for example, this restriction occurs in the first node of the element, the curvature at this node can be computed as:

also, noting that ${\textstyle ^{o}\mathbf {t} }$ is constant, the curvature variation results:

Free or simple supported edge: this is a natural boundary condition. The most straightforward approximation is to assume a null curvature at this point.

4 Curved and non-smooth beams/shells

In the above developments it is assumed a straight beam/shell, thus, the original angle ${\textstyle \phi }$ between two elements is zero in the reference configuration. This is implicit in the way the in-plane curvatures where computed at each element node (3.3).

Consider now two elements that in the original configuration are not aligned. Assume that for the reference element the node under consideration is the first node. Using only the nodes position, it is not possible from the tangent vectors ${\textstyle \mathbf {t} }$ (or the normals ${\textstyle \mathbf {n} }$ ), to classify the node neighborhood as a curved line or a non-smooth line (a kink). This information must be explicitly supplied. Anyway, the handling of both cases is the same from practical purposes, because the effects due to the original curvatures are disregarded in the present formulation (see expression (23)).

Figure 5: Initial configuration for a non-smooth surface.

Consider the configuration shown in Figure 5 where there exists an angle ${\textstyle ^{o}\phi ^{1}}$ between the original tangents at each element defined by:

$\cos ^{o}\phi ^{1}=^{o}\mathbf {t} \cdot ^{o}\mathbf {t} _{i}=^{o}\mathbf {n} \cdot ^{o}\mathbf {n} _{l}$	(69.a)
$\sin ^{o}\phi ^{1}=-^{o}\mathbf {t} \cdot ^{o}\mathbf {n} _{l}=^{o}\mathbf {n} \cdot ^{o}\mathbf {t} _{l}$	(69.b)

In the deformed configurations the tangent vectors ${\textstyle \mathbf {t} }$ and ${\textstyle \mathbf {t} _{l}}$ define a different angle ${\textstyle \phi ^{1}}$ ,

$\cos \phi ^{1}=\mathbf {t} \cdot \mathbf {t} _{l}$	(70.a)
$\sin \phi ^{1}=\mathbf {n} \cdot \mathbf {t} _{l}$	(70.b)

that represents a change with respect to the original configuration:

(71)

The angle change so defined allows to compute the curvature changes at each of the two adjacent elements applying the same expressions derived in a previous section, i.e.: use Eqs. (54) to compute the relative rotation ${\textstyle \gamma ^{\left(n\right)}}$

$\gamma ^{\left(1\right)}=c_{l}\Delta \phi ^{1}$	(72.a)
$\gamma ^{\left(2\right)}=c_{r}\Delta \phi ^{2}$	(72.b)

and then use (48-49) to evaluate the curvature change at each element node

(73)

5 More than two elements sharing a node (branching)

In a general case there can be ${\textstyle n}$ elements sharing a node. The procedure for computing the in-plane curvature in terms of the nodal coordinates of the elements adjacents to the branching node is described below. Indeed the formulation of just two elements sharing a node described in the previous section is a particular case of the general approach here presented next.

Figure 6: Branching in beams.

To alleviate the notation, it will be assumed that the intersection node (denoted by ${\textstyle 0}$ ) is the first node of each element, while the second node of each beam ${\textstyle i}$ will also be defined (arbitrarily) as node ${\textstyle i}$ (see Figure 6). In the original configuration, each element ( ${\textstyle i}$ ) configuration is completely defined by its tangent vector:

(74)

As usual the element forms with the Cartesian axis ${\textstyle \mathbf {e} _{1}}$ an angle ${\textstyle ^{o}\alpha _{i}}$ .

In the deformed configuration, the direction of element ${\textstyle i}$ is computed as

(75)

defining with the Cartesian axis ${\textstyle x_{1}}$ an angle ${\textstyle \alpha _{i}}$ . The angle rotated by each element with respect to the original configuration is:

(76)

Using the nodal equilibrium as above the rotation of the node is defined as the weighted average:

(77)

Denoting again by ${\textstyle \gamma _{i}^{\left(0\right)}}$ twice the difference between the rotation of the node ${\textstyle 0}$ and the rotation of the element ${\textstyle i}$

(78)

it is possible to write the in-plane change of curvature at the first node ( ${\textstyle 0}$ ) of each element as:

(79)

The variation of this in-plane curvature at the first node of each element gives:

(80)

where for each element sharing the node

(81)

and with this

(82)

finally

+2{\frac {\lambda _{si}}{^{o}L_{i}}}\left[-{\frac {1}{L_{i}}}\mathbf {n} _{i}\cdot \left(\delta \mathbf {\varphi } ^{i}-\delta \mathbf {\varphi } ^{0}\right)+{\frac {1}{\sum _{k=1}^{n}R_{k}}}\sum _{j=1}^{n}{\frac {R_{j}}{L_{j}}}\mathbf {n} _{j}\cdot \left(\delta \mathbf {\varphi } ^{j}-\delta \mathbf {\varphi } ^{0}\right)\right]

(83)

This allows to evaluate the variations of the in-plane curvature at the different elements sharing a branching node in order to compute the weak form of the momentum equations. It must be noted that the curvature variation depends on the nodal displacement variations of the nodes of all the elements sharing the branching node. Thus, the stiffness matrix of each element depends on all the nodes of the patch.

Note that above derivations hold both for beam and axisymmetric shell elements. The remaining of the strains in both elements are computed in the standard manner.

6 Numerical examples

In this section several examples are presented to assess the element described. This will be denoted as B2 (2-node Beam). Initially three linear examples are considered, then three geometrically non-linear (elastic) examples are shown, and finally a couple of examples including large strain plasticity (in smooth surfaces) are studied. All units are in the International System. The present element performance is compared with two transverse shear deformable two-node elements with 3 nodal degrees of freedom (two displacements and one rotation) and one integration point. Namely element B21 present in program ABAQUS [20], denoted here by Abq, and the element of Ref [17] denoted here by Sh, that should give identical results at least in the small strain elastic range. Note in the comparisons that using the shear deformable elements implies 50% more DOFs than when using element B2 for the same mesh.

6.1 Linear examples

The examples in the linear range are studied to assess the convergence properties of the B2 element. Comparisons with exact and converged numerical solutions are presented.

6.1.1 Simple supported beam under uniform load

The first example is a very simple one. A simple supported straight beam simply supported at both ends under a uniform transverse load in all its length. Due to symmetry only the left half of the beam is discretized. This allows to test both boundary conditions defined before for the normal rotation. Figure 7.a depicts the whole geometry.

(a)

Figure 7: Simple supported beam under uniform load. (a) geometry; (b) displacements; (c) bending moment.

Different meshes have been considered to assess convergence, namely with 1, 2 and 4 elements. Nodes are uniformly distributed in the beam. In Figure 7.c the normalized displacement (respect to the exact maximum) has been plotted at the center of each element for the three meshes used, while Figure 7.d shows the normalized bending moment at the same points. Results obtained with element Abq are also shown for comparison.

It can be seen that when the mesh is refined convergence to the correct results is obtained. Four elements (8 for the whole beam) are enough to get accurate enough results. Convergence is faster for moments than for displacements, due probably to the fact that the element is non-conforming. For values computed at the element centers better results are obtained with element B2 than with the shear deformable elements using the same meshes.

6.1.2 Simple rigid-joined frame

The second examples considered is a two-beam frame with a sharp change of angle between two adjacent elements (see Figure 8.a). The horizontal member has a larger moment of inertia than the vertical member. An horizontal point load is applied as shown in the figure.


(a)	(b)

(c)
Figure 8: Simple frame under point load. (a) geometry (b) horizontal displacements on the vertical member. (c) Deformed configurations in the non-linear range.

Figure 8.b shows the horizontal displacement over the vertical beam for the different meshes considered of 3, 6, and 12 elements discretizing each member. Results are normalized with respect to the exact displacement at joint A. Results obtained with the shear deformable element Sh with the same number of DOFs (2, 4 and 8 elements respectively) has also been plotted for reference. In this case it can be seen again that six elements at each member are needed to obtain accurate enough results from the engineering point of view, while four shear deformable elements would suffice for this case with a simple linear variation of the bending moment.

6.1.3 Branching axisymmetric shell

The third example considers a branching shell to assess the element formulation when more than two elements are joined at a node. Figure 9.a shows the geometry of the problem. The three parts of the shell have different thicknesses. The material is linear, elastic and isotropic with ${\textstyle E=10^{7}}$ and ${\textstyle \nu =0.3}$ . The spherical shell and the lower part of the cylinder are both subjected to internal pressure ${\textstyle P=1000}$ . This pressure is globally equilibrated by equal boundary forces at the ends of the cylinder.

Two meshes where considered in the discretization, one with 42 elements(18 for the lower cylinder, 16 for the sphere and 8 for the upper cylinder) and a second mesh with half the number of elements. The normal displacement along the cylinder has been plotted in Figure 9.b. Results obtained with the shear deformable element using the mesh of 42 elements and a converged finite element solution [18] are included for comparison (an analytical solution is also possible). For the finer mesh considered the results presented show an excellent agreement.


(a)	(b)

(c)
Figure 9: Krauss' branching shell, (a) geometry $R=20$ , $L_{1}=20$ , $L_{2}=10$ , $h_{1}=0.3$ , $h_{2}=0.4$ , and $h_{3}=0.5$ (b) normal displacement along the cylinder. (c) Normal displacements at points 1 and 2 versus the internal pressure.

6.2 Examples including geometric non-linearity

In this section examples including large displacements and rotations are considered. Results obtained with present element are compared with results obtained with the shear deformable element. Two examples from the previous section, including sharp angles and branching, are now studied in the non-linear range. Also a well studied benchmark for non-smooth shells is analyzed to assess the convergence properties in the non-linear range.

6.2.1 Simple rigid-joined frame

The first example in this section is the simple two-bar frame considered before. A rather fine mesh of 12 elements on each bar has been considered. Figure 8.c shows the displacement profiles of the vertical bar for different values of the load ${\textstyle P}$ . For the mesh considered the results are quite similar with those obtained using the shear deformable element (mesh with 8 elements in each bar and the same number of DOFs).

6.2.2 Branching axisymmetric shell

The second non-linear example considered includes a branching. This allows to assess convergence when more than two elements share a node in the non-linear range. Figure 9.c plots the normal displacements of points denoted as “1” and “2” in Figure 9.a as the internal pressure is increased.

The same finer mesh considered in the linear case (42 elements) was used. Results are compared with those obtained with the shear deformable element for the same mesh. The normal displacement of node 1 (branching node) is almost linear with the pressure and both elements give indistinguishable results. The displacements of point 2 (upper end of the cylinder) show some differences probably due the larger flexibility of the shear deformable element in this part of the shell where the thickness is maximum.

6.2.3 Elastic large deflection response of a Z-shaped cantilever

This third example was extracted from a set of benchmarks for geometric non-linear behavior of structures [19]. This benchmark is intended to assess membrane and bending actions, tension stiffening and change in sign of bending moments under large rotations and large displacements. Figure 10.a shows the cross section properties, the material parameters and the initial geometry of a Z-shaped cantilever under a conservative load at the free end. Figure 10.b shows the deformed configurations for three different values of the point load. These load values are defined as target values in Ref. [19] to assess the element formulation.

In Figure 10.c and d the results for meshes with 3, 6 and 12 elements per zone are shown along all the process. In Figure 10.c the load versus the tip displacement has been plotted while the bending moment at A versus the load is shown in Figure 10.d. The bending moment plotted is obtained as the average of the bending moments computed at the integration points adjacent to point A. The results obtained with the B2 element show a very good agreement with the three expected values for the mesh with 6 elements per zone.




Figure 10: Z-shaped cantilever. (a) Geometry and material; (b) Normalized displacements and moments as the mesh is refined; (c) Load versus tip displacement; (d) Bending moment at A versus load (e) and (f) convergence of the results in the non-linear range as the mesh is refined for P = 104.5 and $P=1263$ respectively.

To assess convergence as the mesh is refined in the non-linear range Figures 10.e and f plots for ${\textstyle P=104.5}$ and ${\textstyle P=1263}$ the normalized tip displacement and normalized moment at point A as a function of the number of DOFs in the mesh. The target values from Ref. [19] where used for normalization. It can be seen that the element has good convergence properties for both displacements and moments, while the element B21 of program ABAQUS [20] shows excellent convergence properties for the displacement but a very slow convergence for the moments.

6.3 Examples including large elastic-plastic strains

In this section a couple of examples including large elastic-plastic strains are studied. In contrast with previous examples, the initial geometry of both cases is smooth. These simulations are performed with the explicitly dynamic program STAMPACK [21], oriented to the simulation of sheet metal forming processes.

6.3.1 Stretching of a circular sheet with an hemispherical punch

First a well studied benchmark is considered [22]. Figure 11.a shows the initial configuration of the stretching of circular plate with an hemispherical punch. The material of the blank is a mild steel with ${\textstyle E=6.9\times {10}^{10}}$ , ${\textstyle \nu =0.3}$ and tensile strength ${\textstyle \sigma _{y}=5.89\times {10}^{8}(0.0001+e^{p})^{0.216}}$ .

The tools are assumed to be rigid and a frictionless interface was used. Contact between the tools and the sheet was considered via a penalty formulation. Two different meshes of 28 and 56 elements where used in the simulations with no practical differences in the results. Punch force versus punch travel is plotted in Figure 11.b. Results obtained with the B2 element have been compared with those obtained using two solid element: a non-conforming triangle [23] (TR) and a standard quadrilateral with constant pressure [24] (QUAD). The results show a excellent agreement.

Figures 11.c and 11.d show the profiles of thickness ratio and equivalent plastic strain along the original radius for different stages of the simulation. These values compare quite well with the results obtained with the solid elements (not plotted for clarity).



Figure 11: Punch stretching simulation of a thin metal sheet. (a)geometry; (b)punch force versus punch travel; (c)thickness ratio profiles as a function of the radius for different punch travels. (d)effective plastic strain (thickness averaged) profiles as a function of the radius for different punch travels.

6.3.2 Deep drawing and springback of a flat strip

The last example is a benchmark of NUMISHEET'93 [25]. The aim of this example is to assess the forming process (deep drawing) of a thin sheet and the effects of the elastic springback after the tools are removed. Figure 12.a shows the setting of the experiment. The geometrical parameters to be computed and compared with experimental results are shown in Figure 12.b. Namely the angles ${\textstyle \theta _{1}}$ and ${\textstyle \theta _{2}\,}$ and the radius of curvature ${\textstyle \rho }$ of the lateral wall, obtained through the circle defined by points A, B and C. Due to the friction with the tools the shell is assumed to work in plane strain conditions (direction e ${\textstyle _{3}}$ ) during the forming process. This assumption is kept for the springback step although the restriction of the tools is removed. Two different materials have been used for the simulations: aluminum and a mild steel. A complete definition of these materials can be found in [25]. The material is assumed elastic-plastic. An isotropic behavior is used for the elastic component and an orthotropic yield function [26] with exponential isotropic hardening is considered for the plastic component. For the simulations a low blank holder force ( ${\textstyle 2.45}$ kN) has been used.

The deformed shape in Figure 12.b corresponds to the simulation with aluminum. Figures 12.c and d (obtained for aluminum and mild steel, respectively) show the geometrical parameters defined above obtained by different experimental groups (each abscissa corresponds to a different group). The average of the experimental values and those obtained with the B2 element are also shown. The numerical simulation agrees reasonably well between the rather scattered experimental values.



Figure 12: 2-D Draw bending. (a) Geometry. (b) Geometrical parameters to measure the springback. (c) Comparison with experimental values (aluminum); (d) Comparison with experimental values (mild steel).

7 Conclusions

A finite element for non-linear, large strain, analysis of axisymmetric shells, beams has been presented. Smooth, kinked and branching situations have been considered. The main characteristic of the element formulation is that it does not include rotational degrees of freedom (rotation-free element). The curvature is computed resorting to the positions of the neighbor elements. The element can be considered as an “assumed strain element”, with a linear interpolation of the curvature. The numerical examples show that the element converges to the correct solution in all cases, and that the formulation can handle non-smooth surfaces and branching situations. Two integration points are necessary. Compared with the standard two-node shear deformable element this is a slight disadvantage when an explicit code is used. The element has given excellent results in problems including plasticity with large strains, contact with friction, different boundary conditions and loads.

Acknowledgments

The first author is a member of the scientific staff of the Science Research Council of Argentina (CONICET). The support provided by a grant of CONICET is gratefully acknowledged. Financial support from CIMNE and support from Quantech ATZ (Barcelona-Spain) for providing Stampack, are also acknowledged.

References

[1] A.C. Ugural. Stresses in Plates and Shells, McGraw Hill, New York, 1981.

[2] D. Bushnell and B.O. Almroth, “Finite difference energy method for non linear shell analysis”, J. Computers and Structures, Vol. 1, 361, 1971.

[3] D. Bushnell, “Computerized analysis of shells–-governing equations”, J. Computers and Structures, Vol. 18(3), 471-536, 1984.

[4] R.A. Nay and S. Utku. “An alternative to the finite element method”. Variational Methods Eng., Vol. 1, 1972.

[5] M.R. Barnes. Form finding and analysis of tension space structures by dynamic relaxation. Ph.D. Tesis, Department of Civil Engineering, The City University, London, 1977.

[6] J.K. Hampshire, B.H.V. Topping and H.C. Chan. “Three node triangular elements with one degree of freedom per node”. Engng. Comput. Vol. 9, pp. 49–62, 1992.

[7] R. Phaal and C.R. Calladine. “A simple class of finite elements for plate and shell problems. II: An element for thin shells with only translational degrees of freedom”. Int. J. Num. Meth. Engng., Vol. 5, pp. 979–996, 1992.

[8] E. Oñate and M. Cervera. “Derivation of thin plate bending elements with one degree of freedom per node”. Engineering Computations, Vol. 10, pp 553–561, 1993.

[9] 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 Technology, Ed. M.J.M. Barata Marques, 18th IDDRG Biennial Congress, Lisbon, 1994.

[10] 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.

[11] 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, 813-833, 2001.

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

[13] E. Oñate, P. Cendoya and J. Miquel. Non linear explicit dynamic analysis of shells using the BST rotation-free triangle. Engineering Computations, vol.9(6), 662–706, 2002.

[14] F.G. Flores and E. Oñate. Improvements in the membrane behaviors of the three node rotation-free BST shell triangle using an assumed strain approach. Computer Methods in Applied Mechanics and Engineering, Vol. 194, pp 907-932, 2005.

[15] R. Phaal and C.R. Calladine. A simple class of finite elements for plate and shell problems. I: Elements for beams and thin plates. Int. J. Num. Meth. Engng., Vol. 5, pp. 955–977, 1992.

[16] J. Jovicevic and E. Oñate. Analysis of beams and shells using a rotation-free finite element-finite volume formulation, Monograph 43, CIMNE, Barcelona, 1999.

[17] F.G. Flores. Two Dimensional Shell Element for Nonlinear Analysis. Applied Mechanics Reviews, Vol. 48(11 part.2), pp.30-35. 1995.

[18] F.G. Flores y L.A. Godoy. Finite Element Applications to the Internal Pressure Loadings on Spherical Shells and other Shells of Revolution, Finite Element Analysis of Thin Walled Structures (Editor John Bull), Barking (England), Cap.9 pp. 259-296, 1990.

[19] Prinja NK and Clegg RA , A review of benchmark problems for geometric non-linear behaviors of 3-D beams and shells, NAFEMS R0024, 1993.

[20] Hibbit, Karlson and Sorensen Inc., ABAQUS, version 5.2.1, Pawtucket, USA.

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

[22] J.K. Lee, R. Wagoner, and E. Nakamachi. A benchmark test for sheet metal forming analysis, Technical Report, Ohio State University, 1990.

[23] F.G. Flores. A two-dimensional linear assumed strain triangular element for finite deformations analysis. Journal of Applied Mechanics (submitted) 2005.

[24] C.G. García Garino. A Numerical Model for Large Strain Analysis of Elasto-plastic Solids, (in Spanish). Ph.D. Thesis, Universidad Politécnica de Cataluña, 1993.

[25] A. Makinouchi et al (Ed), Proceedings of NUMISHEET 93 conference, Isehari, Japan, 1993.

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

APPENDIX A

The variation ${\textstyle \delta \phi ^{1}}$ can be obtained noting that

$\cos \phi ^{1}=\mathbf {t} \cdot \mathbf {t} _{l}=\mathbf {n} \cdot \mathbf {n} _{l}$	(84.a)
$\sin \phi ^{1}=\mathbf {n} \cdot \mathbf {t} _{l}=-\mathbf {t} \cdot \mathbf {n} _{l}$	(84.b)

Taking the variation of the Eq.(84.a)

(85)

In Eq.(85) ${\textstyle \delta \mathbf {t} }$ has been obtained observing that t and n are unit orthogonal vectors

(86)

then ${\textstyle \delta \mathbf {t} }$ has components only along the direction ${\textstyle \mathbf {n} }$ :

(87)

Finally, from Eq.(85)

(88)

Similarly for node 2 where the angle ${\textstyle \phi ^{2}}$ is defined from ${\textstyle \mathbf {t} _{r}}$ to ${\textstyle \mathbf {t} }$ in such a way that:

$\cos \phi ^{2}=\mathbf {t} \cdot \mathbf {t} _{r}$	(89.a)
$\sin \phi ^{2}=-\mathbf {n} \cdot \mathbf {t} _{r}$	(89.b)