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

$$\mathbf{x}\_1 = \mathbf{x}\_1(\mathbf{X}\_1, \mathbf{X}\_2),\tag{17}$$

$$\mathbf{x}\_2 = \mathbf{x}\_2(\mathbf{X}\_1, \mathbf{X}\_2),\tag{18}$$

$$\mathbf{x}\_3 = \mathbf{X}\_3. \tag{19}$$

Consequently, the associated deformation gradient takes the following form

$$\mathbf{F} = \begin{bmatrix} F\_{11} & F\_{12} & \mathbf{0} \\ F\_{21} & F\_{22} & \mathbf{0} \\ \mathbf{0} & \mathbf{0} & \mathbf{1} \end{bmatrix}, \mathbf{F}\_{2 \times 2} \begin{bmatrix} F\_{11} & F\_{12} \\ F\_{21} & F\_{22} \end{bmatrix}. \tag{20}$$

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

$$\begin{aligned} N\_1 &= \frac{1}{4} (\mathbf{1} - \boldsymbol{\zeta\_1}) (\mathbf{1} - \boldsymbol{\zeta\_2}), N\_2 = \frac{1}{4} (\mathbf{1} + \boldsymbol{\zeta\_1}) (\mathbf{1} - \boldsymbol{\zeta\_2}) \\\\ N\_3 &= \frac{1}{4} (\mathbf{1} + \boldsymbol{\zeta\_1}) (\mathbf{1} + \boldsymbol{\zeta\_2}), N\_4 = \frac{1}{4} (\mathbf{1} - \boldsymbol{\zeta\_1}) (\mathbf{1} + \boldsymbol{\zeta\_2}) \end{aligned} \tag{21}$$

being *<sup>ζ</sup>* <sup>¼</sup> *<sup>ζ</sup>*1*ζ*<sup>2</sup> ½ �<sup>T</sup> 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

$$\mathbf{X}(\zeta\_1, \zeta\_2) = N\_1 \mathbf{X}\_1 + N\_2 \mathbf{X}\_2 + N\_3 \mathbf{X}\_3 + N\_4 \mathbf{X}\_4. \tag{22}$$

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

#### **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).*

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 *xi*ð Þ *i* ¼ 1 … 4 and numeric vectors *Xi*ð Þ *i* ¼ 1 … 4 . The deformation gradient can be evaluated using always Eq. (11) where now the operators are

$$\begin{aligned} \frac{\partial \mathbf{x}}{\partial \zeta} &= \left[ \left( -\frac{(\mathbf{1} - \zeta\_2)}{4} \mathbf{x}\_1 + \frac{(\mathbf{1} - \zeta\_2)}{4} \mathbf{x}\_2 + \frac{(\mathbf{1} + \zeta\_2)}{4} \mathbf{x}\_3 - \frac{(\mathbf{1} + \zeta\_2)}{4} \mathbf{x}\_4 \right) \\\\ &- \left( -\frac{(\mathbf{1} - \zeta\_1)}{4} \mathbf{x}\_1 - \frac{(\mathbf{1} + \zeta\_1)}{4} \mathbf{x}\_2 + \frac{(\mathbf{1} + \zeta\_1)}{4} \mathbf{x}\_3 + \frac{(\mathbf{1} - \zeta\_1)}{4} \mathbf{x}\_4 \right) \right] \end{aligned} \tag{24}$$

and

$$\begin{split} \frac{\partial \mathbf{X}}{\partial \boldsymbol{\zeta}} &= \left[ \left( -\frac{(\mathbf{1} - \boldsymbol{\zeta}\_{2})}{4} \mathbf{X}\_{1} + \frac{(\mathbf{1} - \boldsymbol{\zeta}\_{2})}{4} \mathbf{X}\_{2} + \frac{(\mathbf{1} + \boldsymbol{\zeta}\_{2})}{4} \mathbf{X}\_{3} - \frac{(\mathbf{1} + \boldsymbol{\zeta}\_{2})}{4} \mathbf{X}\_{4} \right) \right. \\ & \left. \left( -\frac{(\mathbf{1} - \boldsymbol{\zeta}\_{1})}{4} \mathbf{X}\_{1} - \frac{(\mathbf{1} + \boldsymbol{\zeta}\_{1})}{4} \mathbf{X}\_{2} + \frac{(\mathbf{1} + \boldsymbol{\zeta}\_{1})}{4} \mathbf{X}\_{3} + \frac{(\mathbf{1} - \boldsymbol{\zeta}\_{1})}{4} \mathbf{X}\_{4} \right) \right] \end{split} \tag{25}$$

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)

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

