# Problem Solving Steady-State Euler Equations in 2D Using FEM

I am working with gas flows, and my Reynolds number is very large (>10^6). Still as a toy exercise, I solve the steady-state Navier Stokes equations for my volume and it seems to work ok for Reynolds number <~100. My code:

Needs["NDSolveFEM"];
h[x_] := 0.5 (Erf[0.64 x/0.5642] + 1);
rng = 12;
goal = 5;
mcm = 0.016;
solnRegn = ImplicitRegion[{y > h[x]}, {{x, -rng, rng}, {y, 0, rng}}];
nsOp = {Inactive[
u[x, y], {x, y}]), {x,
y}] + {{u[x, y], v[x, y]}}.Inactive[Grad][u[x, y], {x, y}] +
\!$$\*SuperscriptBox[\(p$$,
TagBox[
RowBox[{"(",
RowBox[{"1", ",", "0"}], ")"}],
Derivative],
MultilineFunction->None]\)[x, y],
Inactive[
v[x, y], {x, y}]), {x,
y}] + {{u[x, y], v[x, y]}}.Inactive[Grad][v[x, y], {x, y}] +
\!$$\*SuperscriptBox[\(p$$,
TagBox[
RowBox[{"(",
RowBox[{"0", ",", "1"}], ")"}],
Derivative],
MultilineFunction->None]\)[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]} /. {a -> 200.};
,
(1,0))[x,y]+(v^(0,1))[x,y]};*)
pde = nsOp == {0, 0, 0};
y0 = 10^-4;
bcs = {DirichletCondition[u[x, y] == 0, y == h[x]],
DirichletCondition[v[x, y] == 0, y == h[x]],
DirichletCondition[u[x, y] == 0.1 Log[1 + (y - h[x])/y0],
x == -rng],
DirichletCondition[u[x, y] == 0.1 Log[1 + (y - h[x])/y0],
x == rng],
DirichletCondition[u[x, y] == 0.1 Log[1 + (y - h[x])/y0],
y == rng], DirichletCondition[v[x, y] == 0, x == -rng],
DirichletCondition[v[x, y] == 0, x == rng],
DirichletCondition[v[x, y] == 0, y == rng],
DirichletCondition[p[x, y] == 0., x == rng]};
refinementRegion =
ImplicitRegion[{y > h[x], y < h[x] + 2.}, {{x, -3, 6}, {y, 0, 3}}];
mrf = With[{rmf = RegionMember[refinementRegion]},
Function[{vertices, area},
Block[{x, y}, {x, y} = Mean[vertices];
If[rmf[{x, y}], area > 0.001, area > 0.02]]]];
soln = NDSolveValue[{pde, bcs}, {u, v, p}, {x, y} \[Element]
solnRegn,
Method -> {"FiniteElement",
"InterpolationOrder" -> {u -> 2, v -> 2, p -> 1},
"MeshOptions" -> {{"MaxCellMeasure" -> mcm},
AccuracyGoal -> goal, PrecisionGoal -> goal,
MeshRefinementFunction -> mrf}}];
Print[soln];


I get into trouble if I try to increase my Reynolds number >~100, receiving the following error message:

I kind-of understand that as Re increases, the viscosity term in the NS equations gets smaller and at some point, it will be too small compared with numerical approximations, and will fail. I guess this is what happens here at Re>~100 so I can never get anywhere near 10^6. I understand that this is not the way to proceed.

So, instead, I want to try to solve for just the "free flow" away from the boundary (and I will look to some specific turbulence/boundary models near the boundary as a later step). Thus, I completely turn-off the viscosity term by setting viscosity=0, and this leads me to the steady-state Euler equations in 2D (I plan to go to 3D once I get things working well in 2D).

So, I simply change the differential operator to the Euler one:

eulOp = {{{u[x, y], v[x, y]}}.Inactive[Grad][u[x, y], {x, y}] +
\!$$\*SuperscriptBox[\(p$$,
TagBox[
RowBox[{"(",
RowBox[{"1", ",", "0"}], ")"}],
Derivative],
MultilineFunction->None]\)[x, y],
{{u[x, y], v[x, y]}}.Inactive[Grad][v[x, y], {x, y}] +
\!$$\*SuperscriptBox[\(p$$,
TagBox[
RowBox[{"(",
RowBox[{"0", ",", "1"}], ")"}],
Derivative],
MultilineFunction->None]\)[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]};
pde = eulOp == {0, 0, 0};


I change nothing else and run again.

The code fails with the following two errors.

FindRoot::nosol: Linear equation encountered that has no solution.

