# Solutions for Stokes flow in 3D

I am trying to extend Stokes flow example to 3D but I get an error. Not sure what's wrong.

For example we define the region and this works:

Ω =
ImplicitRegion[
-0.5 <= z <= 0.5 && 0 <= x <= 2 && 0 <= y <= 0.5 &&
!(x >= 1 && y <= 0.1) && !(x >= 1 && y >= 0.4),
{x, y, z}];

RegionPlot3D[Ω]


Then I define the Stokes flow operator, but I am not sure this is right:

stokesFlowOperator =
{Inactive[Div][{{-1, 0}, {0, -1}}.Inactive[Grad][u[x, y, z], {x, y, z}], {x, y, z}] +
Derivative[1, 0][w][x, y],
Inactive[Div][{{-1, 0}, {0, -1}}.Inactive[Grad][v[x, y, z], {x, y, z}], {x, y, z}] +
Derivative[0, 1][w][x, y, z],
Inactive[Div][{{-1, 0}, {0, -1}}.Inactive[Grad][v[x, y, z], {x, y, z}], {x, y, z}] +
Derivative[0, 1, 1][w][x, y, z],
Derivative[0, 1, 0][v][x, y, z] + Derivative[1, 0, 1][u][x, y, z]};

Subscript[Γ, D] =
{DirichletCondition[u[x, y, z] == z + 4*0.3*y*((0.5 - y)/0.41^2), {x == 0.}],
DirichletCondition[{u[x, y, z] == 0., v[x, y, z] == 0.}, 0 < x < 2],
DirichletCondition[w[x, y, z] == 0., x == 2]}; This lines fails:

{xVel, yVel, zVel, pressure} =
NDSolveValue[
{stokesFlowOperator == {0, 0, 0}, Subscript[Γ, D]},
{u, v, w}, {x, y, z} ∈ Ω,
Method ->
{"FiniteElement",
"InterpolationOrder" -> {u -> 2, v -> 2, w -> 1}}]


And I get this error:

NDSolveValue::dsvar: 0.35 cannot be used as a variable.
Set::shape: Lists {xVel, yVel, zVel, pressure} and NDSolveValue[{False, {DirichletCondition[u[x, y, 0.35] == 0.35 + 7.13861 Plus[<<2>>] y, {x == 0.}], DirichletCondition[{u[x, y, 0.35] == 0., v[x, y, 0.35] ==0.}, 0 < x < 2], DirichletCondition[w[x, y, 0.35] == 0., x == 2]}}, {u, v, w}, {x, y, 0.35} ∈ ImplicitRegion[-0.5 <= z <= 0.5 && 0 <= x <=2 && 0 <= y <= 0.5 && !(x >= 1 && y <= 0.1) && !(x >= 1 && y >= 0.4), {x, y, z}], Method -> {FiniteElement, InterpolationOrder -> {u -> 2, v -> 2, w -> 1}}] are not the same shape.

I am not sure what I am doing wrong when trying to extend to 3D.

These are not the Stokes equations. You must have built in several errors while translating the 2D equations into 3D. The following works although it complains because there are no boundary conditions for the pressure that would make it unique. Also, I don't know at all whether these are the boundary conditions that you would like to apply.

a = IdentityMatrix[3];
stokesFlowOperator = {
Inactive[Div][a.Inactive[Grad][u[x, y, z], {x, y, z}], {x, y, z}] - D[p[x, y, z], x],
Inactive[Div][a.Inactive[Grad][v[x, y, z], {x, y, z}], {x, y, z}] - D[p[x, y, z], y],
Inactive[Div][a.Inactive[Grad][w[x, y, z], {x, y, z}], {x, y, z}] - D[p[x, y, z], z],
Div[{u[x, y, z], v[x, y, z], w[x, y, z]}, {x, y, z}]
};
ΓD = {
DirichletCondition[u[x, y, z] == z + 4*0.3*y*(0.5 - y)/(0.41)^2, x == 0.],
DirichletCondition[u[x, y, z] == 0., 0 < x < 2],
DirichletCondition[v[x, y, z] == 0., 0 < x < 2],
DirichletCondition[w[x, y, z] == 0., x == 2]
};

{xVel, yVel, zVel, pressure} = NDSolveValue[
{
stokesFlowOperator == {0, 0, 0, 0},
ΓD
},
{u, v, w, p},
{x, y, z} ∈ Ω,
Method -> {"FiniteElement",
"InterpolationOrder" -> {u -> 2, v -> 2, w -> 2, p -> 1}}
]

• How can I visualize the returned array of InterpolatingFunction? Assume I added DirichletCondition[p[x, y, z] == 2., 0 < x < 2] to make it unique.
– 0x90
Commented Oct 3, 2018 at 14:22
• The values can be obtained with xVel[[4]];. The corresponding mesh nodes are these: xVel[[3, 1, 1]];. Notice that these nodes can also be located on edge migpoints if interpolation order 2 is used. (The underlying mesh is xVel[[3, 1]], plottable with xVel[[3, 1]]["Wireframe"].) Commented Oct 3, 2018 at 14:30

An example with a Poiseuille profile at the inlet.

H = 1/2; L = 2;
reg = ImplicitRegion[0 <= z <= H && 0 <= x <= L && 0 <= y <= H &&  !(x >= 1 && y <= 0.1) &&  !(x >= 1 && y >= 0.4), {x, y, z}];
RegionPlot3D[reg]
Um = 45/100; nu = 1;
U0[y_, z_] := 16*Um*y*z*(H - y)*((H - z)/H^4)
{U, V, W, P} = NDSolveValue[{{(-nu)*Laplacian[u[x, y, z], {x, y, z}] + Derivative[1, 0, 0][p][x, y, z], (-nu)*Laplacian[v[x, y, z], {x, y, z}] +
Derivative[0, 1, 0][p][x, y, z], (-nu)*Laplacian[w[x, y, z], {x, y, z}] + Derivative[0, 0, 1][p][x, y, z],
Derivative[1, 0, 0][u][x, y, z] + Derivative[0, 1, 0][v][x, y, z] + Derivative[0, 0, 1][w][x, y, z]} == {0, 0, 0, 0},
{DirichletCondition[{u[x, y, z] == U0[y, z], v[x, y, z] == 0, w[x, y, z] == 0}, x == 0], DirichletCondition[{u[x, y, z] == 0, v[x, y, z] == 0, w[x, y, z] == 0},
0 < x < L], DirichletCondition[p[x, y, z] == 0, x == L]}}, {u, v, w, p}, Element[{x, y, z}, reg],
Method -> {"FiniteElement", "InterpolationOrder" -> {u -> 2, v -> 2, w -> 2, p -> 1}, "MeshOptions" -> {"MaxCellMeasure" -> 0.0001}}];
ContourPlot[U[x, y, H/2], {x, 0, L}, {y, 0, H}, AspectRatio -> Automatic, Contours -> 20, PlotPoints -> 50]


• Your answer is great. You may need the browser plug-in in this link to edit the code format. Commented Feb 17, 2020 at 23:20
• @PleaseCorrectGrammarMistakes Thank you. I edited the post using RawInputForm` Commented Feb 17, 2020 at 23:53