Welcome to the World of Modelling and Simulation

What is Modelling?

This blog is all about system dynamics modelling and simulation applied in the engineering field, especially mechanical, electrical, and ...

Stress and Strain Analysis of Linear Finite Elements

In this tutorial, we will use principle of virtual work method to determine the stress and displacement of the nodes of linear one dimensional finite elements. Following figure 1 represents a cantilever beam which is stressed axially. Let's assume the parameters L, A, E, c have the same value of 1 unit for simplicity. Our goal is to compare the finite element results with the exact solutions of a typical cantilever problem.
For this problem, we need to consider two quadratic one dimensional element. So, the displacement is interpolated among the three nodes and has the following form:

 Fig.1: Showing one dimensional axially load bar given in the assignment problem.

 Fig.2: Showing one dimensional piecewise linear interpolation.

Mathematica Program for Element 1

Clear[x, a0, a1, a2, u1, u2, u3, L]

(*Definition of quadratic form of displacement*)
u = a0 + a1*x + a2*x^2;

eq1 = (u /. x -> 0) - u1;
eq2 = (u /. x -> L/4) - u2;
eq3 = (u /. x -> L/2) - u3;

sol = Solve[{eq1 == 0, eq2 == 0, eq3 == 0}, {a0, a1, a2}];
u = u /. sol[[1]]

N1 = Coefficient[u, u1];
N2 = Coefficient[u, u2];
N3 = Coefficient[u, u3];

(*Shape Functions*)
Nn = {N1, N2, N3}

p = c*x;

(*Nodal Force Vectors*)
fe = Integrate[Nn*p, {x, 0, L/2}];
c = 1;
L = 1;
fe // MatrixForm

Program Outputs

Displacement Vector for Element 1:

Shape Functions for Element 1:

Nodal Force Vector for Element 1:

Mathematica Program for Element 2

Clear[x, a0, a1, a2, u1, u2, u3, L]

(*Definition of quadratic form of displacement*)
u = a0 + a1*x + a2*x^2;

eq1 = (u /. x -> 0) - u1;
eq2 = (u /. x -> L/4) - u2;
eq3 = (u /. x -> L/2) - u3;

sol = Solve[{eq1 == 0, eq2 == 0, eq3 == 0}, {a0, a1, a2}];
u = u /. sol[[1]]

N1 = Coefficient[u, u1];
N2 = Coefficient[u, u2];
N3 = Coefficient[u, u3];

(*Shape Functions*)
Nn = {N1, N2, N3}

p = c*(x + L/2);

(*Nodal Force Vectors*)
fe = Integrate[Nn*p, {x, 0, L/2}];
c = 1;
L = 1;
fe // MatrixForm

Program Outputs

Displacement Vector for Element 2

Shape Functions for Element 2

Local Stiffness Matrix for Element 1

Clear[x, a0, a1, a2, u1, u2, u3, L, EA]

(*Definition of quadratic form of displacement*)
u = a0 + a1*x + a2*x^2;

eq1 = (u /. x -> 0) - u1;
eq2 = (u /. x -> L/4) - u2;
eq3 = (u /. x -> L/2) - u3;

sol = Solve[{eq1 == 0, eq2 == 0, eq3 == 0}, {a0, a1, a2}];
u = u /. sol[[1]];

N1 = Coefficient[u, u1];
N2 = Coefficient[u, u2];
N3 = Coefficient[u, u3];

B = {{D[N1, x], D[N2, x], D[N3, x]}};

(*Local Stiffness Matrix*)
K = Integrate[EA*Transpose[B].B, {x, 0, L/2}];

L = 1;
EA = 1;
K // MatrixForm

Comparison between the Exact and Solution using Four Elements

We know that the exact solution for the axially loaded bar shown in Fig. 1 may be written as,

Where, for the consecutive four elements, we may write the displacement for each element as,

From the assignment problem, the parameter values are,

From the virtual work method, we may write,

The internal virtual work (IVW) may be expressed as,

And, the external virtual work may be presented as,

MATHEMATICA PROGRAM

The following program is written to determine the nodal displacement using the finite element method considering the principle of virtual work.

Clear[x, a0, a1, a2, u1, u2, u3, L]

