Avoiding artificial diffusion and minimize changes to code

Question

I am currently working through Solving Partial Differential Equations with Finite Elements specifically the fluid flow problems. I took the Stokes flow problem and replaced it with Euler's equations as well as the Navier-Stokes equations but the Euler equations are causing issues.

Ω = RegionUnion[Rectangle[{0, 0}, {1, 1/2}], Rectangle[{1, 1/10}, {2, 2/5}]];
RegionPlot[Ω, AspectRatio -> Automatic]
eqnsEuler = {ρ (u[x, y] D[u[x, y], x] + v[x, y] D[u[x, y], y]) + 
     D[p[x, y], x] - μ (D[u[x, y], x, x] + 
        D[u[x, y], y, y]), ρ (u[x, y] D[v[x, y], x] + 
        v[x, y] D[v[x, y], y]) + 
     D[p[x, y], y] - μ (D[v[x, y], x, x] + D[v[x, y], y, y]), 
    D[u[x, y], x] + D[v[x, y], y]} /. {μ -> 0, ρ -> 1};
pde = eqnsEuler == {0, 0, 0};
bcs = {DirichletCondition[{u[x, y] == 40.3y*(0.5 - y)/(0.41)^2, 
     v[x, y] == 0.}, x == 0.], 
   DirichletCondition[{v[x, y] == 0.}, 0 < x < 2], 
   DirichletCondition[p[x, y] == 0., x == 2]};
{xVel, yVel, pressure} = 
  NDSolveValue[{pde, bcs}, {u, v, p}, {x, y} ∈ Ω, Method -> {"FiniteElement", 
  "InterpolationOrder" -> {u -> 2, v -> 2, p -> 1}, 
  "MeshOptions" -> {"MaxCellMeasure" -> 0.0005}}];

Which gives a few errors

InitializePDECoefficients::femcscd: The PDE is convection dominated and the result may not be stable. Adding artificial diffusion may help.

FindRoot::stfail: The method AffineCovariantNewton failed to compute the next step.

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.

So referring to a few posts here on Stack Exchange as well as the documentation on Finite Element best practices and ElementMesh creation I added a mesh.

mesh = ToElementMesh[
  RegionUnion[Rectangle[{0, 0}, {1, 1/2}], 
   Rectangle[{1, 1/10}, {2, 2/5}]], 
  "MaxBoundaryCellMeasure" -> 0.0005, "MaxCellMeasure" -> 0.0005]

which I made really fine keeping in mind that it may cause a long running time. That didn't help.

I also tried defining a refinement region and then using MeshRefinementFunction in NDSolveValue with no luck as well.

With minimal changes to the code I've been successful with both Stokes and Navier-Stokes on different regions, including flow around cylinders etc. However with convection driven flow I am experiencing difficulties. The Stabilization of Convection-Dominated Equations in the above links doesn't help as I:

1. Do not want to add artificial diffusion
2. Am not having luck with different meshes (Maybe it's my choice of mesh? Runtime is not much of an issue, I am aiming for minimal changes to original code to compare the methods for the 3 cases)

Are there small changes that will allow the Euler equations to be solvable or do more complex changes need to be made? I am trying to avoid methods similar to what's described here or here

Answer 1

To be realistic, we can make boundary conditions as for inviscid flow and save small $\mu =10^{-3}$ to use FEM and NDSolve as it is, we have

Needs["NDSolve`FEM`"]
mesh = ToElementMesh[
   RegionUnion[Rectangle[{0, 0}, {1, 1/2}], 
    Rectangle[{1, 1/10}, {2, 2/5}]], "MaxBoundaryCellMeasure" -> .01, 
   "MaxCellMeasure" -> 5  10^-4];
mesh["Wireframe"]
eqnsEuler = {[Rho]  (u[x, y]  D[u[x, y], x] + 
        v[x, y]  D[u[x, y], y]) + 
     D[p[x, y], 
      x] - [Mu]  (D[u[x, y], x, x] + 
        D[u[x, y], y, y]), [Rho]  (u[x, y]  D[v[x, y], x] + 
        v[x, y]  D[v[x, y], y]) + 
     D[p[x, y], y] - [Mu]  (D[v[x, y], x, x] + D[v[x, y], y, y]), 
    D[u[x, y], x] + D[v[x, y], y]} /. {[Mu] -> 10^-3, [Rho] -> 1};
pde = eqnsEuler == {0, 0, 0};
bcs = {DirichletCondition[{u[x, y] == 1, v[x, y] == 0.}, x == 0.], 
   DirichletCondition[{v[x, y] == 0.}, 0 < x < 1 || 1 < x < 2], 
   DirichletCondition[{u[x, y] == 0.}, 
    x == 1 && 0 <= y <= 1/10 || 2/5 <= y <= 1/2], 
   DirichletCondition[{p[x, y] == 0., v[x, y] == 0}, x == 2]};
{xVel, yVel, pressure} = 
  NDSolveValue[{pde, bcs}, {u, v, p}, {x, y} [Element] mesh, 
   Method -> {"FiniteElement", 
     "InterpolationOrder" -> {u -> 2, v -> 2, p -> 1}}];

Visualization

Show[DensityPlot[
  Norm[{xVel[x, y], yVel[x, y]}], {x, y} \[Element] mesh, 
  ColorFunction -> "AvocadoColors", AspectRatio -> Automatic, 
  PlotPoints -> 50, PlotLegends -> Automatic], 
 StreamPlot[{xVel[x, y], yVel[x, y]}, {x, y} \[Element] mesh, 
  StreamColorFunction -> Hue]]

It is interesting, that flow is nonsymmetric around axis y=0.25. The reason is not clear. Probably mesh is not symmetric.