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

Improve this question
$\endgroup$

2 Answers 2

Reset to default
14
$\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$
2
  • 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
    Commented Oct 3, 2018 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
    Commented Oct 3, 2018 at 14:30
11
$\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$
2
  • $\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
    Commented Feb 17, 2020 at 23:20
  • $\begingroup$ @PleaseCorrectGrammarMistakes Thank you. I edited the post using RawInputForm $\endgroup$
    – Alex Trounev
    Commented Feb 17, 2020 at 23:53

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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