Chapter 8

Theory and Implementation of the Finite Element Method

 

 

The derivation and implementation of the finite element method outlined in the previous chapter is simple and easy to follow, but it gives the misleading impression that the finite element method relies on the principle of minimum potential energy in linear elasticity, and so is applicable only to linear elastic solids.  This is not the case, of course $–$ the finite element method can solve problems involving very complex materials and large shape changes.

 

This chapter contains

 

1. A more rigorous derivation of the finite element equations, based on the principle of virtual work, which was derived in Section 2.4.

 

2. A more sophisticated implementation of finite element method for static linear elasticity. This includes more accurate element interpolation schemes, and also extends the finite element method to three dimensions.

 

3. A discussion of time integration schemes that are used in finite element simulations of dynamic problems, and a discussion of modal techniques for dynamic linear elasticity problems;

 

4. An extension of the finite element method to nonlinear materials, using the hypoelastic material model described in Section 3.2 as a representative nonlinear material;

 

5. An extension of the finite element method to account for large shape changes, using finite strain elasticity as a representative example;

 

6. A discussion of finite element procedures for history dependent solids, using small strain viscoplasticity as a representative example.

 

7. A discussion of the phenomenon of `locking’ that can cause the standard finite element method to fail in certain circumstances.   Several techniques for avoiding locking are presented.

 

8. Finite element methods for structural analysis, including trusses, beams and plates;

 

9. An outline of special procedures in finite element analysis, including use of penalty methods and Lagrange multipliers to enforce constraints; modeling interface fracture using cohesive zones; and contact.

 

In addition, a set of sample finite element codes (implemented in MATLAB) are provided to illustrate the how the various finite element procedures are implemented in practice.  Sample codes may be downloaded from https://github.com/albower/Applied_Mechanics_of_Solids.

 

 

 

 

 

8.1 Generalized FEM for static linear elasticity

 

This section gives a more general derivation and implementation of the finite element method for static linear elastic solids than the energy-based derivation given in Chapter 7.

 

 

 

8.1.1 Review of the principle of virtual work

 

Governing equations: We begin by summarizing the usual governing equations of linear elasticity, which must be solved by the FEA code.   The problem to be solved is illustrated in the figure. We are given:

 

1. The shape of the solid in its unloaded condition $R_0\approx R$

 

2. The initial stress field in the solid (we will take this to be zero in setting up our FEM code)

 

3. The elastic constants for the solid $C_{ijkl}$

 

4. The thermal expansion coefficients for the solid, and temperature distribution (we will take this to be zero for our FEM code, for simplicity)

 

5. A body force distribution $b$ acting on the solid (Note that in this section we will use b to denote force per unit volume rather than force per unit mass, to avoid having to write out the mass density all the time)

 

6. Boundary conditions, specifying displacements $u^{*}(x)$ on a portion ${\partial }_{1}R$ or tractions on a portion ${\partial }_{2}R$ of the boundary of R

 

We then wish to calculate displacements, strains and stresses $u_i,{\epsilon }_{ij},{\sigma }_{ij}$ satisfying the governing equations of static linear elasticity

 

1. The strain-displacement equation ${\epsilon }_{ij}={\displaystyle \frac{1}{2}(\partial u_i/\partial x_j+\partial u_j/\partial x_i)}$

 

2. The elastic stress-strain law ${\sigma }_{ij}=C_{ijkl}{\epsilon }_{kl}$

 

3. The equation of static equilibrium for stresses $\partial {\sigma }_{ij}/\partial x_i+b_j=0$

 

4. The boundary conditions on displacement and stress $u_i=u_i^{*}\text{on}{\partial }_{1}R{\sigma }_{ij}{n}_i=t_j^{*}\text{on}{\partial }_{2}R$

 

 

The principle of virtual work: As we discussed in section 2.4, the principle of virtual work can be used to replace the stress equilibrium equations.   To express the principle, we define a kinematically admissible virtual displacement field $\delta v(x)$, satisfying  $\delta v=0$ on ${\partial }_{1}R$.  You can visualize this field as a small change in the displacement of the solid, if you like, but it is really just an arbitrary differentiable vector field.  The term `kinematically admissible’ is just a complicated way of saying that $\delta v=0$ on ${\partial }_{1}R$ - that is to say, if you perturb the displacement slightly, the boundary conditions on displacement are still satisfied.

 

In addition, we define an associated virtual strain field

$\delta {\epsilon }_{ij}={\displaystyle \frac{1}{2}\left({\displaystyle \frac{\partial \delta v_i}{\partial x_j}+{\displaystyle \frac{\partial \delta v_j}{\partial x_i}}}\right)}$ 

The principle of virtual work states that if the stress field ${\sigma }_{ij}$ satisfies

${\displaystyle \underset{R}{\int }{\sigma }_{ij}\delta {\epsilon }_{ij}d{V}_0-{\displaystyle \underset{R}{\int }{b}_i\delta v_{i}d{V}_0-{\displaystyle \underset{{\partial }_{2}R}{\int }{t}_i\delta v_i}}}dA=0$

for all possible virtual displacement fields and corresponding virtual strains, it will automatically satisfy the equation of stress equilibrium $\partial {\sigma }_{ij}/\partial x_i+b_j=0$, and also the traction boundary condition ${\sigma }_{ij}{n}_i=t_j^{*}\text{on}{\partial }_{2}R$

 

 

 

8.1.2 Integral (weak) form of the governing equations of linear elasticity

 

