Open Educational Resources

Finite Element Analysis: One and Two Dimensional Isoparametric Elements and Gauss Integration

In the previous section, the basic elements that are used for discretizing a plane domain were introduced. These elements can be used to mesh material bodies with regular geometries. For irregular geometries, quadrilateral and triangle elements will not necessarily follow the shapes of the basic elements described (a quadrilateral element will not necessarily be a square or a rectangle, and the triangle element will not necessarily have a right angle, as described before). Isoparametric elements enable meshing an irregular domain with triangular or quadratic elements that do not maintain orthogonality between the sides of the element. The purpose of the isoparametric formulation is to create shape functions that would ensure the compatibility of the displacement between neighboring elements while maintaining the requirements for shape functions mentioned in the previous section. When the elements’ sides do not necessarily have right angles, the isoparametric formulation is used to map those irregular elements into elements with regular shapes. The mapping functions used are the same shape functions that are used for the displacements, and hence, the name isoparametric (same parameters). After the mapping, integration is performed in the coordinate system of the element with the regular shape. The integration is performed numerically rather than using exact integrals. Numerical integration is preferred over exact integration for two reasons. The first is to reduce the computational resources required, and the second because the irregular geometries render radical (fraction) quantities that could be problematic for exact integrals. Many numerical integration techniques are available but the technique that is widely used for finite element analysis is the Gauss numerical integration technique.

Isoparametric Elements

One Dimensional 3-Node Quadratic Isoparametric Element

the one dimensional quadratic element described previously was defined with the middle node equidistant from the end nodes. However, the isoparametric formulation allows defining a general one dimensional quadratic element with nodal coordinates x_1, x_2, and x_3 such that the only requirement is that x_1<x_2<x_3. Another coordinate system \xi can be defined for the element such that the coordinates of the nodes 1, 2, and 3 are \xi_1=-1, \xi_2=0, and \xi_3=1, respectively. A mapping between the general coordinate system x and the coordinate system \xi can be defined (Figure 1)using the same shape functions that define the displacement in the coordinate system \xi as follows:

    \[ x=N_1x_1+N_2x_2+N_3x_3 \]

where the shape functions have the form:

    \[\begin{split} N_1&=\frac{-\xi(1-\xi)}{2}\\ N_2&=(1-\xi)(1+\xi)\\ N_3&=\frac{\xi(1+\xi)}{2} \end{split} \]

The gradient of this mapping has the following form:

    \[ \frac{\mathrm{d}x}{\mathrm{d}\xi}=\sum_{i=1}^3\frac{\mathrm{d}N_i}{\mathrm{d}\xi}x_i \]

Since x_1<x_2<x_3, and since \xi\in[-1,1], then

    \[ \frac{\mathrm{d}x}{\mathrm{d}\xi}>0 \]

Therefore, the inverse of the gradient is well defined as:

    \[ \frac{\mathrm{d}\xi}{\mathrm{d}x}=\frac{1}{\frac{\mathrm{d}x}{\mathrm{d}\xi}}=\frac{1}{\sum_{i=1}^3\frac{\mathrm{d}N_i}{\mathrm{d}\xi}x_i} \]

If u_1, u_2, and u_3 are the horizontal displacements of the nodes 1, 2, and 3, respectively, then, the horizontal displacement function is defined in the \xi coordinate system as follows:

    \[ u=N_1u_1+N_2u_2+N_3u_3 \]

The normal strain \varepsilon_{11} is defined in the x coordinate system as follows:

    \[ \varepsilon_{11}=\frac{\mathrm{d}u}{\mathrm{d}x}=\frac{\mathrm{d}u}{\mathrm{d}\xi}\frac{\mathrm{d}\xi}{\mathrm{d}x}=\frac{\mathrm{d}\xi}{\mathrm{d}x}\left<\begin{array}{ccc} \frac{\mathrm{d}N_1}{\mathrm{d}\xi} & \frac{\mathrm{d}N_2}{\mathrm{d}\xi} & \frac{\mathrm{d}N_3}{\mathrm{d}\xi} \end{array}\right>\left(\begin{array}{c}u_1\\u_2\\u_3\end{array}\right) \]


Figure 1. One Dimensional Quadratic Isoparametric Mapping

Figure 1. One Dimensional Quadratic Isoparametric Mapping

The virtual displacement u^* and the associated virtual strain \varepsilon_{11}^* can be chosen as follows:

    \[\begin{split} u^*&=N_1u_1^*+N_2u_2^*+N_3u_3^*\\ \varepsilon_{11}^*&=\frac{\mathrm{d}\xi}{\mathrm{d}x}\left<\begin{array}{ccc} u_1^*&u_2^*&u_3^*\end{array}\right>\left(\begin{array}{c} \frac{\mathrm{d}N_1}{\mathrm{d}\xi} \\ \frac{\mathrm{d}N_2}{\mathrm{d}\xi}\\ \frac{\mathrm{d}N_3}{\mathrm{d}\xi} \end{array}\right) \end{split} \]

Assuming p as a distributed load per unit length of the one dimensional element, then, the contribution to the virtual work equation of this element are as follows:

    \[\begin{split} IVW_e&=A\int_{x_1}^{x_3}\! \sigma_{11}\varepsilon_{11}^*\,\mathrm{d}x\\ EVW_e&=\int_{x_1}^{x_3}\! pu^*\,\mathrm{d}x \end{split} \]

To utilize the numerical Gauss integration scheme (as described later) Integration is performed in the \xi coordinate system. In general, for a function g(x) we have:

    \[ \int_{x_1}^{x_3}\! g(x)\,\mathrm{d}x=\int_{-1}^{1}\! g(x(\xi))\,\frac{\mathrm{d}x}{\mathrm{d}\xi}\mathrm{d}\xi \]

