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, 2, 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, 7, 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, 10, 11, 12, 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 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
1. properties (SetAccess = private)
2. % symbolic properties
3. % numeric properties
4. end
5. methods
6. function E = Element()
7. end
8. function E = Initialize (E,D, i)
9. end
10. function E = Compute(E)
11. end
12. function sig = Stress(E,gx)
13. end
14. 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, 27, 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, 30, 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

∏x≡∏intx+∏extx=stat.E1

x is the global vector of the current positions of the nodal points defining the mesh used to describe the geometry of the solid. ∏intx, 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.

∏intx=∑eΨex,E2

being Ψex the hyperelastic strain energy relative to the generic finite element. ∏extx 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

∏extx=−f⋅x,E3

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

Aegex−f=0.E4

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.

gex=∂Ψex∂x.E5

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

AeJexjxj+1−xj=−Aegexj−f,E6

where x_j and xj+1 are the estimated solutions at j-th and (j + 1)-th iterations and Jex is the Jacobian matrix of the finite element given by

Jex=∂gex∂x.E7

The gradient vector gex and the Jacobian matrix Jex 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

Ψex=12Eε2V,ε=εx=lnlL,E8

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.

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

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)
1. % symbolic properties
2. xa % current coordinates of node a
3. xb % current coordinates of node b
4. xe % current element coordinates
5. ge % gradient g (xe)
6. Je % Jacobian J (xe)
7. eps % strain eps (xe)
8. N % axial force N (xe)
9. % numeric properties
10. a % node a global index
11. b % node b global index
12. Xe % reference element coordinates
13. % …
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.

Listing 1.3. Truss class: function Initialize.

function T = Initialize (T,D, i)
1. % D brings all the problem data and its use
2. % is not shown here
3. T. xa = sym ('xa', [3 l], 'real');
4. T. xb = sym ('xb', [3 1], 'rea1');
5. T. xe = [T. xa; T. xb];
6. % reference configuration
7. L = dot(Xb–Xa, Xb–Xa);
8. L = sqrt(L);
9. % current configuration
10. 1 = dot (T. xb–T. xa, T. xb–T. xa);
11. 1 = sqrt(1);
12. T. eps = log(1/L);
13. Psi = 1/2 * E * T.eps^2 * (L*A);
14. T. ge = gradient(Psi, T.xe);
15. T. Je = jacobian (T.ge, T.xe);
16. % preliminary symbolic evaluation of N
17. 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_e and J_e 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_j, see Eq. (6), by calling the following class method (Listing 1.4).

Listing 1.4. Truss class: function Compute.

function T = Compute(T)
1. T.se = subs (T.ge, T.xe, T.xxe);
2. 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)
1. % extraction of local vector lx from global gx
2. N = subs(T.N, T.xe, lx);
end

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

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.

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

Xζ1ζ2ζ3ζ4=N1X1+N2X2+N3X3+N4X4=ζ1X1+ζ2X2+ζ3X3+ζ4X4,E9

xζ1ζ2ζ3ζ4=N1x1+N2x2+N3x3+N4x4=ζ1x1+ζ2x2+ζ3x3+ζ4x4.E10

The element local coordinates ζ=ζ1ζ2ζ3ζ4T 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

F=∂x∂X=∂x∂ζ∂ζ∂X=∂x∂ζ∂X∂ζ−1=Fx,E11

with the operator

∂x∂ζ=x1x2x3x4E12

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

∂X∂ζ=X1X2X3X4E13

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 ∂X∂ζ−1 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)
1. % symbolic propetries
2. xl % current coordinates of node 1
3. x2 % current coordinates of node 2
4. x3 % current coordinates of node 3
5. x4 % current coordinates of node 4
6. xe % current element coordinates
7. Fe % deformation gradient
8. % numeric properties
9. 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)
1. % D brings all the problem data and its use
2. % is not shown here …
3. % reference configuration
4. dXdzeta = [X1 X2 X3 X4];
5. A = [1 1 1 1; dXdzeta];
6. V = det(A)/6;
7. iA = inv(A);
8. dzetadX = iA(1:4 ,2:4);
9. % current configuration (symbolic)
10. dxdzeta = [T.xl T.x2 T.x3 T.x4];
11. % deformation gradient
12. T.Fe = dxdzeta*dzetadX;
13. % …
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

