Engineering at Alberta Courses » Two Dimensional Solid Elements

Finite Element Analysis: Two Dimensional Solid Elements

Triangular Elements

One of the ways to mesh a domain in finite element analysis is using triangular elements. The advantages of using triangular elements is the ability to develop meshing algorithms that can easily mesh any irregular domain with triangular elements. However, triangular elements provide a crude approximation of the variable being interpolated. This can be overcome by either using a very fine mesh of triangular elements or using the quadratic version of triangular elements.

Linear Triangular Element (Constant Strain Triangle)

The linear triangular element is obtained by assuming that the displacement within a triangular shape is a linear function of the displacements at the three corner nodes. This linear triangular element is also called the constant strain triangle, since as will be shown in the derivation below, the strain across the whole element is always constant (the matrix $B$ has constant values). This is a major disadvantage of such element since it tends to be extremely stiff. Approximating the domain with triangular elements can only be considered within good accuracy if the difference in the strains between neighboring elements is relatively small.

Consider the plane triangle shown in Figure 1. For simplicity, two of the triangle sides are assumed to have unit lengths. The assumed (approximate) linear-displacement function across the element has the form:

$u_{2D}=\left(\begin{array}{c}u\\v\end{array}\right)=\left(\begin{array}{c}a_0+a_1X_1+a_2X_2\\b_0+b_1X_1+b_2X_2\end{array}\right)$

where $a_0$ , $a_1$ , $a_2$ , $b_0$ , $b_1$ , and $b_2$ are six generalized degrees of freedom of the element. The approximate displacement function in terms of the nodal degrees of freedom has the form:

$u_{2D}=\left(\begin{array}{c}u\\v\end{array}\right)=\left(\begin{array}{c}N_1u_1+N_2u_2+N_3u_3\\N_1v_1+N_2v_2+N_3v_3\end{array}\right)$

Figure 1. Linear Triangular Element (Constant Strain Triangle)

There are two ways to find the shape functions $N_i$ . The first way is to replace the generalized degrees of freedom with the nodal degrees of freedom in the first equation. The multipliers of the nodal degrees of freedom would then be the shape function. The following mathematica code does that for you:
View Mathematica Code

u = a0 + a1 * X1 + a2*X2;
v = b0 + b1*X1 + b2*X2;
Coordinates = {{X1 -> 0, X2 -> 0}, {X1 -> 1, X2 -> 0}, {X1 -> 0,X2 -> 1}};
Eq1 = u /. Coordinates[[1]];
Eq2 = u /. Coordinates[[2]];
Eq3 = u /. Coordinates[[3]];
Eq4 = v /. Coordinates[[1]];
Eq5 = v /. Coordinates[[2]];
Eq6 = v /. Coordinates[[3]];
a = Solve[{Eq1 == u1, Eq2 == u2, Eq3 == u3, Eq4 == v1, Eq5 == v2, Eq6 == v3}, {a0, a1, a2, b0, b1, b2}]
u = u /. a[[1]];
v = v /. a[[1]];
N1 = Coefficient[u, u1]
N1 = Coefficient[v, v1]
N2 = Coefficient[u, u2]
N2 = Coefficient[v, v2]
N3 = Coefficient[u, u3]
N3 = Coefficient[v, v3]

An alternate way is directly by realizing that $N_1$ has to satisfy that the condition that at node 1, $N_1=1$ , it varies linearly and is equal to zero on the line $1-X_1-X_2=0$ . Similarly, $N_2$ is equal to 1 at node 2 and is equal to zero on the line $X_1=0$ . Also, $N_3$ is equal to 1 at node 3 and is equal to zero on the line $X_2=0$ . These result in the following expressions for the shape functions:

$\begin{split} N_1&=1-X_1-X_2\\ N_2&=X_1\\ N_3&=X_2 \end{split}$

The distribution of the shape functions on the element are illustrated in Figure 2.

Figure 2. Shape functions in linear triangular elements

The strain associated with the assumed displacement field has the following form:

$\begin{split} \varepsilon=\left(\begin{array}{c}\varepsilon_{11}\\\varepsilon_{22}\\2\varepsilon_{12}\end{array}\right) & =\left(\begin{array}{c}\frac{\partial u}{\partial X_1}\\\frac{\partial v}{\partial X_2}\\\frac{\partial u}{\partial X_2}+\frac{\partial v}{\partial X_1}\end{array}\right)\\ &=\left(\begin{array}{cccccc}\frac{\partial N_1}{\partial X_1}&0&\frac{\partial N_2}{\partial X_1}&0&\frac{\partial N_3}{\partial X_1}&0\\0&\frac{\partial N_1}{\partial X_2}&0&\frac{\partial N_2}{\partial X_2}&0&\frac{\partial N_3}{\partial X_2}\\\frac{\partial N_1}{\partial X_2} & \frac{\partial N_1}{\partial X_1} & \frac{\partial N_2}{\partial X_2} & \frac{\partial N_2}{\partial X_1} & \frac{\partial N_3}{\partial X_2} & \frac{\partial N_3}{\partial X_1} \end{array}\right)\left(\begin{array}{c}u_1\\v_1\\u_2\\v_2\\u_3\\v_3\end{array}\right)\\ &=\left(\begin{array}{cccccc}-1&0&1&0&0&0\\0&-1&0&0&0&1\\-1 & -1 & 0& 1 & 1 & 0 \end{array}\right)\left(\begin{array}{c}u_1\\v_1\\u_2\\v_2\\u_3\\v_3\end{array}\right)\\ &=B u_e \end{split}$

Notice that as mentioned above, the values of the strain components are independent of the location inside the triangle; rather the strain components are constant across the element.

Stiffness Matrix

