Formulation of Second-Order Triangular Finite Element

In this blog post, we will go through the detailed formulation of a second-order triangular finite element used in static structural analysis. This type of element is commonly used in 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 second-order triangular element has 6 nodes and 2 degrees of freedom (DOF) per node, leading to a total of 12 DOF. Therefore, we can already be sure that the stiffness matrix for this element will be a 12×12 matrix. This matrix will also be symmetric and singular.

{\left[k\right]}_{12x12}\;{\left[\delta\right]}_{12x1}\;={\left[f\right]}_{12x1}

The degrees of freedom for an element used in structural analysis are displacements. For second-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 second-order triangular element, the nodes are located at the vertices and midpoints of the sides. The coordinates of the nodes in the xy-plane can be denoted as shown below:

One of the initial steps in formulating any finite element is to define functions that calculate the displacement components at any position within the element. The simplest type of function to work with is a polynomial. To define a polynomial, it is essential to determine the number of coefficients it requires. The method for determining the number of coefficients needed for each polynomial is shown below.

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

Therefore, the displacement function u is as follows:

u(x,\;y)=a_0+a_1\cdot x+a_2\cdot y+a_3\cdot x\cdot y+a_4\cdot x^2+a_5\cdot y^2

That can be written using matrix notation as:

u(x,y)=\begin{bmatrix}1&x&y&x\cdot y&x^2&y^2\end{bmatrix}\cdot\begin{bmatrix}a_0\\a_1\\a_2\\a_3\\a_4\\a_5\end{bmatrix}

or in a more compact manner:

u(x,y)=\lbrack V\rbrack\cdot\{a_u\}

And, the displacement function v is as follows:

v(x,\;y)=a_6+a_7\cdot x+a_8\cdot y+a_9\cdot x\cdot y+a_{10}\cdot x^2+a_{11}\cdot y^2

That can be written using matrix notation as:

v(x,y)=\begin{bmatrix}1&x&y&x\cdot y&x^2&y^2\end{bmatrix}\cdot\begin{bmatrix}a_6\\a_7\\a_8\\a_9\\a_{10}\\a_{11}\end{bmatrix}

or in a more compact manner:

v(x,y)=\lbrack V\rbrack\cdot\{a_v\}

Note that the equations above include second-order (quadratic) terms. This is why this element is referred to as a Second-Order Triangular Finite Element.

To determine the coefficients a_i , the nodal displacements can be used as boundary conditions, leading to the formation of 12 algebraic equations. Since there are 12 coefficients to calculate, this system will have a unique solution. The equations are represented below using matrix notation.

\begin{bmatrix}u_1\\u_2\\u_3\\u_4\\u_5\\u_6\end{bmatrix}=\begin{bmatrix}1&x_1&y_1&x_1\cdot y_1&x_1^2&y_1^2\\1&x_2&y_2&x_2\cdot y_2&x_2^2&y_2^2\\1&x_3&y_3&x_3\cdot y_3&x_3^2&y_3^2\\1&x_4&y_4&x_4\cdot y_4&x_4^2&y_4^2\\1&x_5&y_5&x_5\cdot y_5&x_5^2&y_5^2\\1&x_6&y_6&x_6\cdot y_6&x_6^2&y_6^2\end{bmatrix}\cdot\begin{bmatrix}a_0\\a_1\\a_2\\a_3\\a_4\\a_5\end{bmatrix}

\begin{bmatrix}v_1\\v_2\\v_3\\v_4\\v_5\\v_6\end{bmatrix}=\begin{bmatrix}1&x_1&y_1&x_1\cdot y_1&x_1^2&y_1^2\\1&x_2&y_2&x_2\cdot y_2&x_2^2&y_2^2\\1&x_3&y_3&x_3\cdot y_3&x_3^2&y_3^2\\1&x_4&y_4&x_4\cdot y_4&x_4^2&y_4^2\\1&x_5&y_5&x_5\cdot y_5&x_5^2&y_5^2\\1&x_6&y_6&x_6\cdot y_6&x_6^2&y_6^2\end{bmatrix}\cdot\begin{bmatrix}a_6\\a_7\\a_8\\a_9\\a_{10}\\a_{11}\end{bmatrix}

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

\{u\}=\lbrack P\rbrack\cdot\{a_u\}

\{v\}=\lbrack P\rbrack\cdot\{a_v\}

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

\{a_u\}={\lbrack P\rbrack}^{-1}\cdot\{u\}

\{a_v\}={\lbrack P\rbrack}^{-1}\cdot\{v\}

Thus, the displacement functions can be calculated as follows:

u(x,y)=\lbrack V\rbrack\cdot{\lbrack P\rbrack}^{-1}\cdot\{u\}=\lbrack N\rbrack\cdot\{u\}

v(x,y)=\lbrack V\rbrack\cdot{\lbrack P\rbrack}^{-1}\cdot\{v\}=\lbrack N\rbrack\cdot\{v\}

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

