The development of an **efficient shell element** has gathered a lot of work since
the beginning of structural finite element history. In this article, I will
present the linear quad shell element (improved MITC4 element) that has been
freshly implemented in SesamX.
I will start with the assumptions made as well as the derivation of the element
stiffness matrix. Finally, I will showcase some SesamX data cards and a
**comparison of results with Abaqus S4 element**.

## Quick remark on notations

Before starting, let’s define some notations that are used through this article:

vectors are denoted with an underline $ \underline{u} $ and tensors with a double underline $ \underline{\underline{\varepsilon}} $ ,

matrices are represented with brackets, as well as vector and tensor components $ \lbrack \ \rbrack $ ,

Einstein summation convention is used on repeated indices,

partial derivatives are denoted with the comma notation $ u_{,x} $ ,

infinitesimal strain and stress tensors are represented in column matrix notation $ \lbrack \epsilon \rbrack = \begin{bmatrix} \varepsilon_{11} \\ \varepsilon_{22} \\ \varepsilon_{33} \\ 2 \varepsilon_{12} \\ 2 \varepsilon_{23} \\ 2 \varepsilon_{13} \end{bmatrix} $ and $ \lbrack \sigma \rbrack = \begin{bmatrix} \sigma_{11} \\ \sigma_{22} \\ \sigma_{33} \\ \sigma_{12} \\ \sigma_{23} \\ \sigma_{13} \end{bmatrix} $ , where $ \lbrack \varepsilon \rbrack $ represents the

**engineering strains**.

## Introduction

**Shell element are made to analyse thin structures**
undergoing forces and moments in any direction. The structure can be flat or
curved but it should have **one dimension smaller than the other 2**. For instance,
it is the perfect element to analyse a slab or an aircraft fuselage.

Slab |
Aircraft fuselage |

For multiple reasons, **developing an effective shell element is not an easy task**
(we will cover some of these reasons through this article). Many
shell elements are available in the literature, and they are still being
improved as time goes on. **In SesamX we chose to
implement the MITC shell** (**M**ixed **I**nterpolation of **T**ensorial
**C**omponents) because it recently received new improvements making it **one of
the most practical and efficient shell element**. You will find all the details
in the following resources:

Dvorkin EN, Bathe KJ. A continuum mechanics based four-node shell element for general nonlinear analysis: original paper about the MITC4 element,

Ko Y, Lee PS, Bathe KJ. A new MITC4+ shell element: recent paper exposing new improvements of the MITC4 element,

Ko Y, Lee PS, Bathe KJ. A new 4-node MITC element for analysis of two dimensional solids and its application in shell analyses: other improvements of the MITC4 element,

Ko Y, Lee Y, Lee PS, Bathe KJ. Performance of the MITC3+ and MITC4+ shell elements in widely-used benchmark problems: comparison of the improved MITC4 element against other well known shell elements (including Abaqus S4 element).

Basically, SesamX quad shell element is a concatenation of all these papers leading to the MITC4+ element mentioned in [4].

To fully define the shell geometry, the element needs the following 3 characteristics:

The

**mid-surface geometry**,The

**thickness**at each point of the mid-surface: $ h $ ,The

**director vector**at each point of the mid-surface: $ \underline{v} $

The mid-surface and thickness are what comes naturally in mind when thinking about shells.

**The director vector is a specific MITC feature**. It allows to control how the
“cross section” is tilted relative to the mid-surface.

Theses features make it possible to modelize shells with varying thickness along the mid-surface, as well as arbitrary shells connections.

Varying thickness shell |
Shell connection |

Furthermore, this element is not restricted to flat shell and **can also represent
curved geometries**. The element is meshed on the mid-surface using 4 nodes. At
each node $ I $
the following quantities are available:

the

**position**of the node $ \underline{X^I} $ ,the

**thickness**at the node $ h^I $ ,the

**shell director vector**at the node $ \underline{v^I} $ ,the

**displacement vector**at the node $ \underline{u^I} $ ,and the

**rotation vector**at the node $ \underline{\theta^I} $ .