Ψex=∫ΩeΨCdV=μ2I1−3−μlnJ+λ2lnJ2V.E14

C=Cx=FTF is the right Cauchy strain tensor and I₁ its first invariant, J = 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 )
1. % …
2. C = T. Fe.'*T.Fe;
3. Il = trace (C);
4. J = det (T.Fe);
5. Psi = (mi/2*(Il-3)-mi*log(J) + lam/2*log (J)^2) *V;
6. T.ge = gradient(Psi, T.xe);
7. T.Je = jacobian(T.ge, T.xe);
8. % …
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

S=2∂Ψ∂C=SCE15

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

σ=J−1FSCFT.E16

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)
1. % …
2. % symbolic properties
3. Se % second Piola-Kirchhof stress tensor S(C)
4. Ce % symbolic tensor C from which Se depends
end
function T = Initialize (T,D, i )
1. % …
2. T.Ce = sym ('C', [3, 3], 'real');
3. I1 = trace (T. Ce);
4. I3 = det (T.Ce);
5. Psi = mi/2*(Il – 3)–mi*log ( sqrt (I3))+ …
  1. lam/2*log (sqrt( I3))^2;
6. 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)
1. % extraction of local vector lx from global gx
2. F = subs (T.Fe, T.xe, lx);
3. C = F.'*F;
4. sig = F*subs (T.Se, T.Ce, C)*F.'/det(F);
end

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

5. Plane strain 4 nodes element

The use of finite elements specifically formulated for the analysis of problems which admit a 2D reduction is very common and quadrangular elements play an important role in the case of simple geometries. In this section a 4 nodes quadrangular element subjected to plain strain condition is discussed. The element is very basic but it allows to discuss also the use of the Gauss integration points in the calculation of the required FEM operators. The use of the Gauss integration point is an important cornerstone for all finite element formulations.

Plane strain condition stems from the following assumption on the transformation defining the new configuration of each point of the body

x1=x1X1X2,E17

x2=x2X1X2,E18

x3=X3.E19

Consequently, the associated deformation gradient takes the following form

F=F11F120F21F220001,F2×2F11F12F21F22.E20

Eqs. (18)–(20) allow the dealing with a 2D kinematic description. The stress solution, however, is not strictly plane because Eq. (19) constitutes an internal constraint determining also the presence of the component σ33. This component anyway depends only from the 2D kinematic solution as it will be shown in the following.

The standard shape function of the four nodes plane element are

N1=141−ζ11−ζ2,N2=141+ζ11−ζ2N3=141+ζ11+ζ2,N4=141−ζ11+ζ2E21

being ζ=ζ1ζ2T the element local coordinates used for quadrangular elements, see for example [6]. Shape functions (21) can be properly used to describe, see Figure 3, the reference configuration and the current configuration of the element giving

Figure 3.
Four nodes plane element: definition of the geometric quantities relative to the reference configuration (uppercase letters) and the current configuration (lower-case letters).

Xζ1ζ2=N1X1+N2X2+N3X3+N4X4.E22

xζ1ζ2=N1x1+N2x2+N3x3+N4x4.E23

We have exactly the same pattern of the tetrahedron element, see Eqs. (9) and (10), except for the meaning of the shape function and the 2D dimension of the symbolic vectors xii=1…4 and numeric vectors Xii=1…4. The deformation gradient can be evaluated using always Eq. (11) where now the operators are

∂x∂ζ=−1−ζ24x1+1−ζ24x2+1+ζ24x3−1+ζ24x4−1−ζ14x1−1+ζ14x2+1+ζ14x3+1−ζ14x4E24

and

∂X∂ζ=−1−ζ24X1+1−ζ24X2+1+ζ24X3−1+ζ24X4−1−ζ14X1−1+ζ14X2+1+ζ14X3+1−ζ14X4E25

are 2 × 2 matrices depending on the local coordinates of the element. Then the necessity to use the Gauss integration points in the evaluation of the strain energy of the element and, as a consequence, of the gradient and Jacobian of the element, see Eqs. (5) and (7). In particular four Gauss points are used, their coordinates and weights can be found in any FEM text book and are also shown in the complete listing of the class available in [34].

Previous discussion introduces the implementation details of the MATLAB® class PF4, PF stays for Plane F, whose kinematic properties, see Listing 1.11, are similar to those used for the class Tetra4 plus other properties required for the Gauss integration points. These properties are used to implement Eqs. (11), (24) and (25) which must to be evaluated in each Gauss point (bulky details are not shown but they can be found in the complete listing of the class, see [34]).

