Implementing finite-rate chemistry in a CFD code

[Dave Holland]

Hi all,

 

My colleague and I are currently developing a high-speed, compressible CFD solver which we would like to incorporate finite-rate chemistry using Cantera (We have already implemented equilibrium chemistry using Cantera). To do this we would like to use the Strang splitting method described by Fedkiw et al [1].

 

In this paper, the temperature is updated implicitly using a Newton-Raphson scheme, rather than solving an ODE for temperature along with evolving the mass fractions. To test this methodology, we have written a test piece of c++ code to perform a single cell calculation. Here we use Cantera to calculate the net production rates and gas properties and use Sundials to solve the system of ODE’s. Within this piece of code, we also use a Cantera reactor to solve the same problem to give us a result to compare against. We have found that there is a small difference in the final temperature between these two methods. So far, we are unable to explain the difference in these results. If anyone has any ideas on why this would be they would be greatly appreciated.

 

In looking into the correct implementation of this, we have a few questions that we were wondering if anyone here have any views on.

 

In Fedkiw’s paper it appears clear to us that the chemistry should be performed under the assumption of constant volume within our CFD simulations, however, looking at several other papers they choose to use to use a constant pressure assumption without giving much (if any) reasoning why. Does anyone here know or have a view on why one would choose constant pressure over constant volume?

 

When it comes to implementing this in our CFD code we are trying to decide if we should store an IdealGasMix object for each cell with the correct conditions always set. Alternatively, we could have only one IdealGasMix object that we would reset the conditions for whenever we need to do work for that cell. This is really a CPU usage vs memory usage question. How do other people handle this in their codes?

 

Are there any other sticking points we should be looking out for when implementing finite-rate chemistry into our code?

 

Regards,

 

Dave

 

[1] Ronald P. Fedkiw, Barry Merriman, and Stanley Osher, High Accuracy Numerical Methods for Thermally Perfect Gas Flows with Chemistry, JOURNAL OF COMPUTATIONAL PHYSICS 132, 175–190 (1997)

[Ray Speth]

Dave,

It is certainly correct that you should hold either internal energy or enthalpy constant within the chemistry step. However, expressing the equations in terms of the temperature can still be preferable, assuming you are using an ideal gas model. If you are using the internal energy or enthalpy as the state variable, you don't need to do the iteration yourself to set the state of the Cantera objects -- they provide functions which do these iterations internally, e.g. setState_UV for the combination of internal energy and specific volume. These different options essentially cover the range of models implemented in Cantera's reactor network simulations.

As far as the set up for working with a large mesh, I think the most efficient strategy is to use one set of Cantera objects for each thread. The overhead of initializing and storing an object for each grid point would be prohibitive, and I don't think it would actually provide any performance benefit. If you're interested in seeing an operator-split code which uses Cantera for finite-rate chemistry, you can take a look at Ember, which solves time-dependent 1D flames.

Regards,

Ray

[Dave Holland]

Hi Ray,

 

Thank you for your response.

 

Yes, we are looking at ideal gases and I agree that one needs to keep either the energy or the enthalpy constant during the chemistry step. However, in addition to this we must choose whether to keep the volume constant or the pressure constant as the chemistry is evolved. In terms of Cantera reactors this would be a choice between using IdealGasReactor or IdealGasConstPressureReactor.

 

Thank you for pointing me towards your code, Ember. Having scanned this code it appears to me that you are keeping the pressure constant during the chemistry step (I’m looking at the function SourceSystemCVODE::f in the file sourceSystem.cpp). What was the reasoning behind choosing this over a constant volume formulation?

 

Thank you for your advice on using one set of Cantera objects per thread, this is the strategy we will take.

Finally, if you have thoughts, could you indicate whether you can see any problems down the road in terms of calculating temperature implicitly (as per Fedkiw reference I provided) versus what appears to be the more common and would appear unnecessary solution of a specific temperature ODE?

 

Regards,

Dave

[Ray Speth]

Dave,

I think that the best choice between holding internal energy and enthalpy constant depends on which is most representative of what's happening in your system. For most flames (i.e., not detonations) burning at a given location results in volume expansion and occurs at nearly constant pressure even if the volume of the overall system is constrained as in an internal combustion engine, and it makes the most sense to hold the enthalpy constant. In Ember, the choice is even simpler, since the pressure is constant in space and time.

For the chemistry step, it is true that you can skip integration of an ODE for the temperature and calculate the temperature from the composition and the initial value of the enthalpy. However, in practice it can be harder to construct a good Jacobian for the integrator to use, depending on how you handle issues such as the constraint on the sum of the mass fractions. This is part of the reason that Cantera has specialized reactor models for ideal gases, where we write the energy equation with the temperature as the state variable, and we use the more general form of the energy equation only when it is required by the equation of state.

Regards,

Ray

[Matt Z]

Hi Dave,

I am working on modeling detonations using finite rate chemistry with cantera. I currently have been storing a reactor object per cell and for smaller meshes this works fine. This method is simple to implement but the impacts on performance are huge. When increasing the mesh densities to accurately model detonations this method is extremely prohibitive. 

Also, If you are working with high pressures and temperatures it can be beneficial to use high temperature thermodynamic polynomials in your kinetic model files. The NASA9 polynomials seem to work fine and you should be able to pull those from the CEA database. 

Regards 

Matt