## Assumptions

In essence, the shell element is a simplification of a more complex 3D structure. When designing such an element, the aim is to decrease the computational cost, while making relevant assumptions on how the element behaves. The main kinematic assumptions underlying the MITC element are given here with some explanations.

### Stress assumption

As the shell must be able to sustain transverse and plane loading, the stress state must allow for transverse shear stress, as well as in plane stress. It is then sufficient to enforce:

$$ \boxed{\sigma_{33} = 0} \tag{1} $$

This is **the classical stress assumption** made for shell elements when the “3”
direction is normal to the mid-surface. However, in the case of the MITC element,
**the “3” direction is collinear to the director vector at each point of the shell
mid-surface** (basis and coordinate systems are detailed further ahead in this
article).

### Kinematic assumption

Before detailing the kinematic assumption, we need to define a convenient parametric coordinate system on the shell.

#### Parametric coordinates system and position vector

We rely on the position vector to define the parametric system. Without providing more details, the first 2 coordinates $ r $ and $ s $ can be any of the mid-surface curvilinear coordinates. The last coordinate $ t $ is such that for a given $ r $ and $ s $ , $ t $ evolves when we move along the director vector.

Using these parametric coordinates, we can describe the **position of every
point in the space**:

$$ \boxed{\underline{X}(r,s,t) = \underline{X_0}(r,s) + \cfrac{t}{2} h(r,s) \underline{v}(r,s)} \tag{2} $$

Where **$ \underline{X_0}(r,s) $
describes the mid-surface
position**. To define clearly what $ r $
and
$ s $
are, we need to detail further
$ \underline{X_0}(r,s) $
, $ h(r,s) $
and $ \underline{v}(r,s) $
.

Using the position of the nodes and the shape functions we have:

- $ \underline{X_0}(r,s) = N^I(r,s)\underline{X^I} $
- $ h(r,s) \underline{v}(r,s) = N^I(r,s)h^I \underline{v}^I $

Where the shape functions $ N^I $ are:

$$ \begin{cases} N^1(r,s) = \cfrac{1}{4} (1-r) (1-s) \\ N^2(r,s) = \cfrac{1}{4} (1+r) (1-s) \\ N^3(r,s) = \cfrac{1}{4} (1+r) (1+s) \\ N^4(r,s) = \cfrac{1}{4} (1-r) (1+s) \\ \end{cases} $$

Which yields the full expression of the **position from the parametric
coordinates**:

$$ \boxed{\underline{X}(r,s,t) = N^I(r,s)\underline{X^I} + \cfrac{t}{2} N^I(r,s)h^I \underline{v}^I} \tag{3} $$

Please note the following remarks:

for $ t = 0 $ the domain $ r, s \in [-1, 1] $ describes the shell mid-surface,

$ t $ moving from -1 to 1 corresponds to a point moving from the inferior face to the superior face of the shell,

the full product $ h(r,s) \underline{v}(r,s) $ is linearly interpolated from its values at the nodes (we want to stress that $ h(r,s) $ and $ \underline{v}(r,s) $ are not interpolated on their own but through their product),

even though expression (3) is useful to define the parametric space, we use extensively expression (2) in the following for its conciseness.

#### Displacement vector

Next we express the **displacement vector at each point inside the shell**, relative
to its parametric coordinates. We have:

$$ \underline{u}(r,s,t) = \underline{u_0}(r,s) + \cfrac{t}{2} h(r,s) (\underline{\theta}(r,s) \times \underline{v}(r,s)) $$

Where:

$ \underline{u_0}(r,s) $ is the displacement vector on the shell mid-surface,

$ \underline{\theta}(r,s) $ is the rotation vector on the shell mid-surface.

We see from this formula that, for a given point on the mid-surface, the rotation vector $ \underline{\theta}(r,s) $ describes simply the rotation of the director vector. Hence $ \cfrac{t}{2} h(r,s) (\underline{\theta}(r,s) \times \underline{v}(r,s))$ is simply the displacement resulting from the director vector rotation.