Using the constitutive law \sigma_{11}=E\varepsilon_{11}, the virtual work contributions of the element have the forms:

    \[\begin{split} IVW_e&=\left<\begin{array}{ccc} u_1^* & u_2^* & u_3^*\end{array}\right> EA\int_{-1}^{1}\! \left(\frac{\mathrm{d}\xi}{\mathrm{d}x}\right)^2\frac{\mathrm{d}x}{\mathrm{d}\xi}\left(\begin{array}{c} \frac{\mathrm{d}N_1}{\mathrm{d}\xi}\\ \frac{\mathrm{d}N_2}{\mathrm{d}\xi}\\\frac{\mathrm{d}N_3}{\mathrm{d}\xi} \end{array}\right)\left<\begin{array}{ccc} \frac{\mathrm{d}N_1}{\mathrm{d}\xi}&\frac{\mathrm{d}N_2}{\mathrm{d}\xi}&\frac{\mathrm{d}N_3}{\mathrm{d}\xi} \end{array}\right>\,\mathrm{d}\xi \left(\begin{array}{c}u_1\\u_2\\u_3\end{array}\right) \\ EVW_e&=\left<\begin{array}{ccc} u_1^* & u_2^* & u_3^*\end{array}\right> \int_{-1}^{1}\! p\left(\begin{array}{c}N_1\\N_2\\N_3\end{array}\right) \frac{\mathrm{d}x}{\mathrm{d}\xi} \,\mathrm{d}\xi \end{split} \]

Therefore, the stiffness matrix and the nodal forces vector for this element have the form:

    \[\begin{split} K^e&=EA\int_{-1}^{1}\! \left(\frac{\mathrm{d}\xi}{\mathrm{d}x}\right)^2\frac{\mathrm{d}x}{\mathrm{d}\xi}\left(\begin{array}{c} \frac{\mathrm{d}N_1}{\mathrm{d}\xi}\\ \frac{\mathrm{d}N_2}{\mathrm{d}\xi}\\\frac{\mathrm{d}N_3}{\mathrm{d}\xi} \end{array}\right)\left<\begin{array}{ccc} \frac{\mathrm{d}N_1}{\mathrm{d}\xi}&\frac{\mathrm{d}N_2}{\mathrm{d}\xi}&\frac{\mathrm{d}N_3}{\mathrm{d}\xi} \end{array}\right>\,\mathrm{d}\xi\\ f^e&=\int_{-1}^{1}\! p\left(\begin{array}{c}N_1\\N_2\\N_3\end{array}\right) \frac{\mathrm{d}x}{\mathrm{d}\xi} \,\mathrm{d}\xi \end{split} \]

Example

Let a nonlinear isoparametric one dimensional element be placed such that the nodal coordinates x_1, x_2, and x_3 are equal to 0, 4.5, and 10 units, respectively. Assume that a constant load per unit length p is applied in the direction of x. Find the stiffness matrix and the nodal forces for that element.

Solution

The mapping function between the two coordinate systems has the form:

    \[ x=N_1x_1+N_2x_2+N_3x_3=4.5+(5+0.5\xi)\xi \]

Applying the above equations, the stiffness matrix and the nodal forces for the element have the form:

    \[\begin{split} K^e&=EA\left(\begin{matrix}0.2653 & -0.3006& 0.0353\\-0.3006 & 0.5465 & -0.2459\\0.0353 & -0.2459 & 0.21067\end{matrix}\right)\\ f^e&=p\left(\begin{array}{c}1.333\\6.667\\2.000\end{array}\right) \end{split} \]

The following Mathematica code was utilized.
View Mathematica Code

Clear[x1, x2, x3, x, xi, p]
N1 = -xi (1 - xi)/2;
N2 = (1 - xi) (1 + xi);
N3 = xi (1 + xi)/2;
x = N1*x1 + N2*x2 + N3*x3;
dxxi = FullSimplify[D[x, xi]];
dxix = 1/dxxi;
x1 = 0; x2 = 4.5; x3 = 10;
FullSimplify[x]
B = FullSimplify[{{D[N1, xi], D[N2, xi], D[N3, xi]}}]
K1 = FullSimplify[EA*dxix*Transpose[B].B];
Kb = Integrate[K1, {xi, -1, 1}];
Kb // MatrixForm
a = Integrate[p*N1*dxxi, {xi, -1, 1}];
b = Integrate[p*N2*dxxi, {xi, -1., 1}];
c = Integrate[p*N3*dxxi, {xi, -1, 1}];
ff = {a, b, c};
Chop[ff] // MatrixForm

Two Dimensional 4-node Bilinear Isoparametric Element

Similar to the one dimensional case, an irregular quadrilateral can be defined with four nodes such that the only requirement is that the node numbering is counterclockwise and that the element area is positive. Let (x_1,y_1), (x_2,y_2), (x_3,y_3), and (x_4,y_4) be the two dimensional coordinates of nodes 1, 2, 3, and 4, respectively. The element coordinate system \xi and \eta is chosen such that the coordinates of nodes 1, 2, 3, and 4 are: (-1,-1), (1,-1), (1,1), and (-1,1), respectively (Figure 2). The mapping between the two coordinate systems is performed using the same shape functions for the bilinear quadrilateral as follows:

    \[ \begin{split} x&=N_1x_1+N_2x_2+N_3x_3+N_4x_4\\ y&=N_1y_1+N_2y_2+N_3y_3+N_4y_4 \end{split} \]

where

    \[\begin{split} N_1&=\frac{(1-\xi)(1-\eta)}{4}\\ N_2&=\frac{(1+\xi)(1-\eta)}{4}\\ N_3&=\frac{(1+\xi)(1+\eta)}{4}\\ N_4&=\frac{(1-\xi)(1+\eta)}{4} \end{split} \]

Figure 2. Isoparametric Mapping in the Two Dimensional Case of a 4-node Bilinear Element

Figure 2. Isoparametric Mapping in the Two Dimensional Case of a 4-node Bilinear Element

The requirement that the element area is positive and that the node numbering is counterclockwise ensures that the Jacobian J of the coordinate transformation is invertible at every point within the element. J and its inverse J^{-1} have the forms:

    \[\begin{split} J&=\left(\begin{matrix}\frac{\partial x}{\partial \xi} & \frac{\partial x}{\partial \eta} \\ \frac{\partial y}{\partial \xi} & \frac{\partial y}{\partial \eta}\end{matrix}\right)\\ J^{-1}&=\left(\begin{matrix}\frac{\partial \xi}{\partial x} & \frac{\partial \xi}{\partial y} \\ \frac{\partial \eta}{\partial x} & \frac{\partial \eta}{\partial y}\end{matrix}\right) \end{split} \]