The constitutive relationship of plane linear elastic materials is defined using a $3\times 3$ matrix $C$ that depends on whether the material is in a plane strain or a plane stress state. For the sake of illustration, the plane strain state is chosen here, the material constitutive relationship matrix in that case is:

$\left(\begin{array}{c}\sigma_{11}\\\sigma_{22}\\\sigma_{12}\end{array}\right)=\frac{E}{(1+\nu)(1-2\nu)}\left(\begin{matrix}1-\nu &\nu&0\\\nu & 1-\nu & 0\\0 & 0 & \frac{1-2\nu}{2}\end{matrix}\right)\left(\begin{array}{c}\varepsilon_{11}\\\varepsilon_{22}\\2\varepsilon_{12}\end{array}\right)$

The stiffness matrix of the constant strain triangle can be evaluated by the following integral:

$K^e=\int_e \! B^T C B\,\mathrm{d}X=\int_e \! \left(\begin{array}{ccc}-1&0&-1\\0&-1&-1\\1&0&0\\0&0&1\\0&0&1\\0&1&0\end{array}\right) C \left(\begin{array}{cccccc}-1&0&1&0&0&0\\0&-1&0&0&0&1\\-1 & -1 & 0& 1 & 1 & 0\end{array}\right)\,\mathrm{d}X$

Assuming that the triangle has a constant thickness $t$ , then the differential volume $\mathrm{d}X= t\mathrm{d}X_1\mathrm{d}X_2$ . In addition, all the components of the matrix $B^TCB$ are constant; therefore, the integral can be evaluated by simply multiplying the matrix $B^TCB$ by $\frac{t}{2}$ , which represents the volume of the triangle:

$K^e=\frac{tE}{2(1+\nu)}\left(\begin{matrix}\frac{3-4\nu}{2-4\nu}&\frac{1}{2-4\nu}&\frac{1-\nu}{-1+2\nu}&\frac{-1}{2}&\frac{-1}{2}&\frac{\nu}{-1+2\nu}\\ &\frac{3-4\nu}{2-4\nu}&\frac{\nu}{-1+2\nu}&\frac{-1}{2}&\frac{-1}{2}&\frac{1-\nu}{-1+2\nu}\\ & & \frac{1-\nu}{1-2\nu}& 0 & 0 &\frac{\nu}{1-2\nu}\\ & & &\frac{1}{2} & \frac{1}{2} & 0\\ & \multicolumn{2}{c}{\text{sym.}}& & \frac{1}{2}&0\\ & & &&&\frac{1-\nu}{1-2\nu} \end{matrix}\right)$

View Mathematica Code

B={{-1,0,1,0,0,0},{0,-1,0,0,0,1},{-1,-1,0,1,1,0}};
Cc=Ee/(1+nu)*{{(1-nu)/(1-2nu),nu/(1-2nu),0},{(nu)/(1-2nu),(1-nu)/
(1-2nu),0},{0,0,1/2}};
K=FullSimplify[Transpose[B].Cc.B];
FullSimplify[1/Ee*(1+nu)*K]//MatrixForm

Nodal Forces due to Body Forces

Assuming that the distributed body forces vector per unit mass is constant and is given by $b={b_1, b_2}$ , and that the mass density is $\rho$ , then the nodal forces due to the distributed body forces can be obtained using the integral:

$f^e=\int_e \! N^T \rho b \, \mathrm{d}X=\rho t \int_0^1\int_0^{1-X_2} \! \left(\begin{array}{cc}N_1&0\\0&N_1\\N_2&0\\0&N_2\\N_3&0\\0&N_3\end{array}\right) \left(\begin{array}{c}b_1\\b_2\end{array}\right) \, \mathrm{d}X_1 \mathrm{d}X_2=\frac{\rho t}{6}\left(\begin{array}{c}b_1\\b_2\\b_1\\b_2\\b_1\\b_2\end{array}\right)$

i.e., the constant body forces vector can be represented by lumped equal concentrated loads at the nodes (Figure 3a).
View Mathematica Code

Shapefun=Table[0,{i,1,3}]
Shapefun[[1]]=1-x1-x2;
Shapefun[[2]]=x1;
Shapefun[[3]]=x2;
Nn=Table[0,{i,1,2},{j,1,6}];
Do[Nn[[1,2i-1]]=Nn[[2,2i]]=Shapefun[[i]],{i,1,3}];
rb={rb1,rb2};
Integrate[Transpose[Nn].rb,{x2,0,1},{x1,0,1-x2}]//MatrixForm

Nodal Forces due to Traction Vector on One Side

Assuming that a constant distributed pressure per unit area of $t_n={t)1, t_2}$ is acting on the surface joining nodes 1 and 3, then the nodal forces due to the distributed traction vector can be obtained using the following integral evaluated on the left side $X_1 =0$ :

$f^e=\int_{\partial e} \! N^T t_n \, \mathrm{d}S=t \int_0^1 \! \left(\begin{array}{cc}N_1&0\\0&N_1\\N_2&0\\0&N_2\\N_3&0\\0&N_3\end{array}\right) \left(\begin{array}{c}t_1\\t_2\end{array}\right) \, \mathrm{d}X_2 =\frac{t}{2}\left(\begin{array}{c}t_1\\t_2\\0\\0\\t_1\\t_2\end{array}\right)$

Thus, the distributed constant traction on one side of the triangle can be lumped into equal nodal loads applied on the nodes of that side (Figure 3b).
View Mathematica Code

Shapefun=Table[0,{i,1,3}]
Shapefun[[1]]=1-x1-x2;
Shapefun[[2]]=x1;
Shapefun[[3]]=x2;
Nn=Table[0,{i,1,2},{j,1,6}];
Do[Nn[[1,2i-1]]=Nn[[2,2i]]=Shapefun[[i]],{i,1,3}];
tn={t1,t2};
Integrate[(Transpose[Nn].tn/.x1->0),{x2,0,1}]//MatrixForm