The principle of virtual work can be used to write the governing equation for the displacement field in a linear elastic solid in an integral form (called the `weak form’).  Instead of solving the governing equations listed in the preceding section, the displacements, strains and stresses $u_i,{\epsilon }_{ij},{\sigma }_{ij}$ are calculated as follows. 

 

1. Find a displacement field $u_i(x_j)$ satisfying

${\displaystyle \underset{R}{\int }{C}_{ijkl}{\displaystyle \frac{\partial u_k}{\partial x_l}{\displaystyle \frac{\partial \delta v_i}{\partial x_j}dV-{\displaystyle \underset{R}{\int }{b}_i\delta v_{i}dV-{\displaystyle \underset{{\partial }_{2}R}{\int }{t}_i^{*}\delta v_{i}dA=0,}}}}}{u}_i=u_i^{*}\text{on}{\partial }_{1}R$

for all virtual velocity fields $\delta v_i$ satisfying $\delta v_i=0$ on ${\partial }_{1}R$.

 

2. Compute the strains from the definition ${\epsilon }_{ij}={\displaystyle \frac{1}{2}(\partial u_i/\partial x_j+\partial u_j/\partial x_i)}$

 

3. Compute the stresses from the stress-strain law ${\sigma }_{ij}=C_{ijkl}{\epsilon }_{kl}$

 

 

The stress will automatically satisfy the equilibrium equation and boundary conditions, so all the field equations and boundary conditions will be satisfied.

 

The significance of this result is that it replaces the derivatives in the partial differential equations of equilibrium with an equivalent integral, which is easier to handle numerically.  It is essentially equivalent to replacing the equilibrium equation with the principle of minimum potential energy, but the procedure based on the principle of virtual work is very easily extended to dynamic problems, other stress-strain laws, and even to problems involving large shape changes.

 

 

Derivation: start with the virtual work equation

${\displaystyle \underset{R}{\int }{\sigma }_{ij}\delta {\epsilon }_{ij}d{V}_0-{\displaystyle \underset{R}{\int }{b}_i\delta v_{i}d{V}_0-{\displaystyle \underset{{\partial }_{2}R}{\int }{t}_i\delta v_i}}}dA=0$

Recall that ${\sigma }_{ij}={\sigma }_{ji}$, and that $\delta {\epsilon }_{ij}=(\partial \delta v_i/\partial x_j+\partial \delta v_j/\partial x_i)/2$, so that

${\sigma }_{ij}\delta {\epsilon }_{ij}={\displaystyle \frac{1}{2}\left({\sigma }_{ij}{\displaystyle \frac{\partial \delta v_i}{\partial x_j}+{\sigma }_{ji}{\displaystyle \frac{\partial \delta v_j}{\partial x_i}}}\right)={\sigma }_{ij}{\displaystyle \frac{\partial \delta v_i}{\partial x_j}}}$

Finally, recall that ${\sigma }_{ij}=C_{ijkl}{\epsilon }_{kl}=C_{ijkl}(\partial u_k/\partial x_l+\partial u_l/\partial x_k)/2$ and that the elastic compliances must satisfy $C_{ijkl}=C_{ijlk}$ so that ${\sigma }_{ij}=C_{ijkl}(\partial u_k/\partial x_l)$.  Finally, this shows that

${\sigma }_{ij}\delta {\epsilon }_{ij}=C_{ijkl}{\displaystyle \frac{\partial u_i}{\partial x_j}{\displaystyle \frac{\partial \delta v_i}{\partial x_j}}}$

Substituting into the virtual work equation gives the result we need.

 

 

 

8.1.3 Interpolating the displacement field and the virtual velocity field

 

To solve the integral form of the elasticity equations given in 8.1.2, we discretize the displacement field.  A representative finite element mesh is sketched in the figure. We choose to calculate the displacement field at a set of n discrete points in the solid (called `nodes’ in finite element terminology).  We will denote the coordinates of these special points by $x_i^a$, where the superscript a ranges from 1 to N.  The unknown displacement vector at each nodal point will be denoted by $u_i^a$.

 

The displacement field at an arbitrary point within the solid will be specified by interpolating between nodal values in some convenient way.  An efficient and robust implementation of the finite element method requires a careful choice of interpolation scheme, but for now we will denote the interpolation in a general way as

$u_i(x)={\displaystyle \sum _{a=1}^n{N}^a(x)u_i^a}$

Here, x denotes the coordinates of an arbitrary point in the solid.  The interpolation functions $N^a(x)$ are functions of position only, which must have the property that

$u_i^b={\displaystyle \sum _{a=1}^n{N}^a(x^b)u_i^a}$

for all b=1…N.  (This is to make sure that the displacement field has the correct value at each node).  Recently developed meshless finite element methods use very complex interpolation functions, but the more traditional approach is to choose them so that

