# Curved Dirichlet boundary in FEM

Based on this fascinating blog post (https://blog.wolfram.com/2013/07/09/using-mathematica-to-simulate-and-visualize-fluid-flow-in-a-box/) from a few years ago that shows how to simulate a driven cavity flow in a rectangle, I want to simulate the driven cavity flow in a droplet.

With the following Mathematica code:

reynolds = 0.1;
pde = {
-
\!$$\*SuperscriptBox[\(p$$,
TagBox[
RowBox[{"(",
RowBox[{"1", ",", "0"}], ")"}],
Derivative],
MultilineFunction->None]\)[x, y] + 1/reynolds \!$$\*SubsuperscriptBox[\(\[Del]$$, $${x, y}$$, $$2$$]$$u[x, y]$$\),
-
\!$$\*SuperscriptBox[\(p$$,
TagBox[
RowBox[{"(",
RowBox[{"0", ",", "1"}], ")"}],
Derivative],
MultilineFunction->None]\)[x, y] + 1/reynolds \!$$\*SubsuperscriptBox[\(\[Del]$$, $${x, y}$$, $$2$$]$$v[x, y]$$\),

\!$$\*SuperscriptBox[\(u$$,
TagBox[
RowBox[{"(",
RowBox[{"1", ",", "0"}], ")"}],
Derivative],
MultilineFunction->None]\)[x, y] +
\!$$\*SuperscriptBox[\(v$$,
TagBox[
RowBox[{"(",
RowBox[{"0", ",", "1"}], ")"}],
Derivative],
MultilineFunction->None]\)[x, y]
} == {0, 0, 0};

domain = Disk[{0, 0}, 1, {\[Pi], 2 \[Pi]}];

boundaryConditions = {
DirichletCondition[{u[x, y] == 1, v[x, y] == 0, p[x, y] == 0}, y == 0],
DirichletCondition[{u[x, y] == 0, v[x, y] == 0}, x^2 + y^2 == 1 \[And] y <= 0]
};

(** solve flow velocities with a method that is an order higher than \
pressure for stable solutions **)
{xVel, yVel, pressure} = NDSolveValue[
{pde, boundaryConditions},
{u, v, p},
{x, y} \[Element] domain,
Method -> {"FiniteElement",
"InterpolationOrder" -> {u -> 2, v -> 2, p -> 1}}];

(** visualization **)
flowPlot =
VectorPlot[{xVel[x, y], yVel[x, y]}, {x, y} \[Element] domain,
PlotRangePadding -> None, AspectRatio -> Automatic,
StreamPoints -> 6, StreamColorFunction -> "TemperatureMap",
StreamColorFunctionScaling -> False,
Epilog -> Circle[{0, 0}, 1, {\[Pi], 2 \[Pi]}],
FrameLabel -> {"x", "y"}, FrameStyle -> Directive[Black, 18]]


I get a nice result:

Here is what is giving me trouble: In the previous example I had no-slip boundary conditions, meaning that the horizontal and vertical velocity components vanish at the circular boundary. Now I want to see what happens, when the drop swims on a viscous liquid. In this case I don't have no-slip boundary conditions. Instead only the velocity normal to the boundary vanishes and the parallel velocity remains intact.

My strategy to implement the vanishing normal velocity on the circular boundary is this: First I express the normal (=radial) velocity $$v_r$$ as a function of the horizontal velocity $$u$$ and vertical velocity $$v$$:

$$v_r = (u,v) \cdot (\cos\phi,\sin\phi)$$

The polar angle $$\phi$$ can be expressed in $$x$$ and $$y$$ as $$\arctan y/x$$. This results in the radial velocity:

$$v_r = \frac{1}{\sqrt{1+(y/x)^2}}u + \frac{y/x}{\sqrt{1+(y/x)^2}}v$$

The problem arises when I try to enforce $$v_r$$ to be zero with a Dirichlet boundary condition:

reynolds = 0.1;
pde = {
-
\!$$\*SuperscriptBox[\(p$$,
TagBox[
RowBox[{"(",
RowBox[{"1", ",", "0"}], ")"}],
Derivative],
MultilineFunction->None]\)[x, y] + 1/reynolds \!$$\*SubsuperscriptBox[\(\[Del]$$, $${x, y}$$, $$2$$]$$u[x, y]$$\) == 0,
-
\!$$\*SuperscriptBox[\(p$$,
TagBox[
RowBox[{"(",
RowBox[{"0", ",", "1"}], ")"}],
Derivative],
MultilineFunction->None]\)[x, y] + 1/reynolds \!$$\*SubsuperscriptBox[\(\[Del]$$, $${x, y}$$, $$2$$]$$v[x, y]$$\) == 0,

\!$$\*SuperscriptBox[\(u$$,
TagBox[
RowBox[{"(",
RowBox[{"1", ",", "0"}], ")"}],
Derivative],
MultilineFunction->None]\)[x, y] +
\!$$\*SuperscriptBox[\(v$$,
TagBox[
RowBox[{"(",
RowBox[{"0", ",", "1"}], ")"}],
Derivative],
MultilineFunction->None]\)[x, y] == 0
};

domain = Disk[{0, 0}, 1, {\[Pi], 2 \[Pi]}];

boundaryConditions = {
DirichletCondition[{u[x, y] == 1, v[x, y] == 0, p[x, y] == 0}, y == 0],
DirichletCondition[{(1/Sqrt[1 + (y/x)^2]) u[x, y] + (y/x/Sqrt[1 + (y/x)^2]) v[x, y] == 0},  x^2 + y^2 == 1 \[And] y <= 0]
};

(** solve flow velocities with a method that is an order higher than \
pressure for stable solutions **)
{xVel, yVel, pressure} = NDSolveValue[
{pde, boundaryConditions},
{u, v, p},
{x, y} \[Element] domain,
Method -> {"FiniteElement",
"InterpolationOrder" -> {u -> 2, v -> 2, p -> 1}}];


Instead of a solution, I get this error: NDSolveValue::fembdcc: Cross-coupling of dependent variables in DirichletCondition[u/Sqrt[1+y^2/x^2]+(v y)/(x Sqrt[1+y^2/x^2])==0,x^2+y^2==1&&y<=0] is not supported in this version.

Does this mean I can not solve fluid dynamics in a system with curved boundaries other than no-slip?

(Transforming into polar coordinates might work, but this is a simplified toy problem for illustration. In my actual problem, I can not transform into polar coordinates.)

You could workaround the slip BC limitation by modeling a thin lubricated layer. I cannot take credit for the following code and my initial search did not find it in the Wolfram Function Repository.

The following code will create a green liquid drop with a yellow thin lubricated annular region:

Needs["NDSolveFEM"]
(* Code to Join Multiple Boundary Meshes *)
ClearAll[validInputQ]
validInputQ[bm1_, bm2_] :=
BoundaryElementMeshQ[bm1] &&
BoundaryElementMeshQ[
bm2] && (bm1["EmbeddingDimension"] ===
bm2["EmbeddingDimension"]) && (bm1["MeshOrder"] ===
bm2["MeshOrder"] === 1)
BoundaryElementMeshJoin[bm1_, bm2_,
opts : OptionsPattern[ToBoundaryMesh]] /; validInputQ[bm1, bm2] :=
Module[{c1, c2, nc1, newBCEle, newPEle, eleTypes, markers},
c1 = bm1["Coordinates"];
c2 = bm2["Coordinates"];
nc1 = Length[c1];
newBCEle = bm2["BoundaryElements"];
If[ElementMarkersQ[newBCEle], markers = ElementMarkers[newBCEle],
markers = Sequence[]];
newBCEle =
MapThread[#1[##2] &, {eleTypes, ElementIncidents[newBCEle] + nc1,
markers}];
newPEle = bm2["PointElements"];
If[ElementMarkersQ[newPEle], markers = ElementMarkers[newPEle],
markers = Sequence[]];
newPEle =
MapThread[#1[##2] &, {eleTypes, ElementIncidents[newPEle] + nc1,
markers}];
ToBoundaryMesh["Coordinates" -> Join[c1, c2],
"BoundaryElements" -> Flatten[{bm1["BoundaryElements"], newBCEle}],
"PointElements" -> Flatten[{bm1["PointElements"], newPEle}], opts]]
BoundaryElementMeshJoin[r1_, r2_, r3__] :=
BoundaryElementMeshJoin[BoundaryElementMeshJoin[r1, r2], r3];
(* Model Specific Code *)
(* Set up region associations for easy assignment later *)
regs = <|"liquid" -> 10, "air" -> 20|>;
(* Geometry *)
r1 = 1;
thick = r1/20;
rt = r1 + thick;
liquiddomain = Disk[{0, 0}, r1, {\[Pi], 2 \[Pi]}];
symbolicRegion = Disk[{0, 0}, rt, {\[Pi], 2 \[Pi]}];
interface = symbolicRegion~RegionDifference~liquiddomain;
nr = ToNumericalRegion[symbolicRegion];
symbolicBounds = RegionBounds[symbolicRegion];
(bm1 = ToBoundaryMesh[liquiddomain, symbolicBounds])["Wireframe"]
(bm2 = ToBoundaryMesh[interface, symbolicBounds])["Wireframe"]
(bm = BoundaryElementMeshJoin[bm1, bm2])["Wireframe"]
SetNumericalRegionElementMesh[nr, bm];
meshTriangle = ToElementMesh[nr,
"RegionMarker" -> {{{0, -(r1 + thick/2)}, regs["air"], 0.00025},
{{0, -r1/2}, regs["liquid"], 0.012}}
];
meshTriangle[
"Wireframe"[
"MeshElementStyle" -> {FaceForm[Green], FaceForm[Yellow]},
ImageSize -> Large]]