Figure 3. Nodal forces in a linear triangular element due to (a) constant body forces vectors, (b) constant traction vector on one side.

Figure 3. Nodal forces in a linear triangular element with a constant unit thickness due to (a) constant body forces vectors, (b) constant traction vector on one side.

Quadratic Triangular Elements

The quadratic triangular element offers a better approximation to the displacement field within a triangular element by introducing additional nodes on the straight sides of the triangle. The quadratic triangular element is called the linear strain triangle since, as will be shown in the derivation below, the matrix $B$ contains linear expressions in the coordinates $X_1$ and $X_2$ and therefore, this element is capable of modeling linear strains (for example, bending). However, the addition of nodes comes with a higher computational price compared to its linear counterpart.

Consider the plane triangle shown in Figure 4. For simplicity, two of the triangle sides are assumed to have unit lengths. The displacement function on the element can be assumed to have the following form:

$u_{2D}=\left(\begin{array}{c}u\\v\end{array}\right)=\left(\begin{array}{c}a_0+a_1X_1+a_2X_2+a_3X_1X_2+a_4X_1^2+a_5X_2^2\\b_0+b_1X_1+b_2X_2+b_3X_1X_2+b_4X_1^2+b_5X_2^2\end{array}\right)$

where $a_0$ , $a_1$ , $a_2$ , $a_3$ , $a_4$ , $a_5$ , $b_0$ , $b_1$ , $b_2$ , $b_3$ , $b_4$ , and $b_5$ are twelve generalized degrees of freedom of the element. The approximate displacement function in terms of the nodal degrees of freedom has the form:

$u_{2D}=\left(\begin{array}{c}u\\v\end{array}\right)=\left(\begin{array}{c}N_1u_1+N_2u_2+N_3u_3+N_4u_4+N_5u_5+N_6u_6\\N_1v_1+N_2v_2+N_3v_3+N_4v_4+N_5v_5+N_6v_6\end{array}\right)$

Figure 4. Quadratic triangle

u = a0 + a1 * X1 + a2*X2 + a3*X1*X2 + a4*X1^2 + a5*X2^2;
v = b0 + b1*X1 + b2*X2 + b3*X1*X2 + b4*X1^2 + b5*X2^2;
Coordinates = {{X1 -> 0, X2 -> 0}, {X1 -> 1, X2 -> 0}, {X1 -> 0,X2 -> 1}, {X1 -> 1/2, X2 -> 0}, {X1 -> 1/2, X2 -> 1/2}, {X1 -> 0, X2 -> 1/2}};
Eq1 = u /. Coordinates[[1]];
Eq2 = u /. Coordinates[[2]];
Eq3 = u /. Coordinates[[3]];
Eq4 = v /. Coordinates[[1]];
Eq5 = v /. Coordinates[[2]];
Eq6 = v /. Coordinates[[3]];
Eq7 = u /. Coordinates[[4]];
Eq8 = u /. Coordinates[[5]];
Eq9 = u /. Coordinates[[6]];
Eq10 = v /. Coordinates[[4]];
Eq11 = v /. Coordinates[[5]];
Eq12 = v /. Coordinates[[6]];
a = Solve[{Eq1 == u1, Eq2 == u2, Eq3 == u3, Eq4 == v1, Eq5 == v2, Eq6 == v3, Eq7 == u4, Eq8 == u5, Eq9 == u6, Eq10 == v4, Eq11 == v5, Eq12 == v6}, {a0, a1, a2, a3, a4, a5, b0, b1, b2, b3, b4, b5}]
u = u /. a[[1]];
v = v /. a[[1]];
N1 = FullSimplify[Coefficient[u, u1]]
N1 = FullSimplify[Coefficient[v, v1]]
N2 = FullSimplify[Coefficient[u, u2]]
N2 = FullSimplify[Coefficient[v, v2]]
N3 = FullSimplify[Coefficient[u, u3]]
N3 = FullSimplify[Coefficient[v, v3]]
N4 = FullSimplify[Coefficient[u, u4]]
N4 = FullSimplify[Coefficient[v, v4]]
N5 = FullSimplify[Coefficient[u, u5]]
N5 = FullSimplify[Coefficient[v, v5]]
N6 = FullSimplify[Coefficient[u, u6]]
N6 = FullSimplify[Coefficient[v, v6]]

An alternate way is directly by realizing that $N_1$ has to satisfy that the condition that at node 1, $N_1=1$ , and it is equal to zero on the lines $1-X_1-X_2=0$ and $\frac{1}{2}-X_1-X_2=0$ . Similarly, $N_2$ is equal to 1 at node 2 and is equal to zero on the lines $X_1=0$ and $\frac{1}{2}-X_1=0$ . Also, $N_3$ is equal to 1 at node 3 and is equal to zero on the lines $X_2=0$ and $\frac{1}{2}-X_2=0$ . $N_4$ is equal to 1 at node 4 and is equal to zero on the lines $X_1=0$ and $1-X_1-X_2=0$ . $N_5$ is equal to 1 at node 5 and is equal to zero on the lines $X_1=0$ and $X_2=0$ . $N_6$ is equal to 1 at node 6 and is equal to zero on the lines $X_2=0$ and $1-X_1-X_2=0$ . . These result in the following expressions for the shape functions:

$\begin{split} N_1&=(1-X_1-X_2)(1-2X_1-2X_2)\\ N_2&=X_1(2X_1-1)\\ N_3&=X_2(2X_2-1)\\ N_4&=4X_1(1-X_1-X_2)\\ N_5&=4X_1X_2\\ N_6&=4X_2(1-X_1-X_2) \end{split}$