To demonstrate the role of these matrices, let v_{\xi,\eta}=\{\Delta \xi,\Delta \eta\} represent an infinitesimal vector in the element coordinate system \xi, and \eta. The coordinates of this vector in the general coordinate system x, and y are given by v_{x,y}=\{\Delta x,\Delta y\} such that:

    \[ \left(\begin{array}{c}\Delta x\\\Delta y\end{array}\right)=\left(\begin{matrix}\frac{\partial x}{\partial \xi} & \frac{\partial x}{\partial \eta} \\ \frac{\partial y}{\partial \xi} & \frac{\partial y}{\partial \eta}\end{matrix}\right)\left(\begin{array}{c}\Delta \xi\\\Delta \eta\end{array}\right) \]

i.e.,

    \[ v_{x,y}=Jv_{\xi,\eta} \]

Similarly,

    \[ v_{\xi,\eta}=J^{-1}v_{x,y} \]

The determinant of J can be used to relate areas in the coordinate system defined by x and y and the coordinate system defined by \xi and \eta (See the determinant section). From Figure 2, let d_1=\{\mathrm{d}\xi,0\} and d_2=\{0,\mathrm{d}\eta\} represent two vectors inscribing an area \mathrm{d}a=\mathrm{d}\xi\mathrm{d}\eta. Let dx_1=Jd_1, and dx_2=Jd_2 be the same two vectors in the coordinate system defined by x and y. Let the area measure in this coordinate system be \mathrm{d}A. The relationship between the two area measures is given by:

    \[ \mathrm{d}A=\det{J} \mathrm{d}a=\det{J}\mathrm{d}\xi\mathrm{d}\eta \]

This formula can be used to transform integration over the x and y coordinate system to integration over the \xi, and \eta coordinate system. The displacement function in \xi and \eta coordinate system is given by:

    \[\begin{split} \left(\begin{array}{c}u\\v\end{array}\right)&=\left(\begin{array}{cccccccc}N_1&0&N_2&0&N_3&0&N_4&0\\0&N_1&0&N_2&0&N_3&0&N_4\end{array}\right)\left(\begin{array}{c}u_1\\v_1\\u_2\\v_2\\u_3\\v_3\\u_4\\v_4\end{array}\right)\\ &=Nu_e \end{split} \]

To evaluate the strains, the gradients of the displacement functions in the two different coordinate systems need to be evaluated. The gradients of the displacement functions in the \xi and \eta coordinate system have the following form:

    \[ \left(\begin{array}{c}\frac{\partial u}{\partial \xi} \\ \frac{\partial u}{\partial \eta} \\ \frac{\partial v}{\partial \xi} \\ \frac{\partial v}{\partial \eta}\end{array}\right)= \left(\begin{array}{cccccccc} \frac{\partial N_1}{\partial \xi} & 0& \frac{\partial N_2}{\partial \xi} & 0 & \frac{\partial N_3}{\partial \xi} & 0 & \frac{\partial N_4}{\partial \xi} & 0\\ \frac{\partial N_1}{\partial \eta} & 0& \frac{\partial N_2}{\partial \eta} & 0 & \frac{\partial N_3}{\partial \eta} & 0 & \frac{\partial N_4}{\partial \eta}& 0\\ 0&\frac{\partial N_1}{\partial \xi} & 0& \frac{\partial N_2}{\partial \xi} & 0 & \frac{\partial N_3}{\partial \xi} & 0 & \frac{\partial N_4}{\partial \xi} \\ 0 &\frac{\partial N_1}{\partial \eta} & 0& \frac{\partial N_2}{\partial \eta} & 0 & \frac{\partial N_3}{\partial \eta} & 0 & \frac{\partial N_4}{\partial \eta} \end{array}\right) \left(\begin{array}{c}u_1\\v_1\\u_2\\v_2\\u_3\\v_3\\u_4\\v_4\end{array}\right) \]

The gradients of the displacements functions in the x and y coordinate system can be evaluated as functions of the gradients of the displacements in the \xi and \eta coordinate system as follows:

    \[\begin{split} \left(\begin{array}{c}\frac{\partial u}{\partial x} \\ \frac{\partial u}{\partial y} \\ \frac{\partial v}{\partial x} \\ \frac{\partial v}{\partial y}\end{array}\right) & = \left(\begin{array}{cccc}\frac{\partial \xi}{\partial x} & \frac{\partial \eta}{\partial x} & 0 & 0\\ \frac{\partial \xi}{\partial y} & \frac{\partial \eta}{\partial y} & 0 & 0\\ 0& 0&\frac{\partial \xi}{\partial x} & \frac{\partial \eta}{\partial x} \\ 0&0&\frac{\partial \xi}{\partial y} & \frac{\partial \eta}{\partial y}\end{array}\right) \left(\begin{array}{c}\frac{\partial u}{\partial \xi} \\ \frac{\partial u}{\partial \eta} \\ \frac{\partial v}{\partial \xi} \\ \frac{\partial v}{\partial \eta}\end{array}\right)\\ &=\left(\begin{array}{cccc}\frac{\partial \xi}{\partial x} & \frac{\partial \eta}{\partial x} & 0 & 0\\ \frac{\partial \xi}{\partial y} & \frac{\partial \eta}{\partial y} & 0 & 0\\ 0& 0&\frac{\partial \xi}{\partial x} & \frac{\partial \eta}{\partial x} \\ 0&0&\frac{\partial \xi}{\partial y} & \frac{\partial \eta}{\partial y}\end{array}\right) \left(\begin{array}{cccccccc} \frac{\partial N_1}{\partial \xi} & 0 &\frac{\partial N_2}{\partial \xi} & 0 & \frac{\partial N_3}{\partial \xi} & 0 & \frac{\partial N_4}{\partial \xi} & 0\\ \frac{\partial N_1}{\partial \eta} & 0& \frac{\partial N_2}{\partial \eta} & 0 & \frac{\partial N_3}{\partial \eta} & 0 & \frac{\partial N_4}{\partial \eta}& 0\\ 0&\frac{\partial N_1}{\partial \xi} & 0& \frac{\partial N_2}{\partial \xi} & 0 & \frac{\partial N_3}{\partial \xi} & 0 & \frac{\partial N_4}{\partial \xi} \\ 0 &\frac{\partial N_1}{\partial \eta} & 0& \frac{\partial N_2}{\partial \eta} & 0 & \frac{\partial N_3}{\partial \eta} & 0 & \frac{\partial N_4}{\partial \eta} \end{array}\right) \left(\begin{array}{c}u_1\\v_1\\u_2\\v_2\\u_3\\v_3\\u_4\\v_4\end{array}\right) \end{split} \]