However, this comes at a price. In general, if we seek to compute the rotation of
a vector $ \underline{v} $
under a given
rotation vector $ \underline{\theta} $
,
we have to apply Rodrigues’ formula (more information here).
**This formula is non linear on $ \underline{\theta} $
**, but
can be simplified to
$ \underline{\theta} \times \underline{v} $
if we assume that
$ \underline{\theta} $
is small (which is the case here).

As for the position vector, we relate the **displacement inside the shell from the
displacements and rotations at the nodes**:

$$ \underline{u}(r,s,t) = N^I(r,s)\underline{u^I} + \cfrac{t}{2} N^I(r,s) h^I (\underline{\theta^I} \times \underline{v^I}) $$

And we note again that the full term $ h(r,s) (\underline{\theta}(r,s) \times \underline{v}(r,s)) $ is interpolated at once from the values at the nodes. To make things easier, we can introduce the vector $ \underline{\phi} = \underline{\theta} \times \underline{v} $ to finally write:

$$ \boxed{\underline{u}(r,s,t) = \underline{u_0}(r,s) + \cfrac{t}{2} h(r,s) \underline{\phi}(r,s)} \tag{4} $$

And after interpolation from the values at the nodes:

$$ \boxed{\underline{u}(r,s,t) = N^I(r,s)\underline{u^I} + \cfrac{t}{2} N^I(r,s) h^I \underline{\phi^I}} \tag{5} $$

## Basis and coordinate systems

Before going further, let’s define and summarize the different basis and coordinate systems that are used while deriving the stiffness matrix.

### Covariant basis

We have defined previously the parametric coordinate system ($ r $
,
$ s $
, $ t $
). It is time now to
give **the associated covariant basis**. Using (2), this basis is defined as:

$$ \begin{cases} \begin{alignedat}{1} \underline{g_r}(r,s,t) &= \underline{X}_{,r} = \underline{X_0}_{,r} + \cfrac{t}{2} (h \underline{v})_{,r} \\ \underline{g_s}(r,s,t) &= \underline{X}_{,s} = \underline{X_0}_{,s} + \cfrac{t}{2} (h \underline{v})_{,s} \\ \underline{g_t}(r,s) &= \underline{X}_{,t} = \cfrac{1}{2} h \underline{v} \end{alignedat} \end{cases} \tag{6} $$

Remarks:

$ \underline{g_r} $ and $ \underline{g_s} $ depend on $ r $ , $ s $ and $ t $ while $ \underline{g_t} $ depends only on $ r $ and $ s $ ,

if the element has the same director vectors and thickness at the nodes, then $ h \underline{v} $ is constant, and now $ \underline{g_r} $ and $ \underline{g_s} $ depend only on $ r $ and $ s $ ,

you can have a look here for more details about the covariant basis.

In a sense, **the covariant basis “follows” the geometry of the element**.
Therefore, it is used to compute specific strains inside the element.

### Contravariant basis

We will also use the **contravariant basis** associated with the previous covariant
one. It is defined as follow:

$$ \begin{cases} \begin{alignedat}{1} \underline{g^r}(r,s,t) &= \cfrac {\underline{g_s} \times \underline{g_t}} {\underline{g_r} \cdot (\underline{g_s} \times \underline{g_t})} \\ \underline{g^s}(r,s,t) &= \cfrac {\underline{g_t} \times \underline{g_r}} {\underline{g_r} \cdot (\underline{g_s} \times \underline{g_t})} \\ \underline{g^t}(r,s,t) &= \cfrac {\underline{g_r} \times \underline{g_s}} {\underline{g_r} \cdot (\underline{g_s} \times \underline{g_t})} \end{alignedat} \end{cases} \tag{7} $$

This basis is used conjointly with the covariant basis to manipulate the strain tensor. In fact, it has the following properties: $ \underline{g^i} \cdot \underline{g_j} = 0 $ $ \forall i,j \ i = \not j $ and $ \underline{g^i} \cdot \underline{g_i} = 1 $

### Local basis and local Cartesian coordinate system

