Solid shell element on prisms

Table of Contents

Work in progress.


A solid-shell element which does not possess rotational degrees of freedom (DOFs) and which is applicable to thin plate/shell problems is considered. The element approximation is constructed in prisms, where displacements on the upper and lower surfaces are approximated in the global coordinate system. In addition, two other fields are defined in the shell natural (local) coordinate system that represent the components of the displacement vector in both the current shell normal direction and the current shell tangent plane. To each field, an arbitrary order of approximation can be defined, and all fields reproduce a complete and conforming polynomial approximation basis for the solid prism element. It is not necessary to augment the formulation with an assumed natural strain (ANS) field or enhanced assumed strain (EAS) field or to use reduced integration, making the element ideally suited for geometrically and physically nonlinear problems.

Pinched cylinder

In this work quadratic (or higher-order for thick shells) displacement approximation in direction normal to mid-surface can be set to avoid finite element locking. As result the thickness of the element changes after the defamation. Moreover, displacements on bottom and top surface are independently approximated with arbitrary polynomial order. As consequence straight lines normal to the mid-surface in the initial configuration in a current configuration are not normal to the mid-surface. For this reason, the kinematics of this solid-shell element is richer than the kinematic assumed in Kirchhoff-Love plate theory.

This work builds on a substantial body of published work by a number of different authors, most notably [23]. However, this implementation propose new formulation exploiting hierarchical, heterogeneous and anisotropic approximation spaces.

Recipe for construction of approximation space

Nodal basis functions on prism

Constructing approximation space in prism we apply similar procedure to shown in [1] for tetrahedral. Nodal basis functions are expressed in terms of barycentric coordinates as follows

\[ \phi^v = \lambda_v,\quad {}_L\phi^v = \phi^v (1-\zeta),\quad {}_U\phi^v = \phi^v \zeta \]

where left subscript in \({}_L\phi^v\) and \({}_U\phi^v\) indicates shape functions on lower and upper triangles respectively, \(\lambda_v\) denote the barycentric coordinates and \({}v\) here stands for vertices of the triangle. Convective coordinate \(\zeta \in [0,1]\) is coordinate through prism thickness. If \(\zeta = 0\) then point is on lower prism triangle and \(\zeta=1\) point is on upper triangle. The barycentric are defined as follows

\[ \lambda_0 = 1-\xi-\eta, \quad \lambda_1 = 1-\xi,\quad \lambda_2 = 1-\eta \]

where property can be noted that \(\lambda_0+\lambda_1+\lambda_2=1\) on triangle.

Nodal basis functions on prism

A set of \(\{L_l: l=0,1,\dots \}\) constitutes hierarchical basis on interval \([-1,1]\);

\[ L_0(s)=1;\quad L_1(s) = s,\quad L_{l+1} = \frac{2l+1}{l+1}sL_l(s)-\frac{l}{l+1}L_{l-1}(s) \]

Hierarchy of edge basis functions on triangles

The edge basis consist basis

\[ \beta_{0i} = \lambda_0\lambda_i,\quad \phi^{e_t}_l=\beta_{0i} L_l(\lambda_i-\lambda_j), \quad {}_L\phi^{e_t}_l=\phi^{e_t}_l(1-\zeta), \quad {}_U\phi^{e_t}_l=\phi^{e_t}_l\zeta, \]

where \(i\) and \(j\) are nodal indices on triangle and \(e_t\) strand for triangle edges. \(L_l\) is Legendre polynomial of order \(l\). If \(p\) is polynomial order on triangle, thus \(0 \geq l \geq p-2\) and number of DOFs on edge is \(p-1\) basis functions.

Hierarchy of triangle basis functions

The triangle consist basis

\[ \beta_{0ij}=\lambda_0\lambda_i\lambda_j,\quad \phi^t_{l,m}=\beta_{0ij}L_l(\lambda_0-\lambda_i)L_m(\lambda_0-\lambda_j), \quad {}_L\phi^t_{l,m}=\phi^t_{l,m}(1-\zeta), \quad {}_U\phi^t_{l,m}=\phi^t_{l,m}\zeta, \]

Here \(t\) stands for triangel. If \(p\) is polynomial order on triangle, thus \(0 \geq l,m,l+m \geq p-2\) and number of DOFs on triangle is \((p-1)(p-2)/2\) basis functions.

Hierarchy of thickness edge basis functions

The edge through thickness basis function

\[ \beta_{00} = \lambda_0\lambda_0,\quad \phi^{e_q}_l = \beta_{00} \zeta (1-\zeta) L_l(2\zeta-1) \]