The above relations allow the calculation of the two dimensional strains as follows:

    \[ \varepsilon=\left(\begin{array}{c}\varepsilon_{11}\\\varepsilon_{22}\\\gamma_{12}\end{array}\right)=\left(\begin{array}{c}\frac{\partial u}{\partial x} \\ \frac{\partial v}{\partial y} \\ \frac{\partial u}{\partial y} + \frac{\partial v}{\partial x}\end{array}\right) =\left(\begin{array}{cccc}1 &0& 0& 0\\0&0&0&1\\0&1&1&0\end{array}\right)\left(\begin{array}{c}\frac{\partial u}{\partial x} \\ \frac{\partial u}{\partial y} \\ \frac{\partial v}{\partial x} \\ \frac{\partial v}{\partial y}\end{array}\right) \]

Setting:

    \[ B=\left(\begin{array}{cccc}1 &0& 0& 0\\0&0&0&1\\0&1&1&0\end{array}\right)\left(\begin{array}{cccc}\frac{\partial \xi}{\partial x} & \frac{\partial \eta}{\partial x} & 0 & 0\\ \frac{\partial \xi}{\partial y} & \frac{\partial \eta}{\partial y} & 0 & 0\\ 0& 0&\frac{\partial \xi}{\partial x} & \frac{\partial \eta}{\partial x} \\ 0&0&\frac{\partial \xi}{\partial y} & \frac{\partial \eta}{\partial y}\end{array}\right) \left(\begin{array}{cccccccc} \frac{\partial N_1}{\partial \xi} & 0 &\frac{\partial N_2}{\partial \xi} & 0 & \frac{\partial N_3}{\partial \xi} & 0 & \frac{\partial N_4}{\partial \xi} & 0\\ \frac{\partial N_1}{\partial \eta} & 0& \frac{\partial N_2}{\partial \eta} & 0 & \frac{\partial N_3}{\partial \eta} & 0 & \frac{\partial N_4}{\partial \eta}& 0\\ 0&\frac{\partial N_1}{\partial \xi} & 0& \frac{\partial N_2}{\partial \xi} & 0 & \frac{\partial N_3}{\partial \xi} & 0 & \frac{\partial N_4}{\partial \xi} \\ 0 &\frac{\partial N_1}{\partial \eta} & 0& \frac{\partial N_2}{\partial \eta} & 0 & \frac{\partial N_3}{\partial \eta} & 0 & \frac{\partial N_4}{\partial \eta} \end{array}\right) \]

Therefore, the strains can be written in terms of the degrees of freedom as follows:

    \[ \varepsilon=Bu_e \]

Then, the stiffness matrix of the element can be evaluated as follows:

    \[ K^e=t\int_{-1}^{1} \int_{-1}^{1} \! B^TCB \det{J} \, \mathrm{d}\xi\mathrm{d}\eta \]

The element nodal forces can be evaluated as follows:

    \[ f^e=t\int_{-1}^{1} \int_{-1}^{1} \! N^T\rho b \det{J} \, \mathrm{d}\xi\mathrm{d}\eta + \int_{\mbox{element surface}}\!N^T t_n \,\mathrm{d}S \]

where t is the thickness of the element, C is the plane strain or plane stress material constitutive relationship, \rho b is the body forces vector per unit volume, and t_n is the traction vector per unit area on the surface \mathrm{d}S. Notice that the surface element \mathrm{d}S is given in the x and y coordinate system and it should be transformed accordingly to the \xi and \eta coordinate system.

Note that the results of this section can be extended in a straightforward manner to the three dimensional 8 node isoparametric brick element.

Example

Find the stiffness matrix and the nodal forces vector for a 4 node isoparametric plane element whose coordinates are (0,0), (1,0.1), (1.2,1.2), and (0.2,1) with a constant body forces vector \rho b=\{\rho b_1, \rho b_2\} and a unit thickness. Assume plane stress conditions with E=1 units and \nu=0.

Solution

The following Mathematica code can be utilized to obtain the stiffness matrix and the nodal forces due to a constant body-force vector per unit volume for any 4 node isoparametric element as a function of the coordinates of the 4 nodes of the element. Study the code carefully. In particular, notice the following: If the nodes of the element are changed such that they are not defined in a counterclockwise fashion, then \det{J} becomes negative and the stiffness matrix has negative diagonal entries (The stiffness matrix is not positive definite in that case!). Most finite element analysis software would give a warning message that such element is not adequately defined since the area of the element is negative. Also, notice that the stiffness matrix is obtained using the built-in numerical integration function NIntegrate in Mathematica. Try switching to the builtin exact integration function Integrate to notice the difference in the computational resources required. The numerical integration is a much faster technique.

View Mathematica Code
coordinates = {{0, 0}, {1, 0.1}, {1.2, 1.2}, {0.2, 1}};
t = 1;
(*Plane Stress*)
Ee = 1;
nu = 0;
Cc = Ee/(1 + nu)/(1 - nu)*{{1, nu, 0}, {nu, 1, 0}, {0, 0, (1 - nu)/2}};
a = Polygon[{coordinates[[1]], coordinates[[2]], coordinates[[3]], coordinates[[4]], coordinates[[1]]}];
Graphics[a]
Shapefun = Table[0, {i, 1, 4}];
Shapefun[[1]] = (1 - xi) (1 - eta)/4;
Shapefun[[2]] = (1 + xi) (1 - eta)/4;
Shapefun[[3]] = (1 + xi) (1 + eta)/4;
Shapefun[[4]] = (1 - xi) (1 + eta)/4;
Nn = Table[0, {i, 1, 2}, {j, 1, 8}];
Do[Nn[[1, 2 i - 1]] = Nn[[2, 2 i]] = Shapefun[[i]], {i, 1, 4}];
x = Sum[Shapefun[[i]]*coordinates[[i, 1]], {i, 1, 4}];
y = Sum[Shapefun[[i]]*coordinates[[i, 2]], {i, 1, 4}];
J = Table[0, {i, 1, 2}, {j, 1, 2}];
J[[1, 1]] = D[x, xi];
J[[1, 2]] = D[x, eta];
J[[2, 1]] = D[y, xi];
J[[2, 2]] = D[y, eta];
JinvT = Transpose[Inverse[J]];
mat1 = {{1, 0, 0, 0}, {0, 0, 0, 1}, {0, 1, 1, 0}};
mat2 = {{JinvT[[1, 1]], JinvT[[1, 2]], 0, 0}, {JinvT[[2, 1]], JinvT[[2, 2]], 0, 0}, {0, 0, JinvT[[1, 1]], JinvT[[1, 2]]}, {0, 0,JinvT[[2, 1]], JinvT[[2, 2]]}};
mat3 = Table[0, {i, 1, 4}, {j, 1, 8}];
Do[mat3[[1, 2 i - 1]] = mat3[[3, 2 i]] = D[Shapefun[[i]], xi]; 
  mat3[[2, 2 i - 1]] = mat3[[4, 2 i]] = D[Shapefun[[i]], eta], {i, 1, 
   4}];