The workaround would be to keep the wall BCs and use a region dependent Reynolds Number as I show here:

reynolds =
Evaluate[Piecewise[{{0.1, ElementMarker == regs["liquid"]}, {20,
True}}]];
pde = {
-

\!$$\*SuperscriptBox[\(p$$,
TagBox[
RowBox[{"(",
RowBox[{"1", ",", "0"}], ")"}],
Derivative],
MultilineFunction->None]\)[x, y] + 1/reynolds \!$$\*SubsuperscriptBox[\(\[Del]$$, $${x, y}$$, $$2$$]$$u[x, y]$$\),
-

\!$$\*SuperscriptBox[\(p$$,
TagBox[
RowBox[{"(",
RowBox[{"0", ",", "1"}], ")"}],
Derivative],
MultilineFunction->None]\)[x, y] + 1/reynolds \!$$\*SubsuperscriptBox[\(\[Del]$$, $${x, y}$$, $$2$$]$$v[x, y]$$\),

\!$$\*SuperscriptBox[\(u$$,
TagBox[
RowBox[{"(",
RowBox[{"1", ",", "0"}], ")"}],
Derivative],
MultilineFunction->None]\)[x, y] +

\!$$\*SuperscriptBox[\(v$$,
TagBox[
RowBox[{"(",
RowBox[{"0", ",", "1"}], ")"}],
Derivative],
MultilineFunction->None]\)[x, y]
} == {0, 0, 0};