(*Definition of quadratic form of displacement*)
u = a0 + a1*x + a2*x^2;

eq1 = (u /. x -> 0) - u1;
eq2 = (u /. x -> L/4) - u2;
eq3 = (u /. x -> L/2) - u3;

sol = Solve[{eq1 == 0, eq2 == 0, eq3 == 0}, {a0, a1, a2}];
u = u /. sol[[1]];

(*Shape Functions*)
N1 = Coefficient[u, u1];
N2 = Coefficient[u, u2];
N3 = Coefficient[u, u3];
Nn = {N1, N2, N3};

B = {{D[N1, x], D[N2, x], D[N3, x]}};

c = 1;
L = 1;
EA = 1;

(*Axial Force on Element 1*)
p1 = c*x;
fe1 = Integrate[Nn*p1, {x, 0, L/2}];

(*Axial Force on Element 2*)
p2 = c*(x + L/2);
fe2 = Integrate[Nn*p2, {x, 0, L/2}];

K1 = Integrate[EA*Transpose[B].B, {x, 0, L/2}];

K2 = Integrate[EA*Transpose[B].B, {x, 0, L/2}];

Ndof = 5;
Ne = 2;
ID1 = {1, 2, 3};
ID2 = {3, 4, 5};

K = ConstantArray[0, {Ndof, Ndof}];
K[[ID1, ID1]] = K1 + K[[ID1, ID1]];
K[[ID2, ID2]] = K2 + K[[ID2, ID2]];

fe = ConstantArray[0, {Ndof, 1}];
fe[[ID1]] = fe1 + fe[[ID1]];
fe[[ID2]] = fe2 + fe[[ID2]];

IDE = {1};
IDN = {2, 3, 4, 5};

KR = Delete[K, IDE];
KR = Delete[Transpose[KR], IDE];
KR = Transpose[KR];

feR = Delete[fe, IDE];

(*Calculation of Nodal Displacements*)
DisR = LinearSolve[KR, feR];

Dis = ConstantArray[0, {Ndof, 1}];
Dis[[IDN]] = DisR + Dis[[IDN]];
Dis // MatrixForm

Program Outputs

So, we see that the exact and finite element solutions are the same.

From the following two plots, we notice that although the nodal displacement is continuous from the FEA result, it is not true for the stress. The stress distribution is indeed discontinuous.

MATHEMATICA PROGRAM FOR COMPARISON

(*Exact Solution*)
uexact = c*L^2/2/EA*x - c/6/EA*x^3;

(*Finite Element Approximation*)
N1 = Piecewise[{{4 x/L, 0 <= x <= L/4}, {4/L (L/2 - x),
L/4 <= x <= L/2}}, 0];
N2 = Piecewise[{{4/L (x - L/4), L/4 <= x <= L/2}, {4/L (3 L/4 - x),
L/2 <= x <= 3 L/4}}, 0];
N3 = Piecewise[{{4/L (x - L/2), L/2 <= x <= 3 L/4}, {4/L (L - x),
3 L/4 <= x <= L}}, 0];
N4 = Piecewise[{{4/L (x - 3 L/4), 3 L/4 <= x <= L}}, 0];

u1 = c*L^3/EA*(47/384);
u2 = c*L^3/EA*(11/48);
u3 = c*L^3/EA*(39/128);
u4 = c*L^3/EA*(1/3);

(*Calculating Nodal Displacement*)
u = u1*N1 + u2*N2 + u3*N3 + u4*N4;
uexactp = uexact /. {c -> 1, EA -> 1, L -> 1}

sigmaexact = D[uexactp, x]
up = u /. {c -> 1, EA -> 1, L -> 1}
sigmap = D[up, x]

(*Plotting to Show the Comparison between Exact and FEA Solution*)
Plot[{uexactp, up}, {x, 0, 1},
AxesLabel -> {"x", "Nodal Displacement"},
PlotLegends -> {"uexact", "u_Finite Element"}]
Plot[{sigmaexact, sigmap}, {x, 0, 1}, AxesLabel -> {"x", "sigma_11"},
PlotLegends -> {"sigma_exact", "sigma_Finite Element"}]