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

where , , , , , and are six generalized degrees of freedom of the element. The approximate displacement function in terms of the nodal degrees of freedom has the form:

There are two ways to find the shape functions . 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 has to satisfy that the condition that at node 1, , it varies linearly and is equal to zero on the line . Similarly, is equal to 1 at node 2 and is equal to zero on the line . Also, is equal to 1 at node 3 and is equal to zero on the line . These result in the following expressions for the shape functions:

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

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

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

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

Assuming that the triangle has a constant thickness , then the differential volume . In addition, all the components of the matrix are constant; therefore, the integral can be evaluated by simply multiplying the matrix by , which represents the volume of the triangle:

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 , and that the mass density is , then the nodal forces due to the distributed body forces can be obtained using the integral:

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

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)


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 contains linear expressions in the coordinates and 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:

where , , , , , , , , , , , and are twelve generalized degrees of freedom of the element. The approximate displacement function in terms of the nodal degrees of freedom has the form:

There are two ways to find the shape functions . 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 has to satisfy that the condition that at node 1, , and it is equal to zero on the lines and . Similarly, is equal to 1 at node 2 and is equal to zero on the lines and . Also, is equal to 1 at node 3 and is equal to zero on the lines and . is equal to 1 at node 4 and is equal to zero on the lines and . is equal to 1 at node 5 and is equal to zero on the lines and . is equal to 1 at node 6 and is equal to zero on the lines and . . These result in the following expressions for the shape functions:

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

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

The different entries of the matrix are indeed linear functions of the coordinates and 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 matrix 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:

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

The stiffness matrix has the dimensions of , 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 , and that the mass density is , then the nodal forces due to the distributed body forces can be obtained using the integral:

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

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)