B = Chop[Simplify[mat1.mat2.mat3]];
k1 = Simplify[Transpose[B].Cc.B];
K = t*NIntegrate[k1*Det[J], {xi, -1, 1}, {eta, -1, 1}];
rb = {rb1, rb2};
f3 = Integrate[Transpose[Nn].rb*Det[J], {xi, -1, 1}, {eta, -1, 1}];
f3 // MatrixForm
K // MatrixForm

Two Dimensional 8-node Quadratic Isoparametric Element

The 8-node quadratic isoparametric plane element is better suited for curved surfaces since the mapping function between the coordinate system of \xi and \eta and the coordinate system x and y uses the quadratic interpolation functions (Figure 3). The interpolation functions in the \xi and \eta coordinate system are given by:

    \[\begin{split} N_1&=-\frac{(1-\eta)(1-\xi)(1+\xi+\eta)}{4}\\ N_2&=-\frac{(1-\eta)(1+\xi)(1-\xi+\eta)}{4}\\ N_3&=-\frac{(1+\eta)(1+\xi)(1-\xi-\eta)}{4}\\ N_4&=-\frac{(1+\eta)(1-\xi)(1+\xi-\eta)}{4}\\ N_5&=\frac{(1-\eta)(1-\xi)(1+\xi)}{2}\\ N_6&=\frac{(1+\xi)(1-\eta)(1+\eta)}{2}\\ N_7&=\frac{(1+\eta)(1-\xi)(1+\xi)}{2}\\ N_8&=\frac{(1-\xi)(1-\eta)(1+\eta)}{2} \end{split} \]

Figure 3. Isoparametric Mapping in the 8 Node Quadratic Two Dimensional Element

Figure 3. Isoparametric Mapping in the 8 Node Quadratic Two Dimensional Element

The development of the equations is similar to the previous section and is left to the reader. The only difference is the number of degrees of freedom of the element. The extension to the three dimensional case (the 20-node brick element) is also straightforward. The following is the Mathematica code that can be utilized to produce the stiffness matrix and the nodal forces due to a constant body forces vector for the 8-node isoparametric element as a function of the coordinates. Notice that the numerical integration built-in function in Mathematica has a higher accuracy than the traditional Gauss integration scheme (which will be presented later) but requires higher computational resources.

View Mathematica Code
coordinates = {{0, 0}, {1, 0.1}, {1.2, 1.2}, {0.2, 1}, {0.5, 0}, {1.1,0.65}, {0.7, 1.2}, {0.1, 0.5}};
t = 1;
Ee = 1;
nu = 0.3;
(*Plane Stress*)
Cc=Ee/(1+nu)/(1-nu)*{{1,nu,0},{nu,1,0},{0,0,(1-nu)/2}};
a=Polygon[{coordinates[[1]],coordinates[[5]],coordinates[[2]],coordinates[[6]],coordinates[[3]],coordinates[[7]],coordinates[[4]],coordinates[[8]],coordinates[[1]]}];
Graphics[a]
Shapefun=Table[0,{i,1,8}];
Shapefun[[1]]=-(1-eta)(1-xi)(1+xi+eta)/4;
Shapefun[[2]]=-(1-eta)(1+xi)(1-xi+eta)/4;
Shapefun[[3]]=-(1+eta)(1+xi)(1-xi-eta)/4;
Shapefun[[4]]=-(1+eta)(1-xi)(1+xi-eta)/4;
Shapefun[[5]]=(1-eta)(1-xi)(1+xi)/2;
Shapefun[[6]]=(1+xi)(1-eta)(1+eta)/2;
Shapefun[[7]]=(1+eta)(1-xi)(1+xi)/2;
Shapefun[[8]]=(1-xi)(1-eta)(1+eta)/2;
Nn=Table[0,{i,1,2},{j,1,16}];
Do[Nn[[1,2i-1]]=Nn[[2,2i]]=Shapefun[[i]],{i,1,8}];
x=Sum[Shapefun[[i]]*coordinates[[i,1]],{i,1,8}];
y=Sum[Shapefun[[i]]*coordinates[[i,2]],{i,1,8}];
J=Table[0,{i,1,2},{j,1,2}];
J[[1,1]]=D[x,xi];
J[[1,2]]=D[x,eta];
J[[2,1]]=D[y,xi];
J[[2,2]]=D[y,eta];
JinvT=Transpose[Inverse[J]];
mat1={{1,0,0,0},{0,0,0,1},{0,1,1,0}};
mat2={{JinvT[[1,1]],JinvT[[1,2]],0,0},{JinvT[[2,1]],JinvT[[2,2]],0,0},{0,0,JinvT[[1,1]],JinvT[[1,2]]},{0,0,JinvT[[2,1]],JinvT[[2,2]]}};
mat3=Table[0,{i,1,4},{j,1,16}];
Do[mat3[[1,2i-1]]=mat3[[3,2i]]=D[Shapefun[[i]],xi];mat3[[2,2i-1]]=mat3[[4,2i]]=D[Shapefun[[i]],eta],{i,1,8}];
B=Chop[FullSimplify[mat1.mat2.mat3]];
k1=FullSimplify[Transpose[B].Cc.B];
Ka=t*NIntegrate[k1*Det[J],{xi,-1,1},{eta,-1,1},AccuracyGoal->8];
Ka//MatrixForm
rb={rb1,rb2};
f3 = Chop[NIntegrate[Transpose[Nn]*Det[J], {xi, -1, 1}, {eta, -1, 1}, AccuracyGoal -> 8].rb];
f3//MatrixForm

Two Dimensional 3 node and 6 node Triangle Isoparametric Elements