Listing 1.11. PF4 class: kinematic properties and evaluation of F.

properties (Constant)
1. nG = 4;
2. xiG = % Gauss coordinates , values not shown here
3. wG = [1 1 1 1];
end
properties (SetAccess = private)
1. % symbolic properties
2. xl % current coordinates of node 1
3. x2 % current coordinates of node 2
4. x3 % current coordinates of node 3
5. x4 % current coordinates of node 4
6. xe % current element coordinates
7. Fe % deformation gradient F(xe) (nGP times)
8. % numeric properties
9. Xe % reference element coordinates
end
function PF = Initialize (PF,D, i)
1. % D brings all the problem data and its use
2. % is not shown here
3. PF . Fe = sym( zeros (2 ,2 ,PF .nG));
4. for g = l:PF.nG
  1. % dzetadX evaluation in g
  2. % …
  3. % dxdzeta evaluation in g
  4. % …
  5. % F in g
  6. F = dxdzeta * dzetadX;
  7. PF.Fe(:, :, g) = F;
  8. % …
5. end
end

In each Gauss integration point the strain energy, the compressible neo-Hookean form is used again, must to be evaluated by taking into account the simplification determined by the plane form assumed by tensor F, then

J=detF=detF2×2,E26

and by tensor C

C=C11C120C21C220001,C2×2=C11C12C21C22,E27

from which

I1=trC=trC22+1.E28

Then the expression of the strain energy density valid for the plane strain condition is

ΨPF=μ2I1−2−μlnJ+λ2lnJ2,E29

where I₁ and J are calculated on the basis of the plane form of kinematic tensors. The resulting strain energy of the generic element can be then evaluated by using the following formula

Ψex=∫ΩeΨPFdV=∑g=14ΨPFAgthwg=∑g=14Ψg,E30

where Ag=det∂X∂ζg is the part of the reference domain pertaining to the Gauss point, w_g is the Gauss point weight and th is the domain thickness usually assumed unitary under plane strain condition. Using Eq. (30), (5), and (7) the following results are valid for the generic element

ge=∑g=14∂Ψg∂x=∑g=14gg,Je=∑g=14∂gg∂x=∑g=14Jg,E31

where g_g and J_g are the gradient and Jacobian, respectively, pertaining to the generic Gauss point.

The following MATLAB® instructions, Listing 1.12, implements, inside the function Initialize of PF4 class, the operations required by Eq. (31).

Listing 1.12. PF4 class: evaluation of ge and Je.

function PF = Initialize (PF,D, i)
1. % …
2. PF. Je = sym (zeros (8 ,8 , PF .nG));
3. PF. Fe = sym ( zeros (2 ,2 ,PF .nG));
4. for g = l:PF.nG
  1. % …
  2. C = F.'*F;
  3. I1 = trace (C);
  4. J = det (F);
  5. Psi = (mi/2*(Il–2)–mi*log(J) + lam/2*log (J)^2) …
    1. A*PF.wG(g)*th;
  6. PF.ge(: , 1, g) = gradient (Psi, PF. xe);
  7. PF.Je(:, :, g) = jacobian (PF. ge(: , 1, g), PF. xe);
  8. % …
5. end
6. % …
end

During the analysis the main task to be performed by the element is the numerical evaluation of g_e and J_e that now must to be performed, see Eq. (31), on the basis of the following implementation of the class function Compute.

Listing 1.13. PF4 class:function Compute.

function PF = Compute(PF)
1. PF. se = zeros (8, 1);
2. PF. Ke = zeros (8, 8);
3. for g = l:PF.nG
  1. PF.se = PF.se + subs (PF.ge(:, 1, g), PF.xe, PF.xxe);
  2. PF.Ke = PF.Ke + subs (PF.Je(:, :, g), PF.xe, PF.xxe);
4. end
end

The last part of the class to be discussed regards the evaluation of stress solution. As already observed in the beginning of this section, Eq. (19) constitutes an internal constraint determining the presence of also the stress component σ33 to be evaluated together with the plane part of the stress tensor. The plane part can be calculated using Eqs. (15) and (16) where the plane version of C and F must be used starting from the strain energy expression given by Eq. (29). The σ33 component, stems from the plane solution, and is given by