u(x,y)=\lbrack N\rbrack\cdot\{u\}=\begin{bmatrix}N_1&N_2&N_3&N_4&N_5&N_6\end{bmatrix}\cdot\begin{bmatrix}u_1\\u_2\\u_3\\u_4\\u_5\\u_6\end{bmatrix}

u(x,y)=N_1\cdot u_1+N_2\cdot u_2+N_3\cdot u_3+N_4\cdot u_4+N_5\cdot u_5+N_6\cdot u_6

v(x,y)=\lbrack N\rbrack\cdot\{v\}=\begin{bmatrix}N_1&N_2&N_3&N_4&N_5&N_6\end{bmatrix}\cdot\begin{bmatrix}v_1\\v_2\\v_3\\v_4\\v_5\\v_6\end{bmatrix}

v(x,y)=N_1\cdot v_1+N_2\cdot v_2+N_3\cdot v_3+N_4\cdot v_4+N_5\cdot v_5+N_6\cdot v_6

Where N_i are the interpolation functions or shape functions. Notice that:

\lbrack N\rbrack=\lbrack V\rbrack\cdot{\lbrack P\rbrack}^{-1}

Where matrix [V] contains the variables x and y, matrix [P] contains the node coordinates, and therefore its inverse is also composed of node coordinates. Consequently, the shape functions will be defined by the nodal coordinates (established once the mesh is created) and the variables x and y.

N_i(x,y)

Strain-Displacement Matrix (B-matrix)

For plane stress or plane strain conditions, the three components of strain can be calculated as follows:

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

\varepsilon_y=\frac{\partial v}{\partial x}

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

The previously developed equations for the displacement functions can now be used to calculate these derivatives. Note that since the variables x and y only appear in the shape functions N_i , only these functions need to be differentiated.

\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+\frac{\partial N_4}{\partial x}\cdot u_4+\frac{\partial N_5}{\partial x}\cdot u_5+\frac{\partial N_6}{\partial x}\cdot u_6

\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+\frac{\partial N_4}{\partial y}\cdot v_4+\frac{\partial N_5}{\partial y}\cdot v_5+\frac{\partial N_6}{\partial y}\cdot v_6

\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_4}{\partial y}\cdot u_4+\frac{\partial N_5}{\partial y}\cdot u_5+\frac{\partial N_6}{\partial y}\cdot u_6+...

...+\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+\frac{\partial N_4}{\partial x}\cdot v_4+\frac{\partial N_5}{\partial x}\cdot v_5+\frac{\partial N_6}{\partial x}\cdot v_6

The above expressions can be written in matrix form as:

\begin{pmatrix}\varepsilon_x\\\varepsilon_y\\\gamma_{xy}\end{pmatrix}=\lbrack B\rbrack\cdot\begin{pmatrix}u_1\\v_1\\\begin{array}{c}u_2\\v_2\\\vdots\\\begin{array}{c}u_6\\v_6\end{array}\end{array}\end{pmatrix}=\lbrack B\rbrack\cdot\{\delta\}

The B-matrix, which is constructed using the derivatives of the shape functions with respect to x and y, is:

B=\begin{pmatrix}\frac{\partial N_1}{\partial x}&0&\frac{\partial N_2}{\partial x}&0&\cdots&\frac{\partial N_6}{\partial x}&0\\0&\frac{\partial N_1}{\partial y}&0&\frac{\partial N_2}{\partial y}&\cdots&0&\frac{\partial N_6}{\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}&\cdots&\frac{\partial N_6}{\partial y}&\frac{\partial N_6}{\partial x}\end{pmatrix}

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

Understanding how the components of strain vary over the area of an element is crucial for any finite element formulation.

For example, the linear triangular element, also known as the Constant Strain Triangle (CST), has constant strain and stress throughout the entire area of the element. This characteristic is undesirable in regions with high gradients of stress and strain, as it necessitates the use of many elements, significantly increasing computational costs.

In our current formulation, we have the Second-Order Triangular Finite Element. We have developed equations to calculate the strain components as derivatives of the shape functions. However, the shape functions were not explicitly derived; only the expressions used to derive it were shown. Consequently, it is not possible to conclude how strain will vary.

Therefore, we will now calculate the strain components differently to illustrate how they vary over 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+a_3\cdot x\cdot y+a_4\cdot x^2+a_5\cdot y^2

v(x,\;y)=a_6+a_7\cdot x+a_8\cdot y+a_9\cdot x\cdot y+a_{10}\cdot x^2+a_{11}\cdot y^2

By differentiating these functions, we can calculate the strain components as shown below:

\varepsilon_x=\frac{\partial u}{\partial x}=a_1+a_3\cdot y+2\cdot a_4\cdot x

\varepsilon_y=\frac{\partial v}{\partial y}=a_8+a_9\cdot x+2\cdot a_{11}\cdot y

\gamma_{xy}=\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}=a_2+a_3\cdot x+2\cdot a_5\cdot y+a_7+a_9\cdot y+2\cdot a_{10}\cdot x

The critical conclusion from these expressions is that the strain components for the Second-Order Triangular Element vary linearly in the x and y directions. This is a significant improvement compared to the Constant Strain Triangle.