From the covariant basis, we also define a **Cartesian local basis**
($ \underline{e_1}, \underline{e_2}, \underline{e_3} $
)
as follow:

$$ \begin{cases} \begin{alignedat}{1} \underline{e_3}(r,s) &= \cfrac {\underline{g_t}} {\| \underline{g_t} \|} \\ \underline{e_1}(r,s,t) &= \cfrac {\underline{g_s} \times \underline{e_3}} {\| \underline{g_s} \times \underline{e_3} \|} \\ \underline{e_2}(r,s,t) &= \underline{e_3} \times \underline{e_1} \end{alignedat} \end{cases} \tag{8} $$

And from this local basis, we define the **local Cartesian coordinate system
($ x $
, $ y $
, $ z $
)**.
We note that $ \underline{e_3} $
is not perpendicular to
the shell mid-surface but is aligned with $ \underline{v} $
.

The local basis is used to express the stress strain relation for the shell. The stress assumption given in (1) applies in this basis. While the local Cartesian coordinate system is used to integrate the virtual work over the element volume.

Note that when $ \underline{v} $ is not perpendicular to the mid-surface, the plane described by $ \underline{e_1} $ and $ \underline{e_2} $ is not on the mid-surface.

## MITC4 shell element derivation

The biggest challenge regarding shell element design is to **avoid shear and membrane
locking**, namely the tendency of the element to represent too rigidly shear and
membrane motions. Multiple techniques have been employed to alleviate this issue,
each having its pros and cons. Regarding the MITC, this issue is tackled without
adding extra degree of freedom or instability to the element. Which make it
**an efficient and conceptually simple shell element**.

In fact, it relies on assuming specific transverse shear and membrane strains
interpolation inside the element from strain values at tying points. And
these points are specifically chosen to avoid locking.
**The fundamental assumption underlying the MITC element is the use of the
covariant strain components** to perform this interpolation. Roughly speaking,
covariant strain components at the tying points are manipulated and inserted as
covariant strain components at any point (for instance Gauss points). Hence the
covariant / contravariant basis act as a hub connecting the strains from
multiple points.

### Strain tensor

The strain tensor can be expressed in whatever basis we choose. To go on with the MITC procedure, we need to express the strains in the contravariant basis to get the covariant strain components. Then the tying procedure is applied on some of these terms to improve the element behavior.

#### Covariant strain tensor components

From definition, the **covariant strains components** are given by:

$$ \underline{\underline{\varepsilon}} = \widetilde{\varepsilon_{ij}} \underline{g^i} \otimes \underline{g^j} \qquad \forall i, j = r, s, t $$

And it can be shown that:

$$ \widetilde{\varepsilon_{ij}} = \cfrac{1}{2} (\underline{u}_{,i} \cdot \underline{g_j} + \underline{u}_{,j} \cdot \underline{g_i}) $$

Using (4) we can simplify further this expression:

$$ \begin{cases} \begin{alignedat}{1} \widetilde{\varepsilon_{rr}} &= \underline{u_0}_{,r} \cdot \underline{g_r} + \cfrac{t}{2} (h \underline{\phi})_{,r} \cdot \underline{g_r} \\ \widetilde{\varepsilon_{ss}} &= \underline{u_0}_{,s} \cdot \underline{g_s} + \cfrac{t}{2} (h \underline{\phi})_{,s} \cdot \underline{g_s} \\ \widetilde{\varepsilon_{tt}} &= \cfrac{1}{2} h \underline{\phi} \cdot \underline{g_t} = 0 \\ \widetilde{\varepsilon_{rs}} &= \cfrac{1}{2} (\underline{u_0}_{,r} \cdot \underline{g_s} + \underline{u_0}_{,s} \cdot \underline{g_r}) + \cfrac{t}{4} \big ((h \underline{\phi})_{,r} \cdot \underline{g_s} + (h \underline{\phi})_{,s} \cdot \underline{g_r} \big) \\ \widetilde{\varepsilon_{st}} &= \cfrac{1}{2} (\underline{u_0}_{,s} \cdot \underline{g_t}) + \cfrac{1}{4} (h \underline{\phi} \cdot \underline{g_s}) + \cfrac{t}{4} \big ((h \underline{\phi})_{,s} \cdot \underline{g_t} \big) \\ \widetilde{\varepsilon_{rt}} &= \cfrac{1}{2} (\underline{u_0}_{,r} \cdot \underline{g_t}) + \cfrac{1}{4} (h \underline{\phi} \cdot \underline{g_r}) + \cfrac{t}{4} \big ((h \underline{\phi})_{,r} \cdot \underline{g_t} \big) \end{alignedat} \end{cases} \tag{9} $$

