The Finite Element Method Basic Concepts and Applications

by Pepper/Heinrich

FEM Utilities.

Locate[x_, l_List] := 
        Module[{iStart, iMid, iEnd}, 
iStart = 
            1; 
iEnd = Length[l]; 
While[
            iEnd - iStart > 1, 
iMid = 
              Floor[(((iStart + iEnd))/2)]; 
If[
              x < l[[iMid]], 
  
iEnd = iMid, 
iStart = 
                iMid];
]; 
Return[
            iStart];
];

feSol[r_, rData_List, tData_List] := 
        Module[{el}, 
el  = Locate[r, rData];   (* 
            Return the element that this value of r is in. *) 
Return[
tData[[el]] + (((tData[[el + 1]] - tData[[el]]))/((rData[[el + 1]] - 
                        rData[[el]]))) (r - 
                    rData[[el]])
];
];

Page 37

m1 = {{(1/4) ξ (ξ - 1) (2  ξ - 
                  1), (1/2) (1 - (ξ^2)) (2  ξ - 
                  1), (1/4) ξ (ξ + 1) (2  ξ - 
                  1)}, 
{-(ξ^2) (ξ - 
                  1), -2 ξ (1 - (ξ^2)), -(ξ^2) (ξ + 1)}, 
{(1/4) ξ (ξ - 
                  1) (2  ξ + 
                  1), (1/2) (1 - (ξ^2)) (2  ξ + 
                  1), (1/4) ξ (ξ + 1) (2  ξ + 
                  1)}};
MatrixForm[m1]
m2 = {{20 (2  ξ - 1), -80 ξ, 
            20 (2  ξ + 1)}, 
{25 (2  ξ - 
                  1), -100 ξ, 
            25 (2  ξ + 1)}, 
{30 (2  ξ - 
                  1), -120 ξ, 
            30 (2  ξ + 1)}};
MatrixForm[m2]

m1 . m2 // Simplify

Factor[20  Out[5]] // MatrixForm

MatrixForm[(1/(50 (ξ + 5))) Out[24]]

Integrate[
          Out[24], {ξ, -1, 1}] + {{100, 0, 0}, {0, 0, 0}, {0, 0, 0}} // 
      MatrixForm

LinearSolve[{{68, -72, 10}, {-72, 160, -88}, {0, 0, 
          1}}, {2400, 0, 39.18}]

Page 42

LinearSolve[{{5, -1.7, 0}, {-1.7, 10, -2.3}, {0, 0, 
          1}}, {237.54, 235.08, 39.18}]

Page 51

m = {{1, 0, 0, 0}, {1, (h/3), ((h^2)/9), ((h^3)/27)}, {1, 2  (h/3), 
            4  ((h^2)/9), 8  ((h^3)/27)}, {1, h, (h^2), (h^3)}};

Inverse[m] . {T1, T2, T3, T4} // Simplify

Plot[{(1 - 3  x) (1 - 3  (x/2)) (1 - x), 
          9  x (1 - 3  (x/2)) (1 - x), ((-9)/2) 
            x (1 - 3  x) (1 - x), 
          x (1 - 3  x) (1 - (3/2) x)}, {x, 0, 1}];

(3/((1 - 3  (x/h)))) + (3/((1 - 3  (x/(2  h))))) + (1/((1 - (x/h)))) // 
      Together

Page 51 Problem 3.3

te1 = ((8  
                  x)/((Sqrt[2]) π)) (1 - ((2  x)/π)) + ((2  
                  x)/π) (((4  x)/π) - 
                1);
te2 = (1 - ((4 (x - (π/2)))/π)) (1 - ((2 (x - (π/2)))/π)) + ((8 (x - (π/2)))/((Sqrt[2]) π)) (1 - ((2 (x - (π/2)))/π));
T1[x_] = If[x < (π/2), te1, te2];

Simplify[te2]

Simplify[
      2 (1 - (x/π)) (3 - 
              4  (x/π)) + (8/(Sqrt[2])) ((2/π) x - 1) (1 - 
              (x/π))]

Plot[{Sin[x], T1[x]}, {x, 0, π}];

Plot[{Sin[x] - T1[x]}, {x, 0, π}];

Simplify[D[te1, x]]
Simplify[D[te2, x]]

T1Prime[x_] = If[x < (π/2), Out[15], Out[16]]

Plot[{T1Prime[x], Cos[x]}, {x, 0, π}];

Plot[{T1Prime[x] - Cos[x]}, {x, 0, π}];

Page 51 Problem 3.4

T = (Subscript[α, 1]) (1 - 3  (x/h)) (1 - 3  (x/(2  h))) (1 - 
              (x/h)) + ((Subscript[α, 1]) + (Subscript[α, 2]) 
                (h/3) + (Subscript[α, 3]) ((h^2)/9)) (9  (x/h)) (1 - 
              3  (x/(2  h))) (1 - 
              (x/h)) + ((Subscript[α, 1]) + (Subscript[α, 2]) ((2  h)/3) + (Subscript[α, 3]) ((4  (h^2))/9)) (-9 (x/(2  h))) (1 - 
              3  (x/h)) (1 - 
              (x/h)) + ((Subscript[α, 1]) + (Subscript[α, 2]) h + (Subscript[α, 3]) 
                (h^2)) ((x/h)) (1 - 
              3  (x/h)) (1 - ((3  x)/(2  h)))

Coefficient[T, (Subscript[α, 1])]
Coefficient[T, (Subscript[α, 2])]
Coefficient[T, (Subscript[α, 3])]

Page 51 Problem 3.5

Tlinear[x_] = 
        If[x < 5, 
          160 (1 - (x/5)) + (142/5) 
              x, 
142 (1 - (1/5) (x - 
                        5)) + (100/5) (x - 5)];

