A simple truss problem

exact solution

Clear[A, x, F, Em, eps, uexact]
A[x_]  :=  10 - 0.09 x
F = 2*10^4;
Em = 3*10^6;
eps[x_] = F/(Em*A[x])
uexact[x_] =  Integrate[eps[s], {s, 0, x}]

uexact[100]

\(\(g1 = 
        Plot[{uexact[x]}, {x, 0, 
            100}, 
		AxesLabel -> {"x", "Displacement"}, 

          PlotStyle -> {Dashing[{0.03, 0.02}], Dashing[{}]}];\)\)

\(\(g2 = 
        Plot[
		Evaluate[{D[uexact[x], x]}], 
		\ \ \ {x, 0, 
            100}, 
		AxesLabel -> {"x", "Strain"}, 

          PlotStyle -> {Dashing[{0.03, 0.02}], Dashing[{}]}];\)\)

three elements, variable nodes

Enter the location of the nodes. The values of x0 and x3 have to be 0 and 100.

x0 = 0; x1 = 100/3; x2 = 200/3; x3 = 100;

The information about the truss in question is entered and then the stiffnes matrix computed.

\(Clear[K, A, Em, F, Knum]; 
    K := Em*{{A1/L1 + A2/L2, \(-A2\)/L2, 0}, 
			{\(-A2\)/L2, 
            A2/L2 + A3/L3, \(-A3\)/L3}, 
			{0, \(-A3\)/L3, 
            A3/L3}}\)
A[x_]  :=  10 - 0.09*x;
Em := 3*10^6; 
F  :=  20*10^3; 
\(\(Knum = 
        K /. {L1 -> \((x1 - x0)\), L2 -> \((x2 - x1)\), 
            L3 -> \((x3 - x2)\), 
		A1\  -> A[\((x0 + x1)\)/2], 
            A2 -> A[\((x1 + x2)\)/2], A3 -> A[\((x2 + x3)\)/2]};\)\)

The system of linear equations is solved. The Eigenvalues are computed to very that we have in fact a solution with minimal energy.

{u1, u2, u3} = LinearSolve[Knum, {0, 0, F}]
u0 = 0;
Eigenvalues[Knum]

With a piecewise linear interpolation we compute the numerical solution and the generate the plots to compare with the exact solution.

\(\(unum = 
        Interpolation[{{x0, u0}, {x1, u1}, {x2, u2}, {x3, u3}}, 
          InterpolationOrder\  -> 1];\)
\)
\(\(g1 = 
        Plot[{uexact[x], unum[x]}, {x, x0, 
            x3}, 
		AxesLabel -> {"x", "Displacement"}, 

          PlotStyle -> {Dashing[{0.03, 0.02}], Dashing[{}]}];\)
\)
\(\(g2 = 
        Plot[
		Evaluate[{D[uexact[x], x], 
              D[unum[x], x]}], 
		\ \ \ {x, x0, 
            x3}, 
		AxesLabel -> {"x", "Strain"}, 
		PlotRange -> {0, 0.006}, 

          PlotStyle -> {Dashing[{0.03, 0.02}], Dashing[{}]}];\)\)

Variable number of elements, variable nodes

Enter the location of the nodes. The values of x0 and x3 have to be 0 and 100.

nodes = {0, 100/3, 200/3, 100};
nodes = Table[x, {x, 0, 100, 10}]
A[x_]  :=  10 - 0.09*x;
Em       =  3*10^6;
Force =  2*10^4;

Compute all necessary data of the elements.

n = Length[nodes] - 1;
leng = Drop[nodes, 1] - Drop[nodes, -1];
midpoint = Drop[nodes, -1] + leng/2;
area = Map[A, midpoint];

Create the global stiffness matrix and the fill in the elements.

K = Table[0, {n}, {n}];
 K[[1, 1]] += Em*area[[1]]/leng[[1]];
\(For[k = 2, 
      k <= n, \(k++\), 
\ \ K[\([k - 1, k - 1]\)] += 
        Em*area[\([k]\)]/leng[\([k]\)]; 
\ \ K[\([k, k]\)] += 
        Em*area[\([k]\)]/leng[\([k]\)]; 
\ \ K[\([k - 1, k]\)] -= 
        Em*area[\([k]\)]/leng[\([k]\)]; 
\ \ K[\([k, k - 1]\)] -= 
        Em*area[\([k]\)]/leng[\([k]\)];
]\)
rhs = Table[0, {n}];
rhs[[n]] = Force;
displacement = Flatten[{0, LinearSolve[K, rhs]}];

With a piecewise linear interpolation we compute the numerical solution and the generate the plots to compare with the exact solution.

data = Transpose[{nodes, displacement}];
\(\(unum = Interpolation[data, InterpolationOrder\  -> 1];\)
\)
\(\(g1 = 
        Plot[{uexact[x], unum[x]}, {x, x0, 
            x3}, 
		AxesLabel -> {"x", "Displacement"}, 

          PlotStyle -> {Dashing[{0.03, 0.02}], Dashing[{}]}];\)
\)
\(\(g2 = 
        Plot[
		Evaluate[{D[uexact[x], x], 
              D[unum[x], x]}], 
		\ \ \ {x, x0, 
            x3}, 
		AxesLabel -> {"x", "Strain"}, 
		PlotRange -> {0, 0.006}, 

          PlotStyle -> {Dashing[{0.03, 0.02}], Dashing[{}]}];\)\)

The test problem solved with elements of order 2

<< BVP.m
<< Interpol2.m

Clear[A, b, f, x, data]
A = Function[x, -(10 - 0.09*x)*3*10^6];
b = Function[x, 0];
f = Function[x, 0];
n = 5;
uexact[
        x_] =  -Integrate[
          2*10^4/A[s], {s, 0, x}]
x = Table[s, {s, 0, 100, 100/n}];
data = 
        BVP[A, b, f, x, {"D", "N"}, {0, 2*10^4}];

\(\(g1 = 
        Plot[{uexact[xt], Interpol2[xt, data]}, {xt, 0, 100}, 
          PlotRange\  -> 
            All, 	
		AxesLabel -> {"x", "Displacement"}, 

          PlotStyle -> {Dashing[{0.03, 0.02}], Dashing[{}]}];\)\)

Duexact[xt_] = D[uexact[xt], xt];
\(\(g2 = 
        Plot[{Duexact[xt], DInterpol2[xt, data]}, {xt, 0, 100}, 
          PlotRange\  -> {0, 
              0.007}, 	
		AxesLabel -> {"x", "Strain"}, 

          PlotStyle -> {Dashing[{0.03, 0.02}], Dashing[{}]}];\)\)

\(\(g3 = 
        Plot[\((uexact[xt] - Interpol2[xt, data])\)/uexact[100], {xt, 0, 
            100}, 
PlotRange\  -> All];\)\)