Where the displacement and rotation vectors (through $ \underline{\phi} $ ) appear more clearly. We note that $ \widetilde{\varepsilon_{tt}} = 0 $ because by definition $ \underline{\phi} $ is perpendicular to $ \underline{g_t} $ . The relations given here are valid for an arbitrary point inside the element. Next step is about modifying these terms from their values at the tying points.

#### Covariant transverse shear tying procedure

The procedure detailed here can be found in details in [1]. We are interested
in reworking the terms $\widetilde{\varepsilon_{st}} $
and $\widetilde{\varepsilon_{rt}} $
that are responsible
for **transverse shear locking**. The idea is to compute
the term $\widetilde{\varepsilon_{st}} $
at the tying points C and D, and the term
$\widetilde{\varepsilon_{rt}} $
at the
tying points A and B (see picture below).
Such that we have, $ \widetilde{\varepsilon^A_{rt}} $
,
$ \widetilde{\varepsilon^B_{rt}} $
,
$ \widetilde{\varepsilon^C_{st}} $
and
$ \widetilde{\varepsilon^D_{st}} $
.

Note that these **tying points are located on the mid-surface**. Eventually,
elsewhere in the element, $\widetilde{\varepsilon_{rt}} $
and
$\widetilde{\varepsilon_{st}} $
are interpolated from the
values at the tying points. Hence, for an arbitrary point located at
$ (r, s)$
we have:

$$ \begin{cases} \begin{alignedat}{1} \widetilde{\varepsilon_{rt}} = \cfrac{1}{2} (1+s) \widetilde{\varepsilon^A_{rt}} + \cfrac{1}{2} (1-s) \widetilde{\varepsilon^B_{rt}} \\ \widetilde{\varepsilon_{st}} = \cfrac{1}{2} (1+r) \widetilde{\varepsilon^C_{st}} + \cfrac{1}{2} (1-r) \widetilde{\varepsilon^D_{st}} \end{alignedat} \end{cases} $$

#### Covariant membrane strain tying procedure

The full procedure can be found in [2] and [3]. It is more involved
than the transverse shear procedure but quite similar.
The idea is to rework the terms $\widetilde{\varepsilon_{rr}} $
,
$\widetilde{\varepsilon_{ss}} $
and
$\widetilde{\varepsilon_{rs}} $
**such that membrane locking
is suppressed** by using specific tying points inside the element. As for the
shear terms, it yields new relations for the terms
$\widetilde{\varepsilon_{rr}} $
,
$\widetilde{\varepsilon_{ss}} $
and
$\widetilde{\varepsilon_{rs}} $
that linearly depend on
strains at the tying points.

#### Summary

Once the tying procedure is completed, we can rewrite (9) using the covariant strains at the tying points. We then insert the displacement and rotation vectors interpolated from the nodes values (5). This yield the full expression of the covariant strain components that linearly depend on the degrees of freedom at the nodes.

The last step is then to **project the covariant strain components into the local
basis** (8). Expressing the strain tensor in the contravariant basis on one hand
and in the local basis on the other hand we have:

$$ \underline{\underline{\varepsilon}} = \widetilde{\varepsilon_{ij}} \underline{g^i} \otimes \underline{g^j} = \varepsilon_{kl} \underline{e_k} \otimes \underline{e_l} \qquad i, j = r, s, t \quad and \quad k, l = 1, 2, 3 $$

Which relates the covariant strain components to the local strain components:

$$ \varepsilon_{kl} = \widetilde{\varepsilon_{ij}} (\underline{g^i} \cdot \underline{e_k}) (\underline{g^j} \cdot \underline{e_l}) $$

Using the fact that $ \widetilde{\varepsilon_{tt}} = 0 $ and that $ \underline{e_3} $ is collinear to $ \underline{g_t} $ (hence $ \underline{e_3} \cdot \underline{g^r} = \underline{e_3} \cdot \underline{g^s} = 0 $ ) we can write the local strains in matrix representation (using engineering strains):

$$ \lbrack \varepsilon \rbrack = \begin{bmatrix} \varepsilon_{11} \\ \varepsilon_{22} \\ \varepsilon_{33} \\ 2 \varepsilon_{12} \\ 2 \varepsilon_{23} \\ 2 \varepsilon_{13} \end{bmatrix} = \begin{bmatrix} \varepsilon_{11} \\ \varepsilon_{22} \\ 0 \\ 2 \varepsilon_{12} \\ 2 \varepsilon_{23} \\ 2 \varepsilon_{13} \end{bmatrix} \tag{10} $$

### Stress tensor

To relate the stresses to the strains we need to apply **Hooke’s law** for
linear elastic material.

$$ \begin{bmatrix} \varepsilon_{11} \\ \varepsilon_{22} \\ \varepsilon_{33} \\ 2 \varepsilon_{12} \\ 2 \varepsilon_{23} \\ 2 \varepsilon_{13} \end{bmatrix} = \frac{1}{E} \begin{bmatrix} 1 && -\nu && -\nu && 0 && 0 && 0 \\ -\nu && 1 && -\nu && 0 && 0 && 0 \\ -\nu && -\nu && 1 && 0 && 0 && 0 \\ 0 && 0 && 0 && 2+2\nu && 0 && 0 \\ 0 && 0 && 0 && 0 && 2+2\nu && 0 \\ 0 && 0 && 0 && 0 && 0 && 2+2\nu \end{bmatrix} \begin{bmatrix} \sigma_{11} \\ \sigma_{22} \\ \sigma_{33} \\ \sigma_{12} \\ \sigma_{23} \\ \sigma_{13} \end{bmatrix} \\ $$

Note that this material law is applied on **strains and stresses expressed in the
local basis** (where $ \underline{e_3} $
is not necessarily
perpendicular to the mid-surface). Using assumption (1) on the right hand side,
as well as the result we got in (10) on the left hand side, leads to a peculiar
relation:

$$ \underbrace{ \begin{bmatrix} \varepsilon_{11} \\ \varepsilon_{22} \\ 0 \\ \varepsilon_{12} \\ \varepsilon_{23} \\ \varepsilon_{13} \\ \end{bmatrix} }_{\text{kinematic assumption}} = \underbrace{ \frac{1}{E} \begin{bmatrix} \sigma_{11} - \nu \sigma_{22} \\ \sigma_{22} - \nu \sigma_{11} \\ -\nu\sigma_{11} - \nu\sigma_{22} \\ (2+2\nu)\sigma_{12} \\ (2+2\nu)\sigma_{23} \\ (2+2\nu)\sigma_{13} \end{bmatrix} }_{\text{stress assumption}} $$

As for the the truss
element, this relation cannot hold. However using a similar
explanation, we can find a way out.
When making the kinematic assumption we were interested in the macroscopic
behavior of the shell. Whereas the stress assumption relates more to a microscopic
behavior. It turns out that even if
$ \varepsilon_{33} $
is not 0 (**microscopic scale**) its
effect on the displacement (**macroscopic scale**) is negligible compared to what
comes from the other terms. Of course, this explanation becomes questionable as
the thickness of the shell increases.

We can therefore simplify this relation and write:

$$ \begin{bmatrix} \varepsilon_{11} \\ \varepsilon_{22} \\ 0 \\ \varepsilon_{12} \\ \varepsilon_{23} \\ \varepsilon_{13} \end{bmatrix} = \frac{1}{E} \begin{bmatrix} \sigma_{11} - \nu \sigma_{22} \\ \sigma_{22} - \nu \sigma_{11} \\ 0 \\ (2+2\nu)\sigma_{12} \\ (2+2\nu)\sigma_{23} \\ (2+2\nu)\sigma_{13} \end{bmatrix} $$