The distribution of the shape functions on the element are illustrated in Figure 5.

Figure 5. Shape function distribution for the quadratic triangular element

The strain associated with the assumed displacement field has the following form:

$\begin{split} \varepsilon&=\left(\begin{array}{c}\varepsilon_{11}\\\varepsilon_{22}\\2\varepsilon_{12}\end{array}\right) =\left(\begin{array}{c}\frac{\partial u}{\partial X_1}\\\frac{\partial v}{\partial X_2}\\\frac{\partial u}{\partial X_2}+\frac{\partial v}{\partial X_1}\end{array}\right)\\ &=\left(\begin{array}{cccccccccccc}\frac{\partial N_1}{\partial X_1}&0&\frac{\partial N_2}{\partial X_1}&0&\frac{\partial N_3}{\partial X_1}&0&\frac{\partial N_4}{\partial X_1}&0&\frac{\partial N_5}{\partial X_1}&0&\frac{\partial N_6}{\partial X_1}&0\\ 0&\frac{\partial N_1}{\partial X_2}&0&\frac{\partial N_2}{\partial X_2}&0&\frac{\partial N_3}{\partial X_2}&0&\frac{\partial N_4}{\partial X_2}&0&\frac{\partial N_5}{\partial X_2}&0&\frac{\partial N_6}{\partial X_2}\\ \frac{\partial N_1}{\partial X_2} & \frac{\partial N_1}{\partial X_1} & \frac{\partial N_2}{\partial X_2} & \frac{\partial N_2}{\partial X_1} & \frac{\partial N_3}{\partial X_2} & \frac{\partial N_3}{\partial X_1} & \frac{\partial N_4}{\partial X_2} & \frac{\partial N_4}{\partial X_1} & \frac{\partial N_5}{\partial X_2} & \frac{\partial N_5}{\partial X_1} & \frac{\partial N_6}{\partial X_2} & \frac{\partial N_6}{\partial X_1} \end{array}\right)\left(\begin{array}{c}u_1\\v_1\\u_2\\v_2\\u_3\\v_3\\u_4\\v_4\\u_5\\v_5\\u_6\\v_6\end{array}\right)\\ &=B u_e \end{split}$

The different entries of the matrix $B$ are indeed linear functions of the coordinates $X_1$ and $X_2$ inside the domain. Thus, this element can be used to accurately model a domain with linear strain across the element.

Stiffness Matrix

The stiffness matrix of the linear strain triangle can be evaluated using the following integral (assuming a constant thickness $t$ ):

$K^e=t\int_0^1\int_0^{1-X_2} \! B^T C B\,\mathrm{d}X_1\mathrm{d}X_2$

The stiffness matrix has the dimensions of $12 \times 12$ , and the following Mathematica code can be utilized to view its components:
View Mathematica Code

Shapefun=Table[0,{i,1,6}];
Shapefun[[1]]=2(1-x1-x2)(1/2-x1-x2);
Shapefun[[2]]=2x1 (x1-1/2);
Shapefun[[3]]=2x2 (x2-1/2);
Shapefun[[4]]=4x1 (1-x1-x2);
Shapefun[[5]]=4x1*x2;
Shapefun[[6]]=4x2*(1-x1-x2);
B=Table[0,{i,1,3},{j,1,12}];
Do[B[[1,2i-1]]=B[[3,2i]]=D[Shapefun[[i]],x1];B[[2,2i]]=B[[3,2i-1]]=D[Shapefun[[i]],x2],{i,1,6}];
Cc=Ee/(1+nu)*{{(1-nu)/(1-2nu),nu/(1-2nu),0},{(nu)/(1-2nu),(1-nu)/(1-2nu),0},{0,0,1/2}};
Kbeforeintegration=t*FullSimplify[Transpose[B].Cc.B];
K=Integrate[Kbeforeintegration,{x2,0,1},{x1,0,1-x2}];
K//MatrixForm

Notice that the integration for this element requires high computational time (Try using Mathematica). If a structure is composed of hundreds of elements, the construction of the stiffness matrix would require high computational resources.

Nodal Forces due to Body Forces

$f^e=\int_e \! N^T \rho b \, \mathrm{d}X=\rho t \int_0^1\int_0^{1-X_2} \! \left(\begin{array}{cc}N_1&0\\0&N_1\\N_2&0\\0&N_2\\N_3&0\\0&N_3\\N_4&0\\0&N_4\\N_5&0\\0&N_5\\N_6&0\\0&N_6\end{array}\right) \left(\begin{array}{c}b_1\\b_2\end{array}\right) \, \mathrm{d}X_1 \mathrm{d}X_2=\frac{\rho t}{6}\left(\begin{array}{c}0\\0\\0\\0\\0\\0\\b_1\\b_2\\b_1\\b_2\\b_1\\b_2\end{array}\right)$

i.e., a constant distributed body forces vector can be lumped into equal loads applied on the mid-side nodes and zero loads applied on the corner nodes (Figure 6a). The following is the Mathematica code used for the above calculations:
View Mathematica Code

Shapefun=Table[0,{i,1,6}];
Shapefun[[1]]=2(1-x1-x2)(1/2-x1-x2);
Shapefun[[2]]=2x1 (x1-1/2);
Shapefun[[3]]=2x2 (x2-1/2);
Shapefun[[4]]=4x1 (1-x1-x2);
Shapefun[[5]]=4x1*x2;
Shapefun[[6]]=4x2*(1-x1-x2);
Nn=Table[0,{i,1,2},{j,1,12}];
Do[Nn[[1,2i-1]]=Nn[[2,2i]]=Shapefun[[i]],{i,1,6}];
rb={rb1,rb2};
fe=Integrate[Transpose[Nn].rb,{x2,0,1},{x1,0,1-x2}];
fe//MatrixForm

