Cylindrical coordinates in FEM

Question

I am trying to solve the Stokes equation for fluid flow in a 3d cylinder. All boundaries are no-slip, apart from the top boundary, which enforces flow in the x-direction.

My problem is that I can't enforce periodic boundary conditions at -pi and pi in the azimuthal direction for the pressure. Instead of a solution, I get the errors:

NDSolve: DirichletCondition can not be present on the target boundary of a PeriodicBoundaryConditon. NDSolve: The boundary condition discretization failed.

When I omit the periodic pressure condition, NDSolve finishes, but the solution has a problem around the origin. Also the flows should be mirror-symmetric across the x-axis due to the top boundary condition, but they are not as can be seen in the x-component of the flow field.

I already incorporated the trick of extending the domain in azimuthal direction from here: Solve Laplace equation in Cylindrical - Polar Coordinates. But that did not seem to help.

What can I do to get a good solution out of NDSolve?

Below is a minimal working example.

(** PDE **)
cs = "Cylindrical";
stokesEqns = {
   Simplify[
      Laplacian[{ur[r, \[Phi], z], u\[Phi][r, \[Phi], z], 
        uz[r, \[Phi], z]}, {r, \[Phi], z}, cs]] - 
     Simplify[Grad[pp[r, \[Phi], z], {r, \[Phi], z}, cs]] == {0, 0, 0},
   Simplify[
     Div[{ur[r, \[Phi], z], u\[Phi][r, \[Phi], z], 
       uz[r, \[Phi], z]}, {r, \[Phi], z}, cs]] == 0
   };

(** boundary conditions **)
{u0r, u0\[Phi], u0z} = 
  TransformedField[
    "Cartesian" -> cs, {1, 0, 0}, {xx, yy, zz} -> {r, \[Phi], z}] /. 
   z -> 1;
boundaryConditions = {
   DirichletCondition[{ur[r, \[Phi], z] == u0r, 
     u\[Phi][r, \[Phi], z] == u0\[Phi], uz[r, \[Phi], z] == u0z}, 
    z == 1 \[And] -\[Pi] < \[Phi] < \[Pi]],
   DirichletCondition[{ur[r, \[Phi], z] == 0, 
     u\[Phi][r, \[Phi], z] == 0, uz[r, \[Phi], z] == 0, 
     pp[r, \[Phi], z] == 0}, z == -1 \[And] -\[Pi] < \[Phi] < \[Pi]],
   DirichletCondition[{ur[r, \[Phi], z] == 0, 
     u\[Phi][r, \[Phi], z] == 0, uz[r, \[Phi], z] == 0, 
     pp[r, \[Phi], z] == 0}, r == 1 \[And] -\[Pi] < \[Phi] < \[Pi]],
   PeriodicBoundaryCondition[ur[r, \[Phi], z], \[Phi] == -\[Pi], 
    TranslationTransform[{0, 2 \[Pi], 0}]],
   PeriodicBoundaryCondition[u\[Phi][r, \[Phi], z], \[Phi] == -\[Pi], 
    TranslationTransform[{0, 2 \[Pi], 0}]],
   PeriodicBoundaryCondition[uz[r, \[Phi], z], \[Phi] == -\[Pi], 
    TranslationTransform[{0, 2 \[Pi], 0}]],
   PeriodicBoundaryCondition[pp[r, \[Phi], z], \[Phi] == -\[Pi], 
    TranslationTransform[{0, 2 \[Pi], 0}]]
   };

(** solve **)
AbsoluteTiming[
  solFEM = 
    NDSolve[{stokesEqns, boundaryConditions}, {ur, u\[Phi], uz, 
       pp}, {r, 0, 1}, {\[Phi], -\[Pi], \[Pi] + \[Pi]/4}, {z, -1, 1}, 
      Method -> {"FiniteElement", 
        "InterpolationOrder" -> {ur -> 2, u\[Phi] -> 2, uz -> 2, 
          pp -> 1}}][[1]];
  ][[1]]


(** plot **)
field[xx_, yy_, zz_] = 
  TransformedField[
   cs -> "Cartesian", {ur[r, \[Phi], z], u\[Phi][r, \[Phi], z], 
     uz[r, \[Phi], z]} /. solFEM, {r, \[Phi], z} -> {xx, yy, zz}];
ppCart[xx_, yy_, zz_] = 
  TransformedField[cs -> "Cartesian", 
   pp[r, \[Phi], z] /. solFEM, {r, \[Phi], z} -> {xx, yy, zz}];
