Message Boards Message Boards

3D stable fluids transition simulation algorithm: laminar to turbulent flow

This algorithm is 3D extension of our 2D algorithm published on this page and here. We suppose that with this code we can simulate transition from laminar to turbulent flow. In this example we compute viscous flow around cuboid at Reynolds number $Re=6250$. Note, that code has been tested up to $Re=10^6$.

dif = 1/6250; pec = .72; U0 = 1.; V0 = 0.; W0 = 0.; dn0 = 1.; n = 80;
n1 = n + 1; sm = 500; r = 20; den = ConstantArray[dn0 , {n1, n1, n1}];
u0 = ConstantArray[U0, {n1, n1, n1}];
v0 = ConstantArray[V0, {n1, n1, n1}]; w0 = 
 ConstantArray[W0, {n1, n1, n1}]; Do[u0[[i, j, k]] = 0; 
 v0[[i, j, k]] = 0; 
 w0[[i, j, k]] = 0;, {i, 10, 20, 1}, {j, Round[n/2 - 5], 
  Round[n/2 + 5], 1}, {k, Round[n/2 - 5], Round[n/2 + 5], 1}];

periodic[n_, up_, ud_, ul_, ur_, ub_] := 
  Module[{bd = ub}, 
   Do[bd[[n + 1, i, j]] = bd[[2, i, j]]; 
    bd[[1, i, j]] = bd[[n, i, j]];, {i, 2, n}, {j, 2, n}];
   Do[bd[[i, 1, j]] = ud;
    bd[[i, n + 1, j]] = up; bd[[i, j, 1]] = ul;
    bd[[i, j, n + 1]] = ur;, {i, 1, n + 1}, {j, 1, n + 1}];
   bd];

diffuse[n_, r_, a_, c_, c0_] := 
  Module[{c1 = c}, 
   Do[Do[Do[
       Do[c1[[i, j, 
            k]] = (c0[[i, j, k]] + 
              a (c1[[i - 1, j, k]] + c1[[i + 1, j, k]] + 
                 c1[[i, j - 1, k]] + c1[[i, j + 1, k]] + 
                 c1[[i, j, k - 1]] + c1[[i, j, k + 1]]))/(1 + 
              6 a);, {k, 2, n}];, {j, 2, n}];, {i, 2, n}];
    Do[c1[[i, j, k]] = 0;, {i, 10, 20, 1}, {j, Round[n/2 - 5], 
      Round[n/2 + 5], 1}, {k, Round[n/2 - 5], Round[n/2 + 5], 
      1}];, {k1, 0, r}];
   c1];

advect[n_, d0_, u_, v_, w_, dt_] := 
  Module[{x, y, z, d1, dt0, i0, i1, j0, j1, k0, k1, s0, s1, t0, t1, 
    p1, p0, d000, d100, d010, d001, d110, d011, d101, d111}, 
   d1 = ConstantArray[0, {n + 1, n + 1, n + 1}]; dt0 = dt n;
   Do[Do[Do[x = i - dt0 u[[i, j, k]]; y = j - dt0 v[[i, j, k]]; 
        z = k - dt0 w[[i, j, k]];
        i0 = Which[x <= 1, 1, 1 < x < n, Floor[x], True, n];
        i1 = i0 + 1;
        j0 = Which[y <= 1, 1, 1 < y < n, Floor[y], True, n];
        j1 = j0 + 1; 
        k0 = Which[z <= 1, 1, 1 < z < n, Floor[z], True, n];
        k1 = k0 + 1; 
        d000 = (d0[[i0, j0, 
            k0]] + (x - i0) (d0[[i1, j0, k0]] - 
              d0[[i0, j0, k0]]) + (y - j0) (d0[[i0, j1, k0]] - 
              d0[[i0, j0, k0]]) + (z - k0) (d0[[i0, j0, k1]] - 
              d0[[i0, j0, k0]])); 
        d100 = (d0[[i1, j0, 
            k0]] + (x - i1) (d0[[i1, j0, k0]] - 
              d0[[i0, j0, k0]]) + (y - j0) (d0[[i1, j1, k0]] - 
              d0[[i1, j0, k0]]) + (z - k0) (d0[[i1, j0, k1]] - 
              d0[[i1, j0, k0]])); 
        d010 = (d0[[i0, j1, 
            k0]] + (x - i0) (d0[[i1, j1, k0]] - 
              d0[[i0, j1, k0]]) + (y - j1) (d0[[i0, j1, k0]] - 
              d0[[i0, j0, k0]]) + (z - k0) (d0[[i0, j1, k1]] - 
              d0[[i0, j1, k0]])); 
        d001 = (d0[[i0, j0, 
            k1]] + (x - i0) (d0[[i1, j0, k1]] - 
              d0[[i0, j0, k1]]) + (y - j0) (d0[[i0, j1, k1]] - 
              d0[[i0, j0, k1]]) + (z - k1) (d0[[i0, j0, k1]] - 
              d0[[i0, j0, k0]])); 
        d110 = (d0[[i1, j1, 
            k0]] + (x - i1) (d0[[i1, j1, k0]] - 
              d0[[i0, j1, k0]]) + (y - j1) (d0[[i1, j1, k0]] - 
              d0[[i1, j0, k0]]) + (z - k0) (d0[[i1, j1, k1]] - 
              d0[[i1, j1, k0]])); 
        d011 = (d0[[i0, j1, 
            k1]] + (x - i0) (d0[[i1, j1, k1]] - 
              d0[[i0, j1, k1]]) + (y - j1) (d0[[i0, j1, k1]] - 
              d0[[i0, j0, k1]]) + (z - k1) (d0[[i0, j1, k1]] - 
              d0[[i0, j1, k0]])); 
        d101 = (d0[[i1, j0, 
            k1]] + (x - i1) (d0[[i1, j0, k1]] - 
              d0[[i0, j0, k1]]) + (y - j0) (d0[[i1, j1, k1]] - 
              d0[[i1, j0, k1]]) + (z - k1) (d0[[i1, j0, k1]] - 
              d0[[i1, j0, k0]])); 
        d111 = (d0[[i1, j1, 
            k1]] + (x - i1) (d0[[i1, j1, k1]] - 
              d0[[i0, j1, k1]]) + (y - j1) (d0[[i1, j1, k1]] - 
              d0[[i1, j0, k1]]) + (z - k1) (d0[[i1, j1, k1]] - 
              d0[[i1, j1, k0]]));
        d1[[i, j, 
          k]] = (d000 + d100 + d010 + d001 + d110 + d011 + d101 + 
            d111)/8;, {k, 2, n}];, {j, 2, n}];, {i, 1, n + 1}]; d1];