Tquad[x_] = 
        160 (1 - 2  (x/10)) (1 - (x/10)) + 
          4  (x/10) (1 - (x/10)) 142 + (x/10) (2  (x/10) - 1) 
            100;

Plot[{Tlinear[x], Tquad[x]}, {x, 0, 10}];

Plot[{Tlinear[x] - Tquad[x]}, {x, 0, 10}];

Page 52 Problem 3.8

m1 = {{140.12, -130.12, 0, 0}, 
{-130.12, 
            280.24, -150.12, 0}, 
{0, -150.12, 
            320.27, -170.15}, 
{0, 0, 0, 
            1}};
MatrixForm[m1]

rad1 = { .2,  .233333,  .266666,  .3};
t1 = LinearSolve[m1, {4000, 0, 0, 39.18}]

m2 = {{220, -210, 0, 0}, 
{-210, 
            366.66, -156.66, 0}, 
{0, -150.66, 
            266.66, -110}, 
{0, 0, 0, 
            1}};
MatrixForm[m2]

rad2 = { .2,  .22,  .25,  .3};
t2 = LinearSolve[m2, {4000, 0, 0, 39.18}]

T[
          r_] := ((h r1 (Tl - TInf))/(K - h r1 Log[(r1/r2)])) 
            Log[(r/r2)] + Tl;

rl = {r1 ->  .2, r2 ->  .3, K -> 20, TInf -> 400, 
          Tl -> 39.18, h -> 50};

exactGr = 
        Plot[T[r] /. rl, {r, r1 /. rl, 
            r2 /. rl}, 
PlotStyle -> {RGBColor[0, 1, 
                0]}];

feES = 
        Plot[feSol[r, rad1, t1], {r, r1 /. rl, r2 /. rl}];

feUES = 
        Plot[feSol[r, rad2, t2], {r, r1 /. rl, r2 /. rl}];

Show[{exactGr, feES, feUES}];

I looks like to me that the Finite Element method with unequally space elements did much worse than the Finite Element method with equally spaced elements.

Page 53 Problem 3.13

a = {{1, 0, 0}, {((-k)/L1), 
          (k/L1) + (k/L2), ((-k)/L2)}, {0, ((-k)/L2), 
          (k/L2) + h}}
b = {Tw, 0, h TInft}
LinearSolve[a, b]

Out[9] /. {k -> 75, h -> 15, TInft -> 60, 
        L1 ->  .04, L2 ->  .06, Tw -> 150}

Page 53 Problem 3.15

<< LinearAlgebra`Tridiagonal`

a = {-25, -27, -29, -31, -33, -35};   
b = {26, 52, 56, 60, 64, 68, 1};   
c = {-25, -27, -29, -31, -33, 
          0};
t3 = 
      TridiagonalSolve[c, b, a, {400, 0, 0, 0, 0, 0. , 39.18}]

x3 = Table[ .(1/6) i, {i, 0, 6}]

sixLinElt = 
        Plot[feSol[x, x3, t3], {x, 
            0,  .1}, 
PlotStyle -> 
            RGBColor[1, 0, 0]];

Example 3.1

x1 = {0. ,  .05,  .1};
t1 = {99.822, 66.469, 39.18};
twoLinElt = 
        Plot[feSol[x, x1, t1], {x, 
            0,  .1}, 
PlotStyle -> 
            RGBColor[0, 1, 0]];

Example 3.5

oneQuadEltSol = 
        99.992 (((((x -  .05) (x -  .1)))/((- .05) (- .1)))) + 
          66.546 ((((x (x -  .1)))/(( .05) (- .05)))) + 39.18 ((((x (x -  .05)))/(( .1) ( .05))));
oneQuadElt = 
        Plot[oneQuadEltSol, {x, 0,  .1}, 
PlotStyle -> 
            RGBColor[0, 0, 1]];

Comparison of the three approaches

Show[{sixLinElt, twoLinElt, oneQuadElt}];

Page 53 Problem 3.16

n = {(1/2) ξ (ξ - 1), 
          1 - (ξ^2), (1/2) ξ (ξ + 1)};

m = Outer[Times, n, n]

Integrate[m, {ξ, -1, 1}] // MatrixForm

Continuing the example of 3.6 with Backwards implicit scheme requres θ=1.

mDt = {{2 +  .03 Δt, 1 -  .027 Δt, 
            0}, 
{1 -  .027 Δt, 
            4 +  .06  Δt, 
            1 -  .033 Δt}, 
{0, 0, 
            1}};
m = {{2, 1, 0}, 
{1, 4, 1}, 
{0, 
            0, 0}};
b = {1.2 Δt, 0, 39.18};

Δt = (10^(-3));

f1[x_] := Inverse[mDt] . (m . x + b)
f2[x_List] := LinearSolve[mDt, m . x + b];

Nest[f1, {39.18, 39.18, 39.18}, 50]
Nest[f2, {39.18, 39.18, 39.18}, 50]

Page 58

Solve[{(Subscript[ϕ, 1]) == (Subscript[α, 1]) + (Subscript[α, 2]) 
              (Subscript[x, 1]) + (Subscript[α, 3]) 
              (Subscript[y, 1]), 
(Subscript[ϕ, 2]) == (Subscript[α, 1]) + (Subscript[α, 2]) (Subscript[x, 2]) + (Subscript[α, 3]) 
              (Subscript[y, 2]), 
(Subscript[ϕ, 3]) == (Subscript[α, 1]) + (Subscript[α, 2]) (Subscript[x, 3]) + (Subscript[α, 3]) 
              (Subscript[y, 3])}, {(Subscript[α, 1]), (Subscript[α, 2]), (Subscript[α, 3])}]