Cylindrical coordinates in FEM
Mathematica Asked by Oscillon on May 30, 2021
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"]
2 Answers
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.
Correct answer by Tim Laska on May 30, 2021
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]
Answered by user21 on May 30, 2021
Add your own answers!
Ask a Question
Get help from others!
Recent Questions
- How can I transform graph image into a tikzpicture LaTeX code?
- How Do I Get The Ifruit App Off Of Gta 5 / Grand Theft Auto 5
- Iv’e designed a space elevator using a series of lasers. do you know anybody i could submit the designs too that could manufacture the concept and put it to use
- Need help finding a book. Female OP protagonist, magic
- Why is the WWF pending games (“Your turn”) area replaced w/ a column of “Bonus & Reward”gift boxes?
Recent Answers
- haakon.io on Why fry rice before boiling?
- Joshua Engel on Why fry rice before boiling?
- Peter Machado on Why fry rice before boiling?
- Lex on Does Google Analytics track 404 page responses as valid page views?
- Jon Church on Why fry rice before boiling?