Open Educational Resources

Two Dimensional Solid Elements: Triangular 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} = \begin{pmatrix} u \\ v \\ \end{pmatrix} = \begin{pmatrix} a_0 + a_1X_1 + a_2X_2 \\ b_0 + b_1X_1 + b_2X_2 \\ \end{pmatrix}\]

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} = \begin{pmatrix} u \\ v \\ \end{pmatrix} = \begin{pmatrix} N_1u_1 + N_2u_2 + N_3u_3 \\ N_1v_1 + N_2v_2 + N_3v_3 \\ \end{pmatrix}\]

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]

View Python Code

import sympy as sp
from sympy import solve
a0,a1,a2,x1,x2 = sp.symbols("a_0 a_1 a_2 x_1 x_2")
b0,b1,b2= sp.symbols("b_0 b_1 b_2")
u1,u2,u3,v1,v2,v3 = sp.symbols("u_1 u_2 u_3 v_1 v_2 v_3")
u = a0+a1*x1+a2*x2
v = b0+b1*x1+b2*x2
display("u =",u,"v =",v)
Coordinates = [{x1:0,x2:0},{x1:1,x2:0},{x1:0,x2:1}]
display("Coordinates: ",Coordinates)
Eq1 = u.subs(Coordinates[0])
Eq2 = u.subs(Coordinates[1])
Eq3 = u.subs(Coordinates[2])
Eq4 = v.subs(Coordinates[0])
Eq5 = v.subs(Coordinates[1])
Eq6 = v.subs(Coordinates[2])
s = solve((Eq1-u1,Eq2-u2,Eq3-u3,Eq4-v1,Eq5-v2,Eq6-v3),
          (a0,a1,a2,b0,b1,b2))
display("Solve: ",s)
u = u.subs(s).expand()
v = v.subs(s).expand()
display("u =",u,"v =",v)
N1_u = u.coeff(u1)
N1_v = v.coeff(v1)
N2_u = u.coeff(u2)
N2_v = v.coeff(v2)
N3_u = u.coeff(u3)
N3_v = v.coeff(v3)
display("N1u1 =",N1_u,"N1v1 =",N1_v,"N2u2 =",N2_u,
        "N2v2 =",N2_v,"N3u3 =",N3_u,"N3v3 =",N3_v)

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:

    \[N_1 = 1 - X_1 - X_2\]

    \[N_2 = X_1\]

    \[N_3 = X_2\]

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:

    \[\varepsilon = \begin{pmatrix} \varepsilon_{11} \\ \varepsilon_{22} \\ 2\varepsilon_{12} \\ \end{pmatrix} = \begin{pmatrix} \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{pmatrix}\]

    \[= \begin{pmatrix} \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{pmatrix} \begin{pmatrix} u_1 \\ v_1 \\ u_2 \\ v_2 \\ u_3 \\ v_3 \\ \end{pmatrix}\]

    \[= \begin{pmatrix} -1 & 0 & 1 & 0 & 0 & 0 \\ 0 & -1 & 0 & 0 & 0 & 1 \\ -1 & -1 & 0 & 1 & 1 & 0 \\ \end{pmatrix} \begin{pmatrix} u_1 \\ v_1 \\ u_2 \\ v_2 \\ u_3 \\ v_3 \\ \end{pmatrix}\]

    \[=Bu_e\]

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:

    \[\begin{pmatrix} \sigma_{11} \\ \sigma_{22} \\ \sigma_{12} \\ \end{pmatrix} = \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} \begin{pmatrix} \varepsilon_{11} \\ \varepsilon_{22} \\ 2\varepsilon_{12} \\ \end{pmatrix}\]

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

    \[K^e = \int_e B^T CB dX = \int_e \begin{pmatrix} -1 & 0 & -1 \\ 0 & -1 & -1 \\ 1 & 0 & 0 \\ 0 & 0 & 1 \\ 0 & 0 & 1 \\ 0 & 1 & 0 \\ \end{pmatrix} C \begin{pmatrix} -1 & 0 & 1 & 0 & 0 & 0 \\ 0 & -1 & 0 & 0 & 0 & 1 \\ -1 & -1 & 0 & 1 & 1 & 0 \\ \end{pmatrix} dX\]

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)} \begin{pmatrix} \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{1}{2-4\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} & \frac{\nu}{-1+2\nu} & \frac{1-\nu}{-1+2\nu} & 0 & 0 & \frac{\nu}{-1+2\nu} \\ \frac{-1}{2} & \frac{-1}{2} & 0 & \frac{1}{2} & \frac{1}{2} & 0 \\ \frac{-1}{2} & \frac{-1}{2} & 0 & \frac{1}{2} & \frac{1}{2} & 0 \\ \frac{\nu}{-1+2\nu} & \frac{1-\nu}{-1+2\nu} & \frac{\nu}{-1+2\nu} & 0 & 0 & \frac{1-\nu}{1-2\nu} \\ \end{pmatrix}\]

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

