Cantera / cantera

Chemical kinetics, thermodynamics, and transport tool suite
https://cantera.org
Other
580 stars 341 forks source link

Selecting `None` as `transport_model` for a 1D flame gives no error but selects default #1680

Closed Naikless closed 2 weeks ago

Naikless commented 3 months ago

Problem description

When setting the transport_model property of a FlameBase or derived object to None, the underlying Phase object's default transport model is selected without any information about what happened. This seems undesirable, because the user might get the impression they disabled transport, even though this is not the case. In contrast, providing a phase without transport model to the FlameBase constructor gives an error.

Steps to reproduce

import cantera as ct

gas = ct.Solution('h2o2.yaml')
gas.TPX = 300, ct.one_atm, 'H2:1.1, O2:1, AR:5'
f = ct.FreeFlame(gas, width=0.03)
f.transport_model = None
f.transport_model
>>> 'mixture-averaged'

whereas

gas = ct.Solution('h2o2.yaml', transport_model=None)
f = ct.FreeFlame(gas, width=0.03)
>>> CanteraError: 
*******************************************************************************
CanteraError thrown by StFlow::StFlow:
An appropriate transport model
should be set when instantiating the Solution ('gas') object.
*******************************************************************************

Behavior

I suppose, at the moment, an error should be triggered or at least a warning that the default model was selected instead.

Curious question though: Would it even be possible to deactivate transport for the one-dimensional solvers? Of course this doesn't make much sense physically in most cases, because diffusion is basically the underlying process for all these flames. But for educational purposes or if you want to compare with purely convectional models, I could see some potential application.

Lastly, if the selection of a transport model is mandatory, this should probably also be reflected in the documentation of the transport_model property, with ideally a list of the available options.

System information

ischoegl commented 3 months ago

Thanks for reporting. From my perspective, disabling the transport model via

f.transport_model = None

should throw an error and the current behavior would qualify as a bug.

The behavior of selecting the default transport model when instantiating Solution is, however, expected. This follows the definition in h2o2.yaml, where by default the first phase definition is picked which has the specifications:

phases:
- name: ohmech
  thermo: ideal-gas
  elements: [O, H, Ar, N]
  species: [H2, H, O, O2, OH, H2O, HO2, H2O2, AR, N2]
  kinetics: gas
  transport: mixture-averaged
  state: {T: 300.0, P: 1 atm}

So setting transport to mixture-average is the correct behavior.

PS: it is correct that a transport model of None makes no physical sense.

Naikless commented 3 months ago

Just to clarify, I was never concerned with automatically selecting "mixture-averaged" when instantiating the Phase without specifying the transport model.

I only brought this up in the context of setting f.transport_model = None. So I guess the check and error message should then be coming from the corresponding setter in the FlameBase object?

While a flame without diffusion doesn't make much physical sense indeed, there are many occasions where numerical models omit diffusion altogether. So in theory, setting the diffusion coefficients to zero should at least still allow for a mathematically solvable system. However, the way I understand it, disabling transport would also remove any information about the viscosity, which might indeed make the whole thing unsolvable.

ischoegl commented 3 months ago

I only brought this up in the context of setting f.transport_model = None. So I guess the check and error message should then be coming from the corresponding setter in the FlameBase object?

Correct - a PR should be relatively straight-forward (although the error should probably be thrown in the C++ layer).

While a flame without diffusion doesn't make much physical sense indeed, there are many occasions where numerical models omit diffusion altogether. So in theory, setting the diffusion coefficients to zero should at least still allow for a mathematically solvable system. However, the way I understand it, disabling transport would also remove any information about the viscosity, which might indeed make the whole thing unsolvable.

Laminar flame propagation physics depends on upstream diffusion of heat / species, so while you can write down equations, disabling transport isn't something that the model is designed to handle. Within 1D formulations handled by Cantera, viscosity only appears for the radial momentum equation. For theory, please refer to Cantera Science or Kee's textbook.

speth commented 3 months ago

Laminar flame propagation physics depends on upstream diffusion of heat / species, so while you can write down equations, disabling transport isn't something that the model is designed to handle.

Neglecting the diffusive terms is the approximation which leads to the plug flow reactor model. If you tried to solve this problem using the 1D (flame) solver, one problem you'd run into is that the finite difference scheme introduces some numerical diffusion. While this has a relatively small effect compared to real diffusion on a sufficiently well-resolved grid, that won't hold if numerical diffusion is the only diffusive term.

In any case, I agree that the useful fix here is simply to raise an exception if the user tries to set the transport model to None.

Naikless commented 3 months ago

Neglecting the diffusive terms is the approximation which leads to the plug flow reactor model.

Thanks, that was also my assumption. Incidentally, this is also why I stumbled upon this in the first place. I was trying to find the conditions, for which the PFR and 1D approach deliver comparable solutions. This would thus be one of the "educative" examples I had in mind. But I assume adjusting the solver isn't worth it.

I'll see if I can provide a PR for adding an error message.

ischoegl commented 3 months ago

Unless I'm mistaken, a PFR model will never be able to provide a solution for flame speed (one of the canonical problems encountered in 1D flow). An "educative" example would really only be good for something like a burner-stabilized flame or other cases where the mass flux itself is not part of the solution.

I'll see if I can provide a PR for adding an error message.

:tada: great - let us know if you run into any problems!

Naikless commented 3 months ago

While fiddling around with this, I believe I found a more serious version of this bug: If you change the transport model of the phase object after initializing the FlameBase, it crashes the python interpreter:

import cantera as ct
gas = ct.Solution('gri30.yaml')
f = ct.FreeFlame(gas, width=0.1)
gas.transport_model = None

I am not sure how best to deal with this, if changing the transport model of an existing phase should remain allowed.

ischoegl commented 3 months ago

Interesting. I am confirming a hard crash on my machine also.

A crash is certainly something we'd like to avoid. For context: the underlying Solution objects are handled as shared objects, so changing one should be propagated to others. As an aside: there's one open issue that may be somewhat related - #1457. I do, however, believe that it should be possible to handle 1-D cases properly.

PS: I can take a look at this, but likely only after #1670 is merged ... the legacy MATLAB toolbox is still using some implementations that make a clean resolution more difficult

Naikless commented 3 months ago

Ah, I see, so the infrastructure of tracking whether a certain phase object is already linked to existing flame objects is already there.

With this, I should be able to copy the behavior for the transport setter.

ischoegl commented 3 months ago

Ah, I see, so the infrastructure of tracking whether a certain phase object is already linked to existing flame objects is already there.

Yes, there is some infrastructure there, but much of this is only available in the Python API. I believe that this should be handled in C++. Let us know if you need any assistance.

PS: I just created #1686, which implements 'locks' in C++. While it should resolve #1457, it does not address the issue with the transport model that is reported here.