DensityPlot3D[
 field[x, y, z][[1]]
 , {x, -1, 1}, {y, -1, 1}, {z, -1, 1}
 , PlotRange -> All, PlotLegends -> Automatic, 
 AxesLabel -> {"x", "y", "z"}, PlotLabel -> "x-component of flow"]

Seen that you transform back to Cartesian coordinates why not just solve the problem in Cartesian coordinates? — user21, Jun 04 '20 at 06:48
The issue with DirichletCondition present at the target can be very subtle. Make sure that the specification for the DirichletCondition predicate and the target do not overlap at any point. Initially try to reduce the position of the DirichletCondition as much as possible (even to a smaller space then what you finally want) and then to get the PBCs working then enlarge to find where the problem is. — user21, Jun 04 '20 at 06:51

Tim Laska · Accepted Answer · 2020-06-04T14:33:41.513

This appears to be a lid driven flow problem. I am in agreement with @user21's perspective that you should solve this in Cartesian Coordinates. It should simplify the boundary condition specification. Since the system is closed, you will need to define pressure at a node. I used OpenCascade to build the half cylinder. Here is the workflow.

(* Load Required Packages *)
Needs["OpenCascadeLink`"]
Needs["NDSolve`FEM`"]
(* Use OpenCascade To Make Half Sym Geometry *)
pp = Polygon[{{0, 0, -1}, {0, 0, 1}, {1, 0, 1}, {1, 0, -1}}];
shape = OpenCascadeShape[pp];
axis = {{0, 0, 0}, {0, 0, 1}};
sweep = OpenCascadeShapeRotationalSweep[shape, axis, -Pi];
(* Create Mesh *)
bmesh = OpenCascadeShapeSurfaceMeshToBoundaryMesh[sweep];
mesh = ToElementMesh[bmesh, MaxCellMeasure -> {"Length" -> .075}, 
   "IncludePoints" -> {{0, 0.5, -1}}];
