**1. Introduction**

The developing of finite element formulations, standard or new ones, requires a lengthy process which involves several steps.

The typical starting point is the formulation of a mathematical model where the main physical or real world phenomena to be described are established. At present time the definition of a mathematical model is at the basis of any serious attempt to obtain previsions in any engineering application [1–3], but not only in the engineering field [4, 5].

The subsequent step is the introduction of a numerical approximation technique. The most popular technique is the Finite Element Method (FEM), see [6–8], but now the number of computational approximation techniques is very large and a synthetic summary can be tried only by citing some of these less conventional methods: mixed finite element methods [9–13]; partition of unity-based discontinuous finite elements [14, 15]; meshless methods [16]; discontinuous Galerkin methods [17]. This operation leads to the identification of the needed discrete operators which define the computational model. For example, in the case of the analysis of solid mechanics problems by FEM, important discrete operators are the mechanical response vector of the finite element and its tangent stiffness matrix. This phase is characterised by the evaluation and the analysis of these operators

through an algebraic manipulator such as MATLAB® [18]. MATLAB® is certainly one of the state-of-the-art mathematical softwares available for performing numeric or symbolic analyses, but it is not the only one and a quite long list, see [19], of packages offering very similar features is available.

The discrete model so defined is usually inserted into a prototype code, often by using again an algebraic manipulator but the use of compiled programming languages is also possible if not common. This prototype code allows to perform basic tests with the aim to check the effectiveness of the adopted model with respect to well known situations and to check for the presence of bugs. Often this phase can highlight also flaws in the mathematical model or in the discretization technique. In any case it is necessary to go back and to repeat the process just described.

The additional last step can be the production of an executable by using compiled programming languages such as C/C++ or fortran. This makes possible to extend the validation of the conceived numerical model by performing the analysis of larger sized problems.

As already said, the previously described work-flow is lengthy and it is often characterised by a gap between the theoretical formulation and its implementation in a numerical code. However some solutions capable to assist the developer in this process already exist ant it is worth to mention some ones. Such solutions typically refer to a specific context by keeping fixed the physical problem but letting open the specific instance of discretization technique which, however, is fixed too. This is the case of open source FEM libraries or commercial packages listed in [20]. In both cases the user must define the procedures or functions needed in order to assign the desired new finite element to be used within the analysis framework already available in the library. Other packages instead solves a generic system of Partial Differential Equations (PDEs) subjected to boundary and initial conditions. Inside this generic form the specific differential problem to be solved must to be fitted by the user, see for example [21, 22], but usually with no control over the discretization technique used by the solver.

The present work, quite far from being an alternative to the hugely developed and rich packages previously cited, proposes a basic approach for the quick early phase developing of solid mechanics finite elements formulation. Its intent is to show how to use MATLAB®, in particular by exploiting the capabilities of the Symbolic Math Toolbox™ [23], to produce numerical approximations of a given solid mechanics problem in a way that the usual gap between the theoretical formulation and its actual implementation in a code is not perceived. This result is obtained by condensing the development process going from the mathematical formulation to the prototype code implementation in a few lines of MATLAB® symbolic instructions. The applicative context is the nonlinear analysis of solids and structures, see [24, 25], by showing the formulation and the subsequent MATLAB® coding of some typical structural and solid finite elements. The mechanical formulation is based on the kinematics of finite deformation and, for the description of the material behaviour, on isotropic hyperelasticity, i.e. the stress solution is found as a derivative of some potential energy function. This allows to express the mechanical problem at hand in terms of the stationary condition of the Total Potential Energy (TPE). The stationary condition, assuming a fem discretization of the given domain, is then easily translated into a nonlinear algebraic problem whose unknowns are the position vectors of the nodes of the mesh used in the discretization.

The following finite elements are discussed: a truss element, a 3D tetrahedral element with four nodes and a four nodes quadrangular element subjected to plane strain condition. The choice allows to discuss gradually the main ingredients present in a finite element formulation and how these can be framed inside the proposed

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

MATLAB® approach. The latter is based on the definition of MATLAB® classes which share the same structuring and which differ only for the particular mechanical response to be implemented. In particular the generic class is structured as shown by the following instructions (Listing 1.1).

