Formulation of First-Order Triangular Finite Element

In this blog post, we will go through the detailed formulation of a first-order triangular finite element used in structural analysis. This type of element is commonly used in the finite element method (FEM) for solving engineering problems.

We will consider the following assumptions:

• Plane stress or Plane strain condition
• Linear elastic, isotropic, and homogeneous material.
• The element is subjected to loading only in the xy plane.
• Small displacements.

Introduction

A first-order triangular element has 3 nodes and 2 degrees of freedom (DOF) per node, leading to a total of 6 DOF. Therefore, we can already be sure that the stiffness matrix for this element will be a 6×6 matrix. This matrix will also be symmetric and singular.

{\left[k\right]}{_{6x6}}{.}{\left[\delta\right]}_{6x1}={\left[f\right]}_{6x1}

The degrees of freedom for an element used in structural analysis are displacements. For first-order triangular elements, there are two possible displacements for each node: one in the x direction and one in the y direction, as shown below:

The material properties required are only the modulus of elasticity (E) and Poisson’s ratio (ν) due to the assumptions related to the material.

Element Geometry and Shape Functions

For a first-order triangular element, the nodes are located at the vertices. The coordinates of the nodes in the xy-plane can be denoted as shown below:

Defining functions to calculate the displacement components at any position within the element is an essential step in formulating any finite element. Polynomials serve as the simplest and most practical functions for this purpose. To establish a polynomial, we must determine the number of coefficients it requires. The process for calculating the necessary number of coefficients for each polynomial is explained below.

To assist in constructing the polynomial, we will use Pascal’s Triangle.

Therefore, the displacement function u is as follows:

u\left(x,\right.y)={a}_0+{a}_1.x+{a}_2.y

That can be written using matrix notation as:

u\left(x,\right.y)=\left[
\begin{array}{llllllllll}
	1 & x & y
\end{array}\right].\left[
\begin{array}{llllllllll}
	{}a_0 \\
	{a}_1 \\
	{a}_2
\end{array}\right]

In a more compact manner:

u\left(x,\right.y)={\left[V\right].}\left\{a_u\right\}

And, the displacement function v is:

v\left(x,\right.y)={a}_3+{a}_4.x+{a}_5.y

That can be written using matrix notation as:

v\left(x,\right.y)=\left[
\begin{array}{llllllllll}
	1 & x & y
\end{array}\right].\left[
\begin{array}{llllllllll}
	{a}_3 \\
	{a}_4 \\
	{a}_5
\end{array}\right]

or in a more compact way:

v\left(x,\right.y)={\left[V\right].}\left\{a_v\right\}

To calculate the coefficients a_{i} , we use the nodal displacements as boundary conditions, which generates a system of 6 algebraic equations. Since there are 6 coefficients to determine, this system provides a unique solution. The equations appear below in matrix notation.

\left[
\begin{array}{llllllllll}
	{u}_1 \\
	{u}_2 \\
	{u}_3
\end{array}\right]=\left[
\begin{array}{llllllllll}
	1 & {x}_1 & {y}_1 \\
	1 & {x}_2 & {y}_2 \\
	1 & {x}_3 & {y}_3
\end{array}\right].\left[
\begin{array}{llllllllll}
	{a}_0 \\
	{a}_1 \\
	{a}_2
\end{array}\right]

\left[
\begin{array}{llllllllll}
	v_1   \\
	v_2   \\
	{v}_3
\end{array}\right]=\left[
\begin{array}{llllllllll}
	1 & {x}_1 & {y}_1 \\
	1 & {x}_2 & {y}_2 \\
	1 & {x}_3 & {y}_3
\end{array}\right].\left[
\begin{array}{llllllllll}
	{a}_3 \\
	{a}_4 \\
	{a}_5
\end{array}\right]

The matrix equations above can be expressed more concisely as follows:

\left\{u\right\}={\left[P\right].}\left\{a_u\right\}

\left\{v\right\}={\left[P\right].}\left\{a_v\right\}

To calculate the coefficients  a_{i} , the column matrix {a} should be isolated in the equations above. Thus:

\left\{{a}_u\right\}={\left[P\right]}^{-1}.\left\{u\right\}

\left\{{a}_v\right\}={\left[P\right]}^{-1}.\left\{v\right\}

Thus, the displacement functions can be calculated as follows:

u\left(x,y\right)={\left[V\right].\left[P\right]}^{-1}.\left\{u\right\}=\left[N\right].\left\{u\right\}