```
nG = 4;
       xiG = % Gauss coordinates , values not shown here
       wG = [1 1 1 1];
end
properties (SetAccess = private)
       % symbolic properties
       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 F(xe) (nGP times)
       % numeric properties
       Xe % reference element coordinates
end
function PF = Initialize (PF,D, i)
       % D brings all the problem data and its use
       % is not shown here
       PF . Fe = sym( zeros (2 ,2 ,PF .nG));
       for g = l:PF.nG
         % dzetadX evaluation in g
         % …
         % dxdzeta evaluation in g
         % …
         % F in g
         F = dxdzeta * dzetadX;
         PF.Fe(:, :, g) = F;
         % …
       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 = \det \mathbf{F} = \det \mathbf{F}\_{2 \times 2},\tag{26}$$

and by tensor *C*

$$\mathbf{C} = \begin{bmatrix} \mathbf{C}\_{11} & \mathbf{C}\_{12} & \mathbf{0} \\\\ \mathbf{C}\_{21} & \mathbf{C}\_{22} & \mathbf{0} \\\\ \mathbf{0} & \mathbf{0} & \mathbf{1} \end{bmatrix}, \quad \mathbf{C}\_{2 \times 2} = \begin{bmatrix} \mathbf{C}\_{11} & \mathbf{C}\_{12} \\\\ \mathbf{C}\_{21} & \mathbf{C}\_{22} \end{bmatrix}, \tag{27}$$

from which

$$I\_1 = \text{tr}\mathbf{C} = \text{tr}\mathbf{C}\_{22} + \mathbf{1}.\tag{28}$$

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

$$
\Psi\_{\rm PF} = \frac{\mu}{2}(I\_1 - 2) - \mu \ln f + \frac{\lambda}{2}(\ln f)^2,\tag{29}
$$

where *I*<sup>1</sup> 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

$$\Psi\_{\epsilon}(\mathbf{x}) = \int\_{\Omega\_{\epsilon}} \Psi\_{\text{PF}} \, dV = \sum\_{\mathbf{g}=1}^{4} [\Psi\_{\text{PF}}] A\_{\mathbf{g}} \, th \, w\_{\mathbf{g}} = \sum\_{\mathbf{g}=1}^{4} \Psi\_{\text{g}},\tag{30}$$

where *Ag* <sup>¼</sup> det *<sup>∂</sup><sup>X</sup> ∂ζ* h i *g* is the part of the reference domain pertaining to the Gauss

point, *wg* 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

$$\mathbf{g}\_{\epsilon} = \sum\_{\mathbf{g}=1}^{4} \frac{\partial \Psi\_{\mathbf{g}}}{\partial \mathbf{x}} = \sum\_{\mathbf{g}=1}^{4} \mathbf{g}\_{\mathbf{g}}, \quad J\_{\epsilon} = \sum\_{\mathbf{g}=1}^{4} \frac{\partial \mathbf{g}\_{\mathbf{g}}}{\partial \mathbf{x}} = \sum\_{\mathbf{g}=1}^{4} J\_{\mathcal{S}}, \tag{31}$$

where *g<sup>g</sup>* and *J<sup>g</sup>* 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)
        % …
        PF. Je = sym (zeros (8 ,8 , PF .nG));
        PF. Fe = sym ( zeros (2 ,2 ,PF .nG));
        for g = l:PF.nG
        % …
        C = F.'*F;
        I1 = trace (C);
        J = det (F);
        Psi = (mi/2*(Il–2)–mi*log(J) + lam/2*log (J)^2) …
        A*PF.wG(g)*th;
        PF.ge(: , 1, g) = gradient (Psi, PF. xe);
        PF.Je(:, :, g) = jacobian (PF. ge(: , 1, g), PF. xe);
        % …
        end
        % …
```
During the analysis the main task to be performed by the element is the numerical evaluation of *g<sup>e</sup>* and *J<sup>e</sup>* that now must to be performed, see Eq. (31), on the basis of the following implementation of the class function Compute.

end

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

**Listing 1.13**. PF4 class:function Compute. function PF = Compute(PF) PF. se = zeros (8, 1); PF. Ke = zeros (8, 8); for g = l:PF.nG PF.se = PF.se + subs (PF.ge(:, 1, g), PF.xe, PF.xxe); PF.Ke = PF.Ke + subs (PF.Je(:, :, g), PF.xe, PF.xxe); 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 *σ*<sup>33</sup> 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 *σ*<sup>33</sup> component, stems from the plane solution, and is given by

$$
\sigma\_{33} = J^{-1} \mathbf{S}\_{33} = J^{-1} \frac{\lambda}{2} \ln \left( \det \mathbf{C} \right) \tag{32}
$$

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)
```

```
% retrieve local vector lx from global
        % solution gx
        sig = zeros(3, 3, PF.nG);
        for g = l:PF.nG
          F = subs (PF.Fe(:, :, g), PF.xe, lx);
          C = F.'*F;
          sig(l:2, l:2, g) = F*subs (PF. Se, PF.Ce,C)*F.'/det(F);
          sig (3, 3, g) = subs (PF.Se33, PF.Ce, C)/det (F);
        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].