View Python Code

import sympy as sp
from sympy import Matrix,simplify
Ee,nu,t = sp.symbols("E nu t")
B = Matrix([[-1,0,1,0,0,0],[0,-1,0,0,0,1],[-1,-1,0,1,1,0]])

Cc = Ee*t/2/(1 + nu)*Matrix([[(1-nu)/(1-2*nu),nu/(1-2*nu),0], 
                         [nu/(1-2*nu),(1-nu)/(1-2*nu), 0], 
                         [0, 0, 1/2]])
K = simplify(B.T*Cc*B)
display("K^e =",K)
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 dX = \rho t \int_0^1 \int_0^{1-X_2} \begin{pmatrix} N_1 & 0 \\ 0 & N_1 \\ N_2 & 0 \\ 0 & N_2 \\ N_3 & 0 \\ 0 & N_3 \\ \end{pmatrix} \begin{pmatrix} b_1 \\ b_2 \\ \end{pmatrix} dX_1dX_2 = \frac{\rho t}{6} \begin{pmatrix} b_1 \\ b_2 \\ b_1 \\ b_2 \\ b_1 \\ b_2 \\ \end{pmatrix}\]

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

View Python Code

import sympy as sp
from sympy import Matrix,integrate
x1,x2,rb1,rb2 = sp.symbols("x_1 x_2 rho_1 rho_2")
Shapefun = Matrix([[1-x1-x2],[x1],[x2]])
Nn = Matrix([[0 for x in range(6)] for y in range(2)])
for i in range(3):
    Nn[0,2*i] = Nn[1,2*i+1] = Shapefun[i]
rb = Matrix([[rb1],[rb2]])
f = Matrix([[integrate(i,(x1,0,1-x2),(x2,0,1))]for i in Nn.T*rb])
display("Nodal Force due to Body Force:")
display("f^e =",f)
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 dS = t \int_0^1 \begin{pmatrix} N_1 & 0 \\ 0 & N_1 \\ N_2 & 0 \\ 0 & N_2 \\ N_3 & 0 \\ 0 & N_3 \\ \end{pmatrix} \begin{pmatrix} t_1 \\ t_2 \\ \end{pmatrix} dX_2 = \frac{t}{2} \begin{pmatrix} t_1 \\ t_2 \\ 0 \\ 0 \\ t_1 \\ t_2 \\ \end{pmatrix}\]

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

View Python Code

import sympy as sp
from sympy import Matrix,integrate
x1,x2,t1,t2 = sp.symbols("x_1 x_2 t_1 t_2")
Shapefun = Matrix([[1-x1-x2],[x1],[x2]])
Nn = Matrix([[0 for x in range(6)] for y in range(2)])
for i in range(3):
    Nn[0,2*i] = Nn[1,2*i+1] = Shapefun[i]
