**4. Tetrahedral element**

The discussion of the implementation of a tetrahedral element, in particular a 4 nodes tetrahedron, allows to introduce an important ingredient of all finite element formulations: the interpolation chosen for the kinematic description. Standard approaches are hinged on the interpolation of the displacement field, in the present approach the focus is on the interpolation of the element coordinates in the reference configuration and in the current one. In the previous section regarding the truss element, this aspect remained hidden because the element elongation is easily formulated with respect to the element nodal coordinates.

Another important aspect which the tetrahedral element bring into play is the use of the continuum mechanics instruments, see [24, 25], and how these can be smoothly framed inside the proposed MATLAB® implementation.

Let us consider the geometry of the 4 node tetrahedron as illustrated in **Figure 2**. The description of the reference and current configurations of the tetrahedron are as follows.

$$\begin{aligned} \mathbf{X}(\zeta\_1, \zeta\_2, \zeta\_3, \zeta\_4) &= N\_1 \mathbf{X}\_1 + N\_2 \mathbf{X}\_2 + N\_3 \mathbf{X}\_3 + N\_4 \mathbf{X}\_4 \\ &= \zeta\_1 \mathbf{X}\_1 + \zeta\_2 \mathbf{X}\_2 + \zeta\_3 \mathbf{X}\_3 + \zeta\_4 \mathbf{X}\_4, \end{aligned} \tag{9}$$

$$\mathbf{x}(\zeta\_1, \zeta\_2, \zeta\_3, \zeta\_4) = N\_1 \mathbf{x}\_1 + N\_2 \mathbf{x}\_2 + N\_3 \mathbf{x}\_3 + N\_4 \mathbf{x}\_4 \tag{10}$$

$$= \zeta\_1 \mathbf{x}\_1 + \zeta\_2 \mathbf{x}\_2 + \zeta\_3 \mathbf{x}\_3 + \zeta\_4 \mathbf{x}\_4. \tag{11}$$

The element local coordinates *<sup>ζ</sup>* <sup>¼</sup> *<sup>ζ</sup>*<sup>1</sup> *<sup>ζ</sup>*<sup>2</sup> *<sup>ζ</sup>*<sup>3</sup> *<sup>ζ</sup>*<sup>4</sup> ½ �<sup>T</sup> are the standard tetrahedral coordinates whose definition can be found in several resources, for example [6, 8]. On this basis the description of the deformation gradient over the tetrahedron can be formulated as follows

#### **Figure 2.**

*Tetrahedral element: definition of the geometric quantities relative to the reference configuration (upper-case letters) and the current configuration (lower-case letters).*

$$F = \frac{\partial \mathbf{x}}{\partial \mathbf{X}} = \frac{\partial \mathbf{x}}{\partial \boldsymbol{\xi}} \frac{\partial \boldsymbol{\xi}}{\partial \mathbf{X}} = \frac{\partial \mathbf{x}}{\partial \boldsymbol{\xi}} \left(\frac{\partial \mathbf{X}}{\partial \boldsymbol{\xi}}\right)^{-1} = F(\mathbf{x}),\tag{11}$$

with the operator

$$\frac{\partial \mathfrak{x}}{\partial \mathfrak{f}} = \begin{bmatrix} \mathfrak{x}\_1 \mathfrak{x}\_2 \mathfrak{x}\_3 \mathfrak{x}\_4 \end{bmatrix} \tag{12}$$

containing in its columns the four coordinate vectors relative to the current configuration, unknown vectors to be expressed in MATLAB® symbolic format, and the operator

$$\frac{\partial \mathbf{X}}{\partial \boldsymbol{\xi}} = [\mathbf{X}\_1 \mathbf{X}\_2 \mathbf{X}\_3 \mathbf{X}\_4] \tag{13}$$

containing the four coordinate vectors relative to the reference configuration to be evaluated numerically for each tetrahedron of the mesh. The apparent problem represented by the evaluation of inverse *<sup>∂</sup> <sup>X</sup> ∂ ζ* �<sup>1</sup> starting from a 3 � 4 matrix is a standard matter in FEM procedures, see for example [6], and it can be easily calculated as shown in Appendix A.

The geometric formulation described above is directly inserted inside the Tetra4 class which can be implemented by following the same scheme already adopted for the class Truss. In particular the geometrical properties of the class are listed below (Listing 1.6).

**Listing 1.6**. Tetra4 class: properties.

properties (SetAccess = private)

% symbolic propetries

xl % current coordinates of node 1

x2 % current coordinates of node 2

x3 % current coordinates of node 3

x4 % current coordinates of node 4

xe % current element coordinates

Fe % deformation gradient

*A MATLAB-Based Symbolic Approach for the Quick Developing of Nonlinear Solid Mechanics… DOI: http://dx.doi.org/10.5772/intechopen.94869*

> % numeric properties Xe % reference element coordinates

end

On this basis the evaluation of the deformation gradient can be performed during the initialisation of the generic element by carrying out the following instructions (Listing 1.7).

