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

Question

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["NDSolve`FEM`"];
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[
       Div][({{-1/a, 0}, {0, -1/a}}.Inactive[Grad][
         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[
       Div][({{-1/a, 0}, {0, -1/a}}.Inactive[Grad][
         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.};
(*nsOp ={{{u[x,y],v[x,y]}}.Inactive[Grad][u[x,y],{x,y}]+(p^(1,0))[x,y]\
, 
     {{u[x,y],v[x,y]}}.Inactive[Grad][v[x,y],{x,y}]+(p^(0,1))[x,y],(u^\
(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.

Answer 1

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.