groups = mesh["BoundaryElementMarkerUnion"];
temp = Most[Range[0, 1, 1/(Length[groups])]];
colors = ColorData["BrightBands"][#] & /@ temp;
mesh["Wireframe"["MeshElementStyle" -> FaceForm /@ colors]]
(* Create PDE System *)
ClearAll[μ]
op = {Inactive[
              Div][({{-μ, 0, 0}, {0, -μ, 0}, {0, 
                    0, -μ}}.Inactive[Grad][
         u[x, y, z], {x, y, z}]), {x, y, 
              z}] + 
     D[p[x, y, z], x], 
        Inactive[
              Div][({{-μ, 0, 0}, {0, -μ, 0}, {0, 
                    0, -μ}}.Inactive[Grad][
         v[x, y, z], {x, y, z}]), {x, y, 
              z}] + 
     D[p[x, y, z], y],
        Inactive[
              Div][({{-μ, 0, 0}, {0, -μ, 0}, {0, 
                    0, -μ}}.Inactive[Grad][
         w[x, y, z], {x, y, z}]), {x, y, 
              z}] + 
     D[p[x, y, z], z], 
    D[u[x, y, z], x] + 
     D[v[x, y, z], y] + 
     D[w[x, y, z], z]} /. μ -> 1;
pde = op == {0, 0, 0, 0};
bcs = {DirichletCondition[
        {u[x, y, z] == 1, v[x, y, z] == 0., w[x, y, z] == 0.}, 
    z == 1.],
      DirichletCondition[
        {u[x, y, z] == 0, v[x, y, z] == 0., w[x, y, z] == 0.}, 
        z == -1. || (x^2 + y^2) > 0.99], 
      DirichletCondition[v[x, y, z] == 0., y > -0.001],
      DirichletCondition[p[x, y, z] == 0., 
        x == 0. && z == -1.](*pressure Point Condition*)};
(* Solve PDE *)
{xVel, yVel, zVel, pressure} = 
    NDSolveValue[{pde, bcs}, {u, v, w, p}, {x, y, z} ∈ mesh, 
      Method -> {"FiniteElement", 
          "InterpolationOrder" -> {u -> 2, v -> 2, w -> 2, p -> 1}}];
(* Visualize Solution *)
surf = {{"YStackedPlanes", {0}}, {"ZStackedPlanes", {-1, 1}}};
Show[SliceContourPlot3D[
    Norm@{xVel[x, y, z], yVel[x, y, z], zVel[x, y, z]}, 
    surf, {x, y, z} ∈ mesh, PlotPoints -> 50, 
    BoxRatios -> Automatic, ColorFunction -> "TemperatureMap"], 
  ImageSize -> Medium, ViewPoint -> Front]
DensityPlot3D[
  Norm[{xVel[x, y, z], yVel[x, y, z], zVel[x, y, z]}], {x, y, 
      z} ∈ mesh, BoxRatios -> Automatic, 
  ColorFunction -> "TemperatureMap", ViewAngle -> 0.3669386546105606`,
  ViewPoint -> {3.7435513617679828`, 1.2106476957796874`, 
   0.9258298223054351`}, 
 ViewVertical -> {0.27079048490259205`, 0.14735018657087556`, 
   0.9512940848148628`}]
SliceVectorPlot3D[{xVel[x, y, z], yVel[x, y, z], 
    zVel[x, y, z]}, surf, {x, y, z} ∈ mesh, 
 VectorPoints -> 20,
   VectorColorFunction -> "BrightBands", BoxRatios -> Automatic, 
 ViewPoint -> Front]

Qualitatively, it agrees with the COMSOL model I threw together.

This is awesome! Could you clarify two points? 1) For the pressure point condition you specfiy x and z-coordinates, but not y. I added it in and it does not seem to change anything. Is it implicitly set to y=0? 2) This is basically a test problem I wanted to understand before continuing with a more complicated two-phase flow problem in spherical coordinates. What I take away from your answer is that with FEM I should stay in cartesian coordinates, b/c curvilinear coordinates introduce messy artifacts. Is that right? — Oscillon, Jun 14 '20 at 23:28
@Oscillon 1) Picking a node with conditional statement is not necessarily a robust way to pick a pressure point, because there is no guarantee the mesher will put a point there. I chose the bottom of the tank because the flow is relatively static (very small pressure gradients), so there is a chance multiple points were grabbed, but had little effect. A more robust approach would be to explicitly define the pressure point in ToElementMesh . 2) I would only use curvilinear coordinates if I could take advantage of symmetry and reduce the dimension of the problem I am solving. Continued- — Tim Laska, Jun 15 '20 at 11:57
@Oscillon If the dimension stays the same, you may as well let the mesher do the heavy lifting of meshing the shape you want to model. Regarding two-phase flow, I hope that you have a defined fluid-fluid interface. Implementing free surfaces in FEM codes is not trivial and the solutions are expensive. — Tim Laska, Jun 15 '20 at 12:04

score 5 · Answer 2 · answered Jun 04 '20 at 07:14

Here is a version in Cartesian coordinates to get you started:

reg = Cylinder[{{0, 0, 0}, {0, 0, 1}}, 1];

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}]};
\[CapitalGamma]D = {
   DirichletCondition[{u[x, y, z] == 1., v[x, y, z] == 0., 
     w[x, y, z] == 0.}, x == 1], 
   DirichletCondition[{u[x, y, z] == 0., v[x, y, z] == 0., 
     w[x, y, z] == 0.}, x < 1], 
   DirichletCondition[p[x, y, z] == 0, x == -1 && y == 0 && z == 1]};

Needs["NDSolve`FEM`"]
mesh = ToElementMesh[reg];

{xVel, yVel, zVel, pressure} = 
  NDSolveValue[{stokesFlowOperator == {0, 0, 0, 
      0}, \[CapitalGamma]D}, {u, v, w, p}, {x, y, z} \[Element] mesh, 
   Method -> {"FiniteElement", 
     "InterpolationOrder" -> {u -> 2, v -> 2, w -> 2, p -> 1}}];

You'd need to think more about the boundary conditions, especially the pressure condition.

rmf = RegionMember[MeshRegion[mesh]];
Quiet[VectorPlot3D[{xVel[x, y, z], yVel[x, y, z], zVel[x, y, z]}, 
  Evaluate[Sequence @@ Join[{{x}, {y}, {z}}, mesh["Bounds"]*1.01, 2]],
   VectorStyle -> "Arrow3D", VectorColorFunction -> "TemperatureMap", 
  VectorScale -> {Tiny, Scaled[0.4], None}, VectorPoints -> {9, 9, 9},
   Axes -> None, Boxed -> False, 
  RegionFunction -> (rmf[{#1, #2, #3}] &)], 
 InterpolatingFunction::femdmval]

Cylindrical coordinates in FEM

2 Answers2