An finally we invert it and omit $ \varepsilon_{33} $ and $ \sigma_{33} $ to get the stresses from the strains:

$$ \begin{bmatrix} \sigma_{11} \\ \sigma_{22} \\ \sigma_{12} \\ \sigma_{23} \\ \sigma_{13} \end{bmatrix} = \underbrace{\begin{bmatrix} \cfrac{E}{1-\nu^2} && \cfrac{\nu E}{1-\nu^2} && 0 && 0 && 0 \\ \cfrac{\nu E}{1-\nu^2} && \cfrac{E}{1-\nu^2} && 0 && 0 && 0 \\ 0 && 0 && \cfrac{E}{2(1+\nu)} && 0 && 0 \\ 0 && 0 && 0 && \cfrac{k E}{2(1+\nu)} && 0 \\ 0 && 0 && 0 && 0 && \cfrac{k E}{2(1+\nu)} \\ \end{bmatrix}}_{\lbrack D \rbrack} \begin{bmatrix} \varepsilon_{11} \\ \varepsilon_{22} \\ \varepsilon_{12} \\ \varepsilon_{23} \\ \varepsilon_{13} \end{bmatrix} $$

Where the transverse shear correction factor $ k = \cfrac{5}{6} $ was introduced.

## Element stiffness derivation

The element stiffness matrix is obtained through the expression of the virtual work. Denoting the virtual strains as $ \lbrack \overline{\varepsilon} \rbrack $ we have for the virtual work over the element volume:

$$ \overline{W} = \iiint_V \lbrack \overline{\varepsilon} \rbrack^T \lbrack \sigma \rbrack dV \tag{11} $$

Where $ dV = dxdydz $
refers to the volume element in the
local Cartesian coordinate system. However, the formula previously derived for
the strains depends on the parametric coordinate system
$ (r, s, t) $
. We can relate this volume element to the
parametric coordinates using the **determinant of the Jacobian matrix**:

$$ dV = |det(\lbrack J \rbrack)| drdsdt $$

It can be shown that the Jacobian matrix is directly given by the **covariant
vectors components expressed in the local basis**:

$$ \lbrack J \rbrack = \big \lbrack \lbrack g_r \rbrack , \lbrack g_s \rbrack , \lbrack g_t \rbrack \big \rbrack $$

Where $ \lbrack g_i \rbrack $ represents the components of the vector $ \underline{g_i} $ in the local basis.

Besides, it can be shown that as long as the $ \lbrack g_i \rbrack $ are components of vectors in an orthonormal basis, we have:

$$ det(\lbrack J \rbrack) = \underline{g_t} \cdot (\underline{g_r} \times \underline{g_s}) $$

Hence after expanding (11) we get the following expression:

$$ \overline{W} = \iiint_V \lbrack \overline{\varepsilon} \rbrack^T \lbrack D \rbrack \lbrack \varepsilon \rbrack |\underline{g_t} \cdot (\underline{g_r} \times \underline{g_s})| drdsdt $$

Unlike classical shell or plate element, **we cannot in general “pre-integrate” along
$ t $
** this volume integral because the dependence on $ t $
is arbitrary. We simply apply a numerical integration over $ r $
,
$ s $
and $ t $
. The 2 x 2 x 2
Gauss rule is enough to integrate exactly along $ r $
and
$ s $
, and has proven to be accurate enough to integrate
along $ t $
for most applications
(even though we cannot ensure an exact integration).

Remark:In the specific case where $ h \underline{v} $ is constant over the element, we have seen that $ \underline{g_r} $ and $ \underline{g_s} $ no longer depend on $ t $ . This makes it possible to integrate exactly along $ t $ without requiring numerical integration.

Eventually, using (10), we can relate each of the terms obtained to the nodes displacements and rotations and recognize the expression of the stiffness matrix.

## SesamX input cards