Nodal Forces due to Traction Vector on One Side

$f^e=\int_{\partial e} \! N^T t_n \, \mathrm{d}S=t \int_0^1 \! \left(\begin{array}{cc}N_1&0\\0&N_1\\N_2&0\\0&N_2\\N_3&0\\0&N_3\\N_4&0\\0&N_4\\N_5&0\\0&N_5\\N_6&0\\0&N_6\end{array}\right) \left(\begin{array}{c}t_1\\t_2\end{array}\right) \, \mathrm{d}X_2 =\frac{t}{6}\left(\begin{array}{c}t_1\\t_2\\0\\0\\t_1\\t_2\\0\\0\\0\\0\\4t_1\\4t_2\end{array}\right)$

The appropriate nodal loads used to lump a constant traction vector on the left side of the quadratic triangular element are shown in Figure 6b. The following Mathematica code was used for the above calculations:

View Mathematica Code

Shapefun=Table[0,{i,1,3}]
Shapefun[[1]]=1-x1-x2;
Shapefun[[2]]=x1;
Shapefun[[3]]=x2;
Nn=Table[0,{i,1,2},{j,1,6}];
Do[Nn[[1,2i-1]]=Nn[[2,2i]]=Shapefun[[i]],{i,1,3}];
tn={t1,t2};
fe=Integrate[(Transpose[Nn].tn/.x1->0),{x2,0,1}];
fe//MatrixForm

Figure 6. Nodal forces in a quadratic triangular element with a constant unit thickness due to (a) constant body forces vectors, (b) constant traction vector on one side.

Quadrilateral Elements

Linear Quadrilateral Elements

It is usually easier to discretize a domain by using triangular elements; however, when the domain to be discretized has a regular shape, quadrilateral elements can be used to offer a more regular displacement discretization, which, in some cases, might offer a better approximation to the displacement shape. The bilinear quadrilateral element offers a bilinear displacement approximation where the displacement within the element is assumed to vary bilinearly (linear in two directions
within the element). This element, however, behaves in a stiff manner when it is used to model linear strains (bending strains) since, as will be shown in the derivation, it is not possible to model pure bending (pure linear normal strain components in one direction) without introducing additional shear strains (usually termed parasitic shear strains).

Consider the plane quadrilateral shown in Figure 7. The bilinear displacement function for this element is assumed to have the following form:

(1) $\begin{equation*} u_{2D}=\left(\begin{array}{c}u\\v\end{array}\right)=\left(\begin{array}{c}a_0+a_1X_1+a_2X_2+a_3X_1X_2\\b_0+b_1X_1+b_2X_2+b_3X_1X_2\end{array}\right) \end{equation*}$

where $a_0$ , $a_1$ , $a_2$ , $a_3$ , $b_0$ , $b_1$ , $b_2$ , and $b_3$ are eight generalized degrees of freedom of the element. The displacement function in terms of the nodal degrees of freedom has the form:

$u_{2D}=\left(\begin{array}{c}u\\v\end{array}\right)=\left(\begin{array}{c}N_1u_1+N_2u_2+N_3u_3+N_4u_4\\N_1v_1+N_2v_2+N_3v_3+N_4v_4\end{array}\right)$

Figure 7. Bilinear quadrilateral

Clear[a, b]
u = a0 + a1*X1 + a2*X2 + a3*X1*X2;
v = b0 + b1*X1 + b2*X2 + b3*X1*X2;
Coordinates = {{X1 -> -a, X2 -> -b}, {X1 -> a, X2 -> -b}, {X1 -> a, X2 -> b}, {X1 -> -a, X2 -> b}};
Eq1 = u /. Coordinates[[1]];
Eq2 = u /. Coordinates[[2]];
Eq3 = u /. Coordinates[[3]];
Eq4 = u /. Coordinates[[4]];
Eq5 = v /. Coordinates[[1]];
Eq6 = v /. Coordinates[[2]];
Eq7 = v /. Coordinates[[3]];
Eq8 = v /. Coordinates[[4]];
sol = Solve[{Eq1 == u1, Eq2 == u2, Eq3 == u3, Eq4 == u4, Eq5 == v1, Eq6 == v2, Eq7 == v3, Eq8 == v4}, {a0, a1, a2, a3, b0, b1, b2, b3}]
u = u /. sol[[1]];
v = v /. sol[[1]];
N1 = FullSimplify[Coefficient[u, u1]]
N1 = FullSimplify[Coefficient[v, v1]]
N2 = FullSimplify[Coefficient[u, u2]]
N2 = FullSimplify[Coefficient[v, v2]]
N3 = FullSimplify[Coefficient[u, u3]]
N3 = FullSimplify[Coefficient[v, v3]]
N4 = FullSimplify[Coefficient[u, u4]]
N4 = FullSimplify[Coefficient[v, v4]]

An alternate way is directly by realizing that $N_1$ has to satisfy that the condition that at node 1, $N_1=1$ , and is equal to zero on the lines $a-X_1=0$ and $b-X_2=0$ . The same can be applied for $N_2$ and $N_3$ . These result in the following expressions for the shape functions:

$\begin{split} N_1&=\frac{(b-X_2)(a-X_1)}{4ab}\\ N_2&=\frac{(b-X_2)(a+X_1)}{4ab}\\ N_3&=\frac{(b+X_2)(a+X_1)}{4ab}\\ N_4&=\frac{(b+X_2)(a-X_1)}{4ab} \end{split}$

The distribution of the shape functions on the element are illustrated in Figure 8.

Figure 8. Shape functions distribution in the bilinear quadrilateral element

The strain associated with the assumed displacement field has the following form:

$\begin{split} \varepsilon&=\left(\begin{array}{c}\varepsilon_{11}\\\varepsilon_{22}\\2\varepsilon_{12}\end{array}\right) =\left(\begin{array}{c}\frac{\partial u}{\partial X_1}\\\frac{\partial v}{\partial X_2}\\\frac{\partial u}{\partial X_2}+\frac{\partial v}{\partial X_1}\end{array}\right)\\ &=\left(\begin{array}{cccccccc}\frac{\partial N_1}{\partial X_1}&0&\frac{\partial N_2}{\partial X_1}&0&\frac{\partial N_3}{\partial X_1}&0&\frac{\partial N_4}{\partial X_1}&0\\ 0&\frac{\partial N_1}{\partial X_2}&0&\frac{\partial N_2}{\partial X_2}&0&\frac{\partial N_3}{\partial X_2}&0&\frac{\partial N_4}{\partial X_2}\\ \frac{\partial N_1}{\partial X_2} & \frac{\partial N_1}{\partial X_1} & \frac{\partial N_2}{\partial X_2} & \frac{\partial N_2}{\partial X_1} & \frac{\partial N_3}{\partial X_2} & \frac{\partial N_3}{\partial X_1} & \frac{\partial N_4}{\partial X_2} & \frac{\partial N_4}{\partial X_1} \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)\\ &=B u_e \end{split}$

Where, the explicit representation of the matrix $B$ is:

$B=\left(\begin{array}{cccccccc}-\frac{b-X_2}{4ab}&0&\frac{b-X_2}{4ab}&0&\frac{b+X_2}{4ab}&0&-\frac{b+X_2}{4ab}&0\\ 0&-\frac{a-X_1}{4ab}&0&-\frac{a+X_1}{4ab}&0&\frac{a+X_1}{4ab}&0&\frac{a-X_1}{4ab}\\ -\frac{a-X_1}{4ab}& -\frac{b-X_2}{4ab}& -\frac{a+X_1}{4ab}& \frac{b-X_2}{4ab} &\frac{a+X_1}{4ab}& \frac{b+X_2}{4ab} & \frac{a-X_1}{4ab}& -\frac{b+X_2}{4ab} \end{array}\right)$

Shear Locking

A quick glance at the matrix $B$ shows that the bilinear quadrilateral can model bending (linear strains) since the strains $\varepsilon_{11}$ and $\varepsilon_{22}$ contain linear expressions in $X_1$ and $X_2$ . However, it will be shown here that to model linear strains, the element predicts an associated shear strain $\gamma_{12}$ as well (termed parasitic shear strains). For simplicity, we will calculate the strains based on Equation 1. The strains have the form:

$\begin{split} \varepsilon_{11}&=a_1+a_3X_2\\ \varepsilon_{22}&=b_2+b_3X_1\\ \gamma_{12}&=a_2+a_3X_1+b_1+b_3X_2 \end{split}$

If we now assume pure bending state of strain with $a_3\neq 0$ . Then, the shear strain cannot be equal to zero since $a_3$ appears in the expression for $\gamma_{12}$ ! This behaviour is termed shear locking and causes the element to be too stiff under pure bending.

Stiffness Matrix

The stiffness matrix of the linear strain triangle can be evaluated using the following integral (assuming a constant thickness $t$ ):

$K^e=t\int_{-a}^{a}\int_{-b}^{b} \! B^T C B\,\mathrm{d}X_1\mathrm{d}X_2$

The stiffness matrix has the dimensions of $8 \times 8$ , and the following Mathematica code can be utilized to view its components:
View Mathematica Code

Shapefun=Table[0,{i,1,4}];
Shapefun[[1]]=(b-x2)(a-x1)/4/a/b;
Shapefun[[2]]=(b-x2)(a+x1)/4/a/b;
Shapefun[[3]]=(b+x2)(a+x1)/4/a/b;
Shapefun[[4]]=(b+x2)(a-x1)/4/a/b;
B=Table[0,{i,1,3},{j,1,8}];
Do[B[[1,2i-1]]=B[[3,2i]]=D[Shapefun[[i]],x1];B[[2,2i]]=B[[3,2i-1]]=D[Shapefun[[i]],x2],{i,1,4}];
Cc=Ee/(1+nu)*{{(1-nu)/(1-2nu),nu/(1-2nu),0},{(nu)/(1-2nu),(1-nu)/(1-2nu),0},{0,0,1/2}};
Kbeforeintegration=t*FullSimplify[Transpose[B].Cc.B];
K=Integrate[Kbeforeintegration,{x2,-b,b},{x1,-a,a}];
K//MatrixForm

Nodal Loads

The nodal loads due to distributed body forces and traction forces on the boundaries are equally distributed among the nodes. The same procedure followed for the triangular elements can be repeated to obtain the distribution shown in Figure 9. The following Mathematica code can be utilized for the calculations producing the distributions in Figure 9.
View Mathematica Code

Shapefun=Table[0,{i,1,4}];
Shapefun[[1]]=(b-x2)(a-x1)/4/a/b;
Shapefun[[2]]=(b-x2)(a+x1)/4/a/b;
Shapefun[[3]]=(b+x2)(a+x1)/4/a/b;
Shapefun[[4]]=(b+x2)(a-x1)/4/a/b;
Nn=Table[0,{i,1,2},{j,1,8}];
Do[Nn[[1,2i-1]]=Nn[[2,2i]]=Shapefun[[i]],{i,1,4}];
tn={t1,t2};
rb={rb1,rb2};
fetraction=Integrate[(Transpose[Nn].tn/.x1->-a),{x2,-b,b}]//MatrixForm
febodyforces=Integrate[(Transpose[Nn].rb),{x2,-b,b},{x1,-a,a}]//MatrixForm

Figure 9. Nodal forces in a bilinear quadrilateral element with a constant unit thickness due to (a) constant body forces vectors, (b) constant traction vector on one side.

