Poiseuille Flow in a Rectangular Channel
Author's Data: Housam Binous
Department of Chemical Engineering
National Institute of Applied Sciences and Technology
Tunis, TUNISIA
Email: binoushousam@yahoo.com

LU Decomposition leads to fill-ins that makes the solution by elimination not manageable.

Np = 30; Np^2

t = Table[a[i, j], {i, 1, 900}, {j, 1, 900}];

Δx = 1/29.

For[i = 1, 
      i < 901, {For[j = 1, 
          j < 901, {a[i, j] = 0, j++}], i++}]

For[k = 31, 
      k < 871, {If[
          Mod[k, 30] ≠ 1   &&  Mod[k, 30] ≠  0,  
          a[k, k] = 4/Δx^2], k++}]

For[k = 31, 
      k < 871, {If[Mod[k, 30] == 0  ||  Mod[k, 30] ==  1,  
          a[k, k] = 1], k++}]

For[k = 1, k < 31, {a[k, k] = 1, k++}]

For[k = 871, k < 901, {a[k, k] = 1, k++}]

For[k = 30, 
      k < 871, {If[
          Mod[k, 30] ≠ 1  &&  
            Mod[k, 30] ≠ 
              0, {a[k, k - 1] = -1/Δx^2, 
            a[k, k + 1] = -1/Δx^2}], k++}]

For[k = 30, 
      k < 871, {If[
          Mod[k, 30] ≠ 1   &&  
            Mod[k, 30] ≠ 
              0, {a[k, k - 30] = -1. /Δx^2, 
            a[k, k + 30] = -1/Δx^2}], k++}]

ArrayPlot[t, 
        ColorRules -> {0 -> White, 
            0.  -> White, _ -> Red}];

{LU, P, cond} = LUDecomposition[t];

U = 
        Transpose[LU] Table[
            If[i ≤ j, 1, 0], {i, Length[LU]}, {j, 
              Length[LU]}];
ArrayPlot[U, 
        ColorRules -> {0 -> White, 
            0.  -> White, _ -> Blue}];

L = 
        Transpose[LU] - U + 
          IdentityMatrix[Length[LU]];
ArrayPlot[L, 
        ColorRules -> {0 -> White, 
            0.  -> White, _ -> Cyan}];

ArrayPlot[Part[IdentityMatrix[900], P], 
        ColorRules -> {0 -> White, 
            0.  -> White, _ -> RGBColor[0, 0.8, 0]}];

k = 0;

For[i = 1, 
      i < 901, {For[j = 1, 
          j < 901, {If[
              L[[i, j]] ≠ 0, k++], j++}], i++}]; 
    Print[k]

Count[Flatten[L], _?Positive] + 
      Count[Flatten[L], _?Negative]

k = 0;

For[i = 1, 
      i < 901, {For[j = 1, 
          j < 901, {If[
              U[[i, j]] ≠ 0, k++], j++}], i++}]; 
    Print[k]

Count[Flatten[U], _?Positive] + 
      Count[Flatten[U], _?Negative]

Count[Flatten[t], _?Positive] + 
      Count[Flatten[t], _?Negative]

(Part[t, P] - L . U // Chop)

Count[(Part[t, P] - L . U // Chop), _?Positive] + 
      Count[(Part[t, P] - L . U // Chop), _?Negative]

Solution to Poiseuille Flow in a rectangular channel using finite differences agree with the analytical result
obtained using separation of variable.

c = Table[b[i, j], {i, 1, 30}, {j, 1, 30}];

For[i = 1, 
      i < 31, {For[j = 1, 
          j < 31, {b[i, j] = 5, j++}], i++}]

For[k = 1, k < 31, {b[1, k] = 0, k++}]

For[k = 1, k < 31, {b[30, k] = 0, k++}]

For[i = 1, 
      i < 31, {For[j = 1, 
          j < 31, {If[Mod[j, 30] ==  1  ||  Mod[j, 30] ==  0, 
              b[i, j] = 0], j++}], i++}]

c

d = Inverse[t] . Flatten[c] // Chop

do = 
      Table[Table[d[[i + (j - 1) 30]], {i, 1, 30}], {j, 1, 
          30}]

ListPlot3D[do]

Plot3D[-5/2 (x - 1) x + 
        Sum[(-4 1^2*5*
                Cosh[(n Pi (1 - 2 y))/(2 1)] Sech[
                  1 n Pi/(2 1)]*
                Sin[n Pi/2]^2 Sin[n Pi x/1])/(n^3 Pi^3), {n, 1, 
            200}], {x, 0, 1}, {y, 0, 1}]

tbl1 = 
        Table[-5/2 (x - 1) x + 
              Sum[(-4 1^2*5*
                      Cosh[(n Pi (1 - 2 y))/(2 1)] Sech[
                        1 n Pi/(2 1)]*
                      Sin[n Pi/2]^2 Sin[n Pi x/1])/(n^3 Pi^3), {n, 
                  1, 200}] /. x ->  4 Δx, {y, 0, 
            1, Δx}];

i = 1 + 4

tbl2 = Table[do[[5]]];

Length[tbl2]

Length[tbl1]

Transpose[{{tbl1}, {tbl2}}] // TableForm

plt1 = 
      ListPlot[Table[{y, -5/2 (x - 1) x + 
                Sum[(-4 1^2*5*
                        Cosh[(n Pi (1 - 2 y))/(2 1)] Sech[
                          1 n Pi/(2 1)]*
                        Sin[n Pi/2]^2 Sin[
                          n Pi x/1])/(n^3 Pi^3), {n, 1, 200}] /. 
              x ->  4 Δx}, {y, 0, 1, Δx}], 
        PlotStyle ->  RGBColor[1, 0, 0], 
        PlotJoined ->  True]

plt2 = 
      ListPlot[Table[{(i - 1) Δx, do[[5, i]]}, {i, 1, 
            30}], PlotStyle ->  RGBColor[0, 0, 1], 
        PlotJoined ->  True]

Show[plt1, plt2]