v\left(x,y\right)={\left[V\right].\left[P\right]}^{-1}.\left\{v\right\}=\left[N\right].\left\{v\right\}

The matrix operation above can be performed yielding the following expressions for the displacement functions:

u\left(x,y\right)=\left[N\right].\left\{u\right\}=\left[
\begin{array}{llllllllll}
	{N}_1 & {N}_2 & {N}_3
\end{array}\right].\left[
\begin{array}{llllllllll}
	{u}_1 \\
	{u}_2 \\
	{u}_3
\end{array}\right]

u\left(x,y\right)={{{N}}}_1.{u}_1+{N}_2.{u}_2+{N}_3.{u}_3

v\left(x,y\right)=\left[N\right].\left\{v\right\}=\left[
\begin{array}{llllllllll}
	{N}_1 & {N}_2 & {N}_3
\end{array}\right].\left[
\begin{array}{llllllllll}
	{v}_1 \\
	{v}_2 \\
	{v}_3
\end{array}\right]

v\left(x,y\right)={{{N}}}_1.{v}_1+{N}_2.{v}_2+{N}_3.{v}_3

Where  N_{i} are the interpolation functions or shape functions. Notice that:

\left[N\right]=\left[V\right].{\left[P\right]}^{-1}

The matrix [V] includes the variables x and y, while the matrix [P] represents the node coordinates. Consequently, its inverse also consists of the nodal coordinates. Thus, the nodal coordinates (defined after generating the mesh) and the variables x and y determine the shape functions.

The matrix operation above was performed, yielding the following result:

N_1(x,\;y)=\frac1{2A}\lbrack(x_2y_3-x_3y_2)+(y_2-y_3)x+(x_3-x_2)y\rbrack=\frac1{2A}(\alpha_1+\beta_1x+\gamma_1y)

N_2(x,\;y)=\frac1{2A}\lbrack(x_3y_1-x_1y_3)+(y_3-y_1)x+(x_1-x_3)y\rbrack=\frac1{2A}(\alpha_2+\beta_2x+\gamma_2y)

N_3(x,\;y)=\frac1{2A}\lbrack(x_1y_2-x_2y_1)+(y_1-y_2)x+(x_2-x_1)y\rbrack=\frac1{2A}(\alpha_3+\beta_3x+\gamma_3y)

Strain-Displacement Matrix (B-matrix)

Under plane stress or plane strain conditions, we calculate the three strain components as follows:

 \varepsilon _{x}= \frac{\partial u}{\partial x}

 \varepsilon _{y}= \frac{\partial v}{\partial y}

\gamma _{xy}= \frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}

These expressions for strain calculation are valid only under the assumption of small deformations. For large deformations, the Green strain tensor should be used. We have a blog article on this topic—click the following link to read it: Green Deformation Calculation for Large Deformations.

The previously derived equations for the displacement functions can now be applied to calculate these derivatives. It is important to note that, since the variables x and y appear exclusively in the shape functions N_{i} , only these functions require differentiation.

\varepsilon _{x}= \frac{\partial u}{\partial x}= \frac{\partial N_{1}}{\partial x}\cdot u_{1}+\frac{\partial N_{2}}{\partial x}\cdot u_{2}+\frac{\partial N_{3}}{\partial x}\cdot u_{3}

\varepsilon _{y}= \frac{\partial v}{\partial y}= \frac{\partial N_{1}}{\partial y}\cdot v_{1}+\frac{\partial N_{2}}{\partial y}\cdot v_{2}+\frac{\partial N_{3}}{\partial y}\cdot v_{3}

\gamma _{xy}= \frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}= \frac{\partial N_{1}}{\partial y}\cdot u_{1}+\frac{\partial N_{2}}{\partial y}\cdot u_{2}+\frac{\partial N_{3}}{\partial y}\cdot u_{3}+\frac{\partial N_{1}}{\partial x}\cdot v_{1}+\frac{\partial N_{2}}{\partial x}\cdot v_{2}+\frac{\partial N_{3}}{\partial x}\cdot v_{3}

We can write the above expressions in matrix form as:

\left\{\begin{array}{l}\varepsilon_x\\\varepsilon_y\\\gamma_{xy}\end{array}\right\}=\left[B\right]\cdot\begin{pmatrix}u_1\\v_1\\u_2\\v_2\\u_3\\v_3\end{pmatrix}=\left[B\right]\cdot\left\{\delta\right\}

We construct the B-matrix using the derivatives of the shape functions with respect to x and y:

