In this article I detail the triangular shell element that is implemented in SesamX.
This element is known as the **MITC3+ element** and is quite similar to the MITC4+
element used for quadrangular shell (presented in this article).
Hence, **most of the element derivation is identical to that of the MITC4+** and
is not detailed again here. This article focuses instead on the specific features
of the MITC3+, such as the presence of an **internal node** and the **stiffness
matrix condensation**.

Finally some SesamX data cards are provided and a **comparison of results with the
Abaqus S3 element is presented**.

## 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

When meshing a complex shell surface, quadrangular and triangular elements are required. Hence, to be practical, a finite element solver must provide both elements and they must be reliable.

A great amount of work has been dedicated to quadrangular shell element, which have led
to the MITC4+, one of the most practical and efficient shell element. Triangular
shell elements have not gathered the same focus. Nevertheless, a similar **M**ixed
**I**nterpolation of **T**ensorial **C**omponents (MITC) can be applied to the
triangular element. And **recent development resulted in one of the best triangular shell
element so far**: the MITC3+ element that we are discussing here. You can find
every detail about this element in the following research papers:

Lee PS, Bathe KJ. Development of MITC isotropic triangular shell finite elements: original paper about the MITC technique applied to triangular element (MITC3),

Lee Y, Lee PS, Bathe KJ. The MITC3+ shell element and its performance: recent paper exposing new improvement to the MITC3 element.

The MITC3+ is conceptually close to the MITC4+ element. However the main
difference between them is that **the MITC3+ element use an additional internal
node in order to improve its bending behavior**. This additional node does not come
from the mesh and is automatically taken into account by the MITC3+ formulation.

As for the MITC4+, the MITC3+ relies on the same shell characteristics and each of the following quantities are available at the nodes:

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} $ .

The director vector and thickness at the internal node are simply assumed to be the mean values from the corner nodes.

## Assumptions

The MITC3+ element relies on assumptions that make it relevant to modelize shell structure. Even though these assumptions are similar to the MITC4+ assumptions, they are specific to the MITC3+ because of the internal node.

### Stress assumptions

The stress assumption mentioned for the MITC4+ element still hold for the MITC3+. In fact, it is the classical stress assumption made for shell and it does not depend on the shape of the element. We assume:

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

We recall that in the case of the MITC element, **the “3” direction is collinear
to the director vector at each point of the shell mid-surface**.

## Kinematic assumption

The kinematic assumptions must take into account the presence of the internal node. Therefore, these kinematic assumptions are specific to the MITC3+ element.

### Position vector

Using the parametric coordinates introduced in the MITC4+ article,
we have for the **position of every point in the element**:

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

In the case of the MITC3+, the parametric coordinates are showcased in the following figure.

However, **the MITC3+ does not rely only on the classical 2D shape functions
of the triangular element**. Indeed, the full expression of the position from the
parametric coordinates is given by:

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

Where the $ N^I $ shape functions are the classical 2D shape functions:

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

And **the $ F^1, F^2, F^3 $
shape functions are the classical
2D shape functions, augmented with the cubic bubble shape function
$ F^4 $
**:

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

This cubic bubble function is associated with the internal node located at the center of the element ($ r = \cfrac{1}{3}, s = \cfrac{1}{3} $ ), and it is such that it vanishes on the element edges.

We note from the formula given in (2) that the internal node appears only in the second part of the relation.

### Displacement vector

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

$$ \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)) $$

Applying the interpolation from the nodes values, we have finally:

$$ \boxed{\underline{u}(r,s,t) = \displaystyle\sum_{I=1}^3 N^I(r,s)\underline{u^I} + \displaystyle\sum_{I=1}^4 \cfrac{t}{2} F^I(r,s) h^I (\underline{\theta^I} \times \underline{v^I})} \tag{3} $$

**Obviously, the internal node appears only through its rotation degrees of freedom
$ \underline{\theta^4} $
.**

It is important to remark here that **these new degrees of freedom are treaded
the same as the other degrees of freedom**. Hence, they do not add complexity when
deriving the stiffness matrix.

However, they are peculiar because they are not
attached to a node belonging to the mesh. Therefore, these
**internal degrees of freedom are not shared between adjacent elements**,
which can have useful consequences when assembling the
stiffness matrix (as discussed a bit further).