where \(\lambda_0\) is barycentric coordinate for the node of the triangle to which edge through thickness is adjacent. If \(p\) is polynomial order in prism, thus \(0 \geq l \geq p-2\) and number of DOFs on edge is \(p-1\) basis functions. The quadrilateral through thickness basis function is

Hierarchy of thickness quadrilateral basis functions

The quadrilateral through thickness basis function

\[ \phi^q_{l,m} = \beta_{0i} \zeta (1-\zeta) L_l(\lambda_0-\lambda_i)L_m(2\zeta-1) \]

where \(0\) and \(j\) indicate nodes on opposite corner nodes of quadrilateral in its own canonical numbering. If \(p\) is polynomial order on prism, thus \(0 \geq l,m,l+m \geq p-4 \) the number of degrees of freedom on quadrilateral is \((p-3)(p-2)/2\).

Hierarchy of bubble thickness prism basis functions

The bubble prism basis functions are given by

\[ \phi^p_{l,m,k} = \beta_{0ij} \zeta(1-\zeta) L_l(\lambda_0-\lambda_i)L_m(\lambda_0-\lambda_j)L_k(2\zeta-1) \]

where \(0,i,j\) are indices of nodes on the triangle. If \(p\) is polynomials order in prism, thus \(0 \geq l,m,k,l+m+k \geq p-5\) and number of DOFs in prism is \((P-5)(P-4)(P-3)/6\).

Linear independence

Similar reasoning presented below can be found in [1]. Suppose that

\[ \sum_{L,U} \left( \sum_v w^v + \sum_{e_t} w^{e_t} + \sum_t w^t \right) + \sum_{e_q} w^{e_q} + \sum_q w^q + w^p = 0 \]

where \(w^v\) is a function which span on vertices, \(w^{e_t}\) span on triangles edges, and similar for remaining terms. To show that base functions are linear independent we need to show that as consequence all terms in above equation are zero. Since base functions for edges, triangles, quads and prism are bubble functions, i.e. vanish on vertices. Consequently that barycentric coordinates are linear independent, first term has to be zero. Analogically, since Lagrange polynomials are linearly independent and all remaining terms are zero on given edge, w^{e_t} vanish in above equation. Repeating the same reasoning we can conclude that all terms are zero, thus base is linearly independent.

See implementation of element here MoFEM::FatPrismElementForcesAndSourcesCore

Geometry approximation in reference configuration

The practicality of approach demands that input surface mesh is generated by a standard mesh generator. In case of non-planar surface, geometry could be defined by 6-noded or higher-order triangles. In this formulation shell geometry is described consistently by hierarchal approximation basis, directly inferred from input triangular mesh. For given input surface triangular mesh no other data are needed with exception to the shell thickness.

Reference position of points on shell mid-surface in the global Cartesian coordinate system is given on triangular mesh by

\[ {^MZ_r}(\xi,\eta) = \sum_v \phi_v(\xi,\eta) \underline{Z}_r^v + \sum_{e_t} \sum_l \phi^{e_t}_l(\xi,\eta) (\underline{Z}_r)^{e_t}_l+ \sum_{l,m} \phi^t_{lm}(\xi,\eta) (\underline{Z}_r)^t_{lm} \]

where \(\underline{Z}^{v,e_t,t}\) are DOFs on vertices, edges and triangles, right subscript \({}(\cdot)_r\) indicate coefficient of position vector in the global coordinate system. The left superscript \({}^M(\cdot)\) indicates that field is defined only on mid-surface. Above expression can be given in vector-matrix form,

\[ ^M\mathbf{Z}(\xi,\eta) = {}^3\boldsymbol{\phi}^v(\xi,\eta) \underline{\mathbf{Z}}^v+ {}^3\boldsymbol{\phi}^{e_t}(\xi,\eta) \underline{\mathbf{Z}}^{e_t}+ {}^3\boldsymbol{\phi}^{t}(\xi,\eta) \underline{\mathbf{Z}}^{t}= \sum_{g=v,e_t,t} {}^3\boldsymbol{\phi}^g(\xi,\eta) \underline{\mathbf{Z}}^{g} \]

where \(\underline{Z}^{v,e_t,t}\) are DOFs on vertices, edges and triangles. Left superscript \({}^3\boldsymbol{\phi}\) indicate that it is matrix base approximation functions for vector field.

In addition we define on mid-surface a field of unit length director vectors (see [40]), as follows

\[ ^M\mathbf{V}(\xi,\eta)= \sum_{g=v,e_t,t} {}^3\boldsymbol{\phi}^g(\xi,\eta) \underline{\mathbf{V}}^{g} \]

