large Div(V) for incompressible flow around corner

[bigfooted]

The problem:

As discussed in the online SU2 developer meetings, here is the problem:
In the tutorial for incompressible flow with species transport here (which can be converged up to machine precision):
https://su2code.github.io/tutorials/Inc_Species_Transport/
It was noticed that when the flow goes around the corner, the value of the species coming from the top was not conserved. This is due to a nonzero velocity divergence term Div(V). One solution would be to solve the species transport equations and explicitly put Y*Div(V) to zero in the species equations.

However, the root cause is the nonzero Div(V) term coming from the Navier Stokes equations. This term does not seem to be physical. Actually, switching on the MUSCL scheme for the flow equations reduces the maximum divergence by a factor of 2.

My question is what the best way forward would be to reduce Div(V) or its impact on species/turbulence transport without resorting to implementing a pressure based (PISO, SIMPLE) solver (or possibly using much finer meshes around corners).

[EDIT: it is clear now that Div(V) does not have to be, and in fact is not, zero in the incompressible flow cases we are considering]

Any thoughts?

[pcarruscag]

The flow solver does not solve for Div(V) = 0, it solves for Div(F) = 0, where F is the mass flux across an edge.
This flux contains dissipation which is why Div(V) != 0, you cannot avoid it, you need some numerical dissipation or the results will be even more unphysical.
Pressure-velocity coupling strategies will not give you Div(V) = 0, in coupling pressure with velocity there is some dissipation introduced.
If you want a divergence free flux to discretize the scalar equations you need to take what the flow solver computes, see #721 and #726

[bigfooted]

Thanks for the feedback! What I understand is that part of the divergence observed around the corner node is numerical and part is physical, but wouldn't a low dissipative scheme help in reducing the numerical part? Why would that lead to unphysical results?

I'm also wondering how you would conserve the mass fractions when you actually have large velocity divergence in a compressible flow.

Other solvers seem to do a better job in conserving mass fractions around the corner in the presence of nonzero flow divergence, although I mainly have experience with pressure based solvers (where Div(rho*V)=0 is explicitly imposed but as you said, it is not zero in the numerical solution).

Regarding #721 and #726, The solution you propose to is the implementation of terms like this?

  [+]a0 = max(0.0, MassFlux);
  [+]a1 = min(0.0, MassFlux);

[bigfooted]

I ran this case in ANSYS Fluent, staying as close as possible to the SU2 setup. Fluent also shows the nonzero velocity divergence around the corner, but the species field looks perfect. There is some influence when switching from density-based to pressure-based, and between first order and third order, but the overall Div(V) field looks similar to the above image.

[pcarruscag]

