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}
(*Applied Axial Load*)
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}
(*Applied Axial Load*)
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
Nodal Force Vector 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
Program Outputs
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"}]
For further information on finite element analysis and Mathematica, readers are encouraged to visit this page.
#FiniteElement #StressStrain #Mathematica #Abaqus #FEA #Blog #Blogger
No comments:
Post a Comment