σ33=J−1S33=J−1λ2lndetCE32

A simple derivation of this expression through MATLAB® is reported in Appendix A. The implementation of the operations required for the evaluation of the stress solution are reported below, Listing 1.14.

Listing 1.14. PF4 class: function Stress.

function sig = Stress (T, gx)
1. % retrieve local vector lx from global
2. % solution gx
3. sig = zeros(3, 3, PF.nG);
4. for g = l:PF.nG
  1. F = subs (PF.Fe(:, :, g), PF.xe, lx);
  2. C = F.'*F;
  3. sig(l:2, l:2, g) = F*subs (PF. Se, PF.Ce,C)*F.'/det(F);
  4. sig (3, 3, g) = subs (PF.Se33, PF.Ce, C)/det (F);
5. end
end

Listing 1.14 shows the use of the symbolic properties PF.Se which is initialised in way similar to the property T.Se shown in Listing 1.9 for the tetrahedral element. Anyway the complete listing of the class can be found in [34].

6. Conclusions

The early phase developing of finite elements can be a lengthy and error prone processes involving the use of different tools. The MATLAB® symbolic approach here presented can be effectively used to test a produce new finite element formulation reducing a lot the distance between the formulation and its actual implementation. In order to be more illustrative the presentation regarded basic solid mechanics finite elements, a truss, tetrahedral and plane quadrangular element, but the developing of finite elements for more specific engineering applications is an objective worth to be pursued and it is the subject of the author’s current work.

The weakness of the proposed approach is the low performance of the final codes making difficult the analysis of real sized problems by using common hardware resources which, however, are adequate if small but significative test cases are chosen. A workaround, already tested by the author but not presented here, is the generation and storing on files of MATLAB® functions for the evaluation of the element operators. This must happens before, and one time for all, the execution of the analysis. The MATLAB® functions so obtained can be called during the analysis for evaluating the required finite element operators avoiding the calls to time-consuming function subs. Anyway the tuning of this operation is less automatic because the generation of the required MATLAB® functions can be, depending on the size of the operator to be translated into a MATLAB® function, time consuming, specially if the optimization flag is active. Then techniques quite common in the field of the symbolic and /or algorithmic differentiation should be exploited for the most intricate cases.

A.

A.1 Tetrahedron reference configuration operator inversion

The problem of the evaluation of the inverse of matrix ∂X∂ζ present in Eq. (11) is circumvented by evaluating the Jacobian of the following system of equations

1=ζ1+ζ2+ζ3+ζ4Xζ1ζ2ζ3ζ4=ζ1X1+ζ2X2+ζ3X3+ζ4X4E33

whose linearisation gives

0dX=1111X1X2X3X4dζ.E34

By inverting this relationship, i. e.

dζ=1111X1X2X3X4−10dX=−∂ζ1∂X1∂ζ1∂X2∂ζ1∂X3−∂ζ2∂X1∂ζ2∂X2∂ζ2∂X3−∂ζ3∂X1∂ζ3∂X2∂ζ3∂X3−∂ζ4∂X1∂ζ4∂X2∂ζ4∂X30dX.E35

the evaluation of the desired 4×3 matrix, ∂X∂ζ−1, is obtained. Moreover the volume of the tetrahedron in its reference configuration is an additional result thanks to relationship

6V=det1111X1X2X3X4.E36

A.2 Out-of-plane normal component for the plane strain condition

The following MATLAB® instructions allow to find an explicit expression of the S₃₃ component, from which σ33=J−1S33.

syms lam mi 'real'
syms C [3 3] 'real'
C(l, 3) = 0; C(3, l) = 0; C(3, 2) = 0; C(2, 3) = 0;
I1 = trace(C); I3 = det(C);
Psi = mi/2*(I1–3)–mi*log(sqrt (I3))+ …
1. lam/2*log (sqrt (13))^2;
S33 = simplify (2*diff (Psi ,C(3, 3)));
S33 = subs(S33, C(3, 3), 1);

Stephen Itoga

Chapters authored

Related collaborators

Abdelrahman Elfaki

Christine Chan

Ahmad Shahi

Claire Laudy

Rodziah Atan

Saravanan Muthaiyah

Chin-Kuan Ho

Somnuk Phon-Amnuaisuk

Simon Suigen Guo

Robert Harrison