FindRoot::sszero: The step size in the search has become less than the tolerance prescribed by the PrecisionGoal option, but the function value is still greater than the tolerance prescribed by the AccuracyGoal option.

I can't see what the problem is. I tried a few different boundary conditions, including simple constants, but I always get the above errors. Could anyone help with this issue? Thanks in advance.

• Is something like shock wave supposed to appear in the solution? If so, turning to Euler equations only makes the problem harder, I'm afraid. Related: mathematica.stackexchange.com/q/11748/1871 Mar 4, 2020 at 1:25
• I am not sure you BCs are correct. You have a pressure 0 at rng and also set u and v on that edge. Also using y>h[x] in solnRegn is probably an issue. Use y>=h[x]. Typically it is possible to get a Reynolds number of ~1000+ for FEM when the setup is right. You can typically reach a bit more if you use and low Reynolds number as an InitialSeeding. Anything beyond that will probably need a different model, like a k-espilon model (never done that). Having said that, solving the Euler equations with FEM is not a good idea since those are convection dominant. You'd be better of .... Mar 4, 2020 at 7:54
• .... with a different method like Discontinous Galerking. But again, I have not done that myself so I can not guarantee that these will work. Start by looking at your BCs once more, there needs to be a place were the fluid can leave the domain. Mar 4, 2020 at 7:56
• @xzczd. Thanks for your comment. I hope there is no shock wave involved here. This is supposed to be a rather gentle physical process. I followed the link to the post you mentioned. It has variable density, but I am working in the incompressible regime (v<<c => rho=const). Mar 4, 2020 at 10:21
• @user21. Thanks for your suggestions. I will start by adjusting some BCs to see if I can get it to work. But note that these same BCs worked fine for the NS equations with largish Re, so physically, the Euler situation is similar. E.g. the fluid was able to escape fine in that case along the edge at x=rng. I have no understanding of why convection dominance would make FEM an inappropriate method. Mar 4, 2020 at 10:26

I am giving an answer to my own question, in case it is of some help to others who find themselves with a similar problem. I am new to this area of physics (fluid dynamics), so have been learning as I go along. Therefore, if I say something wrong, I hope someone will correct me about this.

First, to recall that in my OP, I tried solving two different partial differential equations, the Navier-Stokes ones at moderate Reynolds number and the (simpler) Euler ones, both 2D steady-state incompressible versions. The one I really wanted to solve was the Euler problem, so I will concentrate on that in this answer.

Having tried most of the suggestions posted as comments, nothing I tried worked. So, I began to ask myself if I was somehow trying to do something "impossible" or mathematically "ill-conditioned". This is what I have concluded for the following reasons.

For the Euler Equations, in particular, the 2D, steady-state incompressible case is a special case. It is a case in which it is possible to define a "stream-function". (For the experts, my case is not irrotational - the vorticity is non-zero from my boundary conditions - so that one cannot define a velocity potential). The point of this is that both the u[x,y] and v[x,y] are partial derivatives of psi, and are thus not independent. My three Euler equations (two independent momentum equations plus the continuity equation) are therefore too many for the number of degrees of freedom, and the problem would appear to be over-determined. This is what I think was wrong.

I am hoping to find an equation satisfied by the stream function, so that I can find it. Then the complete solution will follow easily from the definition of the velocity components, u[x,y] and v[x,y], as its derivatives, and finally I can find the pressure from just one of the momentum equations (or equivalently from the Bernoulli equation).

While I am not certain that this fully explains what was wrong, it seems to fit the facts. I still have some work to do, to understand exactly what equation the stream function satisfies.

• Is solving using Laplace's equation any use to you? Irrotational fluid flow reduces to Laplace's equation. This certainly works in Mathematica see here.
– Hugh
Mar 6, 2020 at 11:18
• @Hugh. Thanks for your suggestion. That was the problem I solved first and it worked fine, but is not sufficiently realistic for my needs. The issue here is that my fluid is not irrotational. You can see this from my boundary conditions, where the vorticity is non-zero. In fact, the stream-function here satisfies a more general equation, Poisson's equation, but the source term on the RHS is the vorticity itself, which I know only on the boundaries, not in general. Mar 6, 2020 at 14:56
• You seem to be on top of the issues. If there is no viscosity then the vorticity will "swim" with the irrotational flow. As a suggestion you could solve for the irrotational flow and then let the vorticity be advected by the irrotational flow. You still have to model how the vorticity interacts with itself but if you are lucky it will have left your domain after a brief time so you don't have to work too hard.
– Hugh
Mar 6, 2020 at 18:34