The development of the isoparametric triangle elements is also similar to the previous section. The isoparametric mapping between the coordinate system of \xi and \eta and the coordinate system of x and y is shown in Figure 4 for the two types of triangular elements studied (Linear and quadratic). The development of the equations follows the previous section. The shape functions in the \xi and \eta coordinate system for the triangular elements are similar to the ones obtained previously.



Figure 4. Two Dimensional Isoparametric Mapping for the 3 node and 6 node Triangle Elements

Figure 4. Two Dimensional Isoparametric Mapping for the 3 node and 6 node Triangle Elements

Gauss Integration

The Gauss integration scheme is a very efficient method to perform numerical integration over regular domains. In fact, if the function to be integrated is a polynomial of an appropriate degree, then the Gauss integration scheme produces exact results. The Gauss integration scheme has been implemented in almost every finite element analysis software due to its simplicity and computational efficiency. This section outlines the basic principles behind the Gauss integration scheme, along with its application in the finite element analysis method. In particular, this section will discuss the development of the reduced-integration elements and the benefits and drawbacks of utilizing them in finite element analysis models.

Motivation

The isoparametric elements developed above imply that all the integrated functions are defined such that f:[-1,1]\rightarrow\mathbb{R}. In general, numerical integration of such a function has the form:

    \[ \int_{-1}^{1} \! f\,\mathrm{d}x=\sum_{i=1}^n \Delta x f(x_i)=\sum_{i=1}^n w_i f(x_i) \]

For the Gauss integration method, x_i is called an integration point and w_i is called the associated weight. Now, if f is always affine, i.e., f(x)=ax+b then:

    \[ \int_{-1}^{1} \! ax+b\,\mathrm{d}x=2b+a\frac{x^2}{2}\bigg|_{-1}^1=2b=2f(0) \]

So, for functions that are very close to being affine, a numerical integration scheme with 1 integration point that is x_1=0 with an associated weight of 2 can be employed.

Gauss Integration over One Dimensional Domains

Let f:[-1,1]\rightarrow \mathbb{R}. Then, the definite integral of f over the domain x\in[-1,1] can be approximated by the following formula:

    \[ \int_{-1}^{1} \! f\,\mathrm{d}x=\sum_{i=1}^n w_i f(x_i) \]

where \forall i:x_i\in[-1,1] and w_i\in\mathbb{R} is a weight factor associated with each x_i. x_i are termed the integration points. The number of integration points required to obtain sufficient accuracy depend on the shape or form of the function f. If f is a linear function, then one integration point is enough to exactly integrate f. If f is a cubic polynomial function, then two integration points are sufficient to obtain the exact solution. Figure 5 shows a graphical representation of the Gauss numerical integration scheme for integrating polynomials up to the fifth degree. Following is a table of the number of integration points required to obtain the exact solution for a polynomial function f with different polynomial degrees:

Isotable

Figure 5. Gauss Integration for Polynomials of the 1st, 3rd and 5th Degrees

Figure 5. Gauss Integration for Polynomials of the 1st, 3rd and 5th Degrees

Example

Verify that the exact integral of a general cubic polynomial on the interval [-1, 1] can be obtained by the two integration points x_1=\frac{1}{\sqrt{3}} and x_2=\frac{-1}{\sqrt{3}} with the two weight factors w_1=w_2=1.

Solution

Let f=a_0+a_1x+a_2x^2+a_3x^3. The exact integral is equal to:

    \[ \int_{-1}^1\! f\,\mathrm{d}x=2a_0+\frac{2a_2}{3} \]

On the other hand, Gauss integration with two integration points is equal to:

    \[ \int_{-1}^1\! f\,\mathrm{d}x= w_1 (a_0+a_1x_1+a_2x_1^2+a_3x_1^3)+w_2 (a_0+a_1x_2+a_2x_2^2+a_3x_2^3) \]

Equating the Gauss integration with the exact solution yields the following four equations:

    \[ \begin{split} w_1+w_2&=2\\ w_1x_1+w_2x_2&=0\\ w_1x_1^2+w_2x_2^2&=\frac{2}{3}\\ w_1x_1^3+w_2x_2^3&=0 \end{split} \]

These four equations are obtained by noticing that the exact solution and the Gauss integration solution should be equal given any values for the coefficients a_0, a_1, a_2, and a_3. The solution to the above four equations yields x_1=\frac{1}{\sqrt{3}}, x_2=\frac{-1}{\sqrt{3}}, w_1=w_2=1. The Mathematica code to solve the four equations is:

View Mathematica Code
Solve[{w1 + w2 == 2, w1*x1^3 + w2*x2^3 == 0, w1*x1^2 + w2*x2^2 == 2/3,w1*x1 + w2*x2 == 0}, {w1, w2, x1, x2}]

Note that this example can be extended to show that Gauss integration in one dimension with n points can exactly integrate a polynomial of order 2n-1.

Gauss Integration Over Two Dimensional Domains

Let f:[-1,1]\times[-1,1]\rightarrow \mathbb{R}. Then, the integral of f over the domain can be approximated by the following formula:

    \[ \int_{-1}^1\int_{-1}^1 \! f\,\mathrm{d}x\mathrm{d}y=\sum_{j=1}^n\sum_{i=1}^n w_iw_j f(x_i,y_j) \]

where \forall i,j\leq n:(x_i,y_j) is an integration point while w_i and w_j are weight factors associated with x_i and y_j respectively. It can be easily seen that this is a direct consequence of the one dimensional Gauss integration formula as follows:

    \[ \int_{-1}^1\int_{-1}^1 \! f\,\mathrm{d}x\mathrm{d}y=\int_{-1}^1 \! \sum_{i=1}^nw_if(x_i,y)\,\mathrm{d}y==\sum_{j=1}^n\sum_{i=1}^n w_iw_j f(x_i,y_j) \]

Gauss Integration for Two Dimensional Quadrilateral Isoparametric Elements

As described for the two dimensional linear quadrilateral element, the strains of the 4 node elements are linear, and thus, the entries in the matrix B^TCB are quadratic functions of position. Thus, a 2\times 2 numerical integration scheme is required to produce accurate integrals for the isoparametric 4-node elements (Figure 6). On the other hand, the strains in the two dimensional quadratic quadrilateral are quadratic functions of position, and thus, the entries in the matrix B^TCB are fourth-order polynomials of position. Thus, a 3\times 3 numerical integration scheme is required to produce accurate integrals for the 8 node isoparametric element (Figure 6). It should be noted that this can extended the three dimensional case in a straightforward manner. The 8 node brick element has 8 integration points, while the 20 node brick element has 27 integration points.