$N^a(x^b)=\left\{\begin{array}{c}1a=b\\ 0a\ne b\end{array}\right.$

The simple constant strain triangle elements introduced in 7.1 are one example of this type of interpolation scheme.  We will define more complicated interpolation functions shortly.

 

We can obviously interpolate the virtual velocity field in exactly the same way (since the principle of virtual work must be satisfied for all virtual velocities, it must certainly be satisfied for an interpolated velocity field…) so that

$\delta v_i(x)={\displaystyle \sum _{a=1}^n{N}^a(x)\delta v_i^a}$

where $\delta v_i^a$ are arbitrary nodal values of virtual velocity.

 

 

 

8.1.4 Finite element equations

 

Substituting the interpolated fields into the virtual work equation, we find that

${\displaystyle \underset{R}{\int }{C}_{ijkl}{\displaystyle \frac{\partial N^b(x)}{\partial x_l}{u}_k^b{\displaystyle \frac{\partial N^a(x)}{\partial x_j}\delta v_i^{a}dV}}}-{\displaystyle \underset{R}{\int }{b}_i{N}^a(x)\delta v_i^a}dV-{\displaystyle \underset{\partial 2R}{\int }{t}_i^{*}{N}^a(x)\delta v_i^a}dA=0$

where summation on a and b is implied, in addition to the usual summation on i,j,k,l. 

 

Note that the interpolation functions are known functions of position. We can therefore re-write the virtual work equation in matrix form as

$\left(K_{aibk}{u}_k^b-F_i^a\right)\delta v_i^a=0$

where

$K_{aibk}={\displaystyle \underset{R}{\int }{C}_{ijkl}{\displaystyle \frac{\partial N^a(x)}{\partial x_j}{\displaystyle \frac{\partial N^b(x)}{\partial x_l}dV}}}{F}_i^a={\displaystyle \underset{R}{\int }{b}_i{N}^a(x)}dV+{\displaystyle \underset{\partial 2R}{\int }{t}_i^{*}{N}^a(x)}dA$

Here K is known as the `stiffness matrix’ and F is known as the force vector.  K is a function only of the elastic properties of the solid, its geometry, and the interpolation functions and nodal positions. It is therefore a known matrix.  Similarly, F is a function only of the known boundary loading and body force field, and the interpolation scheme and nodal positions.  Observe that the symmetry of the elasticity tensor implies that K also has some symmetry $–$ specifically $K_{aibk}=K_{bkai}$.

 

The virtual work equation must be satisfied for all possible sets of $\delta v_i^a$ with $\delta v_i^a=0$ for nodes a that lie on ${\partial }_{1}R$.   At these nodes, the displacements must satisfy $u_i^a=u_i^{*}$.   Evidently, this requires

$\begin{array}{l}{K}_{aibk}{u}_k^b=F_i^a\forall \{a,i\}:x_k^a\text{not on }{\partial }_{\text{1}}R\\ u_i^a=u_i^{*}(x_i^a)\forall \{a,i\}:x_k^a\text{on }{\partial }_{\text{1}}R\end{array}$

This is a system of n linear equations for the n nodal displacements.

 

 

 

 

8.1.5 Simple 1D Implementation of the finite element method

 

Before describing a fully general 3D implementation of the finite element method, we will illustrate all the key ideas using a simple 1-D example. Consider a long linear elastic bar, as shown in the figure. Assume

 

1. The bar has shear modulus $\mu$ and Poisson’s ratio $\nu$

 

2. The bar has cross section $h\times h$ and length L

 

3. It is constrained on all its sides so that $u_2=u_3=0$

 

4. The bar is subjected to body force $b=b(x_1)e_1$,

 

5. The bar is either loaded or constrained at its ends, so that the boundary conditions are either  $t_1(0)=t^{*}(0),t_1(L)=t^{*}(L)$ or displacement $u_1(0)=u^{*}(0)u_1(L)=u^{*}(L)$ at x=0 and x=L. 

 

 

For this 1-D example, then, the finite element equations reduce to

$K_{ab}{u}_1^b=F_{}^a$

where

$\begin{array}{l}{K}_{ab}=h^2{\displaystyle \underset{0}{\overset{L}{\int }}{\displaystyle \frac{2\mu (1-\nu )}{1-2\nu }{\displaystyle \frac{\partial N^a(x_1)}{\partial x_1}{\displaystyle \frac{\partial N^b(x_1)}{\partial x_1}d{x}_1}}}}\\ F_{}^a=h^2{\displaystyle \underset{0}{\overset{L}{\int }}{b}_{}{N}^a(x_1)d{x}_1}+h^2{t}_1^{*}(0)N^a(0)+h^2{t}_{}^{*}(L)N^a(L)\end{array}$

We could obviously choose any interpolation scheme, evaluate the necessary integrals and solve the resulting system of equations to compute the solution.  It turns out to be particularly convenient, however, to use a piecewise-Lagrangian interpolation scheme, and to evaluate the integrals numerically using a Gaussian quadrature scheme.

 

To implement the Lagrangian interpolation scheme, we sub-divide the region $0\le x_1\le L$ into a series of elements, as illustrated in the figure. Each element is bounded by two nodal points, and may also contain one or more interior nodes. The displacement field within the element is interpolated between the nodes attached to the element.  So, we would use a linear interpolation between the nodes on a 2-noded element, a quadratic interpolation between the nodes on a 3 noded element, and so on.

 

 

 

 

Generic linear and quadrilateral 1-D elements are illustrated in the table above.  The local nodes on the element are numbered 1 and 2 for the linear element, and 1,2,3 for the quadratic element as shown.  We suppose that the element lies in the region $-1\le {\xi }_1\le 1$.  The displacements within the element are then interpolated as

$u_1({\xi }_1)={\displaystyle \sum _{a=1}^{N_e}{N}^a({\xi }_1)u_1^a}$

where $N_e$ denotes the number of nodes on the element, $u_1^a$ denotes the value of the displacement at each node, and the shape functions are given in the table.                

                             
Of course, the actual nodal coordinates do not lie at $–$1, +1 and 0 for all the elements.  For a general element, we map this special one to the region of interest.  A particularly convenient way to do this is to set

$x_1({\xi }_1)={\displaystyle \sum _{a=1}^{N_e}{N}^a(\xi )x_1^a}$

where $x_1^a$ denotes the coordinate of each node on the element, and $N_e$ is the number of nodes on the element (2 or 3).  Elements that interpolate displacements and position using the same shape functions are called isoparametric elements.

 

Next, we need to devise a way to do the integrals in the expressions for the stiffness matrix and force vector.  We can evidently divide up the integral so as to integrate over each element in turn

$\begin{array}{l}{K}_{ab}={\displaystyle \sum _{l=1}^{N_{lmn}}{h}^2{\displaystyle \underset{x0}{\overset{x1}{\int }}{\displaystyle \frac{2\mu (1-\nu )}{1-2\nu }{\displaystyle \frac{\partial N^a(x_1)}{\partial x_1}{\displaystyle \frac{\partial N^b(x_1)}{\partial x_1}d{x}_1}}}}}\\ F_{}^a={\displaystyle \sum _{l=1}^{N_{lmn}}{h}^2{\displaystyle \underset{x0}{\overset{x1}{\int }}{b}_{}{N}^a(x_1)d{x}_1}}+h^2{t}_i^{*}{N}^a(0)+h^2{t}_i^{*}{N}^a(L)\end{array}$

where $N_{lmn}$ is the total number of elements, and $x0$ and $x1$ denote the coordinates of the ends of the lth element.  We now notice an attractive feature of our interpolation scheme.  The integral over the lth element depends only on the shape functions associated with the nodes on the lth element, since the displacement in this region is completely determined by its values at these nodes.  We can therefore define element stiffness matrix, and element force matrix

$k_{ab}=h^2{\displaystyle \underset{x0}{\overset{x1}{\int }}{\displaystyle \frac{2\mu (1-\nu )}{1-2\nu }{\displaystyle \frac{\partial N^a(x_1)}{\partial x_1}{\displaystyle \frac{\partial N^b(x_1)}{\partial x_1}d{x}_1}}}}{f}_{}^a=h^2{\displaystyle \underset{x0}{\overset{x1}{\int }}{b}_{}{N}^a(x_1)d{x}_1}$

for each element, which depend on the geometry, interpolation functions and material properties of the element.  The first and last elements have additional contributions to the element force vector from the boundary terms $h^2{t}_i^{*}{N}^a(0),h^2{t}_i^{*}{N}^a(L)$.  The global stiffness matrix is computed by summing all the element stiffness matrices

$K_{ab}={\displaystyle \sum _{l=1}^{N_{lmn}}{k}_{ab}}{F}_{}^a={\displaystyle \sum _{l=1}^{N_{lmn}}{f}^a}+h^2{t}_i^{*}{N}^a(0)+h^2{t}_i^{*}{N}^a(L)$

 

Finally we need to devise a way to compute the integrals for each element stiffness matrix.  It is convenient to map the domain of integration to [-1,+1] and integrate with respect to the normalized coordinate $\xi$ - thus

$k_{ab}=h^2{\displaystyle \underset{-1}{\overset{+1}{\int }}{\displaystyle \frac{2\mu (1-\nu )}{1-2\nu }{\displaystyle \frac{\partial N^a(x)}{\partial x}{\displaystyle \frac{\partial N^b(x)}{\partial x}Jd\xi }}}}{f}_{}^a=h^2{\displaystyle \underset{-1}{\overset{+1}{\int }}{b}_{}{N}^a(x_1)Jd\xi }$

where $J=\left|\partial x/d\xi \right|$ is the Jacobian associated with the mapping, which may be computed as

${\displaystyle \frac{\partial x}{\partial \xi }={\displaystyle \frac{\partial }{\partial \xi }{\displaystyle \sum _{a=1}^{N_e}{N}^a({\xi }_{})x_{}^a}={\displaystyle \sum _{a=1}^{N_e}{\displaystyle \frac{\partial N^a}{\partial \xi }{x}_{}^a}}}}$

Note that the mapping also enables us to calculate the shape function derivatives in the element stiffness matrix as

${\displaystyle \frac{\partial N^a}{\partial x}={\displaystyle \frac{\partial N^a}{\partial \xi }{\displaystyle \frac{\partial \xi }{\partial x}={\displaystyle \frac{\partial N^a}{\partial \xi }{\left[{\displaystyle \frac{\partial x}{\partial \xi }}\right]}^{-1}}}}}$

 

Finally, note that integrals may be computed numerically using a quadrature formula, as follows

${\displaystyle \underset{-1}{\overset{+1}{\int }}g(\xi )d\xi }\approx {\displaystyle \sum _{I=1}^M{w}_{I}g({\xi }^{(I)})}$

where  ${\xi }^{(I)}$  I=1…M denotes a set of integration points in the region [-1,+1], and $w_I$ is a set of integration weights, which are chosen so as to make the approximation as accurate as possible. Values are given in the table below for $M=1,2$ and 3. Higher order integration schemes exist but are required only for higher order elements.  For the linear 1-D element described earlier a single integration point is sufficient to evaluate the stiffness exactly. Similarly, for the quaratic element, two integration points will suffice.

 

 

 

 

8.1.6 Summary of the 1D finite element procedure

 

To summarize, then, the finite element solution requires the following steps:

 

1. For each element, compute the element stiffness matrix as follows:

$k_{ab}=h^2{\displaystyle \sum _{I=1}^M{w}_I{\displaystyle \frac{2\mu (1-\nu )}{1-2\nu }{\displaystyle \frac{\partial N^a}{\partial x}{\displaystyle \frac{\partial N^b}{\partial x}J({\xi }_I)}}}}$

where

$J=\left|{\displaystyle \frac{\partial x}{\partial \xi }}\right|=\left|{\displaystyle \sum _{a=1}^{N_e}{\displaystyle \frac{\partial N^a({\xi }_I)}{\partial \xi }{x}_{}^a}}\right|$            ${\displaystyle \frac{\partial N^a}{\partial x}={\displaystyle \frac{\partial N^a({\xi }_I)}{\partial \xi }{\displaystyle \frac{\partial \xi }{\partial x}={\displaystyle \frac{\partial N^a({\xi }_I)}{\partial \xi }{\left[{\displaystyle \frac{\partial x}{\partial \xi }}\right]}^{-1}}}}}$

and the integration points ${\xi }_I,w_{I}I=\mathrm{1...}M$ are listed in the table above, and shape functions $N^{(a)}(\xi ),a=\mathrm{1...}{N}_e$ were listed in Section 8.1.5.

 

2. Assemble the contribution from each element to the global stiffness

$K_{ab}={\displaystyle \sum _{l=1}^{N_{lmn}}{k}_{ab}}$

 

3. Similarly, if there is a non-zero body force, then compute for each element

$f_{}^a=h^2{\displaystyle \sum _{I=1}^M{w}_I{b}_{}{N}^a({\xi }_I)J}({\xi }_I)$

and assemble the global force vector

$F_{}^a={\displaystyle \sum _{l=1}^{N_{lmn}}{f}^a}$

 

4. Add contributions to the force vector from prescribed traction boundary conditions at $x=0$ and $x=L$

$F^1=F^1+h^2{t}_i^{*}(0)F^{(L)}=F^{(L)}+h^2{t}_i^{*}(L)$

where the $(L)$ superscript denotes the node that lies at x=L.

 

5. Modify the stiffness matrix to enforce the constraints

$u^1=u^{*}(0)u^{(L)}=u^{*}(L)$

 

6. Solve the system of linear equations

$K_{ab}{u}_1^b=F_{}^a$

for the unknown displacements $u_1^b$

 

 

 

8.1.7 Example FEM Code and solution

 

A simple example MATLAB code for this 1-D example can be found in the file FEM_1D_Static.m at https://github.com/albower/Applied_Mechanics_of_Solids/

 

It is set up to solve for displacements for the bar sketched in the figure, with the following parameters (in arbitrary units):

 

1. Length L=5, and unit cross-sectional area,

 

2. Shear modulus 50, Poisson’s ratio 0.3,

 

3. Uniform body force magnitude 10,

 

4. Displacement u=0 at x=0

 

5. Traction t=2 at x=L.

 

The code computes the (1D) displacement distribution in the bar. The predicted displacement field is plotted in the figure.

 

Of course, in general we want to calculate more than just displacements $–$ usually we want the stress field too. We can calculate the stress field anywhere within an element by differentiating the displacements to calculate strains, and then substituting into the constitutive relation.  This gives

${\sigma }_{11}={\displaystyle \frac{2\mu (1-\nu )}{(1-2\nu )}{\displaystyle \sum _{a=1}^n{\displaystyle \frac{\partial N^a}{\partial x}}}{u}^a}$

This works well for a uniform body force with quadratic (3 noded elements) as the figure below shows.

 

 

 

However, if we switch to linear elements, the stress results are not so good (displacements are still calculated exactly).  In this case, the stress must be uniform in each element (because strains are constant for a linear displacement field), so the stress plot looks like the graph to the right.

Notice that the stresses are most accurate near the center of each element (at the integration point).  For this reason, FEM codes generally output stress and strain data at integration points. 

 

It is interesting also to examine the stiffness matrix (shown below for 3 linear elements, before addition of the u=0 constraint for the first node)

$\left[\begin{array}{cccc}105& -105& 0& 0\\ -105& 210& -105& 0\\ 0& -105& 210& -105\\ 0& 0& -105& 105\end{array}\right]$

Notice that stiffness is symmetric, as expected, and also banded.  A large FEM matrix is sparse $–$ most of the elements are zero.  This allows the matrix to be stored in compact form $–$ for very large matrices indexed storage (where only the nonzero elements together with their indices are stored) is the best approach; for smaller problems skyline storage or band storage (where only the central, mostly nonzero, band of the matrix is stored) may be preferable.  In this case equation numbers need to be assigned to each degree of freedom so as to minimize the bandwidth of the stiffness matrix.

 

 

 

8.1.8 Extending the 1D finite element method to 2 and 3 dimensions

 

It is straightforward to extend the 1-D case to more general problems.  All the basic ideas remain the same.   Specifically

 

1. In both 2D and 3D we divide up our solid of interest into a number of elements, shown schematically for a 2D region in the figure.

 

2. We define interpolation functions $N^a({\xi }_i)$ for each element in terms of a local, dimensionless, coordinate system within the element.  The coordinates satisfy $-1\le {\xi }_i\le +1$.  The displacement field and the position of a point inside an element are computed in terms of the interpolation functions as

$u_i={\displaystyle \sum _{a=1}^{N_e}{N}^a({\xi }_j)u_i^a}{x}_i={\displaystyle \sum _{a=1}^{N_e}{N}^a({\xi }_j)x_i^a}$

where $N^a({\xi }_j)$ denote the shape functions, $u_i^a,x_i^a$ denote the displacement values and coordinates of the nodes on the element, and $N_e$ is the number of nodes on the element.

 

3. We introduce an element stiffness matrix for each element by defining

$k_{aibk}^{(l)}={\displaystyle \underset{V_e^{(l)}}{\int }{C}_{ijkl}{\displaystyle \frac{\partial N^a(x)}{\partial x_j}{\displaystyle \frac{\partial N^b(x)}{\partial x_l}dV}}}{f}_i^a{}^{(l)}={\displaystyle \underset{V_e^{(l)}}{\int }{b}_i{N}^a(x)}dV+{\displaystyle \underset{{\partial }_2{V}_e^{(l)}}{\int }{t}_i^{*}{N}^a(x)}dA$

where $k_{aibk}^{(l)}$ denotes the element stiffness matrix for the (lth) element, and $V_e^{(l)}$ denotes the volume (in 3D) or area (in 2D) of the (lth) element, while ${\partial }_2{V}_e^{(l)}$ denotes the surface of the (lth) element

 

4. The volume integrals over each element are calculated by expressing the volume or surface integral in terms of the dimensionless coordinates $-1\le {\xi }_i\le +1$, and then evaluating the integrals numerically, using a quadrature formula of the form

${\displaystyle \underset{V_e}{\int }f({\xi }_i)d{V}_{\xi }}={\displaystyle \sum _{I=1}^{N_I}{w}_{I}f({\xi }_i^I)}$

Here, $w_I$ are a set of $I=\mathrm{1...}{N}_I$ integration weights (just numbers), and ${\xi }_i^{(I)}$ are a set of coordinates that are selected to make the integration scheme as accurate as possible (also just numbers).

 

5. The global stiffness matrix

$K_{aibk}={\displaystyle \underset{R}{\int }{C}_{ijkl}{\displaystyle \frac{\partial N^a(x)}{\partial x_j}{\displaystyle \frac{\partial N^b(x)}{\partial x_l}dV}}}{F}_i^a={\displaystyle \underset{R}{\int }{b}_i{N}^a(x)}dV+{\displaystyle \underset{\partial 2R}{\int }{t}_i^{*}{N}^a(x)}dA$

is then computed by summing the contribution from each element as

$K_{aibk}={\displaystyle \sum _{l=1}^{N_{lmn}}{k}_{aibk}^{(l)}}\text{}\text{}{F}_i^a={\displaystyle \sum _{l=1}^{N_{lmn}}{f}_i{{}^a}^{(l)}}$

 

6. The stiffness matrix is modified to enforce any prescribed displacements

 

7. The system of equations

$\begin{array}{l}{K}_{aibk}{u}_k^b=F_i^a\forall \{a,i\}:x_k^a\text{not on }{\partial }_{\text{1}}R\\ u_i^a=u_i^{*}(x_i^a)\forall \{a,i\}:x_k^a\text{on }{\partial }_{\text{1}}R\end{array}$

is solved for the unknown nodal displacements.

 

8. The stresses and strains within each element are then deduced.

 

To implement this procedure, we must (a) Define the element interpolation functions; (b) Express the integrals for the element stiffness matrices and force vectors in terms of normalized coordinates; (c) Formulate a numerical integration scheme to evaluate the element stiffness matrices and force vectors.

These details are addressed in the sections to follow.

 

 

 

8.1.9 Interpolation functions for 2D elements

 

The interpolation functions for 2D elements are listed in Table 8.3. They are defined for the region $0\le {\xi }_1\le 10\le {\xi }_2\le 1$ for triangular elements, and $-1\le {\xi }_1\le 1-1\le {\xi }_2\le 1$ for quadrilateral elements. The numbers shown inside the element show the convention used to number the element faces.

 

 

 

 

 

 

 8.1.10 Interpolation Functions for 3D elements

 

Interpolation functions for 3D elements are listed in the table below. The interpolation functions are defined for the region $-1\le {\xi }_i\le 1$ unless the figures show otherwise.                                      

    

 

 

The element faces are numbered as listed in the table below.

 

 

 

 

 

 

8.1.11 Volume integrals for stiffness and force in terms of normalized coordinates

 

In this section we outline the procedure that is used to re-write the integrals for the element stiffness and force in terms of the normalized coordinates ${\xi }_i$.  The integrals are

$k_{aibk}^{(l)}={\displaystyle \underset{V_e^{(l)}}{\int }{C}_{ijkl}{\displaystyle \frac{\partial N^a(x)}{\partial x_j}{\displaystyle \frac{\partial N^b(x)}{\partial x_l}dV}}}{f}_i^a{}^{(l)}={\displaystyle \underset{V_e^{(l)}}{\int }{b}_i{N}^a(x)}dV+{\displaystyle \underset{{\partial }_2{V}_e^{(l)}}{\int }{t}_i^{*}{N}^a(x)}dA$

To evaluate them, we need to

 

1. Find a way to calculate the derivatives of the shape functions in terms of ${\xi }_i$

 

2. Map the volume (or area) integral to the region $-1\le {\xi }_i\le +1$

 

 

Calculating the shape function derivatives. The shape function derivatives can be evaluated by writing

${\displaystyle \frac{\partial N^a}{\partial x_j}={\displaystyle \frac{\partial N^a}{\partial {\xi }_i}{\displaystyle \frac{\partial {\xi }_i}{\partial x_j}}}}$

where the derivatives $\partial N^a/\partial {\xi }_i$ are easy to compute (just differentiate the expressions given earlier…).  To compute $\partial {\xi }_i/\partial x_j$ recall that the coordinates of a point at position ${\xi }_j$ within an element can be determined as

$x_i={\displaystyle \sum _{a=1}^{N_e}{N}^a({\xi }_j)x_i^a}$

where $N_e$ denotes the number of nodes on the element. Therefore

${\displaystyle \frac{\partial x_i}{\partial {\xi }_j}={\displaystyle \frac{\partial }{\partial {\xi }_j}\left({\displaystyle \sum _{a=1}^{N_e}{N}^a(\xi )x_i^a}\right)={\displaystyle \sum _{a=1}^{N_e}{\displaystyle \frac{\partial N^a(\xi )}{\partial {\xi }_j}{x}_i^a}}}}$

Note that $\partial x_i/\partial {\xi }_j$ is a 2x2 matrix (in 2D) or a 3x3 matrix (in 3D).  Finally, $\partial {\xi }_i/\partial x_j$ follows as the inverse of this matrix

${\displaystyle \frac{\partial {\xi }_j}{\partial x_i}={\left({\displaystyle \frac{\partial x}{\partial \xi }}\right)}_{ji}^{-1}}$

 

 

Mapping the volume integral: To map the region of integration we define

$J=\det \left({\displaystyle \frac{\partial x_i}{\partial {\xi }_j}}\right)$

where the matrix $\partial x_i/\partial {\xi }_j$ was defined earlier.  Then the integral with respect to x is mapped into an integral with respect to $\xi$ by setting

$k_{aibk}={\displaystyle \underset{\Omega }{\int }{C}_{ijkl}{\displaystyle \frac{\partial N^a(x)}{\partial x_j}{\displaystyle \frac{\partial N^b(x)}{\partial x_l}Jd{V}_{\xi }}}}{f}_i^a={\displaystyle \underset{\Omega }{\int }{b}_i{N}^a(x)}Jd{V}_{\xi }+{\displaystyle \underset{\partial \Omega }{\int }{t}_i^{*}{N}^a(x)}\stackrel{⌢}{J}d{A}_{\widehat{\xi }}$

We note in passing that the boundary integral in the element force vector can be regarded as a 1-D line integral for 2D elements and a 2D surface integral for 3D elements.  So the procedures we developed in 8.1.5 (1D elements) can be used to evaluate the surface integral for a 2D element.  Similarly, the procedures we develop to integrate stiffness matrices for 2D elements can be used to evaluate the surface integral for a 3D element.

 

 

 

8.1.12 Numerical integration schemes for 2D and 3D elements

 

Finally, to evaluate the integrals, we once again adopt a quadrature scheme, so that

${\displaystyle \underset{\Omega }{\int }f({\xi }_i)d{V}_{\xi }}={\displaystyle \sum _{I=1}^{N_I}{w}_{I}f({\xi }_i^I)}$

The integration points ${\xi }_j^I$ and weights $w_I$ depend on the element geometry, and are listed in Tables 8.5-8.7 for several element types.

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Choosing the number of integration points: There are two considerations.  If too many integration points are used, time is wasted without gaining any accuracy.  If too few integration points are used, the stiffness matrix may be singular, or else the rate of convergence to the exact solution with mesh refinement will be reduced.  The schemes listed in the table below will avoid both.

 

There are situations where it is preferable to use fewer integration points and purposely make the stiffness singular.   These are discussed in more detail in Section 8.5.

 

 

 

 

 

 

8.1.13 Computer Implementation

 

With these definitions, then, we write the element stiffness matrix as

$\begin{array}{l}{k}_{aibk}={\displaystyle \sum _{I=1}^{N_I}{w}_I{C}_{ijkl}{\displaystyle \frac{\partial N^a({\xi }_i^I)}{\partial {\xi }_p}{\displaystyle \frac{\partial {\xi }_p}{\partial x_j}{\displaystyle \frac{\partial N^b({\xi }_i^I)}{\partial x_q}{\displaystyle \frac{\partial {\xi }_q}{\partial x_l}J({\xi }_i^I)}}}}}\\ f_i^a={\displaystyle \sum _{I=1}^N{w}_I{b}_i{N}^a({\xi }_i^I)J}+{\displaystyle \sum _{I=1}^{\widehat{N}}{w}_I{t}_i^{*}{N}^a({\xi }_i^I)\stackrel{⌢}{J}}\end{array}$

where

${\displaystyle \frac{\partial x_i}{\partial {\xi }_j}={\displaystyle \sum _{a=1}^{N_e}{\displaystyle \frac{\partial N^a(\xi )}{\partial {\xi }_j}{x}_i^a}}}$            $J=\det \left({\displaystyle \frac{\partial x_i}{\partial {\xi }_j}}\right)$            ${\displaystyle \frac{\partial {\xi }_j}{\partial x_i}={\left({\displaystyle \frac{\partial x}{\partial \xi }}\right)}_{ji}^{-1}}$

 

The stiffness matrix is usually computed by re-writing these expressions in matrix form.  To this end, the coordinates of nodes on the element are and shape functions derivatives (at the current integration point) are stored as matrices

$x=\left\{\begin{array}{c}\left[\begin{array}{ccc}{x}_1^1& x_1^2& x_1^3\\ x_2^1& x_2^2& \cdots \end{array}\right](2\text{D)}\\ \left[\begin{array}{ccc}{x}_1^1& x_1^2& x_1^3\\ x_2^1& x_2^2& \cdots \\ x_3^1& & \end{array}\right](3\text{D)}\end{array}\right.$

${\displaystyle \frac{dN}{d\xi }=\left\{\begin{array}{c}\left[\begin{array}{cc}\partial N^1/\partial {\xi }_1& \partial N^1/\partial {\xi }_2\\ \partial N^2/\partial {\xi }_1& \partial N^2/\partial {\xi }_2\\ \partial N^3/\partial {\xi }_1& \partial N^3/\partial {\xi }_2\\ ⋮& \end{array}\right](2\text{D)}\\ \left[\begin{array}{ccc}\partial N^1/\partial {\xi }_1& \partial N^1/\partial {\xi }_2& \partial N^1/\partial {\xi }_3\\ \partial N^2/\partial {\xi }_1& \partial N^2/\partial {\xi }_2& \\ \partial N^3/\partial {\xi }_1& & \ddots \end{array}\right](3\text{D)}\end{array}\right.}$

 

The element stiffness matrix is then computed as follows:

 

Loop over the integration points.   For the current integration point:

 

1.  Calculate the Jacobian matrix, its determinant and the spatial shape function derivatives

${\displaystyle \frac{dx}{d\xi }=x{\displaystyle \frac{dN}{d\xi }J=\det \left({\displaystyle \frac{dx}{d\xi }}\right){\displaystyle \frac{dN}{dx}={\displaystyle \frac{dN}{d\xi }{\left({\displaystyle \frac{dx}{d\xi }}\right)}^{-1}}}}}$

 

2.  Assemble a matrix B that maps the element displacement vector onto a strain vector as defined below

$\epsilon =BU$

$U=\left\{\begin{array}{c}{\left[\begin{array}{ccccc}{u}_1^1& u_2^1& u_1^2& u_2^2& \cdots \end{array}\right]}^T\text{(2D)}\\ {\left[\begin{array}{ccccc}{u}_1^1& u_2^1& u_3^2& u_1^2& \cdots \end{array}\right]}^T\text{(3D)}\end{array}\right.$

$\epsilon =\left\{\begin{array}{c}{\left[\begin{array}{ccc}{\epsilon }_{11}& {\epsilon }_{22}& 2{\epsilon }_{12}\end{array}\right]}^T(2\text{D)}\\ {\left[\begin{array}{cccccc}{\epsilon }_{11}& {\epsilon }_{22}& {\epsilon }_{33}& 2{\epsilon }_{12}& 2{\epsilon }_{13}& 2{\epsilon }_{23}\end{array}\right]}^T(3\text{D)}\end{array}\right.$

$B=\left\{\begin{array}{c}\left[\begin{array}{ccccccc}\partial N^1/\partial x_1& 0& \partial N^2/\partial x_1& 0& \partial N^3/\partial x_1& & \\ 0& \partial N^1/\partial x_2& 0& \partial N^2/\partial x_2& 0& \dots & \\ \partial N^1/\partial x_2& \partial N^1/\partial x_1& \partial N^2/\partial x_2& \partial N^2/\partial x_1& & & \end{array}\right](\text{2D})\\ \left[\begin{array}{ccccccc}\partial N^1/\partial x_1& 0& 0& \partial N^2/\partial x_1& 0& 0& \\ 0& \partial N^1/\partial x_2& 0& 0& \partial N^2/\partial x_2& 0& \\ 0& 0& \partial N^1/\partial x_3& 0& 0& \partial N^2/\partial x_3& \cdots \\ \partial N^1/\partial x_2& \partial N^1/\partial x_1& 0& \partial N^2/\partial x_2& \partial N^2/\partial x_1& & \\ \partial N^1/\partial x_3& 0& \partial N^1/\partial x_1& \partial N^2/\partial x_3& 0& \partial N^2/\partial x_1& \\ 0& \partial N^1/\partial x_3& \partial N^1/\partial x_2& 0& \partial N^2/\partial x_3& \partial N^2/\partial x_2& \end{array}\right](\text{3D})\end{array}\right.$

 

3.  Assemble a matrix D of elastic constants.  For an isotropic material in 3D:

$D={\displaystyle \frac{E}{(1+\nu )(1-2\nu )}\left[\begin{array}{cccccc}1-\nu & \nu & \nu & 0& 0& 0\\ \nu & 1-\nu & \nu & 0& 0& 0\\ \nu & \nu & 1-\nu & 0& 0& 0\\ 0& 0& 0& {\displaystyle \frac{\left(1-2\nu \right)}{2}}& 0& 0\\ 0& 0& 0& 0& {\displaystyle \frac{\left(1-2\nu \right)}{2}}& 0\\ 0& 0& 0& 0& 0& {\displaystyle \frac{\left(1-2\nu \right)}{2}}\end{array}\right]}$

For 2D plane strain

$D={\displaystyle \frac{E}{(1+\nu )(1-2\nu )}\left[\begin{array}{ccc}1-\nu & \nu & 0\\ \nu & 1-\nu & 0\\ 0& 0& {\displaystyle \frac{1-2\nu }{2}}\end{array}\right]}$

while for plane stress

$D={\displaystyle \frac{E}{(1-{\nu }^2)}\left[\begin{array}{ccc}1& \nu & 0\\ \nu & 1& 0\\ 0& 0& (1-\nu )/2\end{array}\right]}$

Matrices for anisotropic materials can be found in Section 3.2.

 

4.  Add the contribution from the current integration point to the sum

$k_{el}={\displaystyle \sum _{I=1}^{N_I}{B}^{T}DBJ{w}_I}$

 

As we shall see in later sections, nearly every finite element stiffness matrix can be rearranged into the form shown in Step 4.    Usually, only the matrices D and B need to be modified to implement new element types.

 

 

 

8.1.14 Sample 2D/3D linear elastostatic FEM code

 

You can find a MATLAB implementations of a simple 2D/3D static linear elasticity code on Github at

https://github.com/albower/Applied_Mechanics_of_Solids.

 

 The code is in the file FEM_2Dor3D_linelast_standard.m.

 

 

The code reads an input file.  Several examples are provided :

 

1. Linear_elastic_triangles.txt: Simple 2D plane strain problem with two triangular elements

 

2. Linear_elastic_quad4.txt: Simple 2D plane strain problem with eight 4 noded quadrilateral elements

 

3. Linear_elastic_quad8.txt: Simple 2D plane strain problem with two 8 noded quadrilateral elements

 

4. Linear_elastic_brick4.txt: Simple 3D problem with 8 noded brick elements

 

 

As an example, we show how to run the program with the first input file. The file sets up the problem illustrated in the figure.  The elements are linear elastic plane strain with $E=10,\nu =0.3$. 

 

The program input file is listed in the figure below. Here is a brief explanation of the data in the file

 

1. The first part of the input file specifies material properties. 

 

2. The second part specifies properties and coordinates of the nodes.  For a 2D problem each node has 2 coordinates and 2 DOF; for a 3D problem each node has 3 coordinates and 3DOF.  Then enter nodal coordinates for each node.

 

3. The third part lists the element properties.  Here, you must specify the number of elements, and the maximum number of nodes on any one element (you can mix element types if you like).  Then you must specify the nodes connected to each element (known as element connectivity).   For each element, you must specify the number of nodes attached to the element; an identifier that specifies the element type (Entering 1 will run a plane strain simulation; 0 will run plane stress), then enter the nodes on each element following the convention shown earlier.

 

4. The fourth part of the file specifies boundary constraints.  For any constrained displacements, enter the node number, the displacement component to be prescribed, and its value.

 

5. The last part of the file specifies distributed loading acting on the element faces.  The loading is assumed to be uniform.  For each loaded boundary, you should specify the element number, the face of the element (the face numbering convention was described in section 7.2.9 and 7.2.10 $–$ note that you must be consistent in numbering nodes and faces on each element), and the components of traction acting on the element face, as a vector with 2 or 3 components.

 

 

Note that the program performs absolutely no error checking on the input file.  If you put in a typo, you will get some bizzarre error message from Matlab $–$ often during element stiffness assembly.

 

For the input file shown, the program plots the deformed mesh, and produces the output file shown below

 

 

 

The code prints the displacements at each node in the mesh, and also the strains and stresses at each integration point (where these quantities are most accurate)  for each element.

 

To run the code, you must complete the following steps

 

1. Open the Matlab executable file;

 

2. Edit the code to insert the full path for the input and output files near the top of the code

 

3. Run the code in the usual way. You will see the code plot the undeformed and deformed finite element mesh at the end.  The stresses and strains in the elements are printed to the output file.