in addition we define on mid-surface plate thickness field as follows

\[ ^Ma(\xi,\eta)= \sum_{g=v,e_t,t} {}^1\boldsymbol{\phi}^g(\xi,\eta) \underline{\mathbf{a}}^{g} \]

where left subscript \({}^1\boldsymbol{\phi}\) indicate that is matrix of base functions for scalar field.

Finally, reference position vector in prism element is given by in global Cartesian coordinate system

\[ \mathbf{Z}(\xi,\eta,\zeta) = {}^M\mathbf{Z}(\xi,\eta)+ {}^Ma(\xi,\eta){}^M\mathbf{V}(\xi,\eta)(\zeta-\frac{1}{2}), \]

where equivalently expressed on prism element

\[ \mathbf{Z}(\xi,\eta,\zeta) = \sum_{s=L,U}\sum_{g=v,e_t,t} {}^3\boldsymbol{\phi}^g_s(\xi,\eta,\zeta) \underline{\mathbf{Z}}^{g}_{s}. \]

In the following derivations for simplicity we restrict ourselves only to the case that shell thickness is constant, i.e. \({}^Ma(\xi,\eta)=\textrm{const}=a\). The Jacobian is defined as

\[ \mathbf{J}(\xi,\eta,\zeta) = \{ J_k^J \} = {}^M\mathbf{J}(\xi,\eta) + {}^V\mathbf{J}(\xi,\eta,\zeta) \]


\[ {}^M\mathbf{J}(\xi,\eta,\zeta) = \left[ \begin{array}{ccc} \frac{\partial {}^M_0Z}{\partial \xi} & \frac{\partial {}^M_0Z}{\partial \eta} & 0 \\ \frac{\partial {}^M_1Z}{\partial \xi} & \frac{\partial {}^M_1Z}{\partial \eta} & 0 \\ 0 & 0 & 0 \end{array} \right] \]


\[ {}^V\mathbf{J}(\xi,\eta,\zeta) = a \left[ \begin{array}{ccc} 0 & 0 & {}^M_0V \\ 0 & 0 & {}^M_1V \\ 0 & 0 & {}^M_2V \end{array} \right] + a(\zeta-\frac{1}{2}) \left[ \begin{array}{ccc} \frac{\partial {}^M_0V}{\partial \xi} & \frac{\partial {}^M_0V}{\partial \eta} & 0 \\ \frac{\partial {}^M_1V}{\partial \xi} & \frac{\partial {}^M_1V}{\partial \eta} & 0 \\ \frac{\partial {}^M_2V}{\partial \xi} & \frac{\partial {}^M_2V}{\partial \eta} & 0 \end{array} \right] \]

Note that in the following derivations we will be using transpose inverse of Jacobian, given by

\[ J^k_I = \left\{\mathbf{J}^{-\textrm{T}}\right\}^k_I. \]

Curvilinear systems in reference and current configuration

The vector coefficients in the Cartesian coordinate system are \(Z^I\) and \(z^i\), for reference and current configurations respectively. For simplicity we use both coordinate systems with the same base vectors and origin. In addition we use a local curvilinear coordinate base in reference configuration, where base vectors of that coordinate system are \(\mathbf{E}_\alpha\) and covariant coefficients in that base are \(X^A\). The field of \(\mathbf{E}_\alpha\) are approximated consequently using hierarchical approximation as follows

\[ ^M\mathbf{E}_0(\xi,\eta) = E^0_I \mathbf{I}^I = \frac{\partial X^0}{\partial Z^I} \mathbf{I}^I = \left\{ \sum_{g=v,e_t,t} {}^3\boldsymbol{\phi}^g(\xi,\eta) \underline{\mathbf{E_0}}^{g}, \right\} \mathbf{I}^I \quad ^M\mathbf{E}_1(\xi,\eta) = E^1_I \mathbf{I}^I = \frac{\partial X^1}{\partial Z^I} \mathbf{I}^I = \left\{ \sum_{g=v,e_t,t} {}^3\boldsymbol{\phi}^g(\xi,\eta) \underline{\mathbf{E_1}}^{g} \right\} \mathbf{I}^I \]

and base vector in thickness direction is given by

\[ ^M\mathbf{E}^2(\xi,\eta) = E^2_I \mathbf{I}^I = \left\{ \textrm{Spin}[^M\mathbf{E}_0(\xi,\eta)]^M\mathbf{E}_1(\xi,\eta) \right\} \mathbf{I}^I \]

where \(\textrm{Spin}[\cdot]\) is a spin operator acting as a vector product.

The local curvilinear coordinate base in current configuration, convected by shell mid-surface motion is given by