Figure 6. Gauss Integration Points for integration (a) 4 node isoparametric elements .... with the weight factors...... (b) 8 node isoparametric elements ..... with the weight factors associated with the coordinates respectively.

Figure 6. Gauss integration points for integrating (a) 4 node isoparametric elements \xi,\eta\in\{\frac{1}{\sqrt{3}},\frac{-1}{\sqrt{3}}\} with the weight factors w_i=w_j=1 (b) 8 node isoparametric elements \xi,\eta\in\{-\sqrt{0.6},0,\sqrt{0.6}\} with the weight factors w_i of \frac{5}{9},\frac{8}{9},\frac{5}{9} associated with the coordinates -\sqrt{0.6},0,\sqrt{0.6} respectively.

Reduced Gauss Integration (Under Integration) for Two Dimensional Quadrilateral Isoparametric Elements

Elements that are integrated as per the previous section are termed “full integration” elements. For example, to integrate the 8 \times 8 stiffness matrix of the 4-node plane isoparametric quadrilateral element, each entry will be evaluated at the four locations (2\times 2) of the integration points in the isoparametric element. Thus, the number of computations (excluding the addition) for the integration of the element are 8 \times 8 \times 4 = 256 computations. The number of computations for the 16 \times 16 stiffness matrix of the 8-node plane isoparametric quadrilateral element is 2,304 computations since it uses a 3 \times 3 Gauss integration points scheme! To save computational resources, a “reduced integration” option for each of those elements is available in finite element analysis software. The reduced integration version of the 4-node plane isoparametric quadrilateral element uses only one integration point in the center of the element (Figure 7), and the number of computations is reduced from 256 to 64 computations. However, this is accompanied by a reduction in accuracy. Figure 7 shows two nonlinear functions to be integrated using one integration point vs. two integration points. In the first case, the function has a small variation along the domain and the integration seems to be accurate. However, in the second case, the function has a large variation and is equal to zero in the centre of the domain. The reduced integration will predict zero stiffness!

For the 8-node plane isoparametric quadrilateral, a 2 \times 2 Gauss integration points scheme is used instead of the 3 \times 3 Gauss integration points scheme, which reduces the number of computations from 2,304 to 1,024 computations (Figure 8). The benefit of using reduced integration is not only in saving computational resources but also in balancing out the over-stiffness introduced by assuming a certain deformation field within an element. In fact, reduced integration can lead to better performance in most cases. However, while reduced integration can dramatically reduce the computational time, it can lead to a structure that is too flexible. In fact, the displacements are over predicted when a structure is discretized using reduced integration elements. The reason, as will be shown in the following section, is because some information is lost during sampling of the stiffness at the fewer integration points of the reduced integration elements.


Figure 8

Figure 7. Full vs. reduced integration in the 4-node quadrilateral element




Figure 8. Full vs. reduced integration in the 8-node quadrilateral element

Figure 8. Full vs. reduced integration in the 8-node quadrilateral element

Spurious Modes Associated with Reduced Integration Elements

Spurious modes, or sometimes called zero energy modes, are combinations of displacements that in general produce strains; however, such strains cannot be detected at the integration or sampling points. So, the internal energy calculated at the integration points is equal to zero for such mode. In other words, a very small amount of load would be enough to produce a very large displacement. If a spurious mode exists for the whole structure, it becomes unstable. Usually spurious modes exist in under-integrated element (reduced integration). This is because the reduced number of integration points cannot detect a spurious strain mode. The two very common spurious modes are those for the bilinear quadrilateral (4-node plane element) and the quadratic displacement quadrilateral (8-node plane element). In general, any bending mode is a spurious mode in a 4-node reduced-integration linear element. The strain modes shown cannot be detected at the integration point; the strain at the center of the element (at the location of the single integration point) is equal to zero (Figure 9). On the other hand, the spurious mode of an 8 node reduced integration quadratic element is called “the hour glass” mode (Figure 10) and is produced with the following configuration of nodal displacements:

    \[ u_1=u_7=v_5=v_7=-1\qquad u_3=u_5=v_1=v_3=1\qquad v_6=u_8=\frac{1}{2}=-u_4=-v_2 \]



Figure 9. Spurious modes of the 4 node reduced integration plane quadrilateral element

Figure 9. Spurious modes of the 4 node reduced integration plane quadrilateral element




Figure 10. Spurious mode of an 8 node reduced integration plane quadrilateral element (Hour glass mode).

Figure 10. Spurious mode of an 8 node reduced integration plane quadrilateral element (Hour glass mode).

The spurious modes will be illustrated using the following problem.

Example

Find the internal energy associated with the shown displacement for the 4-node plane stress quadrilateral element shown below, assuming a bilinear displacement interpolation function.
Use exact integration, full integration (2 \times 2 integration points), and reduced integration (1 integration point).
Iso7

Solution

The shape functions of the shown element are:

    \[\begin{split} N_1&=\frac{(1-\xi)(1-\eta)}{4}\\ N_2&=\frac{(1+\xi)(1-\eta)}{4}\\ N_3&=\frac{(1+\xi)(1+\eta)}{4}\\ N_4&=\frac{(1-\xi)(1+\eta)}{4} \end{split} \]

The displacement functions u and v in the directions \xi and \eta, respectively, are:

    \[ u=u_sN_1-u_sN_2+u_sN_3-u_sN_4\qquad v=0 \]

Therefore, the strains have the form:

    \[ \left(\begin{array}{c}\varepsilon_{\xi\xi}\\\varepsilon_{\eta\eta}\\\gamma_{\xi\eta}\end{array}\right)=\left(\begin{array}{c} \frac{\partial N_1}{\partial \xi}-\frac{\partial N_2}{\partial \xi}+\frac{\partial N_3}{\partial \xi}-\frac{\partial N_4}{\partial \xi} \\0\\ \frac{\partial N_1}{\partial \eta}-\frac{\partial N_2}{\partial \eta}+\frac{\partial N_3}{\partial \eta}-\frac{\partial N_4}{\partial \eta}\end{array}\right)u_s=Bu_s \]