```
Listing 1.1. Generic class.
classdef Element
        properties (SetAccess = private)
          % symbolic properties
          % numeric properties
        end
        methods
        function E = Element()
        end
        function E = Initialize (E,D, i)
        end
        function E = Compute(E)
        end
        function sig = Stress(E,gx)
        end
        end
end
```
The properties section contains a group of symbolic properties devoted to handle the unknowns and the quantities depending on them used in the description of the element mechanical behaviour. The other group of numeric properties are used to handle quantities that are known and then they can have a numeric value. Beyond the constructor, that must to be present in any class, we have the function *Initialize,* belonging to the *pre-processing phase* of a FEM code, whose main task is the inizialization of the i-th element on the basis of the assigned data structure D. This is the moment also for evaluating the element operators, in a symbolic format, needed to the analysis. In the subsequent *analysis phase,* the function Compute evaluates the numeric instances of the symbolic operators previously prepared. The function Stress is typical of the *post-processing phase* of any FEM code and its task is to compute the stress solution inside the generic element starting from the kinematic global solution represented by vector gx.

Before proceeding with the description of the proposed work, it is noteworthy to observe that the use of MATLAB® to perform mathematical and numeric analyses is not certainly new and several books are dedicated to this subject, see [26–28] just to cite a few. Moreover the already cited book [24] employs MATLAB® for the implementation of a FEM software. However the present less comprehensive work is different because it carry out the formulation of the FEM operators by exploiting the potentiality of the symbolic manipulator and advising an object-oriented programming style.

A last further annotation regards the use of the symbolic approach which, with respect to the expected performance of final codes, represents a weakness. This aspect however is to be considered less important in a work regarding the early phase developing of a FEM formulation. Nevertheless techniques, [29–31], for the automatic generation of efficient and highly compressed code is a research theme

which is attracting increasingly interest, making viable the up-scaling of the proposed approach.

The chapter is organised as follows. Section 2 presents the FEM formulation of the Total Potential Energy for a generic structure or solid, showing also the evaluation of the gradient needed to define the discrete equilibrium equations and the evaluation of the Jacobian necessary for their solution. Sections 3, 4 and 5 describes, respectively, the truss element, the tetrahedral element and plane quadrangular element. The closing section furnishes some additional final comments.

### **2. Total Potential Energy**

An effective description of a generic mechanical problem can be obtained through the stationary condition of its Total Potential Energy (TPE) which, see for example [24], can be expressed as follows

$$\prod(\mathfrak{x}) \equiv \prod\_{\mathrm{int}}(\mathfrak{x}) + \prod\_{\mathrm{ext}}(\mathfrak{x}) = \mathrm{stat.} \tag{1}$$

*x* is the global vector of the current positions of the nodal points defining the mesh used to describe the geometry of the solid. Q intð Þ *x* , excluding dynamic and dissipative effects, is given only by the strain energy obtained by summing all the contribution from all the finite elements, i. e.

$$\prod\_{\text{int}}(\mathfrak{x}) = \sum\_{\text{e}} \Psi\_{\text{e}}(\mathfrak{x}),\tag{2}$$

Q being Ψeð Þ *x* the hyperelastic strain energy relative to the generic finite element. extð Þ *x* is the potential energy of external forces. For simplicity the case of a solid body subjected only to external punctual forces will be considered, in this case the potential energy can be written as

$$\prod\_{\text{ext}}(\mathfrak{x}) = -\mathfrak{f} \cdot \mathfrak{x},\tag{3}$$

where *f* is the global vector of the applied forces in each node of the mesh. *f* has the same length of *x*, however it is mainly composed by null entries.

On this basis the equilibrium equations can be easily formulated with respect the degrees of freedom involved in the FEM description of the body. In particular, by imposing the stationary condition (1), the equilibrium equations can be derived, obtaining

$$\mathbf{A}\_{\mathbf{e}} \mathbf{g}\_{\mathbf{e}}(\mathbf{x}) - \mathbf{f} = \mathbf{0}.\tag{4}$$

where the assembly operator A is used to build up the global response vector by using the gradient vector of each finite element strain energy contribution, i.e.

$$\mathbf{g}\_{\mathbf{e}}(\mathfrak{x}) = \frac{\partial \Psi\_{\mathbf{e}}(\mathfrak{x})}{\partial \mathfrak{x}}.\tag{5}$$

The solution of Eq. (4), a typically nonlinear algebraic system whose unknowns are the components of vector *x*, is based on a Newton–Raphson iteration which can be formulated as follows

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

$$\mathbf{A}\_{\mathbf{e}}(\mathbf{x}\_{j}) \left(\boldsymbol{\omega}\_{j+1} - \boldsymbol{\omega}\_{j}\right) = -\left(\mathbf{A}\_{\mathbf{e}}(\boldsymbol{\omega}\_{j}) - \boldsymbol{\mathcal{f}}\right),\tag{6}$$

