Samer Adeeb

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) \]


Linear Triangular Element (Constant Strain Triangle)

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.


Shape functions in linear triangular elements

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 5.  Quadratic triangle

Figure 4. Quadratic 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 + 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

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

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

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\\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.

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

Figure 7. Bilinear quadrilateral

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

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

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

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

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

Figure 11. Shape functions distribution on the quadratic quadrilateral

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

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

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.

Page Comments

  1. Sandeep says:

    These resources are very helpful . Thank you .

    1. Samer Adeeb says:

      You are very welcome

Leave a Reply

Your email address will not be published. Required fields are marked *