B=\begin{bmatrix}\frac{\partial N_1}{\partial x}&0&\frac{\partial N_2}{\partial x}&0&\frac{\partial N_3}{\partial x}&0\\0&\frac{\partial N_1}{\partial y}&0&\frac{\partial N_2}{\partial y}&0&\frac{\partial N_3}{\partial y}\\\frac{\partial N_1}{\partial y}&\frac{\partial N_1}{\partial x}&\frac{\partial N_2}{\partial y}&\frac{\partial N_2}{\partial x}&\frac{\partial N_3}{\partial y}&\frac{\partial N_3}{\partial x}\end{bmatrix}

Important Conclusion about the Strain-Displacement Matrix (B-matrix)

We will now calculate the strain components using a different approach to demonstrate how they vary across the element’s area.

Using the displacement functions in their initial form, we have the following expressions:

u(x,y)=a_{0}+a_{1}\cdot x+a_{2}\cdot y

v(x,y)=a_{3}+a_{4}\cdot x+a_{5}\cdot y

We can calculate the strain components by differentiating these functions, as shown below:

\varepsilon_x=\frac{\partial u}{\partial x}=a_1

\varepsilon_y=\frac{\partial v}{\partial y}=a_5

\gamma_{xy}=\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}=a_2+a_4

These equations lead to an important conclusion: the strain in a First-Order Triangular Element is constant (this is why this element is also named as the Constant Strain Triangle). In other words, strain does not change across the element. This is not a good structural behavior.

Since strain is constant, all stress components are also constant throughout the element. This is a problem because if you use this element in a region where stress or strain changes a lot, you will need a very fine mesh to capture these variations. This increases the computational cost.

A better option is to use a Second-Order Triangular Element. In this type of element, stress and strain vary across the area, so fewer elements are needed in regions with high stress or strain gradients. This reduces computational cost.

We also have a blog post explaining the formulation of the Second-Order Triangular Element. Click the link below to read it:

Formulation of Second-Order Triangular Finite Element

We could also have differentiated the interpolation functions obtained previously to reach to the same conclusion, as follows:

\frac{\partial N_1}{\partial x}=\frac1{2A}\cdot\beta_1

\frac{\partial N_1}{\partial y}=\frac1{2A}\cdot\gamma_1

\frac{\partial N_2}{\partial x}=\frac1{2A}\cdot\beta_2

\frac{\partial N_2}{\partial y}=\frac1{2A}\cdot\gamma_2

\frac{\partial N_3}{\partial x}=\frac1{2A}\cdot\beta_3

\frac{\partial N_3}{\partial y}=\frac1{2A}\cdot\gamma_3

Where:

\beta_1=y_2-y_3

\gamma_1=x_3-x_2

\beta_2=y_3-y_1

\gamma_2=x_1-x_3

\beta_3=y_1-y_2

\gamma_3=x_2-x_1

The strain components are calculated as shown below.

\varepsilon_x=\frac1{2A}\cdot\beta_1\cdot u_1+\frac1{2A}\cdot\beta_2\cdot u_2+\frac1{2A}\cdot\beta_3\cdot u_3

\varepsilon_y=\frac1{2A}\cdot\gamma_1\cdot v_1+\frac1{2A}\cdot\gamma_2\cdot v_2+\frac1{2A}\cdot\gamma_3\cdot v_3

\gamma_{xy}=\frac1{2A}\cdot\beta_1\cdot v_1+\frac1{2A}\cdot\beta_2\cdot v_2+\frac1{2A}\cdot\beta_3\cdot v_3+\frac1{2A}\cdot\gamma_1\cdot u_1+\frac1{2A}\cdot\gamma_2\cdot u_2+\frac1{2A}\cdot\gamma_3\cdot u_3

\begin{Bmatrix}\varepsilon_x\\\varepsilon_y\\\gamma_{xy}\end{Bmatrix}=\begin{bmatrix}\frac{\beta_1}{2A}&0&\frac{\beta_2}{2A}&0&\frac{\beta_3}{2A}&0\\0&\frac{\gamma_1}{2A}&0&\frac{\gamma_2}{2A}&0&\frac{\gamma_3}{2A}\\\frac{\gamma_1}{2A}&\frac{\beta_1}{2A}&\frac{\gamma_2}{2A}&\frac{\beta_2}{2A}&\frac{\gamma_3}{2A}&\frac{\beta_3}{2A}\end{bmatrix}\begin{Bmatrix}u_1\\v_1\\u_2\\v_2\\u_3\\v_3\end{Bmatrix}