where *x<sup>j</sup>* and *x <sup>j</sup>*þ<sup>1</sup> are the estimated solutions at *j*-th and (*j* + 1)-th iterations and *J*eð Þ *x* is the Jacobian matrix of the finite element given by

$$J\_{\mathbf{e}}(\mathbf{x}) = \frac{\partial \mathbf{g}\_{\mathbf{e}}(\mathbf{x})}{\partial \mathbf{x}}.\tag{7}$$

The gradient vector *g*eð Þ *x* and the Jacobian matrix *J*eð Þ *x* can be used as basic building blocks for the finite element formulation. This is the approach at the basis of the MATLAB® implementations to be described in the following sections.

#### **3. Truss element**

The strain energy of the truss element is defined, see [24], as follows

$$\Psi\_{\mathbf{e}}(\mathbf{x}) = \frac{1}{2} E \varepsilon^2 V, \varepsilon = \varepsilon(\mathbf{x}) = \ln\left(\frac{l}{L}\right),\tag{8}$$

where *E* is the Young modulus, *L* and *V* are the length and the volume of the bar in the reference configuration, *l* is the length of the bar in the current configuration. The geometric quantities just described are depicted in **Figure 1** where the coordinate vectors of the nodal points are also shown.

The implementation of the MATLAB® class Truss can stem from the properties reported in Listing 1.2. Some of them, those describing the reference configuration, can be numeric because are fixed. The other properties, which describe the current configuration, are expressed in symbolic form in order to be used as quantities whose the strain energy of truss element depends on.

```
Listing 1.2. Truss class: properties.
properties (SetAccess = private)
        % symbolic properties
        xa % current coordinates of node a
        xb % current coordinates of node b
        xe % current element coordinates
        ge % gradient g (xe)
        Je % Jacobian J (xe)
        eps % strain eps (xe)
        N % axial force N (xe)
        % numeric properties
        a % node a global index
        b % node b global index
        Xe % reference element coordinates
        % …
end
```
During the pre-processing phase the numeric properties of the class are appropriately assigned and the symbolic properties are evaluated as shown in the following listing (Listing 1.3). This happens inside the function Initialize belonging to the methods of the class.

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

**Listing 1.3**. Truss class: function Initialize.

function T = Initialize (T,D, i) % D brings all the problem data and its use % is not shown here T. xa = sym ('xa', [3 l], 'real'); T. xb = sym ('xb', [3 1], 'rea1'); T. xe = [T. xa; T. xb]; % reference configuration L = dot(Xb–Xa, Xb–Xa); L = sqrt(L); % current configuration 1 = dot (T. xb–T. xa, T. xb–T. xa); 1 = sqrt(1); T. eps = log(1/L); Psi = 1/2 \* E \* T.eps^2 \* (L\*A); T. ge = gradient(Psi, T.xe); T. Je = jacobian (T.ge, T.xe); % preliminary symbolic evaluation of N T .N = E\*A\*T. eps;

end

Listing 1.3 shows that, after the computation of the strain energy using Eq. (8), symbolic properties *g<sup>e</sup>* and *J<sup>e</sup>* are evaluated on the basis of Eqs. (5) and (7), respectively, by simply calling the function gradient and the function jacobian both belonging to the Symbolic Math Toolbox™. This highlights the *short distance* between the formulation and its code implementation.

After having prepared each Truss object in the way described above, it is possible to evaluate, whenever it is needed during the solution of the nonlinear equilibrium equations, the gradient and the Jacobian of the generic element with respect to estimated solution *x<sup>j</sup>* , see Eq. (6), by calling the following class method (Listing 1.4).

**Listing 1.4**. Truss class: function Compute. function T = Compute(T)

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

```
T.se = subs (T.ge, T.xe, T.xxe);
T.Ke = subs (T.Je, T.xe, T.xxe);
```
end

The function Compute uses the numeric property T.xxe previously filled with the current nodal coordinates values and it stores the resulting numeric expressions of the gradient and Jacobian in the class properties se and Ke. The desired result is obtained by calling MATLAB® function subs which substitutes the symbolic variable T.xe with its numeric value T.xxe.

The post-processing phase of any FEM codes certainly comprehends the evaluation of the stress solution. In the case of the truss element, the axial force must to be computed with respect to each vector *x* calculated by the solution of the equilibrium Eqs. (4). Listing 1.5 shows the very simple function implementing the required computation.

```
Listing 1.5. Truss class: function Stress.
function N = Stress (T,gx)
        % extraction of local vector lx from global gx
        N = subs(T.N, T.xe, lx);
end
```
The complete listing of the class can be found in [32].