The stress strain relationship for plane stress is given by:

    \[ \left(\begin{array}{c}\sigma_{\xi\xi}\\\sigma_{\eta\eta}\\\sigma_{\xi\eta}\end{array}\right)=\frac{E}{(1+\nu)(1-\nu)}\left(\begin{matrix}1 &\nu&0\\\nu&1&0\\0&0&\frac{(1-\nu)}{2}\end{matrix}\right)\left(\begin{array}{c}\varepsilon_{\xi\xi}\\\varepsilon_{\eta\eta}\\\gamma_{\xi\eta}\end{array}\right) \]

The stiffness matrix K can be evaluated using the following integral where t is the element thickness:

    \[ K=t\int_{-1}^{1}\int_{-1}^{1}\! B^TCB\,\mathrm{d}\xi\mathrm{d}\eta = \frac{tE}{2(\nu^2-1)}\int_{-1}^{1}\int_{-1}^{1}\! (\nu-1)\xi^2-2\eta^2\,\mathrm{d}\xi\mathrm{d}\eta \]

Notice that the system has only one degree of freedom (the unknown displacement variable u_s), and therefore, the stiffness matrix has the dimensions 1\times 1.

Using Exact Integration:
The stiffness matrix evaluated using exact integration is:

    \[ K=\frac{2E(\nu-3)}{3(\nu^2-1)} \]

The strain energy associated with this displacement mode is equal to:

    \[ \frac{1}{2}Ku_s^2=\frac{2E(\nu-3)}{6(\nu^2-1)}u_s^2 \]

Using Full Integration:
The stiffness matrix evaluated using full integration is:

    \[ K=\frac{tE}{2(\nu^2-1)}\left(\sum_{i,j=1}^2 w_iw_j \left((\nu-1)\xi_i^2-2\eta_j^2\right)\right)=\frac{2E(\nu-3)}{3(\nu^2-1)} \]

where \xi,\eta\in\{\frac{1}{\sqrt{3}},\frac{-1}{\sqrt{3}}\} and w_i=w_j=1 which gives the same results as the exact integration.

Using Reduced Integration:
The stiffness matrix evaluated using reduced integration is:

    \[ K=\frac{tE}{2(\nu^2-1)}\left(\sum_{i,j=1}^1 w_iw_j \left((\nu-1)\xi_i^2-2\eta_j^2\right)\right)=0 \]

where \xi,\eta\in\{0\} and w_i=w_j=2. The strain energy associated with this displacement mode is equal to zero. Thus, the reduced integration scheme predicts that the displacement field shown does not require any external forces applied to the element. Thus, this element is unstable if loads produce the displacement mode shown, i.e., bending! The following is the Mathematica code used for the above calculations:

View Mathematica Code
N1=(1-xi)(1-eta)/4;
N2=(1+xi)(1-eta)/4;
N3=(1+xi)(1+eta)/4;
N4=(1-xi)(1+eta)/4;
B={D[N1-N2+N3-N4,xi],0,D[N1-N2+N3-N4,eta]};
Cc=Ee/(1+nu)/(1-nu)*{{1,nu,0},{nu,1,0},{0,0,(1-nu)/2}};
Kbeforeintegration=FullSimplify[B.Cc.B]
(*Using Exact Integration*)
Kexact=Integrate[Kbeforeintegration,{eta,-1,1},{xi,-1,1}]
Energyexact=us*Kexact*us/2
(*Using Full Integration*)
etaset={1/Sqrt[3],-1/Sqrt[3]};
xiset={1/Sqrt[3],-1/Sqrt[3]};
Kfull=FullSimplify[Sum[(Kbeforeintegration/.{xi->xiset[[i]],eta-
>etaset[[j]]}),{i,1,2},{j,1,2}]]
Energyfull=us*Kfull*us/2
(*Using reduced Integration*)
etaset=xiset={0};
Kreduced=Sum[4(Kbeforeintegration/.{xi->xiset[[i]],eta-
>etaset[[j]]}),{i,1,1},{j,1,1}]
Energyreduced=us*Kreduced*us/2

Gauss Integration for Two Dimensional Triangle Isoparametric Elements

The strains of the 3 node triangle element are constant, and thus, the entries in the matrix B^TCB are constant. Thus, a 1 point numerical integration scheme is required to produce accurate integrals for the isoparametric 3-node elements (Figure 11). For the 6 node triangle elements, a 3 points numerical integration scheme is required to produce accurate integrals (Figure 11). It should be noted that the shown locations for the integration points and their number can be extended to the three-dimensional case in a straightforward manner.


Figure 11.

Figure 11. Gauss integration points for integrating (a) 3 node isoparametric triangle elements \xi,\eta\in\{\frac{1}{3}\} and w=\frac{1}{2}. (b) 6 node isoparametric elements (\xi,\eta)\in\{(\frac{1}{6},\frac{2}{3}),(\frac{1}{6},\frac{1}{6}),(\frac{2}{3},\frac{1}{6})\} with equal weight factors w=\frac{1}{6}.

Final Notes

The use of isoparametric elements and numerical integration dramatically increases the robustness of the finite element analysis method. The techniques described in the previous sections allow meshing of irregular domains with compatible elements, i.e., elements that ensure the continuity of the displacement field across boundaries (How?). However, to ensure the continuity of the displacement field, the elements of the mesh should be of the same type; otherwise, the user has to apply mesh constraints to ensure the compatibility of the displacement between elements. For example, a 4-node bilinear displacement element when connected to an 8-node quadratic displacement element will have incompatible displacements at the side of the connection. This can be overcome by constraining the mid-side node of the quadratic element to always lie on the straight side of the bilinear element. Otherwise, the results at the location of discontinuous displacement fields would be erroneous.

The stiffness matrix and the nodal loads of the structure can be calculated using the Gauss numerical integration scheme, which speeds up the computational time required. However, the price of using the Gauss numerical integration is the loss of accuracy. Gauss numerical integration produces accurate results only when the integrals are exact polynomials. The integrals are exact polynomials only when the elements of the mesh are exactly aligned with the isoparametric mapping into the coordinates of \xi and \eta. When the shape of the elements in the mesh deviates from the shape of the element in the coordinate system of \xi and \eta, the functions to be integrated are no longer exact polynomials but rather radicals; thus, the numerical integration scheme is no longer accurate. Most finite element analysis software give warnings when elements are highly distorted since this will lead to inaccurate computations of the stiffness matrix during the numerical integration process. In addition, for nonlinear problems when the stiffness changes with deformation, it is possible that the element will reach a highly distorted shape that the numerical integration scheme is no longer valid.

Leave a Reply

Your email address will not be published.