"Pressure-based" solvers do not impose Div(rho*V) = 0 explicitly. Unless you are thinking of staggered grids or something like that.
They compute the face fluxes with momentum interpolation and in so doing introduce a dependence of the flux on pressure, aka dissipation.
It is this numerical flux that is divergence free, and therefore the one that should be used for scalars (this is indeed what I did in #726).
For turbulence this does not matter much because the source terms are more important for the overall solution.
But as you guys are finding it, matters for passive scalars. Because if you use a flux that is not divergence free you are effectively introducing a source term.
Yes lower dissipation gives Div(V) closer to 0, which may be why in Fluent you see a different Div(V). In the limit of no dissipation you get checkerboarding and all kinds of numerical crap.

[bigfooted]

OK, I see in the PR that
EdgeMassFluxes[iEdge ] = Res_Conv[0]
is added in the upwind residual. But the variable Res_Conv[0] is now (mostly) obsolete. Is this value now computed by the call to Viscous_Residual? And Res_Conv[0] is now stored in residual[0] (after calling Viscous_Residual)?

[pcarruscag]

You can grab the mass flux from residual[0] after the call to numerics->ComputeResidual(config) (in Centered or UpwindResidual) (at least for the incompressible solver).
The viscous fluxes do not contribute to the continuity equation.

[bigfooted]

OK, I get it now.This is the old result:

This is the new result with the consistent mass fluxes:

[pcarruscag]

Looks reasonable

[bigfooted]

I created WIP PR #1697 so we can all (@Cristopher-Morales @TobiKattmann @danielmayer @EvertBunschoten ) have a look at the code and slowly, but surely, work towards a proper implementation for the species solver.
From what I understood from #721 is that convergence issues were the main reason for not pursuing that PR further?

[pcarruscag]

Yes convergence was slower and it did not seem to make a big difference for turbulent cases.

My theory regarding convergence is that it slowed down because we were no longer using the updated velocities to compute fluxes for the turbulence equations.
To mitigate that problem we can probably do something like:
Flux[iEdge] = solver_mass_flux - rho V.n
Flow solver advances the solution to V*
MassFlux = Flux[iEdge] + rho V*.n
So basically storing a difference between the real flux and the simplified way of computing it, at convergence V* = V and everything is consistent.

[bigfooted]

Some convergence plots for the primitive venturi. The scalar is decoupled from the rest, so only the scalar results should change. The results are minimal. You cannot see it in the figure, but changes are of the order O(0.1) between the residuals for the species.
convergence of the old solver:

convergence of the new solver:

[pcarruscag]

For passive it should not affect convergence, but turbulence is not passive. Do you have a rans problem?

[bigfooted]

Buice & Eaton diffuser, inc. rans, SA

Buice&Eaton diffuser, SA model, using develop:

Buice&Eaton diffuser, SA model, using new model:

Buice&Eaton diffuser, SA model, using new model + div(V) correction of Evert:

Buice&Eaton diffuser, SA model, using old model + div(V) correction of Evert:

[pcarruscag]

Difference between old/new is probably due to the lag in velocity update that I was mentioning.
The old strategy uses new velocities whereas using the massflux of the solver implies using the old velocities. So a block-Jacobi vs block-GS type of difference.

[EvertBunschoten]

I wanted to check whether actually the effects of flow divergence would be negated with the current changes. For this, I set up a 1D, incompressible test case with a constant velocity distribution of 1.0m/s. The initial passive scalar solution was set to 0.5. At the inlet, I set a velocity of 1.1m/s, resulting in a front of nonzero flow divergence to traverse the domain. If the effects of flow divergence would have been negated, the scalar solution would remain constant.

In the figure, you can see the solution after 500 iterations at a CFL of 1.0 (top left). As you can see, the scalar solution shows a bump in the wake of the velocity gradient (right line plot, green curve). I added a set of correction terms according to the scalar advection equation. With these terms enabled, the scalar solution remains constant (bottom left, blue curve), even in the region of nonzero velocity divergence.

[pcarruscag]

I'm not sure I follow, the case is not converged then? Is the domain expanding? Because if it is a rectangular channel, and the flow is incompressible, then the solution is not physical (even for an unsteady flow).

What you are doing with the modification I saw in the PR is effectively a change to the scalar descretization, maybe it is equivalent to a centered scheme?

[EvertBunschoten]

The solution is indeed not converged. Even though the flow solution is non-physical at this time step, that does not mean that the scalar solution has to be non-physical as well. With this modification, monotonicity of the scalar solution is ensured, meaning that the scalar will not change in value, unless nonzero scalar gradients are present. You can see in the image that the green curve returns to the correct value of 0.5 after some time, so the converged solution will be correct. However, if for example the scalar here is the progress variable of an FGM simulation, it could be that the fluctuation of 0.03 is enough to go into the flammability limit, resulting in an artificial spark to form in the domain. This would then result in the wrong solution if the goal is to simulate the cold gas flow.

[pcarruscag]

Ok but then I think that correction should be exposed as a different type of scalar convection (instead of Scalar_Upwind)

[bigfooted]

@EvertBunschoten "without correction" means this new massflux implementation, but excluding your addition?
Can you write down the numerical scheme for the complete "with correction" implementation? That should make it more clear where this new term comes from and what the difference is between your correction that you added in CScalarSolver.inl and the massflux correction that is directly above your addition

[EvertBunschoten]

By "without correction" I indeed mean the upwind scheme with the new mass flux implementation, but without the additional terms I added in my commit. We discussed this yesterday during the developer meeting. The incompressible scalar problem without diffusion is given in equation 1.7. If we want to completely negate the effects of flow divergence on the scalar solution and ensure monotonicity, one arrives at equation 1.10; the scalar advection equation. The second term on the left hand side of 1.10 can be computed with 1.12. Here, the first term on the rhs is the residual term as computed by Scalar_Upwind without the correction term. The second term on the rhs of 1.12 is the correction term I implemented.

During the developer meeting we agreed on enabling these correction terms using a different option for the scalar convective numerical scheme (e.g. SCALAR_CENTERED or SCALAR_ADVECTION instead of SCALAR_UPWIND). I'll also run some turbulent simulations in order to see what effect these correction terms have on the turbulent solution.