Quadratic Quadrilateral Elements

By introducing four additional nodes in the mid-sides of a quadrilateral element, the displacement shape within the element can have a quadratic form, and the parasitic shear stiffness of the bilinear quadrilateral can be avoided. The quadratic quadrilateral offers such advantage with the price of higher computational time due to having additional degrees of freedom. Consider the plane quadrilateral shown in Figure 10. The trial displacement function has the following form:

$u_{2D}=\left(\begin{array}{c}u\\v\end{array}\right)=\left(\begin{array}{c}a_0+a_1X_1+a_2X_2+a_3X_1X_2+a_4X_1^2+a_5X_2^2+a_6X_1X_2^2+a_7X_1^2X_2\\b_0+b_1X_1+b_2X_2+b_3X_1X_2+b_4X_1^2+b_5X_2^2+b_6X_1X_2^2+a_7X_1^2X_2\end{array}\right)$

where $a_0$ , $a_1$ , $a_2$ , $a_3$ , $a_4$ , $a_5$ , $a_6$ , $a_7$ , $b_0$ , $b_1$ , $b_2$ , $b_3$ , $b_4$ , $b_5$ , $b_6$ and $b_7$ are 16 generalized degrees of freedom of the element. The displacement function in terms of the nodal degrees of freedom has the form:

$u_{2D}=\left(\begin{array}{c}u\\v\end{array}\right)=\left(\begin{array}{c}N_1u_1+N_2u_2+N_3u_3+N_4u_4+N_5u_5+N_6u_6+N_7u_7+N_8u_8\\N_1v_1+N_2v_2+N_3v_3+N_4v_4+N_5v_5+N_6v_6+N_7v_7+N_8v_8\end{array}\right)$

Figure 10. Quadratic quadrilateral element

Following the procedures shown in the previous sections, the shape functions can be shown to have the following forms:

$\begin{split} N_1&=-\frac{(b-X_2)(a-X_1)(1+\frac{X_1}{a}+\frac{X_2}{b})}{4ab}\\ N_2&=-\frac{(b-X_2)(a+X_1)(1-\frac{X_1}{a}+\frac{X_2}{b})}{4ab}\\ N_3&=-\frac{(b+X_2)(a+X_1)(1-\frac{X_1}{a}-\frac{X_2}{b})}{4ab}\\ N_4&=-\frac{(b+X_2)(a-X_1)(1+\frac{X_1}{a}-\frac{X_2}{b})}{4ab}\\ N_5&=\frac{(b-X_2)(a-X_1)(a+X_1)}{2a^2b}\\ N_6&=\frac{(a+X_1)(b-X_2)(b+X_2)}{2ab^2}\\ N_7&=\frac{(b+X_2)(a-X_1)(a+X_1)}{2a^2b}\\ N_8&=\frac{(a-X_1)(b-X_2)(b+X_2)}{2ab^2} \end{split}$

The distribution of the shape functions on the element are illustrated in Figure 11.

Figure 11. Shape functions distribution on the quadratic quadrilateral

Stiffness Matrix

The stiffness matrix of the linear strain triangle can be evaluated using the following integral (assuming a constant thickness $t$ ):

$K^e=t\int_{-a}^{a}\int_{-b}^{b} \! B^T C B\,\mathrm{d}X_1\mathrm{d}X_2$

The stiffness matrix has the dimensions of $16 \times 16$ , and the following Mathematica code can be utilized to view its components:
View Mathematica Code

Shapefun=Table[0,{i,1,8}];
Shapefun[[1]]=(b-x2)(a-x1)/4/a/b*-(1+x1/a+x2/b);
Shapefun[[2]]=(b-x2)(a+x1)/4/a/b*-(1-x1/a+x2/b);
Shapefun[[3]]=(b+x2)(a+x1)/4/a/b*-(1-x1/a-x2/b);
Shapefun[[4]]=(b+x2)(a-x1)/4/a/b*-(1+x1/a-x2/b);
Shapefun[[5]]=(b-x2)(a-x1)(a+x1)/2/a^2/b;
Shapefun[[6]]=(a+x1)(b-x2)(b+x2)/2/a/b^2;
Shapefun[[7]]=(b+x2)(a-x1)(a+x1)/2/a^2/b;
Shapefun[[8]]=(a-x1)(b-x2)(b+x2)/2/a/b^2;
B=Table[0,{i,1,3},{j,1,16}];
Do[B[[1,2i-1]]=B[[3,2i]]=D[Shapefun[[i]],x1];B[[2,2i]]=B[[3,2i-1]]=D[Shapefun[[i]],x2],{i,1,8}];
B//MatrixForm
Cc=Ee/(1+nu)*{{(1-nu)/(1-2nu),nu/(1-2nu),0},{(nu)/(1-2nu),(1-nu)/(1-2nu),0},{0,0,1/2}};
Kbeforeintegration=t*FullSimplify[Transpose[B].Cc.B];
K=Integrate[Kbeforeintegration,{x2,-b,b},{x1,-a,a}];
K//MatrixForm

Nodal Forces

Assuming that the distributed body forces vector per unit mass is $b={b_1,b_2}$ , the mass density is $\rho$ , and the traction vector per unit area on the left side is $t_n={t_1,t_2}$ , then, the lumped nodal forces due to the distributed body forces can be obtained as follows:

$f^e=\int_e \! N^T \rho b \, \mathrm{d}X=\frac{\rho t A}{3}\left(\begin{array}{c}-\frac{b_1}{4}\\-\frac{b_2}{4}\\-\frac{b_1}{4}\\-\frac{b_2}{4}\\-\frac{b_1}{4}\\-\frac{b_2}{4}\\-\frac{b_1}{4}\\-\frac{b_2}{4}\\b_1\\b_2\\b_1\\b_2\\b_1\\b_2\\b_1\\b_2\end{array}\right)$