```
Listing 1.7. Tetra4 class:function Initialize (evaluation of the deformation
            gradient).
function T = Initialize (T,D, i)
        % D brings all the problem data and its use
        % is not shown here …
        % reference configuration
        dXdzeta = [X1 X2 X3 X4];
        A = [1 1 1 1; dXdzeta];
        V = det(A)/6;
        iA = inv(A);
        dzetadX = iA(1:4 ,2:4);
        % current configuration (symbolic)
        dxdzeta = [T.xl T.x2 T.x3 T.x4];
        % deformation gradient
        T.Fe = dxdzeta*dzetadX;
        % …
end
```
It is now possibile to discuss the strain energy of the tetrahedral element. The choice is for a compressible neo-Hookean material, see [25], which allows to express the strain energy of the generic tetrahedron as follows

$$\Psi\_{\varepsilon}(\mathbf{x}) = \int\_{\Omega\_{\varepsilon}} \Psi(\mathbf{C})dV = \left(\frac{\mu}{2}(I\_1 - \mathfrak{Z}) - \mu \ln f + \frac{\lambda}{2}(\ln f)^2\right)V. \tag{14}$$

*<sup>C</sup>* <sup>¼</sup> *C x*ð Þ¼ *<sup>F</sup>*<sup>T</sup>*<sup>F</sup>* is the right Cauchy strain tensor and I1 its first invariant, J *<sup>=</sup>* det *F*, *λ* and *μ* are the Lamè parameters of the material. Thanks to the constant pattern of **F** over the element domain Ω*e*, the strain energy of the element is simply given by the product between the strain energy density and the reference volume of the element. Such a evaluation, together with the derivation of the gradient vector and Jacobian matrix is implemented inside the function Initialize as shown in Listing 1.8.

**Listing 1.8**. Tetra4 class: function Initialize (strain energy).

function T = Initialize (T,D, i ) % … C = T. Fe.'\*T.Fe; Il = trace (C); J = det (T.Fe); Psi = (mi/2\*(Il-3)-mi\*log(J) + lam/2\*log (J)^2) \*V; T.ge = gradient(Psi, T.xe);

$$\text{T.Je = jacobian(T.ge, T.xe);} \\ \text{9\% ...}$$

end

The symbolic gradient vector and Jacobian matrix evaluated in the initialization phase are then numerically computed during the analysis using a function identical to the function already presented in Listing 1.4 for the Truss class.

The basic operation of the post-processing is the computation of the Cauchy stress solution which, as **F**, is constant over the element domain. This step requires the evaluation of the second Piola-Kirchhoff stress tensor

$$\mathbf{S} = 2\frac{\partial \Psi}{\partial \mathbf{C}} = \mathbf{S}(\mathbf{C}) \tag{15}$$

and, by applying a push-forward operation to **S** [24, 25], the computation of the Cauchy stress tensor is

$$
\sigma = f^{-1} \mathbf{F} \mathbf{S}(\mathbf{C}) \mathbf{F}^{\mathrm{T}}.\tag{16}
$$

The MATLAB® implementation of Eqs. (15) requires the introduction of a symbolic matrix for **C** to be used to perform another evaluation of the strain energy depending, this time, from the components of **C**. The obtained expression, Ψð Þ *C* , can be derived in order obtain **S**. This step can be performed only one time during the initialisation of the Tetra4 class. Listing 1.9 shows these instructions together with the declaration of the necessary symbolic properties.

**Listing 1.9**. Tetra4 class: function Initialize (second Piola-Kirchhoff stress

```
tensor).
        properties (SetAccess = private)
        % …
        % symbolic properties
        Se % second Piola-Kirchhof stress tensor S(C)
        Ce % symbolic tensor C from which Se depends
end
function T = Initialize (T,D, i )
        % …
        T.Ce = sym ('C', [3, 3], 'real');
        I1 = trace (T. Ce);
        I3 = det (T.Ce);
        Psi = mi/2*(Il – 3)–mi*log ( sqrt (I3))+ …
             lam/2*log (sqrt( I3))^2;
        T.Se = reshape (2* gradient (Psi ,T.Ce (:)) ,3 ,3 );
```
#### end

Eq. (16) is used to compute the stress solution for each tetrahedron with respect to all the solutions *x* found by means of equilibrium equations (6). The class function implementing the required operations is reported in Listing 1.10.

```
Listing 1.10. Tetra4 class: function Stress.
function sig = Stress (T,gx)
        % extraction of local vector lx from global gx
        F = subs (T.Fe, T.xe, lx);
```
*A MATLAB-Based Symbolic Approach for the Quick Developing of Nonlinear Solid Mechanics… DOI: http://dx.doi.org/10.5772/intechopen.94869*

$$\begin{array}{l} \text{C = F."} ^\ast \text{F;}\\ \text{sig = F} ^\ast \text{subs (T.Se, T.Ce, C)} ^\ast \text{F.'} \text{/det(F)};\\ \text{end} \end{array}$$

The complete listing of the class can be found in [33].