This linear variation allows for the use of fewer elements in regions with high gradients of stress and strain, thereby reducing computational costs.

Constitutive Matrix (D-matrix)

The stress components of the Second-Order Triangular Element can be calculated using the following expression:

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

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

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

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

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

It should be noted that the only difference in the formulation of the Second-Order Triangular Element between plane stress and plane strain conditions is in the calculation of the constitutive matrix D.

Another important point to note is that since stress is calculated by multiplying the constitutive matrix D by the strain vector, and the strain components vary linearly in the x and y directions, the stress components will consequently also vary linearly in the x and y directions.

Stiffness Matrix [k]

Now that we have the Strain-Displacement Matrix (B-matrix) and the Constitutive Matrix (D-matrix) defined, we can evaluate the stiffness matrix of the Second-Order Triangular Element. To begin, we need to calculate 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 nodal forces within the element. These forces are illustrated in the image below:

Therefore, W can be calculated as:

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

where

{\{\delta\}}^T=\left[\begin{array}{cccc}u_1&v_1&u_2&v_2\end{array}\;\cdots\;\begin{array}{cc}u_6&v_6\end{array}\right]

\{f\}=\begin{bmatrix}f_{1x}\\f_{1y}\\f_{2x}\\f_{2y}\\\vdots\\f_{6x}\\f_{6y}\end{bmatrix}

Thus, the Total Potential Energy of the element can be calculated as follows:

\Pi=U_e-W=\frac12\cdot\iiint_V{\{\delta\}}^T\cdot{\lbrack B\rbrack}^T\cdot\lbrack D\rbrack\cdot\lbrack B\rbrack\cdot\{\delta\}\operatorname dV-{\{\delta\}}^T\cdot\{f\}

Principle of Minimum Total Potential Energy

The Total Potential Energy of a body is minimized at stable equilibrium. To find this minimum, we differentiate the Total Potential Energy function and set it equal to zero.

For the Second-Order Triangular Element, 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,...,12

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

{\left[k\right]}_{12x12}\;{\left[\delta\right]}_{12x1}\;={\left[f\right]}_{12x1}

Where \lbrack k\rbrack is the stiffness matrix of the Second-Order Triangular Element.

It is important to note that this set of 12 algebraic equations represents the equilibrium condition at each node in both the x and y directions. For instance, the first equation corresponds to the static equilibrium of node 1 in the x direction, while the second equation pertains to the static equilibrium of node 1 in the y direction, and so on.

It can be demonstrated that the stiffness matrix is calculated using the following expression, which is derived from those 12 derivatives.

\lbrack k\rbrack=\iiint_V{\lbrack B\rbrack}^T\cdot\lbrack D\rbrack\cdot\lbrack B\rbrack\cdot\operatorname dV

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

To evaluate this integral, numerical integration is typically used, specifically the Gaussian Quadrature Procedure, which is applied over the element’s area.

Conclusion

This blog post details the formulation of a second-order triangular finite element used in static 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:
  • 6 nodes, 2 degrees of freedom (DOF) per node, totaling 12 DOF.
  • Stiffness matrix: 12x12 , symmetric, and singular.
• Displacement Functions:
  • Defined for x and y directions using polynomials.
  • Functions expressed in matrix form:
    • u(x,y) = [V] \cdot {a_u}
    • v(x,y) = [V] \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} = [B] \cdot {\delta}
• Strain Variation:
  • Second-order elements offer linear variation in strain, improving computational efficiency compared to constant strain triangles.
• Constitutive Matrix (D-matrix):
  • For plane stress and plane strain conditions.
  • Stress components calculated using: {\sigma} = [D] \cdot {\varepsilon}
• Stiffness Matrix ([k]):
  • Calculated via total strain energy or total potential energy.
  • Integral expression: [k] = \iiint_V [B]^T [D] [B] dV
  • Numerical integration, often Gaussian Quadrature, is used.

Isoparametric Formulation

The isoparametric formulation can also be applied to this element. The advantage of this approach is that it uses the same shape functions to interpolate both the geometry of the element and the displacements.

Since the shape functions of the Second-Order Triangular Elements are quadratic functions, using the isoparametric formulation allows the triangle’s sides to be curved, conforming to a second-order polynomial.

This characteristic is particularly useful for representing curved parts of the domain, such as holes. Fewer elements are needed to accurately represent the geometry, thereby decreasing computational cost.

The image below illustrates the difference between an element formulated using the isoparametric formulation and one formulated in the manner described here.

The characteristics of the element remain the same regardless of how it has been formulated. Therefore, the conclusions reached in this text are also valid for the isoparametric formulation.

Learn More!

Check out my Finite Element Method Course by clicking here.

Enhance your understanding of FEA with the following blog posts:

• Stress-Strain Curve Approximation: Ramberg-Osgood Relationship
• Axial Connector in ABAQUS: A Comprehensive Guide
• Bolt Pre-tension Techniques in Abaqus
• Abaqus Surface-Based Coupling Constraint