project[n_, r_, u0_, v0_, w0_, u_, v_, w_] := 
  Module[{ux = u, vy = v, wz = w, div, p}, 
   p = ConstantArray[0, {n + 1, n + 1, n + 1}];
   div = ConstantArray[0, {n + 1, n + 1, n + 1}];
   ux = ConstantArray[0, {n + 1, n + 1, n + 1}];
   vy = ConstantArray[0, {n + 1, n + 1, n + 1}]; 
   wz = ConstantArray[0, {n + 1, n + 1, n + 1}];
   Do[div[[i, j, 
       k]] = .5/
        n (u0[[i + 1, j, k]] - u0[[i - 1, j, k]] + v0[[i, 1 + j, k]] -
          v0[[i, j - 1, k]] + w0[[i, j, k + 1]] - 
         w0[[i, j, k - 1]]);, {i, 2, n}, {j, 2, n}, {k, 2, n}]; 
   Do[div[[i, j, k]] = 0;, {i, 10, 20, 1}, {j, Round[n/2 - 5], 
     Round[n/2 + 5], 1}, {k, Round[n/2 - 5], Round[n/2 + 5], 1}]; 
   div = periodic[n, 0, 0, 0, 0, div];
   Do[Do[Do[
       Do[p[[i, j, 
            k]] = (div[[i, j, 
               k]] + (p[[i - 1, j, k]] + p[[i + 1, j, k]] + 
                p[[i, j - 1, k]] + p[[i, j + 1, k]] + 
                p[[i, j, k - 1]] + p[[i, j, k + 1]]))/6;, {k, 2, 
          n}];, {j, 2, n}], {i, 2, n}];, {k1, 0, r}]; 
   Do[p[[i, j, k]] = 0;, {i, 10, 20, 1}, {j, Round[n/2 - 5], 
     Round[n/2 + 5], 1}, {k, Round[n/2 - 5], Round[n/2 + 5], 1}]; 
   p = periodic[n, 0, 0, 0, 0, p];
   Do[ux[[i, j, k]] = 
     u0[[i, j, k]] + .5 n (p[[i + 1, j, k]] - p[[i - 1, j, k]]);
    vy[[i, j, k]] = 
     v0[[i, j, k]] + .5 n (p[[i, j + 1, k]] - p[[i, j - 1, k]]); 
    wz[[i, j, k]] = 
     w0[[i, j, k]] + .5 n (p[[i, j, k + 1]] - p[[i, j, k - 1]]);, {i, 
     2, n}, {j, 2, n}, {k, 2, n}]; 
   Do[ux[[i, j, k]] = 0; vy[[i, j, k]] = 0; 
    wz[[i, j, k]] = 0;, {i, 10, 20, 1}, {j, Round[n/2 - 5], 
     Round[n/2 + 5], 1}, {k, Round[n/2 - 5], Round[n/2 + 5], 1}]; {ux,
     vy, wz}];

onestep[n_, step_, r_, a_, uin_, vin_, win_, dt_, c_] := 
 Module[{u1, v1, w1, u, v, w, u0, v0, w0},
  u0 = ConstantArray[0., {n + 1, n + 1, n + 1}];
  v0 = ConstantArray[0., {n + 1, n + 1, n + 1}]; 
  w0 = ConstantArray[0., {n + 1, n + 1, n + 1}];
  u = ConstantArray[0., {n + 1, n + 1, n + 1}];
  v = ConstantArray[0., {n + 1, n + 1, n + 1}]; 
  w = ConstantArray[0., {n + 1, n + 1, n + 1}];
  u1 = ConstantArray[0., {n + 1, n + 1, n + 1}];
  v1 = ConstantArray[0., {n + 1, n + 1, n + 1}]; 
  w1 = ConstantArray[0., {n + 1, n + 1, n + 1}]; u0 = uin; v0 = vin; 
  w0 = win; 
  Do[u0[[i, j, k]] = 0; v0[[i, j, k]] = 0; 
   w0[[i, j, k]] = 0;, {i, 10, 20, 1}, {j, Round[n/2 - 5], 
    Round[n/2 + 5], 1}, {k, Round[n/2 - 5], Round[n/2 + 5], 1}];
  u0 = advect[n, u0, u0, v0, w0, dt]; 
  v0 = advect[n, v0, u0, v0, w0, dt]; 
  w0 = advect[n, w0, u0, v0, w0, dt]; 
  Do[u0[[i, j, k]] = 0; 
   v0[[i, j, k]] = 0;, {i, 10, 20, 1}, {j, Round[n/2 - 5], 
    Round[n/2 + 5], 1}, {k, Round[n/2 - 5], Round[n/2 + 5], 1}];
  u0 = periodic[n, 0, 0, 0, 0, u0]; v0 = periodic[n, 0, 0, 0, 0, v0]; 
  w0 = periodic[n, 0, 0, 0, 0, w0];
  u0 = diffuse[n, r, a, c, u0]; v0 = diffuse[n, r, a, c, v0]; 
  w0 = diffuse[n, r, a, c, w0]; 
  Do[u0[[i, j, k]] = 0; v0[[i, j, k]] = 0; 
   w0[[i, j, k]] = 0;, {i, 10, 20, 1}, {j, Round[n/2 - 5], 
    Round[n/2 + 5], 1}, {k, Round[n/2 - 5], Round[n/2 + 5], 1}];
  u0 = periodic[n, 0, 0, 0, 0, u0]; v0 = periodic[n, 0, 0, 0, 0, v0]; 
  w0 = periodic[n, 0, 0, 0, 0, w0];
  {u0, v0, w0} = project[n, r, u0, v0, w0, u, v, w]; 
  Do[u0[[i, j, k]] = 0; v0[[i, j, k]] = 0; 
   w0[[i, j, k]] = 0;, {i, 10, 20, 1}, {j, Round[n/2 - 5], 
    Round[n/2 + 5], 1}, {k, Round[n/2 - 5], Round[n/2 + 5], 1}];
  u0 = periodic[n, 0, 0, 0, 0, u0]; v0 = periodic[n, 0, 0, 0, 0, v0]; 
  w0 = periodic[n, 0, 0, 0, 0, w0]; {u0, v0, w0}]

