Solutions for Stokes flow in 3D

Question

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.

Answer 1

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}}
  ]

Answer 2

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]