\[ \mathbf{e}_a = e^a_i \mathbf{i}_i = {}^MF^i_I {}^ME^{AI} \delta^a_A \mathbf{i}_i \]

where \({}^M\mathbf{F}^i_I(\xi,\eta,\zeta)\) is component of gradient of deformation, i.e. for discretized system gradient of defamation is additively decomposed on part resulting from DOFs through thickness and part resulting from DOFs on triangles. The later part is used to push base vectors.

Displacement and geometry approximation

In following approach, for convenience we use local shell coordinate system to evaluate physical equations. The Cartesian global coordinate system is used to express coefficients of DOFs on entities adjacent to upper and lower triangles. For DOFs on entities thorough shell thickness, i.e. edges through thickness, quads and prism itself, the current curvilinear (convected) shell coordinate system is used to express coefficients of displacement vector.

The position of material point in current configuration is given by

\[ z^a\mathbf{i}_a = Z^J\mathbf{I}_J + u^J \mathbf{I}_J + v e_2^J \mathbf{I}_J \]

where displacements are additively decomposed on global and through thickness component. The displacement \(\mathbf{u}\) is approximated in global coordinate system

\[ u^J(\boldsymbol \xi) = \sum_{s=L,U} \sum_{g=v,e_t,t} {}^3\boldsymbol{\phi}^g_s(\boldsymbol \xi) \underline{\mathbf{u}}^{g,J}_s = \left[{}^3\boldsymbol{\overline{\phi}}(\boldsymbol \xi)\right] \underline{\mathbf{u}}^J. \]

Note that displacement degrees of freedom in vector \(\underline{\mathbf{u}}\) are in global Cartesian coordinate system, thus element could be assemble with other tetrahedral elements without any linking/transfer elements. The displacement \(\mathbf{v}=v e_2^J \mathbf{I}_J\) is approximated in local conservative system following global mid-surface deformation. The vector of displacements through thickness is given by

\[ v(\boldsymbol \xi) = \sum_{g=e_q,p} {}^1\boldsymbol{\phi}^g(\boldsymbol \xi) \underline{\mathbf{v}}^{g,a} = \left[{}^1\boldsymbol{\overline{\phi}}(\boldsymbol \xi)\right]\underline{\mathbf{v}}. \]

Note that vector of values of degrees of freedom \(\underline{\mathbf{v}}\) is in local, curvilinear current (conservative) coordinate system. Those degrees of freedom are associated with edges in volume of prism, quads and volume of prism itself, thus are inside volume of shell structure and are not adjacent to any other DOFs except those of the considered shell.

With above at hand we define gradients of displacements

\[ \nabla_{Z^J} u^i \left( \mathbf{i}_i \otimes \mathbf{I}^J \right) = \left\{ \sum_{s=U,L}\sum_{g=v,e_t,t} \frac{\partial {}^3 \{\boldsymbol \phi^g_s\}^i \mathbf{i}_i }{\partial \xi^k} J_J^k \right\} \underline{\mathbf{z}}^g_s \otimes \mathbf{I}^J,\; \nabla_{X^A} u^a = e_j^a \nabla_{Z^J} u^j E^J_A \]

and gradient of higher-order displacements through thickness

\[ \nabla_{X^A} v^2 \left( \mathbf{e}_a \otimes \mathbf{E}^A \right) = \left\{ \sum_{g=e_q,p} \frac{\partial {}^1 \{\boldsymbol \phi^g\}^a \mathbf{e}_a }{\partial \xi^k} J_J^k E_A^J \right\} \underline{\mathbf{z}}^{g,a} \otimes \mathbf{E}^A. \]

As result gradient of deformation is given by

\[ {}^MF^i_I = \delta^i_I + \nabla_{Z^I} u^i = \delta^i_I + \left[\nabla_{Z^I}{}^3\boldsymbol{\overline{\phi}}(\boldsymbol \xi)\right] \underline{\mathbf{u}}^J, \]


\[ F^a_A = e^a_iE^I_A\delta^i_J + \nabla_{X^A} u^a + \nabla_{X^A} v^a = e^a_i{}^MF^i_I E^I_A + \left[\nabla_{X^a}{}^1\boldsymbol{\overline{\phi}}(\boldsymbol \xi)\right] \underline{\mathbf{v}}^a \]

where \(\mathbf{v} = [ 0,0,v ]\), i.e. component of displacement in normal direction is considered for second higher-order approximation through thickness is considered.

The displacement approximation can be easily further extended by higher order through thickness approximation of displacements which are tangent to mid-surface.

Physical equation