cf = With[{cg = Compile`GetElement, hp = HoldPattern, 
     dv = DownValues}, 
    Hold@Compile[{{u0argu, _Real, 3}, {v0argu, _Real, 
              3}, {w0argu, _Real, 3}, {denargu, _Real, 
              3}, {sm, _Integer}, {n, _Integer}, {r, _Integer}, dif, 
             pec}, Module[{u0 = u0argu, v0 = v0argu, w0 = w0argu, uu, 
              vv, ww, dd, den = denargu, 
              c = Table[0., {n + 1}, {n + 1}, {n + 1}], dt = 40./n^2, 
              a, dnup = den[[1, n + 1, 1]], dnd = den[[1, 1, 1]], 
              dnl = den[[1, 1, 1]], dnr = den[[1, 1, n + 1]]}, 
             a = dt dif n n;

             uu = vv = 
               ww = dd = 
                 Table[0., {sm + 1}, {n + 1}, {n + 1}, {n + 1}];

             Do[{u0, v0, w0} = 
               onestep[n, step, r, a, u0, v0, w0, dt, c];
              uu[[step + 1]] = u0;
              vv[[step + 1]] = v0; ww[[step + 1]] = w0;
              den = diffuse[n, r, a/pec, c, den]; 
              Do[den[[i, j, k]] = 0;, {i, 10, 20, 1}, {j, 

                   Round[n/2 - 5], Round[n/2 + 5], 1}, {k, 
                Round[n/2 - 5], Round[n/2 + 5], 1}];
              den = periodic[n, dnup, dnd, dnl, dnr, den];
              den = advect[n, den, u0, v0, w0, dt]; 
              Do[den[[i, j, k]] = 0;, {i, 10, 20, 1}, {j, 
                Round[n/2 - 5], Round[n/2 + 5], 1}, {k, 
                Round[n/2 - 5], Round[n/2 + 5], 1}];
              den = periodic[n, dnup, dnd, dnl, dnr, den];
              dd[[step + 1]] = den;, {step, 0, sm}]; {uu, vv, ww, 
              dd}], CompilationTarget -> "C", 
            RuntimeOptions -> "Speed"] /. dv@onestep /. 
         Flatten[dv /@ {advect, diffuse, periodic, project}] /. 
        hp@ConstantArray[c_, {i_, j_, kc_}] :> 
         Table[0., {i}, {j}, {kc}] /. hp@Part[a__] :> cg[a] /. 
      hp[cg[a__] = rhs_] :> (Part[a] = rhs) // 
     ReleaseHold]; 



rst = cf[u0, v0, w0, den, sm, n, r, dif, pec]; 

Visualization

Do[lst11[s] = 
  Table[{{(i - 1)/n, (j - 1)/n, (k - 1)/n}, rst[[1, s, i, j, k]]}, {i,
     n1}, {j, n1}, {k, n1}];
 lst12[s] = 
  Table[{{(i - 1)/n, (j - 1)/n, (k - 1)/n}, rst[[2, s, i, j, k]]}, {i,
     n1}, {j, n1}, {k, n1}]; 
 lst13[s] = 
  Table[{{(i - 1)/n, (j - 1)/n, (k - 1)/n}, rst[[3, s, i, j, k]]}, {i,
     n1}, {j, n1}, {k, n1}];, {s, 25, sm, 25}]

Do[su1[i] = 
  Interpolation[Flatten[lst11[i], 2], InterpolationOrder -> 3]; 
 sv2[i] = Interpolation[Flatten[lst12[i], 2], 
   InterpolationOrder -> 3]; 
 sw3[i] = Interpolation[Flatten[lst13[i], 2], 
   InterpolationOrder -> 3];, {i, 25, sm, 25}]

Table[Show[
  DensityPlot3D[
   Norm[{su1[i][x, y, z], sv2[i][x, y, z], sw3[i][x, y, z]}], {x, 0, 
    1}, {y, 0.44, 1}, {z, 0, 1}, BoxRatios -> Automatic, 
   ColorFunction -> "BlueGreenYellow", PlotLabel -> i], 
  VectorPlot3D[{su1[i][x, y, z], sv2[i][x, y, z], 
    sw3[i][x, y, z]}, {x, 0, 1}, {y, 0, .5}, {z, 0, 1}, 
   VectorPoints -> Fine, VectorMarkers -> "Arrow"], 
  Graphics3D[{{Blue, 
     Cuboid[{10, Round[n/2 - 5] - 1/2, 
        Round[n/2 - 5] - 1/2}/(n + 1), {20, Round[n/2 + 5] - 1/2, 
        Round[n/2 + 5] - 1/2}/(n + 1)]}}]], {i, 50, sm, 50}]

enter image description here

Animation

Do[lstu1[k] = 
   Flatten[Table[{{(i - 1)/n, (j - 1)/n}, 
      rst[[1, k, i, Round[(n + 1)/2], j]]}, {i, n1}, {j, n1}], 1]; 
  lstw1[k] = 
   Flatten[Table[{{(i - 1)/n, (j - 1)/n}, 
      rst[[2, k, i, Round[(n + 1)/2], j]]}, {i, n1}, {j, n1}], 
    1];, {k, sm}];
Do[Uvel1[i] = Interpolation[lstu1[i], InterpolationOrder -> 3];, {i, 
   1, sm}];
Do[Wvel1[i] = Interpolation[lstw1[i], InterpolationOrder -> 3];, {i, 
   1, sm}];
frame = Table[
   Show[DensityPlot[
     Norm[{Uvel[m][x, y], Vvel[m][x, y]}], {x, 0, 1}, {y, 0, 1}, 
     PlotRange -> All, ColorFunction -> "BlueGreenYellow", 
     Frame -> False, ImageSize -> Tiny, PlotLabel -> m, 
     PlotPoints -> 50], 
    Graphics[{Gray, 
      Rectangle[{10, Round[n/2 - 5] - 1/2}/(n + 1), {20, 
         Round[n/2 + 5] - 1/2}/(n + 1)]}]], {m, 5, sm, 3}];

enter image description here

The question is about code improvement. How can we define parameter r to solve Laplace and Poison equations in this code? Note, r is number of iterations in diffuse and project module (we use Gauss-Seidel relaxation to solve Laplace and Poison's equations).

Update 1. We can reduce r from r=20 as above to r=5 as below. Globally there is no difference in two animations.

enter image description here

But if we compare velocity in some point like $x=0.65, y=0.5$ then we have good agreement for $r=5$ (red points) and $r=20$ (gray points) in laminar flow at $t<0.5$, but a big difference in the turbulent flow at $t>0.5$.

enter image description here

Update 1. To test the Gauss-Seidel relaxation method itself we can compute viscous flow in a plane channel with external force (gravity) as follows

ClearAll["Global`*"]

dif = 1/80; pec = 1; U0 = 0; V0 = 0; n = 11; n1 = n + 1; dt = 
 2./n^2; sm = 3000; r = 7; n2 = Round[n/2]; a = dt dif n n; c = 
 ConstantArray[0, {n1, n1}]; c0 = ConstantArray[0, {n1, n1}];; u0 = 
 ConstantArray[0, {n1, n1}]; v0 = ConstantArray[0, {n1, n1}]; u = 
 ConstantArray[0, {n1, n1}]; v = ConstantArray[0, {n1, n1}];

periodic[n_, up_, ud_, ub_] := 
  Module[{bd = ub}, 
   Do[bd[[1, i]] = .5 (bd[[n, i]] + bd[[2, i]]); 
    bd[[n + 1, i]] = bd[[1, i]]; bd[[i, 1]] = ud; 
    bd[[i, n + 1]] = up;, {i, 2, n}]; bd[[1, 1]] = ud; 
   bd[[n + 1, n + 1]] = up; bd[[n + 1, 1]] = ud; bd[[1, n + 1]] = up; 
   bd];

diffuse[n_, r_, a_, c_, c0_] := 
 Module[{c1 = c}, c1 = ConstantArray[0, {n + 1, n + 1}]; 
  Do[Do[Do[c1[[i, 
         j]] = (c0[[i, j]] + 
           a (c1[[i - 1, j]] + c1[[i + 1, j]] + c1[[i, j - 1]] + 
              c1[[i, j + 1]]))/(1 + 4 a);, {j, 2, n}];, {i, 2, n}]; 
   c1 = periodic[n, 0, 0, c1];, {k, 0, r}]; c1]; 
advect[n_, b_, d_, d0_, u_, v_, dt_] := 
 Module[{x, y}, d1 = ConstantArray[0, {n + 1, n + 1}]; dt0 = dt n; 
  Do[Do[x = i - dt0 u[[i, j]]; y = j - dt0 v[[i, j]]; 
     i0 = Which[x <= 1, 1, 1 < x < n, Floor[x], x >= n, n]; 
     i1 = i0 + 1; 
     j0 = Which[y <= 1, 1, 1 < y < n, Floor[y], y >= n, n];
     j1 = j0 + 1; s1 = x - i0; s0 = 1 - s1; t1 = y - j0; t0 = 1 - t1; 
     d1[[i, j]] = 
      s0 (t0 d0[[i0, j0]] + t1 d0[[i0, j1]]) + 
       s1 (t0 d0[[i1, j0]] + t1 d0[[i1, j1]]);, {j, 1, n + 1}];, {i, 
    1, n + 1}]; d1]; 
project[n_, r_, u0_, v0_, u_, v_] := 
  Module[{ux = u, vy = v}, p = ConstantArray[0, {n1, n1}]; 
   div = ConstantArray[0, {n1, n1}]; ux = ConstantArray[0, {n1, n1}]; 
   vy = ConstantArray[0, {n1, n1}]; 
   Do[div[[i, 
       j]] = -.5 /
        n (u0[[i + 1, j]] - u0[[i - 1, j]] + v0[[i, 1 + j]] - 
         v0[[i, j - 1]]);, {i, 2, n}, {j, 2, n}]; 
   Do[Do[Do[
       p[[i, j]] = (div[[i, 
             j]] + (p[[i - 1, j]] + p[[i + 1, j]] + p[[i, j - 1]] + 
              p[[i, j + 1]]))/4;, {j, 2, n}], {i, 2, n}];, {k, 0, r}];
    Do[ux[[i, j]] = u0[[i, j]] - .5 n (p[[i + 1, j]] - p[[i - 1, j]]);
     vy[[i, j]] = 
     v0[[i, j]] - .5 n (p[[i, j + 1]] - p[[i, j - 1]]);, {i, 2, 
     n}, {j, 2, n}]; {ux, vy}];   
(*force*)
Fx[t_, x_, y_] := 1/10; Fy[t_, x_, y_] := 0;
f1 = ConstantArray[0, {n1, n1}]; f2 = ConstantArray[0, {n1, n1}];

Do[u0 = advect[n, 0, c, u0, u0, v0, dt]; 
  v0 = advect[n, 0, c, v0, u0, v0, dt];
  u0 = periodic[n, 0, 0, u0]; v0 = periodic[n, 0, 0, v0];
  u0 = diffuse[n, r, a, c, u0]; v0 = diffuse[n, r, a, c, v0];
  u0 = periodic[n, 0, 0, u0]; v0 = periodic[n, 0, 0, v0];
  Do[f1[[i, j]] = Fx[dt (step + .5), (i - 1)/n, (j - 1)/n]; 
   f2[[i, j]] = Fy[dt (step + .5), (i - 1)/n, (j - 1)/n];, {i, 2, 
    n1}, {j, 2, n1}]; u0 += f1 dt; 
  v0 += f2 dt; {u0, v0} = project[n, r, u0, v0, u, v]; 
  u0 = periodic[n, 0, 0, u0]; v0 = periodic[n, 0, 0, v0]; 
  uu[step] = u0; vv[step] = v0;, {step, 0, sm}];

Visualization of flow velocity and absolute error

Do[lstu[k] = 
   Flatten[Table[{{(i - 1)/n, (j - 1)/n}, uu[k][[i, j]]}, {i, n1}, {j,
       n1}], 1]; 
  lstv[k] = 
   Flatten[Table[{{(i - 1)/n, (j - 1)/n}, vv[k][[i, j]]}, {i, n1}, {j,
       n1}], 1];, {k, 0, sm}];
Do[Uvel[i] = Interpolation[lstu[i], InterpolationOrder -> 3];, {i, 1, 
  sm}]
Do[Vvel[i] = Interpolation[lstv[i], InterpolationOrder -> 3];, {i, 1, 
   sm}];

{StreamDensityPlot[{Uvel[sm][x, y], Vvel[sm][x, y]}, {x, 0, 1}, {y, 0,
    1}, ColorFunction -> "RoseColors", Frame -> False, 
  PlotLegends -> Automatic], 
 Plot3D[Uvel[sm][x, y], {x, 0, 1}, {y, 0, 1}, 
  ColorFunction -> "RoseColors"]}

{ListPlot[Table[{i dt, Max[uu[i]]}, {i, sm}], 
  AxesLabel -> {"t", "Umax"}], 
 Plot[{-4 (-y + y^2), Uvel[sm][.5, y]}, {y, 0, 1}, 
  PlotStyle -> {Thick, {Red, Dashed}}, Frame -> True, 
  FrameLabel -> {"y", "u"}, 
  PlotLegends -> {"Exact solution", "Numeric solution"}], 
 Plot[Abs[-Uvel[sm][.5, y] - 4 (-y + y^2)], {y, 0, 1}, Frame -> True, 
  FrameLabel -> {"y", "Error"}]}

enter image description here

Actually for r=7 we have a good convergence to the exact solution on the grid of 12 points with error of $10^{−2}$. For r=5 the error is about 0.2 and for r=10 error is about $10^{−2}$, therefore r=7 is optimal value.

Update 2. In the module advect we can use trilinear interpolation as follows

advect[n_, d0_, u_, v_, w_, dt_] := 
  Module[{x, y, z, d1, dt0, i0, i1, j0, j1, k0, k1, s0, s1, t0, t1, 
    p1, p0, d00, d10, d01, d11, cd0, cd1, xd, yd, zd}, 
   d1 = ConstantArray[0, {n + 1, n + 1, n + 1}]; dt0 = dt n;
   Do[Do[Do[x = i - dt0 u[[i, j, k]]; y = j - dt0 v[[i, j, k]]; 
        z = k - dt0 w[[i, j, k]];
        i0 = Which[x <= 1, 1, 1 < x < n, Floor[x], True, n];
        i1 = i0 + 1;
        j0 = Which[y <= 1, 1, 1 < y < n, Floor[y], True, n];
        j1 = j0 + 1; 
        k0 = Which[z <= 1, 1, 1 < z < n, Floor[z], True, n];
        k1 = k0 + 1;(*Trilinear interpolation*)xd = x - i0; 
        yd = y - j0; zd = z - k0;
        d00 = d0[[i0, j0, k0]] (1 - xd) + d0[[i1, j0, k0]] xd; 
        d01 = d0[[i0, j0, k1]] (1 - xd) + d0[[i1, j0, k1]] xd; 
        d10 = d0[[i0, j1, k0]] (1 - xd) + d0[[i1, j1, k0]] xd; 
        d11 = d0[[i0, j1, k1]] (1 - xd) + d0[[i1, j1, k1]] xd;
        cd0 = d00 (1 - yd) + d10 yd; cd1 = d01 (1 - yd) + d11 yd; 
        d1[[i, j, k]] = cd0 (1 - zd) + cd1 zd;
        , {k, 2, n}];, {j, 2, n}];, {i, 1, n + 1}]; d1];

Update: testing projection step

To test projection step itself we can use Mathematica FEM and exact benchmark solution from well known paper EXACT FULLY 3D NAVIER-STOKES SOLUTIONS FOR BENCHMARKING by C. ROSS ETHIER AND D. A. STEINMAN as follows. First, we consider projection step as an implementation of the predictor-corrector algorithm $$\frac{u-u_n}{\tau}+(u.\nabla)u-\nu\nabla^2 u=0$$ $$\frac{u_{n+1}-u}{\tau}+\nabla p =0$$ here $\tau$ is time step, $u_n, u, u_{n+1}$ if velocity field on previous, intermediate and next step consequently, and $p$ is a pressure. We suppose that $\nabla.u_{n+1}=0$, and therefore $$\nabla^2p-\frac{\nabla.u}{\tau}=0$$ Second, we solve NSE in the unit cuboid with Dirichlet condition using exact solution and FEM, we have code

Clear["Global`*"]
Needs["NDSolve`FEM`"]
reg = Cuboid[];   
mesh = ToElementMesh[reg, MaxCellMeasure -> 0.001];
(*Exact solution*)
U[x_, y_, z_, 
  t_] := -a Exp[-d^2 t] (Exp[a x] Sin[a y + d z] + 
    Exp[a z] Cos[a x + d y]); 
V[x_, y_, z_, 
  t_] := -a Exp[-d^2 t] (Exp[a y] Sin[a z + d x] + 
    Exp[a x] Cos[a y + d z]); 
W[x_, y_, z_, 
  t_] := -a Exp[-d^2 t] (Exp[a z] Sin[a x + d y] + 
    Exp[a y] Cos[a z + d x]);
P[x_, y_, z_, 
  t_] := -a ^2/
   2 Exp[-2 d^2 t] (Exp[2 a x] + Exp[2 a y] + Exp[2 a z] + 
    2 Sin[a x + d y] Exp[a (y + z)] Cos[a z + d x] + 
    2 Sin[a y + d z] Exp[a (x + z)] Cos[a x + d y] + 
    2 Sin[a z + d x] Exp[a (y + x)] Cos[
      a y + d z]); 
(*t0 is time step, nn is number of iterations *)
a = 1; d = 1; t0 = 1/400; nn = 200; \[Nu] = 1;
(*FEM implementation of predictor-corrector algorithm*)
UX[0] = U[x, y, z, 0];
VY[0] = V[x, y, z, 0]; WZ[0] = W[x, y, z, 0];
P0[0] = P[x, y, z, 0];
Do[
   {UX[i], VY[i], WZ[i], P0[i]} = 
     NDSolveValue[{{-\[Nu]*
           Laplacian[
            u[x, y, z], {x, y, z}] + (u[x, y, z] - UX[i - 1])/t0 + 
          UX[i - 1]*D[u[x, y, z], x] + VY[i - 1]*D[u[x, y, z], y] + 
          WZ[i - 1]*D[u[x, y, z], z], -\[Nu]*
           Laplacian[
            v[x, y, z], {x, y, z}] + (v[x, y, z] - VY[i - 1])/t0 + 
          UX[i - 1]*D[v[x, y, z], x] + VY[i - 1]*D[v[x, y, z], y] + 
          WZ[i - 1]*D[v[x, y, z], z], -\[Nu]*
           Laplacian[
            w[x, y, z], {x, y, z}] + (w[x, y, z] - WZ[i - 1])/t0 + 
          UX[i - 1]*D[w[x, y, z], x] + VY[i - 1]*D[w[x, y, z], y] + 
          WZ[i - 1]*D[w[x, y, z], z], 
         Laplacian[p[x, y, z], {x, y, z}] - (
\!\(\*SuperscriptBox[\(u\), 
TagBox[
RowBox[{"(", 
RowBox[{"1", ",", "0", ",", "0"}], ")"}],
Derivative],
MultilineFunction->None]\)[x, y, z] + 
\!\(\*SuperscriptBox[\(v\), 
TagBox[
RowBox[{"(", 
RowBox[{"0", ",", "1", ",", "0"}], ")"}],
Derivative],
MultilineFunction->None]\)[x, y, z] + 
\!\(\*SuperscriptBox[\(w\), 
TagBox[
RowBox[{"(", 
RowBox[{"0", ",", "0", ",", "1"}], ")"}],
Derivative],
MultilineFunction->None]\)[x, y, z])/t0} == {0, 0, 0, 0}, {
        DirichletCondition[{u[x, y, z] == 
           U[x, y, z, i t0 - t0/2] + t0 D[P[x, y, z, i t0], x], 
          v[x, y, z] == 
           V[x, y, z, i t0 - t0/2] + t0 D[P[x, y, z, i t0], y], 
          w[x, y, z] == 
           W[x, y, z, i t0 - t0/2] + t0 D[P[x, y, z, i t0], z], 
          p[x, y, z] == P[x, y, z, i t0 ]}, True]}}, {u[x, y, z] - t0 
\!\(\*SuperscriptBox[\(p\), 
TagBox[
RowBox[{"(", 
RowBox[{"1", ",", "0", ",", "0"}], ")"}],
Derivative],
MultilineFunction->None]\)[x, y, z], v[x, y, z] - t0 
\!\(\*SuperscriptBox[\(p\), 
TagBox[
RowBox[{"(", 
RowBox[{"0", ",", "1", ",", "0"}], ")"}],
Derivative],
MultilineFunction->None]\)[x, y, z], w[x, y, z] - t0 
\!\(\*SuperscriptBox[\(p\), 
TagBox[
RowBox[{"(", 
RowBox[{"0", ",", "0", ",", "1"}], ")"}],
Derivative],
MultilineFunction->None]\)[x, y, z], 
       p[x, y, z]}, {x, y, z} \[Element] mesh, 
      Method -> {"FiniteElement", 
        "InterpolationOrder" -> {u -> 2, v -> 2, w -> 2, 
          p -> 2}}];, {i, 1, nn}]; 

Visualization of error on every 5 steps on the line $x=1/2, z=1/2$ for $u_n=(UX[n],VY[n],WZ[n])$, and $p=P0[n]$

Table[{k t0 // N, 
  Plot[Evaluate[{U[x, y, z, k t0]/UX[k] - 1} /. {x -> 1/2, 
      z -> 1/2}], {y, 0, 1}, PlotRange -> All], 
  Plot[Evaluate[{V[x, y, z, k t0]/VY[k] - 1} /. {x -> 1/2, 
      z -> 1/2}], {y, 0, 1}, PlotRange -> All], 
  Plot[Evaluate[{W[x, y, z, k t0]/WZ[k] - 1} /. {x -> 1/2, 
      z -> 1/2}], {y, 0, 1}, PlotRange -> All], 
  Plot[Evaluate[{P[x, y, z, k t0]/P0[k] - 1} /. {x -> 1/2, 
      z -> 1/2}], {y, 0, 1}, PlotRange -> All]}, {k, 5, nn, 5}]

Please note, that in the picture shown last steps only. The relative error for pressure is about $3\times 10^{-3}$, and for velocity $1.25\times10^{-3}$ on the grid $10\times10\times 10$ on 200 steps in time with time step $\tau =1/400$. Figure 1

Using exact benchmark solution we also can test linear FEM algorithm described on Solver for unsteady flow with the use of Mathematica FEM and extended to 3D below as follows

Needs["NDSolve`FEM`"]
reg = Cuboid[];
mesh = ToElementMesh[reg, MaxCellMeasure -> 0.001]

U[x_, y_, z_, 
  t_] := -a Exp[-d^2 t] (Exp[a x] Sin[a y + d z] + 
    Exp[a z] Cos[a x + d y]); 
V[x_, y_, z_, 
  t_] := -a Exp[-d^2 t] (Exp[a y] Sin[a z + d x] + 
    Exp[a x] Cos[a y + d z]); 
W[x_, y_, z_, 
  t_] := -a Exp[-d^2 t] (Exp[a z] Sin[a x + d y] + 
    Exp[a y] Cos[a z + d x]);
P[x_, y_, z_, 
  t_] := -a ^2/
   2 Exp[-2 d^2 t] (Exp[2 a x] + Exp[2 a y] + Exp[2 a z] + 
    2 Sin[a x + d y] Exp[a (y + z)] Cos[a z + d x] + 
    2 Sin[a y + d z] Exp[a (x + z)] Cos[a x + d y] + 
    2 Sin[a z + d x] Exp[a (y + x)] Cos[
      a y + d z]); a = 1; d = 1; t0 = 1/400; nn = 200; \[Nu] = 1;
UX[0] = U[x, y, z, 0];
VY[0] = V[x, y, z, 0]; WZ[0] = W[x, y, z, 0];
P0[0] = P[x, y, z, 0];
Do[
   {UX[i], VY[i], WZ[i], P0[i]} = 
     NDSolveValue[{{-\[Nu]*
           Laplacian[
            u[x, y, z], {x, y, z}] + (u[x, y, z] - UX[i - 1])/t0 + 
          UX[i - 1]*D[u[x, y, z], x] + VY[i - 1]*D[u[x, y, z], y] + 
          WZ[i - 1]*D[u[x, y, z], z] + 
\!\(\*SuperscriptBox[\(p\), 
TagBox[
RowBox[{"(", 
RowBox[{"1", ",", "0", ",", "0"}], ")"}],
Derivative],
MultilineFunction->None]\)[x, y, z], -\[Nu]*
           Laplacian[
            v[x, y, z], {x, y, z}] + (v[x, y, z] - VY[i - 1])/t0 + 
          UX[i - 1]*D[v[x, y, z], x] + VY[i - 1]*D[v[x, y, z], y] + 
          WZ[i - 1]*D[v[x, y, z], z] + 
\!\(\*SuperscriptBox[\(p\), 
TagBox[
RowBox[{"(", 
RowBox[{"0", ",", "1", ",", "0"}], ")"}],
Derivative],
MultilineFunction->None]\)[x, y, z], -\[Nu]*
           Laplacian[
            w[x, y, z], {x, y, z}] + (w[x, y, z] - WZ[i - 1])/t0 + 
          UX[i - 1]*D[w[x, y, z], x] + VY[i - 1]*D[w[x, y, z], y] + 
          WZ[i - 1]*D[w[x, y, z], z] + 
\!\(\*SuperscriptBox[\(p\), 
TagBox[
RowBox[{"(", 
RowBox[{"0", ",", "0", ",", "1"}], ")"}],
Derivative],
MultilineFunction->None]\)[x, y, z], (
\!\(\*SuperscriptBox[\(u\), 
TagBox[
RowBox[{"(", 
RowBox[{"1", ",", "0", ",", "0"}], ")"}],
Derivative],
MultilineFunction->None]\)[x, y, z] + 
\!\(\*SuperscriptBox[\(v\), 
TagBox[
RowBox[{"(", 
RowBox[{"0", ",", "1", ",", "0"}], ")"}],
Derivative],
MultilineFunction->None]\)[x, y, z] + 
\!\(\*SuperscriptBox[\(w\), 
TagBox[
RowBox[{"(", 
RowBox[{"0", ",", "0", ",", "1"}], ")"}],
Derivative],
MultilineFunction->None]\)[x, y, z])} == {0, 0, 0, 0}, {
        DirichletCondition[{u[x, y, z] == U[x, y, z, i t0], 
          v[x, y, z] == V[x, y, z, i t0], 
          w[x, y, z] == W[x, y, z, i t0], 
          p[x, y, z] == P[x, y, z, i t0]}, True]}}, {u[x, y, z] , 
       v[x, y, z] , w[x, y, z] , p[x, y, z]}, {x, y, z} \[Element] 
       mesh, Method -> {"FiniteElement", 
        "InterpolationOrder" -> {u -> 2, v -> 2, w -> 2, 
          p -> 1}}];, {i, 1, nn}]; 

Visualization of error on every 5 steps on the line $y=1/2,z=1/2$ for $u_n=(UX[n],VY[n],WZ[n])$, and $p=P0[n]$

Table[{k t0 // N, 
  Plot[Evaluate[{U[x, y, z, k t0]/UX[k] - 1} /. {y -> 1/2, 
      z -> 1/2}], {x, 0, 1}, PlotRange -> All], 
  Plot[Evaluate[{V[x, y, z, k t0]/VY[k] - 1} /. {y -> 1/2, 
      z -> 1/2}], {x, 0, 1}, PlotRange -> All], 
  Plot[Evaluate[{W[x, y, z, k t0]/WZ[k] - 1} /. {y -> 1/2, 
      z -> 1/2}], {x, 0, 1}, PlotRange -> All], 
  Plot[Evaluate[{P[x, y, z, k t0]/P0[k] - 1} /. {y -> 1/2, 
      z -> 1/2}], {x, 0, 1}, PlotRange -> All]}, {k, 5, nn, 5}]

Note, that in the picture shown last steps only. The relative error for pressure is about $2×10^{−3}$, and for velocity $1.5×10^{−4}$ on the grid 10×10×10 on 200 steps in time with time step τ=1/400. Figure 1

4 Replies

Congratulations! Your post was highlighted on the Wolfram's official social media channels. Thank you for your contribution. We are looking forward to your future posts.

POSTED BY: Moderation Team

Using exact benchmark solution we also can test linear FEM algorithm described on Solver for unsteady flow with the use of Mathematica FEM and extended to 3D below as follows

Needs["NDSolve`FEM`"]
reg = Cuboid[];
mesh = ToElementMesh[reg, MaxCellMeasure -> 0.001]

U[x_, y_, z_, 
  t_] := -a Exp[-d^2 t] (Exp[a x] Sin[a y + d z] + 
    Exp[a z] Cos[a x + d y]); 
V[x_, y_, z_, 
  t_] := -a Exp[-d^2 t] (Exp[a y] Sin[a z + d x] + 
    Exp[a x] Cos[a y + d z]); 
W[x_, y_, z_, 
  t_] := -a Exp[-d^2 t] (Exp[a z] Sin[a x + d y] + 
    Exp[a y] Cos[a z + d x]);
P[x_, y_, z_, 
  t_] := -a ^2/
   2 Exp[-2 d^2 t] (Exp[2 a x] + Exp[2 a y] + Exp[2 a z] + 
    2 Sin[a x + d y] Exp[a (y + z)] Cos[a z + d x] + 
    2 Sin[a y + d z] Exp[a (x + z)] Cos[a x + d y] + 
    2 Sin[a z + d x] Exp[a (y + x)] Cos[
      a y + d z]); a = 1; d = 1; t0 = 1/400; nn = 200; \[Nu] = 1;
UX[0] = U[x, y, z, 0];
VY[0] = V[x, y, z, 0]; WZ[0] = W[x, y, z, 0];
P0[0] = P[x, y, z, 0];
Do[
   {UX[i], VY[i], WZ[i], P0[i]} = 
     NDSolveValue[{{-\[Nu]*
           Laplacian[
            u[x, y, z], {x, y, z}] + (u[x, y, z] - UX[i - 1])/t0 + 
          UX[i - 1]*D[u[x, y, z], x] + VY[i - 1]*D[u[x, y, z], y] + 
          WZ[i - 1]*D[u[x, y, z], z] + 
\!\(\*SuperscriptBox[\(p\), 
TagBox[
RowBox[{"(", 
RowBox[{"1", ",", "0", ",", "0"}], ")"}],
Derivative],
MultilineFunction->None]\)[x, y, z], -\[Nu]*
           Laplacian[
            v[x, y, z], {x, y, z}] + (v[x, y, z] - VY[i - 1])/t0 + 
          UX[i - 1]*D[v[x, y, z], x] + VY[i - 1]*D[v[x, y, z], y] + 
          WZ[i - 1]*D[v[x, y, z], z] + 
\!\(\*SuperscriptBox[\(p\), 
TagBox[
RowBox[{"(", 
RowBox[{"0", ",", "1", ",", "0"}], ")"}],
Derivative],
MultilineFunction->None]\)[x, y, z], -\[Nu]*
           Laplacian[
            w[x, y, z], {x, y, z}] + (w[x, y, z] - WZ[i - 1])/t0 + 
          UX[i - 1]*D[w[x, y, z], x] + VY[i - 1]*D[w[x, y, z], y] + 
          WZ[i - 1]*D[w[x, y, z], z] + 
\!\(\*SuperscriptBox[\(p\), 
TagBox[
RowBox[{"(", 
RowBox[{"0", ",", "0", ",", "1"}], ")"}],
Derivative],
MultilineFunction->None]\)[x, y, z], (
\!\(\*SuperscriptBox[\(u\), 
TagBox[
RowBox[{"(", 
RowBox[{"1", ",", "0", ",", "0"}], ")"}],
Derivative],
MultilineFunction->None]\)[x, y, z] + 
\!\(\*SuperscriptBox[\(v\), 
TagBox[
RowBox[{"(", 
RowBox[{"0", ",", "1", ",", "0"}], ")"}],
Derivative],
MultilineFunction->None]\)[x, y, z] + 
\!\(\*SuperscriptBox[\(w\), 
TagBox[
RowBox[{"(", 
RowBox[{"0", ",", "0", ",", "1"}], ")"}],
Derivative],
MultilineFunction->None]\)[x, y, z])} == {0, 0, 0, 0}, {
        DirichletCondition[{u[x, y, z] == U[x, y, z, i t0], 
          v[x, y, z] == V[x, y, z, i t0], 
          w[x, y, z] == W[x, y, z, i t0], 
          p[x, y, z] == P[x, y, z, i t0]}, True]}}, {u[x, y, z] , 
       v[x, y, z] , w[x, y, z] , p[x, y, z]}, {x, y, z} \[Element] 
       mesh, Method -> {"FiniteElement", 
        "InterpolationOrder" -> {u -> 2, v -> 2, w -> 2, 
          p -> 1}}];, {i, 1, nn}]; 

Visualization of error on every 5 steps on the line $y=1/2,z=1/2$ for $u_n=(UX[n],VY[n],WZ[n])$, and $p=P0[n]$

Table[{k t0 // N, 
  Plot[Evaluate[{U[x, y, z, k t0]/UX[k] - 1} /. {y -> 1/2, 
      z -> 1/2}], {x, 0, 1}, PlotRange -> All], 
  Plot[Evaluate[{V[x, y, z, k t0]/VY[k] - 1} /. {y -> 1/2, 
      z -> 1/2}], {x, 0, 1}, PlotRange -> All], 
  Plot[Evaluate[{W[x, y, z, k t0]/WZ[k] - 1} /. {y -> 1/2, 
      z -> 1/2}], {x, 0, 1}, PlotRange -> All], 
  Plot[Evaluate[{P[x, y, z, k t0]/P0[k] - 1} /. {y -> 1/2, 
      z -> 1/2}], {x, 0, 1}, PlotRange -> All]}, {k, 5, nn, 5}]

Note, that in the picture shown last steps only. The relative error for pressure is about $2×10^{−3}$, and for velocity $1.5×10^{−4}$ on the grid 10×10×10 on 200 steps in time with time step τ=1/400. Figure 1

To test projection step itself we can use Mathematica FEM and exact benchmark solution from well known paper EXACT FULLY 3D NAVIER-STOKES SOLUTIONS FOR BENCHMARKING by C. ROSS ETHIER AND D. A. STEINMAN as follows. First, we consider projection step as an implementation of the predictor-corrector algorithm $$\frac{u-u_n}{\tau}+(u.\nabla)u-\nu\nabla^2 u=0$$ $$\frac{u_{n+1}-u}{\tau}+\nabla p =0$$ here $\tau$ is time step, $u_n, u, u_{n+1}$ if velocity field on previous, intermediate and next step consequently, and $p$ is a pressure. We suppose that $\nabla.u_{n+1}=0$, and therefore $$\nabla^2p-\frac{\nabla.u}{\tau}=0$$ Second, we solve NSE in the unit cuboid with Dirichlet condition using exact solution and FEM, we have code

Clear["Global`*"]
Needs["NDSolve`FEM`"]
reg = Cuboid[];   
mesh = ToElementMesh[reg, MaxCellMeasure -> 0.001];
(*Exact solution*)
U[x_, y_, z_, 
  t_] := -a Exp[-d^2 t] (Exp[a x] Sin[a y + d z] + 
    Exp[a z] Cos[a x + d y]); 
V[x_, y_, z_, 
  t_] := -a Exp[-d^2 t] (Exp[a y] Sin[a z + d x] + 
    Exp[a x] Cos[a y + d z]); 
W[x_, y_, z_, 
  t_] := -a Exp[-d^2 t] (Exp[a z] Sin[a x + d y] + 
    Exp[a y] Cos[a z + d x]);
P[x_, y_, z_, 
  t_] := -a ^2/
   2 Exp[-2 d^2 t] (Exp[2 a x] + Exp[2 a y] + Exp[2 a z] + 
    2 Sin[a x + d y] Exp[a (y + z)] Cos[a z + d x] + 
    2 Sin[a y + d z] Exp[a (x + z)] Cos[a x + d y] + 
    2 Sin[a z + d x] Exp[a (y + x)] Cos[
      a y + d z]); 
(*t0 is time step, nn is number of iterations *)
a = 1; d = 1; t0 = 1/400; nn = 200; \[Nu] = 1;
(*FEM implementation of predictor-corrector algorithm*)
UX[0] = U[x, y, z, 0];
VY[0] = V[x, y, z, 0]; WZ[0] = W[x, y, z, 0];
P0[0] = P[x, y, z, 0];
Do[
   {UX[i], VY[i], WZ[i], P0[i]} = 
     NDSolveValue[{{-\[Nu]*
           Laplacian[
            u[x, y, z], {x, y, z}] + (u[x, y, z] - UX[i - 1])/t0 + 
          UX[i - 1]*D[u[x, y, z], x] + VY[i - 1]*D[u[x, y, z], y] + 
          WZ[i - 1]*D[u[x, y, z], z], -\[Nu]*
           Laplacian[
            v[x, y, z], {x, y, z}] + (v[x, y, z] - VY[i - 1])/t0 + 
          UX[i - 1]*D[v[x, y, z], x] + VY[i - 1]*D[v[x, y, z], y] + 
          WZ[i - 1]*D[v[x, y, z], z], -\[Nu]*
           Laplacian[
            w[x, y, z], {x, y, z}] + (w[x, y, z] - WZ[i - 1])/t0 + 
          UX[i - 1]*D[w[x, y, z], x] + VY[i - 1]*D[w[x, y, z], y] + 
          WZ[i - 1]*D[w[x, y, z], z], 
         Laplacian[p[x, y, z], {x, y, z}] - (
\!\(\*SuperscriptBox[\(u\), 
TagBox[
RowBox[{"(", 
RowBox[{"1", ",", "0", ",", "0"}], ")"}],
Derivative],
MultilineFunction->None]\)[x, y, z] + 
\!\(\*SuperscriptBox[\(v\), 
TagBox[
RowBox[{"(", 
RowBox[{"0", ",", "1", ",", "0"}], ")"}],
Derivative],
MultilineFunction->None]\)[x, y, z] + 
\!\(\*SuperscriptBox[\(w\), 
TagBox[
RowBox[{"(", 
RowBox[{"0", ",", "0", ",", "1"}], ")"}],
Derivative],
MultilineFunction->None]\)[x, y, z])/t0} == {0, 0, 0, 0}, {
        DirichletCondition[{u[x, y, z] == 
           U[x, y, z, i t0 - t0/2] + t0 D[P[x, y, z, i t0], x], 
          v[x, y, z] == 
           V[x, y, z, i t0 - t0/2] + t0 D[P[x, y, z, i t0], y], 
          w[x, y, z] == 
           W[x, y, z, i t0 - t0/2] + t0 D[P[x, y, z, i t0], z], 
          p[x, y, z] == P[x, y, z, i t0 ]}, True]}}, {u[x, y, z] - t0 
\!\(\*SuperscriptBox[\(p\), 
TagBox[
RowBox[{"(", 
RowBox[{"1", ",", "0", ",", "0"}], ")"}],
Derivative],
MultilineFunction->None]\)[x, y, z], v[x, y, z] - t0 
\!\(\*SuperscriptBox[\(p\), 
TagBox[
RowBox[{"(", 
RowBox[{"0", ",", "1", ",", "0"}], ")"}],
Derivative],
MultilineFunction->None]\)[x, y, z], w[x, y, z] - t0 
\!\(\*SuperscriptBox[\(p\), 
TagBox[
RowBox[{"(", 
RowBox[{"0", ",", "0", ",", "1"}], ")"}],
Derivative],
MultilineFunction->None]\)[x, y, z], 
       p[x, y, z]}, {x, y, z} \[Element] mesh, 
      Method -> {"FiniteElement", 
        "InterpolationOrder" -> {u -> 2, v -> 2, w -> 2, 
          p -> 2}}];, {i, 1, nn}]; 

Visualization of error on every 5 steps on the line $x=1/2, z=1/2$ for $u_n=(UX[n],VY[n],WZ[n])$, and $p=P0[n]$

Table[{k t0 // N, 
  Plot[Evaluate[{U[x, y, z, k t0]/UX[k] - 1} /. {x -> 1/2, 
      z -> 1/2}], {y, 0, 1}, PlotRange -> All], 
  Plot[Evaluate[{V[x, y, z, k t0]/VY[k] - 1} /. {x -> 1/2, 
      z -> 1/2}], {y, 0, 1}, PlotRange -> All], 
  Plot[Evaluate[{W[x, y, z, k t0]/WZ[k] - 1} /. {x -> 1/2, 
      z -> 1/2}], {y, 0, 1}, PlotRange -> All], 
  Plot[Evaluate[{P[x, y, z, k t0]/P0[k] - 1} /. {x -> 1/2, 
      z -> 1/2}], {y, 0, 1}, PlotRange -> All]}, {k, 5, nn, 5}]

Please note, that in the picture shown last steps only. The relative error for pressure is about $3\times 10^{-3}$, and for velocity $1.25\times10^{-3}$ on the grid $10\times10\times 10$ on 200 steps in time with time step $\tau =1/400$. Figure 1

enter image description here -- you have earned Featured Contributor Badge enter image description here Your exceptional post has been selected for our editorial column Staff Picks http://wolfr.am/StaffPicks and Your Profile is now distinguished by a Featured Contributor Badge and is displayed on the Featured Contributor Board. Thank you!

POSTED BY: Moderation Team
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract