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$.
