# Wall boundary condition in inviscid flows with NDSolve

I am trying to duplicate https://www.wolfram.com/language/12/nonlinear-finite-elements/navier-stokes-equation.html?product=language Naiver-Stokes type simulation for a charge current density. If I am not wrong, I can use inviscid flow. How can I write wall boundary conditions where the normal velocity is zero?

\[CapitalOmega]=RegionUnion[Rectangle[{0.02,0},{0.1,0.0101}],Rectangle[{0,0.0049},{0.02,0.0101}]];

Show[RegionPlot[\[CapitalOmega], Axes -> False, Frame -> None,
ImageSize -> Large],
VectorPlot[Evaluate[{1, 0}], {x, y} \[Element] \[CapitalOmega],
Sequence[VectorScale -> Small, VectorStyle -> Red,
VectorMarkers -> Placed["Arrow", "Start"],
VectorPoints -> Table[{0, y}, {y, 0.005, 0.01, 0.001}]]],
AspectRatio -> Automatic]

inflowBC = DirichletCondition[{u[x, y] == 1, v[x, y] == 0}, x == 0];
outflowBC =
DirichletCondition[{u[x, y] == 0.5, v[x, y] == 0}, x == 0.1];
wallBC = DirichletCondition[{u[x, y] == 0, v[x, y] == 0},
0 < x < 0.1]; (*This is wrong condition. How can I convert to \
inviscid flow?*)
bcs = {inflowBC, outflowBC, wallBC};

{uVel, vVel} = NDSolveValue[{
\!$$\*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,
\!$$\*SuperscriptBox[\(u$$,
TagBox[
RowBox[{"(",
RowBox[{"0", ",", "1"}], ")"}],
Derivative],
MultilineFunction->None]\)[x, y] -
\!$$\*SuperscriptBox[\(v$$,
TagBox[
RowBox[{"(",
RowBox[{"1", ",", "0"}], ")"}],
Derivative],
MultilineFunction->None]\)[x, y] == 0, bcs}, {u,
v}, {x, y} \[Element] \[CapitalOmega],
Method -> {"FiniteElement",
"InterpolationOrder" -> {u -> 2, v -> 2},
"MeshOptions" -> {"MaxCellMeasure" -> 0.00000025}}];

• Hi Emrah. You are in 2D. So in that case, v[x,y] is the "normal" velocity along the horizontal wall. Thus half of your wall BC is correct, namely that v[x,y]==0 along the inlet. I guess you don't want u[x,y]=0 there, which is not right for inviscid flow. I suppose you could try setting u[x,y]=1 there, instead? Let me know if that helps. Commented Apr 19, 2020 at 18:03
• Thanks a lot for the point. Yes, it will be useful to try. Indeed, I also want to try circular shapes. I have to think about more general solution. I know comsol can solve similar systems. The problem does not seem that hard, but may be I am missing some physical points. I will look into potential and stream scalar options, and use Laplacian. Commented Apr 19, 2020 at 19:59
• I had a similar problem recently (in fact, am still working on it). I had good success with both potential and streamfunction formulations (I have been learning as I go along). Note that in potential formulation, you can use Neumann boundary conditions for the velocity normal components at boundaries, where you need them. Commented Apr 19, 2020 at 20:53

When I plotted the OP's self answered question, I saw vertical flow versus the left to right flow that I was expecting.

ContourPlot[sol[x, y], {x, y} ∈ Ω]
VectorPlot[
Evaluate@Grad[sol[x, y], {x, y}], {x, y} ∈ Ω]


As @Paul Harrison alluded to in the comments, that simpler formulation may be to use a NeumannValueflux values at the inlet or outlet boundaries. We will define a NeumannValue flux condition at the inlet and an essential DirichletCondition. Also, we will recast the equations in coefficient form so material flows from high potential to low potential. Here is an example workflow.

Ω =
RegionUnion[Rectangle[{0.02, 0}, {0.1, 0.0101}],
Rectangle[{0, 0.0049}, {0.02, 0.0101}]];
pressurePoint = DirichletCondition[phi[x, y] == 0, x == 0.1];
nv = NeumannValue[1, x == 0];
bcs = {pressurePoint};
sol = NDSolveValue[{Inactive[
Div][{{-1, 0}, {0, -1}}.Inactive[Grad][phi[x, y], {x, y}], {x,
y}] == nv, bcs}, phi, {x, y} ∈ Ω,
Method -> {"FiniteElement",
"MeshOptions" -> {"MaxCellMeasure" -> 0.00000025}}];
vel[x_, y_] := Evaluate@Grad[-sol[x, y], {x, y}]
nvel[x_, y_] := {vel[x, y][[2]], -vel[x, y][[1]]}
Show[ContourPlot[sol[x, y], {x, y} ∈ Ω,
Sequence[PlotRange -> All, Frame -> None, Axes -> None,
Contours -> 40, ColorFunction -> "TemperatureMap",
ContourStyle -> Black]], Graphics[RegionBoundary[Ω]],
AspectRatio -> Automatic, ImageSize -> Medium]
VectorPlot[vel[x, y], {x, y} ∈ Ω,
AspectRatio -> Automatic]


Now, the flow is left to right as I initially expected. We can zoom in on the backward facing step to look at more details of the potential and the flows.

roi = RegionIntersection[
Rectangle[{0.015, 0}, {0.025, 0.01}], Ω];
rmfroi = RegionMember[roi];
bounds = RegionBounds[roi];
vals = sol["ValuesOnGrid"];
mesh = sol["ElementMesh"];
coords = mesh["Coordinates"];
subRange = MinMax@Pick[vals, rmfroi[coords]];
legendBar =
BarLegend[{"TemperatureMap", subRange}, 50, LegendLayout -> "Row",
LegendLabel -> Style["[phi]", Opacity[0.6]],
LegendMarkerSize -> 500, ChartingTickSide -> Left];
options = {PlotRange -> subRange,
ColorFunction -> ColorData[{"TemperatureMap", subRange}],
ContourStyle -> Opacity[0.1], ColorFunctionScaling -> False,
Contours -> 30, AspectRatio -> Automatic,
PlotPoints -> 41, FrameLabel -> {"x", "y"},
PlotLabel -> Style["Potential: phi(x,y)", 18], ImageSize -> 500};
Column[{legendBar,
ContourPlot[sol[x, y], {x, y} ∈ roi, Evaluate[options]]},
Spacings -> 0]
VectorPlot[vel[x, y], {x, y} ∈ roi, AspectRatio -> Automatic]


I got the solution by help of a CFD friend. Using stream function makes much more sense. Some background. http://cfd.mace.manchester.ac.uk/twiki/pub/Main/TimCraftNotes_All_Access/fl3-inviscid.pdf

\[CapitalOmega] =
RegionUnion[Rectangle[{0.02, 0}, {0.1, 0.0101}],
Rectangle[{0, 0.0049}, {0.02, 0.0101}]];
inflowBC =
DirichletCondition[phi[x, y] == 5.82688 (y - 0.0049)/0.0101,
x == 0];
outflowBC = DirichletCondition[phi[x, y] == 3 y/0.0101, x == 0.1];
wallBC1 = DirichletCondition[phi[x, y] == 3, y > 0.007];
wallBC2 = DirichletCondition[phi[x, y] == 0, y < 0.007];
bcs = {inflowBC, outflowBC, wallBC1, wallBC2};
sol = NDSolveValue[{D[phi[x, y], x, x] + D[phi[x, y], y, y] == 0,
bcs}, phi, {x, y} \[Element] \[CapitalOmega]]
`