You can see that the derivatives of the shape functions do not depend on x or y. The only terms in this calculation are geometric parameters, such as the element area and node coordinates. These parameters are defined once the mesh is created, before solving the mathematical model, and therefore remain constant.

The Strain-Displacement Matrix is multiplied by the nodal displacement vector to determine the strain components. Keep in mind that once the mathematical model is solved, the nodal displacements are computed and remain unchanged.

As a result, the strain does not vary across the element’s area, confirming our previous conclusion.

Constitutive Matrix (D-matrix)

We calculate the stress components of the First-Order Triangular Element using the following expression:

\begin{bmatrix}\sigma_x\\\sigma_y\\\tau_{xy}\end{bmatrix}=\lbrack D\rbrack\cdot\begin{bmatrix}\begin{array}{c}\varepsilon_x\\\varepsilon_y\\\gamma_{xy}\end{array}\end{bmatrix}

For plane stress conditions, the constitutive matrix D is defined as:

D=\frac{E}{1-\nu^2}\begin{bmatrix}1 & \nu & 0 \\ \nu & 1 & 0 \\ 0 & 0 & \frac{1-\nu}{2}\end{bmatrix}

And for plane strain condition, the constitutive matrix D is defined as:

D=\frac{E}{(1+\nu)(1-2\nu)}\begin{bmatrix}1-\nu & \nu & 0 \\ \nu & 1-\nu & 0 \\ 0 & 0 & \frac{1-2\nu}{2}\end{bmatrix}

Let’s take the plane stress condition as an example to derive the expressions:

\begin{Bmatrix}\sigma_x\\\sigma_y\\\tau_{xy}\end{Bmatrix}=\frac E{1-\nu^2}\begin{bmatrix}1&\nu&0\\\nu&1&0\\0&0&\frac{1-\nu}2\end{bmatrix}\cdot\begin{Bmatrix}\varepsilon_x\\\varepsilon_y\\\gamma_{xy}\end{Bmatrix}

We can substitute the previously calculated strain values into the strain vector. To simplify the equations, we will use the “a” coefficients, but you can also substitute the full strain expressions derived earlier.

\begin{Bmatrix}\sigma_x\\\sigma_y\\\tau_{xy}\end{Bmatrix}=\frac E{1-\nu^2}\begin{bmatrix}1&\nu&0\\\nu&1&0\\0&0&\frac{1-\nu}2\end{bmatrix}\cdot\begin{Bmatrix}a_1\\a_5\\a_2+a_4\end{Bmatrix}

Performing the matrix multiplication:

\sigma_x=\frac E{1-v^2}\left(a_1+v\cdot a_5\right)

\sigma_y=\frac E{1-v^2}\left(a_5+v\cdot a_1\right)

\tau_{xy}=\frac E{1-v^2}\cdot\left(\frac{1-v}2\right)\cdot\left(a_2+a_4\right)

Notice that none of the stress components depend on either x or y. This means that stress does not change across the element’s area, which is expected since the strain components are also constant.

Stiffness Matrix [k]

With the Strain-Displacement Matrix (B-matrix) and the Constitutive Matrix (D-matrix) defined, we can now calculate the stiffness matrix of the First-Order Triangular Element. To start, we first determine the total strain energy of the element.

Total Strain Energy

Calculating the total strain energy is important because it enables us to evaluate the element stiffness matrix using one of these two methods:

1. Castigliano’s First Theorem
2. Principle of Minimum Total Potential Energy

Both of these methods require the total strain energy of the element, which can be calculated using the following expression:

U_e=\frac12\cdot\iiint_V{\{\varepsilon\}}^T\cdot\lbrack D\rbrack\cdot\{\varepsilon\}\cdot\operatorname dV=\frac12\cdot\iiint_V{\{\delta\}}^T\cdot{\lbrack B\rbrack}^T\cdot\lbrack D\rbrack\cdot\lbrack B\rbrack\cdot\{\delta\}\operatorname dV

Total Potential Energy

The Principle of Minimum Total Potential Energy will be used to evaluate the stiffness matrix K. Therefore, we first need to calculate the Total Potential Energy of the element, which can be done using the following expression:

\Pi=U_e-W

Here, W represents the work associated with all possible external nodal forces. These forces are illustrated in the image below:

Therefore, W can be calculated as:

W=\left\{\delta\right\}^T\cdot\left\{f\right\}

where:

\left\{\delta\right\}^T=\begin{bmatrix}u_1&v_1&u_2&v_2&u_3&v_3\end{bmatrix}