To define a shell element in SesamX the **first step is to define a mesh**.
This is done with the `CREATE-SUBMESH`

function:

```
CREATE-SUBMESH
MODEL: SHELL_ASSEMBLY
NODES
- ID: 1 POINT: 0.,0.,0.
- ID: 2 POINT: 10.,0.,0.
- ID: 3 POINT: 20.,1.,0.
- ID: 4 POINT: 0.,10.,0.
- ID: 5 POINT: 10.,10.,0.
- ID: 6 POINT: 20.,10.,0.
ELEMENTS
TYPE: QD4
- ID: 1 NODES: 1,2,5,4
- ID: 1 NODES: 2,3,6,5
```

Here we define 6 nodes and we create 2 quad elements to connect the nodes. Next we apply the shell property on these 2 elements:

```
CREATE-PROPERTIES
MODEL: SHELL_ASSEMBLY
SHELL-STANDARD
BEHAVIOR: LINEAR
NAME: PROP1
TH: 80
MATERIAL: STEEL
ON-ELEMENTS-FROM: ALL_SHELL
```

Here we apply a `SHELL-STANDARD`

property on the elements from
`ALL_SHELL`

providing a material name `STEEL`

(that we defined previously,
not shown here) and a thickness.
`BEHAVIOR: LINEAR`

indicates that we apply a linear elastic material. This
keyword may seem unnecessary at the moment, but it is a provision for future
material implementations (such as hyper-elastic materials).

## Comparison with Abaqus S4 element

Eventually, the last part of this article focuses on the comparison of SesamX
MITC4 quad shell element against Abaqus equivalent **S4 element**. A comparison is
showcased here on a **distorted curved mesh**. This mesh can be qualified as a bad
mesh, however it is useful to **assess the robustness of the element** against the
well known Abaqus S4. If you are looking for a more detailed comparison please
read [4].

### Model

An hyperboloid-like model is studied for this comparison. The mesh is voluntarily distorted and curved (the quads are not planar). As can be seen on the following picture.

Mesh |
Loading and boundary conditions |

The mesh definition can be found in the Abaqus input cards given below. A linear elastic material is applied with $ E = 200 GPa $ and $ \nu = 0.33 $ and a thickness of $ 80 mm $ is applied on all the elements. The circle at $ Z = 0 $ is clamped while a uniformly distributed load along $ X $ is applied at the top circle of the hyperboloid.

This case is then solved with a **linear static resolution**.

### Results

We expect to see some differences between Abaqus S4 and SesamX shell.
In fact, the **element formulations are not identical** (as discussed in [4]).
However, these differences should not be large. The following figure illustrates
the displacement of the model.

And the table below gives the comparison of the nodal displacements and rotations between Abaqus and SesamX.

Node 1 is one of the clamped nodes. Nodes 571, 583, 596 and 608 are located on the
top circle of the hyperboloid. The **error on displacements is less than 5%** while
the **error on the rotation is larger (at most 30%)**.

Remark:Abaqus S4 element seems to have 6 degree of freedom per node, while the MITC4+ element has only 5 (no drilling rotation). Hence the S4 adds inevitably some terms on the rotation degrees of freedom. This may explain why the differences are larger on rotations than on displacements.

Even though the error values may seem large (especially for rotations), **these
results are not alarming**. As mentioned previously, the mesh is quite poor
but was chosen precisely to amplify any difference that may appear. Besides,
Abaqus S4 element formulation is different than the formulation presented in
this article. To conclude, it has been proven that the **MITC4+ element
may behave better than Abaqus S4** when the mesh is distorted ([4]).

Following these links you have access to the SesamX input cards as well as the Abaqus input cards that I used for this comparison. I imported the Abaqus mesh and selections in SesamX (to make sure I have the same model description) and then I defined the properties, loads and boundary conditions.

## Conclusion

I hope you had a pleasant reading. Starting from
**version V2020_02 of SesamX**
you can use the shell element presented in this article. Feel free to share
your thoughts or simply ask for more information!

#### Did you like this content?

Register to our newsletter and get notified of new articles