boundaryConditions = {
DirichletCondition[{u[x, y] == 1, v[x, y] == 0, p[x, y] == 0},
y == 0],
DirichletCondition[{u[x, y] == 0, v[x, y] == 0},
x^2 + y^2 == rt^2 \[And] y <= 0]
};

(** solve flow velocities with a method that is an order higher than \
\

pressure for stable solutions **)
{xVel, yVel, pressure} = NDSolveValue[
{pde, boundaryConditions},
{u, v, p},
{x, y} \[Element] meshTriangle,
Method -> {"FiniteElement",
"InterpolationOrder" -> {u -> 2, v -> 2, p -> 1}}];

(** visualization **)
flowPlotSlip =
VectorPlot[{xVel[x, y], yVel[x, y]}, {x, y} \[Element] meshTriangle,
PlotRangePadding -> None, AspectRatio -> Automatic,
StreamPoints -> 6, StreamColorFunction -> "TemperatureMap",
StreamColorFunctionScaling -> False,
Epilog -> Circle[{0, 0}, 1, {\[Pi], 2 \[Pi]}],
FrameLabel -> {{y, ""}, {x, "Slip Layer"}},
FrameStyle -> Directive[Black, 18], ImageSize -> Large]


The code is far from optimized, but the slip layer does appear to give non-zero velocities at the interface. I was able to decrease the thickness to 1/100, refine the mesh, and it still solves reasonable fast.

## Modeling Slip

In commercial solvers, sometimes special boundary conditions are implemented by the use of shell elements behind the scenes. Since we do not have access to all the internals of the Mathematica solver, I am explicitly modeling a thin "shell" with actual elements to try to approximate the slip. Rest assured that there are modeling errors in this approach, but also note that the slip boundary condition essentially assumes that gas drag on a liquid surface is zero. Since viscous forces increase with decreasing diameter, this may not be a good assumption at small scales.

To understand the initial modeling error, I used a FEM fluid solver that supports a slip BC on a curved surface and compared the results to a modeled lubricated layer as shown below: For a first guess at an interface layer, the center line velocity magnitude has error, but it looks like something that we could iterate on and improve. Using the features of the curve, we could play with parameters (e.g., Interface Thickness, Mesh Resolution, etc.) to see if we either shifted the vortex center to the right or increased the velocity magnitude at $$R$$. A plot along the center line of the Mathematica solution looks similar to the other commercial code.

In my experience, if a feature requires additional development in a commercial code before I can use it, then it is usually months to years away. In that case, I will look for a work around that gets me to an inferentially valid model.

## Error Reduction with Commercial Codes

I was able to reduce the error significantly in two commercial FEM CFD codes (Altair's AcuSolve and COMSOL) by lowering the viscosity of the lubricated layer and capping the ends of the interface layer with a wall. I did not have the same success with Mathematica. My current hypothesis is that I used boundary layer meshing (flat elements parallel to the curved wall as shown below) with both codes.

The lubricated layer approach has potential, but may need the correct meshing strategy to get the require accuracy in Mathematica.

## Reynolds Number as a Proxy for Inverse Viscosity

Since the equations have been nondimensionalized, I am using the $$Re$$ as a proxy for inverse viscosity. If we keep density, diameter, and characteristic velocity constant, then the following should hold.

$$Re = \frac{{\rho VD}}{\mu } \Rightarrow \mu \sim \frac{1}{{Re }}$$

So, a large $$Re$$ would imply a small viscosity. Pressure coupled FEM fluid solvers can tolerate large viscosity changes, but large density changes make them unstable and that is why I kept density constant.

• Could you elaborate a little more on the background of 1) how the intermediate layer equals/approximates a Dirichlet boundary condition and 2) what is the idea behind using different Reynolds numbers? – Oscillon Nov 4 at 19:38
• I updated the post with some additional elaboration. The Reynolds number is a proxy for inverse viscosity in the nondimensionalized form of the PDE. FEM fluid solvers are pretty tolerant to large viscosity changes so I am creating thin layer that cannot transfer much momentum to the wall. – Tim Laska Nov 5 at 13:33