\left\{f\right\}=\begin{bmatrix}f_{1x}\\f_{1y}\\f_{2x}\\f_{2y}\\f_{3x}\\f_{3y}\end{bmatrix}

We calculate the Total Potential Energy of the element using the following expression:

\Pi=U_e-W=\frac12\cdot\iiint_V\left\{\delta\right\}^T\cdot\left[B\right]^T\cdot\left[D\right]\cdot\left[B\right]\cdot\left\{\delta\right\}\operatorname dV-\left\{\delta\right\}^T\cdot\left\{f\right\}

Principle of Minimum Total Potential Energy

At stable equilibrium, the Total Potential Energy of a body reaches its minimum. To find this minimum, we differentiate the Total Potential Energy function and equate it to zero.

We need to differentiate the Total Potential Energy function with respect to all nodal displacements and set these derivatives to zero, as shown below.

\frac{\partial\Pi}{\partial\delta_i}=0\;for\;i=1,...,6

These derivatives will yield 6 algebraic equations, which can be represented in matrix notation as follows:

{\left[k\right]}_{6x6}{\left[\delta\right]}_{6x1}={\left[f\right]}_{6x1}

Where [k] is the stiffness matrix of the Triangular Element.

This set of 6 algebraic equations represents the equilibrium condition at each node in both the x and y directions. For example, the first equation corresponds to the static equilibrium of node 1 in the x direction, while the second equation represents the static equilibrium of node 1 in the y direction, and so forth.

We calculate the stiffness matrix using the following expression, derived from the 6 derivatives.

\left[k\right]=\iiint_v\left[B\right]^T\cdot\left[D\right]\cdot\left[B\right]\cdot\operatorname dV

The stiffness matrix is a 6×6, singular and symmetric matrix, as expected.

Notice that this integral is taken over the volume of the element. Since neither matrix [B] nor [D] depends on x, y, or z, they can be removed from the integral. The result of the integral is simply the volume of the element. Therefore, the stiffness matrix can be computed using the following expression:

\left[k\right]=V^e\cdot\left[B\right]^T\cdot\left[D\right]\cdot\left[B\right]

The volume of the element is obtained by multiplying its area A by the plate’s thickness t. For this calculation, we will use the plane stress matrix D, resulting in the following expression:

\lbrack k\rbrack=\frac{Et}{4A(1-\nu^2)}\begin{bmatrix}\beta_1&0&\gamma_1\\\beta_2&0&\gamma_2\\\beta_3&0&\gamma_3\\0&\gamma_1&\beta_1\\0&\gamma_2&\beta_2\\0&\gamma_3&\beta_3\end{bmatrix}\begin{bmatrix}1&\nu&0\\\nu&1&0\\0&0&\frac{1-\nu}2\end{bmatrix}\begin{bmatrix}\beta_1&\beta_2&\beta_3&0&0&0\\0&0&0&\gamma_1&\gamma_2&\gamma_3\\\gamma_1&\gamma_2&\gamma_3&\beta_1&\beta_2&\beta_3\end{bmatrix}

Conclusion

This blog post details the formulation of a first-order triangular finite element used in structural analysis. The element, used in the finite element method (FEM), assumes plane stress or strain conditions, linear elastic, isotropic, homogeneous material, xy-plane loading, and small displacements.

Key Points

• Element Characteristics:
  • 3 nodes, 2 degrees of freedom (DOF) per node, totaling 6 DOF.
  • Stiffness matrix: 6x6, symmetric, and singular.
• Displacement Functions:
  • Defined for x and y directions using polynomials.
  • Functions expressed in matrix form:
    • u(x,y)=\left[V\right]\cdot a_u
    • v(x,y)=\left[V\right]\cdot a_v
• Strain-Displacement Matrix (B-matrix):
  • Strain components derived from displacement functions:
    • \varepsilon_x,\varepsilon_y,\gamma_{xy}
  • Expressed in matrix form: \varepsilon=\left[B\right]\cdot{\delta}
• Strain Variation:
  • Strain does not vary across the element’s area
• Constitutive Matrix (D-matrix):
  • For plane stress and plane strain conditions.
  • Stress components calculated using: \sigma=\left[D\right]\cdot{\varepsilon}
• Stiffness Matrix ([k]):
  • Calculated via Principle of Total Potential Energy.
  • Stiffness expression: \left[k\right]=\iiint_V\left[B\right]^T\left[D\right]\left[B\right]\operatorname dV