The lumped nodal forces due to the distributed traction vector has the following form:

$f^e=\int_{\partial e} \! N^T t_n \, \mathrm{d}S=\int_{-b}^{b} \! N^T t_n|_{X_1=a} \, \mathrm{d}X_2=\frac{bt}{3}\left(\begin{array}{c}t_1\\t_2\\0\\0\\0\\0\\t_1\\t_2\\0\\0\\0\\0\\0\\0\\4t_1\\4t_2\end{array}\right)$

The nodal loads due to a constant distributed body forces vector and the traction vector on the left side are shown in Figure 12. Notice the surprising result of having small forces applied on the corner nodes opposite in direction to the applied body forces.

View Mathematica Code

Shapefun=Table[0,{i,1,8}];
Shapefun[[1]]=(b-x2)(a-x1)/4/a/b*-(1+x1/a+x2/b);
Shapefun[[2]]=(b-x2)(a+x1)/4/a/b*-(1-x1/a+x2/b);
Shapefun[[3]]=(b+x2)(a+x1)/4/a/b*-(1-x1/a-x2/b);
Shapefun[[4]]=(b+x2)(a-x1)/4/a/b*-(1+x1/a-x2/b);
Shapefun[[5]]=(b-x2)(a-x1)(a+x1)/2/a^2/b;
Shapefun[[6]]=(a+x1)(b-x2)(b+x2)/2/a/b^2;
Shapefun[[7]]=(b+x2)(a-x1)(a+x1)/2/a^2/b;
Shapefun[[8]]=(a-x1)(b-x2)(b+x2)/2/a/b^2;
Nn=Table[0,{i,1,2},{j,1,16}];
Do[Nn[[1,2i-1]]=Nn[[2,2i]]=Shapefun[[i]],{i,1,8}];
tn={t1,t2};
rb={rb1,rb2};
fetraction=Integrate[(Transpose[Nn].tn/.x1->-a),{x2,-b,b}]//MatrixForm
febodyforces=Integrate[(Transpose[Nn].rb),{x2,-b,b},{x1,-a,a}]//MatrixForm

Figure 12. Nodal forces in a quadratic quadrilateral element with a constant unit thickness due to (a) constant body forces vectors, (b) constant traction vector on one side. $l=2b$ and $A=4ab$ .

Extension to 3D

The 2D elements discussed previously can be extended in a straightforward manner to the three dimensional case. In three dimensions, the 4-node tetrahedron is the 3D version of the linear triangular element, while the 10-node tetrahedron is the nonlinear 3D version corresponding to the quadratic triangular element. The 8-node brick element (trilinear quadrilateral) is the 3D version of the bilinear quadrilateral, while the 20-node quadratic brick element is the nonlinear 3D version corresponding to the quadratic quadrilateral. The properties of the two-dimensional elements are inherited by their three-dimensional counterparts. In addition, a wedge elements can be defined in 3D for meshing irregular regions or when combining tetrahedrons with brick elements. See the following table for illustration of possible 3D elements that exist in the variety of available FEA software.

It should be mentioned that most FEA software have built-in meshing capabilities. For 3D solids, the automatic meshing algorithms that produce brick elements might fail in meshing irregular shapes. However, the tetrahedron meshing will most always succeed in meshing an irregular shape. In such cases, if accuracy is sought, either a fine mesh or the nonlinear tetrahedron should be utilized.

2D elements and their extension to 3D

General Requirements for an Element

Finding an appropriate displacement shape function or interpolation function has been a subject of extensive study since the development of the finite element analysis method. While new types of elements are being introduced on a regular basis, some basic requirements characterize a good approximation for the displacement. Among those requirements are isotropy, ability to model rigid body motions and constant strains, and element compatibility. An isotropic interpolation function is a function that would not favor a direction over the other. An interpolation function that has the form $u = a_1 + a_2X_1 + a_3X_2^2$ would clearly produce different results if the structure is rotated such that the coordinates $X_1$ and $X_2$ are switched. However, an interpolation function of the form $u=a_1 + a_2X_1 + a_3X_2$ is isotropic. Another basic requirement is that the displacement interpolation function should be able to model rigid body motion. An interpolation function that has the form $u = a_1X_1 + a_2X_2$ cannot model a constant displacement, and thus, cannot be acceptable. Such interpolation function, however, can model constant strain. Another major requirement for the interpolation functions is ensuring the compatibility between elements. All the elements presented in the previous section ensure that if two elements of the same type share the same nodes along a boundary, then the displacement along that boundary (the shared side) between the elements is continuous. This was ensured by choosing interpolation or shape functions that guarantee that the displacement on the boundaries (sides) of the element is completely determined by the nodes on that side, and thus, the elements that share that side along with the associated nodes share the same displacement across that side. The displacement interpolation functions presented so far are thus truly $C_0$ interpolation functions. In some applications, more (or less) stringent requirements can be considered, and the user should be aware of the capabilities or lack thereof of the elements being used.

Open Educational Resources

Finite Element Analysis: Two Dimensional Solid Elements

Triangular Elements

Linear Triangular Element (Constant Strain Triangle)

Stiffness Matrix

Nodal Forces due to Body Forces

Nodal Forces due to Traction Vector on One Side

Quadratic Triangular Elements

Stiffness Matrix

Nodal Forces due to Body Forces

Nodal Forces due to Traction Vector on One Side

Quadrilateral Elements

Linear Quadrilateral Elements

Shear Locking

Stiffness Matrix

Nodal Loads

Quadratic Quadrilateral Elements

Stiffness Matrix

Nodal Forces

Extension to 3D

General Requirements for an Element

Leave a Reply Cancel reply