Improved combustion solver robustness for coarse grids

[EvertBunschoten]

Simulations of flames in large domains is difficult because volumetric source terms blow up in large cells, causing the species solver to diverge. I implemented a damping term which scales the species source terms in cells with a length scale above a user-defined threshold. This allows for the flame to propagate through regions of the domain with large cells (such as in the far-field) without diverging, while the source terms in refined regions remain unaffected.

To illustrate how this feature improves robustness, I ran a 2D simulation of a circular flame propagating from a refined region in the middle domain towards the edges where the mesh is significantly more coarse.

The image below shows the temperature field after 30 iterations as computed with the previous implementation (left), and the new implementation with a flame length scale specified as 1e-4 m (right). The solution is unaffected while the flame is within the refined region of the mesh.

However, the instance the flame reaches the coarser cells, the solver quickly diverges without damping. With the damping active, the flame propagates through the entire domain without an issue.

PR Checklist

Put an X by all that apply. You can fill this out after submitting the PR. If you have any questions, don't hesitate to ask! We want to help. These are a guide for you to know what the reviewers will be looking for in your contribution.

• I am submitting my contribution to the develop branch.
• My contribution generates no new compiler warnings (try with --warnlevel=3 when using meson).
• My contribution is commented and consistent with SU2 style (https://su2code.github.io/docs_v7/Style-Guide/).
• I used the pre-commit hook to prevent dirty commits and used pre-commit run --all to format old commits.
• I have added a test case that demonstrates my contribution, if necessary.
• I have updated appropriate documentation (Tutorials, Docs Page, config_template.cpp), if necessary.

[pcarruscag]

Hmmm but then you are not solving the same equations anymore.
If the issue is "propagation speed" it sounds like something should be done about the CFL.

[EvertBunschoten]

Hmmm but then you are not solving the same equations anymore. If the issue is "propagation speed" it sounds like something should be done about the CFL.

It is not about accurately modeling the propagation of the flame, it is just to allow the flame to propagate through the coarse regions of the domain without the solver diverging. For cells larger than the flame length scale, we can't accurately resolve the propagation of the flame front anyway because the diffusive time scales and reaction source time scales differ at those length scales.

[pcarruscag]

So is it possible to adjust the diffusion coefficients in those regions?
Numerical diffusion as a way of dealing with instabilities is a more established practice than disabling source terms.

[bigfooted]

I agree with Pedro, it is great that it doesn't blow up, but the solution is not physical anymore. You can look into 'thickened flame' models. In these models, the propagation speed is physically correct, but the flame is artificially thickened by increasing the diffusion coefficient.

[bigfooted]

It's mentioned in Poinsot & Veynante btw.

[EvertBunschoten]

I was not clear when initially describing the problem. The images I showed were of a steady simulation, whereas the thickened flame models described in literature were developed for unsteady problems.
For unsteady flames, adding the correction factor to only the diffusivity seems to improve the propagation of flames in coarse cells very nicely.
The image below shows the comparison between unsteady simulations of an expanding flame front with (right) and without (left) scaling the diffusivity according to the cell volume. The mesh is refined within the purple circle and is very coarse beyond it.

The flame propagates at the same speed for both cases within the refined region, but changes once the flame goes into the coarse region.

With the artificially thickened flame correction (diffusion only), the flame propagates at approximately the same speed through the coarse cells, whereas the flame speed drops significantly without the correction.

[EvertBunschoten]

For the steady case (constant CFL number=10), it seems that suppressing the source term and enhancing diffusivity in the course cells improves the robustness of the simulations. The image below shows the results of a steady simulation on the same grid as the unsteady simulation without any correction term (left) and with diffusion increased and source terms decreased (right).

At around 80 iterations, the flame front reaches the coarse cells, where the uncorrected simulation quickly diverges.

With the correction terms, the flame front propagates rapidly through the coarse cells, which is expected for a steady simulation, given that larger cells have a larger time step.

[bigfooted]

looks good. Note that the thickened flame model of the paper we discussed https://share.google/Ia0hcPhfoOKvHypb3 recovers the flame speed by also modifying the diffusivity and source term in the energy equation, and therefore modifies the flame temperature.
If you plan to silently do this by default, maybe make the lengthscale depend on the local estimate of the flame thickness (in all regions with non-negligible source terms) and do this when the cell size is > 5 times the flame thickness or so. And start giving warnings to the user on screen about it.

[EvertBunschoten]

looks good. Note that the thickened flame model of the paper we discussed https://share.google/Ia0hcPhfoOKvHypb3 recovers the flame speed by also modifying the diffusivity and source term in the energy equation, and therefore modifies the flame temperature. If you plan to silently do this by default, maybe make the lengthscale depend on the local estimate of the flame thickness (in all regions with non-negligible source terms) and do this when the cell size is > 5 times the flame thickness or so. And start giving warnings to the user on screen about it.

It is possible to calculate the flame thickness during simulations with Equation 4 of the paper you linked. However, the flame speed correction model will not work when the correction factor depends on the local value of the flame thickness.
As the flame propagates through coarse cells, its thickness will increase due to the artificially increased diffusivity. Therefore, the flame will grow in thickness relative to the local cell size, which will reduce the value of the correction factor, leading to a reduction in flame propagation speed.

[bigfooted]

is it possible to compute an estimate of the original flame thickness locally?

[EvertBunschoten]

is it possible to compute an estimate of the original flame thickness locally?

I have tried the following method to try to obtain a measure for the thickness of the flame on a sufficiently refined grid.
Each iteration, the flame thickness is calculated with
$$t=\frac{\max(T) - \min(T)}{\max(\nabla_c T)},$$
where
$$\nabla_c T = \frac{\nabla c \cdot \nabla T}{|\nabla c|}.$$
To account for the thickening of the flame when propagating through coarse cells, the flame thickness value used to calculate the correction factor is only updated if the flame thickness decreases.
$$t_{i+1} = \min(t_i, t)$$
This approach seems to work well for adiabatic flames when the flame starts in a refined section of the domain. However, if there are isothermal or CHT BCs, the flame thickness drastically reduces when the flame front propagates through the thermal boundary layer.
The decrease in flame thickness increases the value of the correction factor, causing the flame front to stagnate or the solver to diverge.
I have also tried to only use the gradient of the progress variable to calculate the flame thickness, but I encounter the same issue.

[bigfooted]

So the issue I have is: how much are we going to accommodate a bad setup (a mesh that is too coarse) for a user?
If it is a matter of adding a 'coarse_mesh_correction = ON' and that's it, then that's fine by me. Preferably when it can be on by default and does not do anything when the mesh is fine enough.
If the adiabatic correction reproduces the correct flamespeed (steady and unsteady) without any other magic constants, I think we should definitely put that in.
Since your testcase is non-adiabatic, more correction is needed. Is there some additional correction possible that takes into account the change in enthalpy? In this heat loss region, can you scale the correction to get it to converge, so have an additional correction when the enthalpy has changed compared to a reference enthalpy? If yes, then maybe we can convert this 'magic scaling factor' into an local_enthalpy/adiabatic_enthalpy ratio or something.

[EvertBunschoten]

So the issue I have is: how much are we going to accommodate a bad setup (a mesh that is too coarse) for a user? If it is a matter of adding a 'coarse_mesh_correction = ON' and that's it, then that's fine by me. Preferably when it can be on by default and does not do anything when the mesh is fine enough. If the adiabatic correction reproduces the correct flamespeed (steady and unsteady) without any other magic constants, I think we should definitely put that in. Since your testcase is non-adiabatic, more correction is needed. Is there some additional correction possible that takes into account the change in enthalpy? In this heat loss region, can you scale the correction to get it to converge, so have an additional correction when the enthalpy has changed compared to a reference enthalpy? If yes, then maybe we can convert this 'magic scaling factor' into an local_enthalpy/adiabatic_enthalpy ratio or something.

The purpose of this feature is not to accommodate a bad setup where the mesh in the flame zone is too coarse to accurately resolve the flame front. For complicated 3D laminar test cases and future support for tubulent flames, the flow field needs to develop before igniting the flame. Consequently, the flame front will propagate through all regions of the computational domain downstream of the flame, including regions close to the outlet where the mesh is less refined.
Modifying the correction factor based on the change in the total enthalpy w.r.t. the initial value is something I'm working on, and by using the suggestion of Josh to have it be dependent on the wall distance.

[bigfooted]

ok, so for this PR, do you just want to implement the steady case, where any unsteady propagation is nonphysical anyway? There is still some kind of length scale involved then?

[pcarruscag]

Perhaps we can move forward with a default that does not compromise results, and it is then an option to set an adequate lenght scale?

[EvertBunschoten]

Perhaps we can move forward with a default that does not compromise results, and it is then an option to set an adequate lenght scale?

I developed a new way to approximate the flame thickness:

$t_\mathrm{flame} = \frac{\max(c) - \min(c)}{\max(|\nabla_{Tu}c|)},$
where
$\nabla_{Tu}c = \frac{\nabla_u c\cdot\nabla T}{\nabla T \cdot \nabla T}\nabla T,$
and
$\nabla_u c = \frac{\nabla c\cdot u}{u\cdot u} u.$
This formulation works better in non-adiabatic regions.

The flame thickness value is now also included in the history output when the option is enabled in the configuration so users are better informed on whether the mesh is sufficiently refined.

Finally, the global flame thickness value is now adjusted upwards or downwards depending on whether the flame thickness increases or decreases. This allows for the flame thickness value to become more accurate as the solution converges, while providing sufficient dampening while the flame propagates through coarse cells.

[bigfooted]

Nice work! Is there also no correction when the mesh is sufficiently fine?

[bigfooted]

@EvertBunschoten can you add a testcase?

[pcarruscag]

It's not a residual

[EvertBunschoten]

I have placed the flame thickness under its own category as a user-defined coefficient.

[pcarruscag]

Is this temperature a tunable parameter?

[EvertBunschoten]

Now it is

[pcarruscag]

This does not work with OpenMP because the local variables are also thread-local, you can take a look at how we compute time steps for the pattern you need to use.

[EvertBunschoten]

Nice work! Is there also no correction when the mesh is sufficiently fine?

No correction is applied if the length scale of the cell is lower than the calculated flame length scale.