# Avoiding artificial diffusion and minimize changes to code

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] == 4*0.3*y*(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

• Are you trying to solve this system of equation using NDSolve with Automatic options? There are standard methods for solving incompressible flows, and these methods are still not programmed into NDSolve. Feb 13 at 19:10
• Euler equation (and in irregular domain)? I'm not that familiar with static case, but if the eventual goal is dealing with time-dependent case, it'll be troublesome, and even problem in regular domain is hard enough. Related: mathematica.stackexchange.com/a/267899/1871 Feb 14 at 8:23
• @AlexTrounev I am aware of the standard methods, I was hoping to achieve it with NDSolve. Any idea why it's not programmed into NDSolve? Feb 14 at 14:16
• @xzczd I am eventually moving to time-dependent cases yes. So it appears the only work around is to add artificial viscosity? Feb 14 at 14:19
• @Kendall To be precise, the easiest work-around. More stable but advanced approaches are, as shown in the link above, turning to the packages in other languages in e.g. julia or python using ExternalEvaluate, or coding our own solver. "Any idea why it's not programmed into NDSolve?" Who knows, perhaps it's just too hard. Feb 14 at 14:54

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["NDSolveFEM"]

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.

• Thank you Alex, it seems I cannot avoid using an artificial diffusion? I've actually come across a few of your answers here on Stack Exchange that have been very helpful. How did you learn what you know? Can you recommend some resources please. Feb 14 at 14:21
• It's much, much easier and much less error prone to just use the operator to generate the equations: vars = {{u[x, y], v[x, y], p[x, y]}, {x, y}}; pars = <|"DynamicViscosity" -> 10^-3, "MassDensity" -> 1|>; eqnsEuler = FluidFlowPDEComponent[vars, pars] You'd get the Euler equations for a DynamicViscosity of 0. But we currently can not solve them, unless you add artificial diffusion. Here are some tips for that. Feb 14 at 15:50
• @Kendall We use artificial diffusion even for Euler equations in a case of compressible flow as well, see for example mathematica.stackexchange.com/questions/287067/… Feb 14 at 15:52
• @AlexTrounev The solution to the nonsymmetric part was bugging me so I played around a bit. It's due to how Mathematica is interpreting the boundary condition as well as the mesh.DirichletCondition[{u[x, y] == 0.}, x == 1 && 0 <= y <= 1/10 || 2/5 <= y <= 1/2] Adding brackets resolves the issue and the solution becomes symmetric. DirichletCondition[{u[x, y] == 0.}, x == 1 && (0 <= y <= 1/10 || 2/5 <= y <= 1/2)] I also decreased the "MaxBoundaryCellMeasure" -> .001, "MaxCellMeasure" -> 10^-4 Feb 16 at 16:15
• @Kendall Thank you very much for your efforts. This update works fine. Feb 16 at 18:33