In the spirit of solid shell element, a physical equation is expressed in curvilinear current shell coordinate system. The right Cauchy–Green defamation tensor is given by

\[ C^A_B = g_{ab}G^{AC} F^b_C F^a_B, \]

where \(g_{ab}\) and \(G^{AC}\) are coefficients of metric tensors in current and reference configuration respectively. The Green strain is

\[ E^A_B = \frac{1}{2}(C^A_B-\delta^A_B). \]

For simplicity in this description we restrict ourselves to St. Venant–Kirchhoff material,

\[ S^A_B = \frac{1}{2} \lambda E^A_B + \mu E^A_B E^A_B, \]

where \(S^A_B\) are coefficients of second Piola-Kirchhoff stress tensor in reference shell local coordinate system. Consequently the first Piola-Kirchhoff stress is given by

\[ P^A_a = g_{ac} G^{AC} F^c_B S^B_C, \]

where \(P^A_a\) are coefficients of stress tensor in current local and reference local coordinate system. Finally, alternatively the first Piola-Kirchhoff stress can be given in global coordinate system

\[ P^I_i = E_A^I e^a_i P^A_a. \]

Element formulation

Problem is solved using standard Newton-Raphson method, where residual vector is expanded in Taylor series and truncated after linear term. Since solution can have local instabilities and is strongly non-linear, classical quadratic arc-length method with line searchers is applied.

The vectors of internal forces are expressed by

\[ \mathbf{f}^\textrm{int}_\textrm{u}(\underline{\mathbf{u}},\underline{\mathbf{v}}) = \int_\Omega \left[\nabla_{Z^I}{}^3\boldsymbol{\overline{\phi}} \right] P^I_i(\underline{\mathbf{u}},\underline{\mathbf{v}}) \textrm{d}\Omega, \]


\[ \mathbf{f}^\textrm{int}_\textrm{v}(\underline{\mathbf{u}},\underline{\mathbf{v}}) = \int_\Omega \left[\nabla_{X^A}{}^1\boldsymbol{\overline{\phi}} \right] P^A_a(\underline{\mathbf{u}},\underline{\mathbf{v}}) \textrm{d}\Omega, \]

Consequently tangent matrix is given by

\[ \left[ \begin{array}{cc} \frac{\partial \mathbf{f}^\textrm{int}_\textrm{u}(\underline{\mathbf{u}},\underline{\mathbf{v}})}{\partial\underline{\mathbf{u}}} & \frac{\partial \mathbf{f}^\textrm{int}_\textrm{u}(\underline{\mathbf{u}},\underline{\mathbf{v}})}{\partial\underline{\mathbf{v}}} \\ \frac{\partial \mathbf{f}^\textrm{int}_\textrm{v}(\underline{\mathbf{u}},\underline{\mathbf{v}})}{\partial\underline{\mathbf{u}}} & \frac{\partial \mathbf{f}^\textrm{int}_\textrm{v}(\underline{\mathbf{u}},\underline{\mathbf{v}})}{\partial\underline{\mathbf{v}}} \end{array} \right] \left\{ \begin{array}{c} \delta \underline{\mathbf{u}} \\ \delta \underline{\mathbf{v}} \end{array} \right\} = \left[ \begin{array}{c} \mathbf{f}^\textrm{ext}_\textrm{u} - \mathbf{f}^\textrm{int}_\textrm{u}(\underline{\mathbf{u}},\underline{\mathbf{v}})\\ - \mathbf{f}^\textrm{int}_\textrm{v}(\underline{\mathbf{u}},\underline{\mathbf{v}}) \end{array} \right] \]

where \(\delta \underline{\mathbf{u}}\) and \(\delta \underline{\mathbf{v}}\) are subincrements of displacements in global and current local coordinate system.

In presented here formulation current curvilinear system follows shell deformation given by DOFs on upper and lower triangle. This in some sense lead to co-rotational like formulation [15]. In this formulation we note that the tangent stiffness is non-symmetric. This observation is consistent with that made by Crisfield [15], who adapted co-rotational formulation where coordinate systems is rotated while shell is deformed. However, numerical experiments have shown that for the problems addressed here, the tangent stiffness matrix becomes symmetric as the iterative procedure reaches equilibrium, thus matrix can be symmetrized without deteriorating convergence rate.

Benchmark test

All tests in this document are executed by running script

./pinched_cylinder_convetgence_test.sh 2>&1 | tee log

User have to edit that file to select desired test.

Small strain/displacement pinched cylinder

Pinched cylinder no diaphragm

Pinched cylinder with diaphragm


If you find this document useful and you would like to cite us, pleas follow that link this link..