10
$\begingroup$

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.

$\endgroup$
12
$\begingroup$

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}}
  ]
| improve this answer | |
$\endgroup$
  • 1
    $\begingroup$ 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. $\endgroup$ – 0x90 Oct 3 '18 at 14:22
  • $\begingroup$ 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"].) $\endgroup$ – Henrik Schumacher Oct 3 '18 at 14:30
10
$\begingroup$

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]  

fig1

| improve this answer | |
$\endgroup$
  • $\begingroup$ Your answer is great. You may need the browser plug-in in this link to edit the code format. $\endgroup$ – A little mouse on the pampas Feb 17 at 23:20
  • $\begingroup$ @PleaseCorrectGrammarMistakes Thank you. I edited the post using RawInputForm $\endgroup$ – Alex Trounev Feb 17 at 23:53

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.