tn = Matrix([[t1],[t2]])
f = Matrix([[(integrate(i,(x2,0,1))).subs({x1:0})]for i in Nn.T*tn])
display("Nodal Force due to Traction Vector on One Side:")
display("f^e =",f)
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} = \begin{pmatrix} u \\ v \\ \end{pmatrix} = \begin{pmatrix} 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{pmatrix}\]

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} = \begin{pmatrix} u \\ v \\ \end{pmatrix} = \begin{pmatrix} 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{pmatrix}\]

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

View Python Code

import sympy as sp
from sympy import solve
a0,a1,a2,a3,a4,a5,x1,x2 = sp.symbols("a_0 a_1 a_2 a_3 a_4 a_5 x_1 x_2")
b0,b1,b2,b3,b4,b5 = sp.symbols("b_0 b_1 b_2 b_3 b_4 b_5")
u1,u2,u3,u4,u5,u6 = sp.symbols("u_1 u_2 u_3 u_4 u_5 u_6")
v1,v2,v3,v4,v5,v6 = sp.symbols("v_1 v_2 v_3 v_4 v_5 v_6")
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
display("u =",u,"v =",v)
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}]
display("Coordinates =",Coordinates)
Eq1 = u.subs(Coordinates[0])
Eq2 = u.subs(Coordinates[1])
Eq3 = u.subs(Coordinates[2])
Eq4 = v.subs(Coordinates[0])
Eq5 = v.subs(Coordinates[1])
Eq6 = v.subs(Coordinates[2])
Eq7 = u.subs(Coordinates[3])
Eq8 = u.subs(Coordinates[4])
Eq9 = u.subs(Coordinates[5])
Eq10 = v.subs(Coordinates[3])
Eq11 = v.subs(Coordinates[4])
Eq12 = v.subs(Coordinates[5])
s = 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))
display("Solve: ",s)
u = u.subs(s).expand()
v = v.subs(s).expand()
display("u =",u,"v =",v)
N1_u = u.coeff(u1)
N1_v = v.coeff(v1)
N2_u = u.coeff(u2)
N2_v = v.coeff(v2)
N3_u = u.coeff(u3)
N3_v = v.coeff(v3)
N4_u = u.coeff(u4)
N4_v = v.coeff(v4)
N5_u = u.coeff(u5)
N5_v = v.coeff(v5)
N6_u = u.coeff(u6)
N6_v = v.coeff(v6)
display("N1u1 =",N1_u,"N1v1 =",N1_v,"N2u2 =",N2_u,
        "N2v2 =",N2_v,"N3u3 =",N3_u,"N3v3 =",N3_v,
        "N4u4 =",N4_u,"N4v4 =",N4_v,"N5u5 =",N5_u,
        "N5v5 =",N5_v,"N6u6 =",N6_u,"N6v6 =",N6_v)

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:

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

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:

    \[\varepsilon = \begin{pmatrix} \varepsilon_{11} \\ \varepsilon_{22} \\ 2\varepsilon_{12} \\ \end{pmatrix} = \begin{pmatrix} \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{pmatrix}\]

    \begin{equation*} \setcounter{MaxMatrixCols}{12} \small = \begin{pmatrix} \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{pmatrix} \begin{pmatrix} 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{pmatrix} \end{equation*}

    \[= Bu_e\]

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:

    \[\begin{pmatrix} \sigma_{11} \\ \sigma_{22} \\ \sigma_{12} \\ \end{pmatrix} = \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} \begin{pmatrix} \varepsilon_{11} \\ \varepsilon_{22} \\ 2\varepsilon_{12} \\ \end{pmatrix}\]

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 CB dX_1 dX_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

View Python Code

import sympy as sp
from sympy import simplify,diff,integrate
x1,x2,t,Ee,nu = sp.symbols("x_1 x_2 t Ee nu")
half=sp.sympify('1/2') 
Shapefun = Matrix([[2*(1-x1-x2)*(half-x1-x2)],[2*x1*(x1-half)],
                   [2*x2*(x2-1/2)],[4*x1*(1-x1-x2)],
                   [4*x1*x2],[4*x2*(1-x1-x2)]])
