MITC3 shell element explained

Learn what are the specific features of the MITC3 shell element. Discover how the SesamX element compares to Abaqus S3 element.

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.

2D mesh

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 Mixed Interpolation of Tensorial Components (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:

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

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

Shell element

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.

Shell parametric coordinates

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.

Bubble shape function

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  
    SHAPE: 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.

Hyperboloid mesh Hyperboloid loading
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.

Model displacement

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

SesamX and Abaqus comparison

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