On a kinematic point of view, these additional degrees of freedom enrich the
element displacement capabilities. In fact, **they improve the element
bending behavior**.

## Basis and coordinate systems

The basis and coordinate systems used are the same as those of the MITC4+ element. They are not detailed again here, please refer to the MITC4+ article.

## MITC3+ shell element derivation

Next, the classical MITC procedure is applied to derive the element strain components. These strain components are computed at specific locations (known as tying points) inside the element, and they are interpolated to give the strain values anywhere inside the element (for instance at Gauss points).

If you are looking for the full MITC procedure employed in the MITC3+ element, you can refer to [2]. The procedure is similar to that of the MITC4+ element.

## 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 $$

As for the MITC4+ element, applying the change of variable from the local Cartesian coordinate system to the parametric coordinate system leads to:

$$ \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 $$

Where now $ V $ denotes the element volume in the parametric coordinate system.

We need to apply the **7-point Gauss integration to integrate exactly on the plane**
defined by $ r $
and $ s $
, because
of the cubic bubble function we used to express the displacement. Besides, this
integrand depends arbitrarily on $ t $
and therefore
we cannot “pre-integrate” through the thickness. We use then a **2-point Gauss
integration on $ t $
** that has proven to be accurate
enough for most applications.

Finally we can write:

$$ \boxed{ \overline{W} = \displaystyle\sum_{i=1}^{14} w_i \underbrace{\lbrack \overline{\varepsilon} \rbrack^T \lbrack D \rbrack \lbrack \varepsilon \rbrack |\underline{g_t} \cdot (\underline{g_r} \times \underline{g_s})|} _{\text{evaluated at the point }r_i, s_i, t_i} } \tag{4} $$

Where $ w_i, r_i, s_i, t_i $ are given in the following table (source G. R. Cowper):

Point | $ w_i $ | $ r_i $ | $ s_i $ | $ t_i $ |
---|---|---|---|---|

1 | 0.0629695902724 | 0.1012865073235 | 0.1012865073235 | -0.5773502691896 |

2 | 0.0629695902724 | 0.7974269853531 | 0.1012865073235 | -0.5773502691896 |

3 | 0.0629695902724 | 0.1012865073235 | 0.7974269853531 | -0.5773502691896 |

4 | 0.0661970763943 | 0.4701420641051 | 0.0597158717898 | -0.5773502691896 |

5 | 0.0661970763943 | 0.4701420641051 | 0.4701420641051 | -0.5773502691896 |

6 | 0.0661970763943 | 0.0597158717898 | 0.4701420641051 | -0.5773502691896 |

7 | 0.1125 | 0.3333333333333 | 0.3333333333333 | -0.5773502691896 |

8 | 0.0629695902724 | 0.1012865073235 | 0.1012865073235 | 0.5773502691896 |

9 | 0.0629695902724 | 0.7974269853531 | 0.1012865073235 | 0.5773502691896 |

10 | 0.0629695902724 | 0.1012865073235 | 0.7974269853531 | 0.5773502691896 |

11 | 0.0661970763943 | 0.4701420641051 | 0.0597158717898 | 0.5773502691896 |

12 | 0.0661970763943 | 0.4701420641051 | 0.4701420641051 | 0.5773502691896 |

13 | 0.0661970763943 | 0.0597158717898 | 0.4701420641051 | 0.5773502691896 |

14 | 0.1125 | 0.3333333333333 | 0.3333333333333 | 0.5773502691896 |

## Condensation of the element stiffness matrix

The element stiffness matrix is expressed on the degrees of freedom at the 3 corner nodes plus the internal node. This matrix has $ \underbrace{3 \times 6}_{\text{corner nodes}} + \underbrace{3}_{\text{internal node}} = 21$ rows and columns.

However, under certain circumstances, **we can factorize this stiffness matrix
to express it only on the degrees of freedom at the corner nodes**.
This factorization is known as the stiffness matrix condensation, and is detailed
here.

Regrouping the external degrees of freedom on one hand and the internal ones on the other hand, we can write the element stiffness matrix as follows:

$$ \lbrack K_{elem} \rbrack = \begin{bmatrix} \lbrack k_{ee} \rbrack & \lbrack k_{ei} \rbrack \\ \lbrack k_{ei} \rbrack ^T & \lbrack k_{ii} \rbrack \\ \end{bmatrix} $$