B = Matrix([[0 for x in range(12)] for y in range(3)])
for i in range(6):
    B[0,2*i] = B[2,2*i+1] = diff(Shapefun[i],x1)
    B[1,2*i+1] = B[2,2*i] = diff(Shapefun[i],x2)                                                                                                   
Cc = Ee/(1 + nu)*Matrix([[(1-nu)/(1-2*nu),nu/(1-2*nu),0], 
                         [nu/(1-2*nu),(1-nu)/(1-2*nu), 0], 
                         [0, 0, 1/2]])   
Kbeforeintegration = t*simplify(B.T*Cc*B)
K = Matrix([[simplify(integrate(Kbeforeintegration[i,j],(x1,0,1-x2),(x2,0,1)))for i in range (12)] for j in range (12)])
display("K^e =",K)

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 dX = \rho t \int_0^1 \int_0^{1-X_2} \begin{pmatrix} 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{pmatrix} \begin{pmatrix} b_1 \\ b_2 \\ \end{pmatrix} dX_1dX_2 = \frac{\rho t}{6} \begin{pmatrix} 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ b_1 \\ b_2 \\ b_1 \\ b_2 \\ b_1 \\ b_2 \\ \end{pmatrix}\]

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

View Python Code

import sympy as sp
from sympy import simplify,integrate,N,Matrix
x1,x2,rb1,rb2 = sp.symbols("x_1 x_2 rho_1 rho_2")
half=sp.sympify('1/2') 
Shapefun = Matrix([[2*(1-x1-x2)*(half-x1-x2)],[2*x1*(x1-half)],
                   [2*x2*(x2-1/2)],[4*x1*(1-x1-x2)],
                   [4*x1*x2],[4*x2*(1-x1-x2)]])
Nn = Matrix([[0 for x in range(12)] for y in range(2)])
for i in range(6):
    Nn[0,2*i] = Nn[1,2*i+1] = Shapefun[i]                                                                                                 
rb = Matrix([rb1,rb2])
solve = Nn.T*rb
f = Matrix([[simplify(integrate(solve[i],(x1,0,1-x2),(x2,0,1)))]for i in range (12)])
display("Nodal Force due to Body Force:")
display("f^e =",f)
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 dS = t \int_0^1 \begin{pmatrix} 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{pmatrix} \begin{pmatrix} t_1 \\ t_2 \\ \end{pmatrix} dX_2 = \frac{t}{6} \begin{pmatrix} t_1 \\ t_2 \\ 0 \\ 0 \\ t_1 \\ t_2 \\ 0 \\ 0 \\ 0 \\ 0 \\ 4t_1 \\ 4t_2 \\ \end{pmatrix}\]

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

View Python Code

import sympy as sp
from sympy import Matrix,integrate
x1,x2,t1,t2 = sp.symbols("x_1 x_2 t_1 t_2")
Shapefun = Matrix([[1-x1-x2],[x1],[x2]])
Nn = Matrix([[0 for x in range(6)] for y in range(2)])
for i in range(3):
    Nn[0,2*i] = Nn[1,2*i+1] = Shapefun[i]
tn = Matrix([[t1],[t2]])
f = Matrix([[(integrate(i,(x2,0,1))).subs({x1:0})]for i in Nn.T*tn])
display("Nodal Force due to Traction Vector on One Side:")
display("f^e =",f)
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.

Video:

Page Comments

  1. Isha says:
    March 9, 2022 at 7:38 pm

    Hey, I was wondering where you discuss about doing the integration quickly (at timestamp 12:50)? Is it explained in another video? Thanks.

    Reply
    1. Samer Adeeb says:
      March 9, 2022 at 7:57 pm

      Yes, please review Gauss integration at a later section.

      Reply

Leave a Reply

Your email address will not be published.