Where:

$ \lbrack k_{ee} \rbrack $ relates the external dofs to the external dofs (size $ 18 \times 18 $ ),

$ \lbrack k_{ii} \rbrack $ relates the internal dofs to the internal dofs (size $ 3 \times 3 $ ),

$ \lbrack k_{ei} \rbrack $ relates the external dofs to the internal dofs (size $ 18 \times 3 $ ).

In the same manner, we can write **the element displacement and force column matrices**
(both have 21 components) as:

$$ \lbrack X_{elem} \rbrack = \begin{bmatrix} \lbrack X_e \rbrack \\ \lbrack X_i\rbrack \\ \end{bmatrix} $$

and

$$ \lbrack F_{elem} \rbrack = \begin{bmatrix} \lbrack F_e \rbrack \\ \lbrack F_i\rbrack \\ \end{bmatrix} $$

Let’s assume we are computing the stiffness matrix in order to perform a **static
linear resolution**. For a given displacement vector on the element, the reaction
force vector is given by:

$$ \lbrack F_{elem} \rbrack = \lbrack K_{elem} \rbrack \lbrack X_{elem} \rbrack $$

And after expanding, we have the following relations:

$$ \boxed{ \begin{cases} \lbrack F_e \rbrack = \lbrack k_{ee} \rbrack \lbrack X_e \rbrack + \lbrack k_{ei} \rbrack \lbrack X_i \rbrack \\ \lbrack F_i \rbrack = \lbrack k_{ei} \rbrack^T \lbrack X_e \rbrack + \lbrack k_{ii} \rbrack \lbrack X_i \rbrack \end{cases} }\tag{5} $$

However, as mentioned previously, the internal degrees of freedom are not
attached to a node belonging to the mesh. Hence, **the internal dofs do not carry any load**.
This translates to:

$$ \lbrack F_i \rbrack = \lbrack 0 \rbrack $$

Solving the second relation of (5) for $ \lbrack X_{i} \rbrack $ gives:

$$ \lbrack X_{i} \rbrack = - \lbrack k_{ii} \rbrack^{-1} \lbrack k_{ei} \rbrack^T \lbrack X_e \rbrack $$

And inserting this result in the first relation of (5) yields:

$$ \boxed{ \lbrack F_e \rbrack = \underbrace{\lbrack k_{ee} \rbrack \lbrack X_e \rbrack - \lbrack k_{ei} \rbrack \lbrack k_{ii} \rbrack^{-1} \lbrack k_{ei} \rbrack^T} _{\lbrack K_{elem}' \rbrack} \lbrack X_e \rbrack } $$

Where $ \lbrack X_{i} \rbrack $
does no longer appear. We
finally arrive at the expression of **the condensed stiffness matrix
$ \lbrack K_{elem}' \rbrack $
relating only the external degrees of freedom**.

The main advantage of this condensation is that **it can be performed on the
element level before assembling the stiffness matrix**. Whether we apply this
condensation or not will lead to the same displacement result in the end.

Finally, **in SesamX we chose not to implement this condensation**. Indeed,it works
well in the case of linear resolution, however in non linear analysis it is more
involved and this procedure cannot be employed (because the relation between the
reactions and the degrees of freedom is non linear).

## 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: 1.,0 ,0.
- ID: 3 POINT: 1.,1.,0.
- ID: 4 POINT: 0.,1.,0.
ELEMENTS
TYPE: TR3
- ID: 1 NODES: 1,2,4
- ID: 2 NODES: 2,3,4
```

Here we define 4 nodes and we create 2 triangular 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 S3 element

Eventually, the last part of this article focuses on the comparison of SesamX
MITC3+ triangular shell element against Abaqus equivalent **S3 element**.
A comparison is
showcased here on a **distorted 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 S3.

### Model

An hyperboloid-like model is studied for this comparison.

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 every element. 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 S3 and SesamX shell.
In fact, the **element formulations are not identical**.
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 0.3%** while
the **error on the rotation is larger (at most 70%)**.

Remark:As for the S4 element, the S3 element seems to have 6 degrees of freedom per node, while the MITC3+ element has only 5 (no drilling rotation). Hence the S3 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 on 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. **From the
displacement values we see that the MITC3+ element is close to the S3 element**.

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_04 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