aka majordome.engineering
Module warnings can be controlled throught the following attributes:
toggle_reactor_warnings(**kwargs : Any): -> None:Reverse truth value of warning flags.
toggle_non_key_value : bool | None = Nonetoggle_missing_species_name : bool | None = Nonetoggle_unknown_species : bool | None = NoneLet’s start by creating a standard Cantera solution:
solution = ct.Solution("airish.yaml")
solution.TPY = 273.15, ct.one_atm, "N2: 0.78 O2: 0.21, AR: 0.01"Conversion of composition strings with name filtering is available:
composition_to_dict(
Y : str,
species_names : list[str] = []
) -> dict[str, float]:Convert a Cantera composition string to dictionary.
Y : strspecies_names : list[str] = []composition_to_dict("O2: 1, teste: 1", solution.species_names){'O2': 1.0}It is also possible to set unit composition with species names; if you do not provide a validation list, all are are kept:
composition_to_dict("O2, N2, hello"){'O2': 1.0, 'N2': 1.0, 'hello': 1.0}When working with arrays, take care not to end up in the following situation:
composition_to_array(", teste: 1", solution.species_names)array([0., 0., 0.])There is also a helper for generating inputs for use with tabulate.tabulate:
solution_report(
sol : cantera.composite.Solution | cantera.composite.Quantity,
specific_props : bool = True,
composition_spec : str = 'mass',
selected_species : list[str] = [],
**kwargs : Any
) -> list[tuple[str, str, typing.Any]]:Generate a solution report for tabulation.
sol : cantera.composite.Solution | cantera.composite.Quantityspecific_props : bool = Truecomposition_spec : str = 'mass'mass or mole.
selected_species : list[str] = []show_mass : bool = 'False'sol is a ct.composite.Quantity.
data = solution_report(solution, specific_props=True,
composition_spec="mass", selected_species=[])
print(tabulate(data))---------------------- -------- ------------
Temperature K 273.15
Pressure Pa 101325
Density kg/m³ 1.28735
Specific enthalpy J/(kg.K) -25114.9
Specific heat capacity J/(kg.K) 1004.95
mass: AR - 0.01
mass: N2 - 0.78
mass: O2 - 0.21
---------------------- -------- ------------Because sometimes Cantera lacks hard-copy utilities for certain classes, we provide simple wrappers that create new instances and set the state to the same of the source object. Nothing checked, nothing tested. A first of this kind is copy_solution, illustrated below:
copy_solution(sol : Solution): -> Solution:Makes a hard copy of a Solution object.
sol : Solutionnewairs = copy_solution(solution)
newairs.TPX = 373.15, None, newairs.X
print(tabulate(solution_report(newairs)))---------------------- -------- -------------
Temperature K 373.15
Pressure Pa 101325
Density kg/m³ 0.942356
Specific enthalpy J/(kg.K) 75902.1
Specific heat capacity J/(kg.K) 1015.82
mass: AR - 0.01
mass: N2 - 0.78
mass: O2 - 0.21
---------------------- -------- -------------print(tabulate(solution_report(solution)))---------------------- -------- ------------
Temperature K 273.15
Pressure Pa 101325
Density kg/m³ 1.28735
Specific enthalpy J/(kg.K) -25114.9
Specific heat capacity J/(kg.K) 1004.95
mass: AR - 0.01
mass: N2 - 0.78
mass: O2 - 0.21
---------------------- -------- ------------In addition to this, there is copy_quantity, which proves quite useful in establishing an algebra of mixtures.
copy_quantity(qty : Quantity): -> Quantity:Makes a hard copy of a ct.composite.Quantity object.
qty : Quantityair = copy_solution(solution)
air1 = ct.Quantity(air, mass=1.0)
air1.TPX = 373.15, None, newairs.X
air2 = copy_quantity(air1)
air2.TPX = 273.15, None, air2.X
mixair = air1 + air2
air1.T, air2.T, mixair.T(373.15, 273.15, 323.337485365699)Common daily work activity for the process engineer is to perform mass balances, but wait, …, gas flow rates are generally provided under normal conditions, and compositions may vary, so you need to compute normal densities first… whatever. This class provides a calculator wrapping a Cantera solution object so that your life gets easier.
NormalFlowRate(
mech : str | pathlib.Path,
*,
X : str | dict[str, float] | None = None,
Y : str | dict[str, float] | None = None,
T_ref : float = 273.15,
P_ref : float = 101325.0,
name : str | None = None
) -> None:Compute normal flow rate for a given composition.
This class makes use of the user defined state to create a function object that converts industrial scale flow rates in normal cubic meters per hour to kilograms per second. Nothing more, nothing less, it aims at helping the process engineer in daily life for this quite repetitive need when performing mass balances.
mech : str | pathlib.PathX : str | dict[str, float] | None = NoneX and Y are mutally exclusive keyword arguments.
Y : str | dict[str, float] | None = NoneX and Y are mutally exclusive keyword arguments.
T_ref : float = 273.15P_ref : float = 101325.0name : str | None = NoneIts simples use case is as follows:
nfr = NormalFlowRate("airish.yaml")
print(f"Convert 1000 Nm³/h to {nfr(1000.0):.5f} kg/s")Convert 1000 Nm³/h to 0.35903 kg/sIf the database file default composition does not suit you, no problems:
nfr = NormalFlowRate("airish.yaml", X="N2: 1")
print(f"Convert 1000 Nm³/h to {nfr(1000.0):.5f} kg/s")Convert 1000 Nm³/h to 0.34718 kg/sYou can also print a nice report to inspect the details of internal state. For more, please check its API documentation.
print(nfr.report())|------------------------|----------|--------------|
| Temperature | K | 273.15 |
| Pressure | Pa | 101325 |
| Density | kg/m³ | 1.24985 |
| Specific enthalpy | J/(kg.K) | -25864.9 |
| Specific heat capacity | J/(kg.K) | 1035.52 |
| mass: N2 | - | 1 |PlugFlowAxialSources(
n_reactors : int,
n_species : int
) -> None:Provides a data structure for use with PlugFlowChainCantera.
Helper data class for use with the solution method loop of the plug-flow reactor implementation. It provides the required memory for storage of source terms distributed along the reactor.
n_reactors : intn_species : intQ : NDArray[np.float64] = Nonem : NDArray[np.float64] = Noneh : NDArray[np.float64] = NoneY : NDArray[np.float64, np.float64] = NonePlugFlowChainCantera(
mechanism : str,
phase : str,
z : numpy.ndarray[tuple[typing.Any, ...], numpy.dtype[numpy.float64]],
V : numpy.ndarray[tuple[typing.Any, ...], numpy.dtype[numpy.float64]],
P : float = 101325.0,
K : float = 1.0,
smoot_flux : bool = False,
cantera_steady : bool = True
) -> None:Plug-flow reactor as a chain of 0-D reactors with Cantera.
mechanism : strphase : strz : numpy.ndarray[tuple[typing.Any, ...], numpy.dtype[numpy.float64]]V : numpy.ndarray[tuple[typing.Any, ...], numpy.dtype[numpy.float64]]P : float = 101325.0K : float = 1.0smoot_flux : bool = Falsecantera_steady : bool = Trueadvance_to_steady_state to solve problem; otherwise advance over meaninful time-scale of the problem.
get_reactor_data(pfr : PlugFlowChainCantera): -> PlugFlowAxialSources:Wrapper to allocate properly dimensioned solver data.
pfr : PlugFlowChainCanteraThis module is concerned with providing classes for computation of energy sources; these are mostly related to estimation of fuel consumption or integration to another simulation modules. Hereafter we illustrate the in cascading complexity the available energy sources. Some of these are simple wrapper around features provided by the combustion utilities.
A common need in industry is to evaluate the real heating value from a gas composition reported by the supplier. This is the main goal of CombustionAtmosphereCHON, which has a very simple API provided below:
gas = CombustionAtmosphereCHON("gri30.yaml")
lhv = gas.solution_heating_value("CH4: 1", "O2: 1")
print(f"Solution heating value: {lhv:.2f} MJ/kg")Solution heating value: 50.03 MJ/kgBuilt upon CombustionAtmosphereCHON, CombustionPowerSupply performs the basic calculations for retrieving required flow rates from a given power specification; this proves very handy for the CFD engineer.
supply = CombustionPowerSupply(500.0, 1.0, "CH4: 1", "O2: 1", "gri30.yaml")
print(supply.report())| Property | Unit | Value |
|---------------------------|--------|----------|
| Required power | kW | 500 |
| Lower heating value | MJ/kg | 50.0254 |
| Fuel mass flow rate | kg/h | 35.9817 |
| Oxidizer mass flow rate | kg/h | 143.532 |
| Total mass flow rate | kg/h | 179.514 |
| Fuel volume flow rate | Nm³/h | 50.2707 |
| Oxidizer volume flow rate | Nm³/h | 100.541 |
| Total volume flow rate | Nm³/h | 150.812 |
| Water production | kg/h | 80.8092 |
| Carbon dioxide production | kg/h | 98.7047 |
| Total emissions | kg/h | 179.514 |source = HeatedGasEnergySource("airish.yaml", 500.0, mass_flow_rate=1.0)
print(source.report())| Property | Unit | Value |
|------------------------|----------|-----------------------|
| Source kind | | HeatedGasEnergySource |
| Provided power | kW | 500.0 |
| Source | | airish.yaml |
| Phase | | air |
| Reference area | m² | 1.0 |
| Mass flow rate | kg/s | 1.0 |
| Volume flow rate | m³/s | 2.2051590808893957 |
| Momentum flux | kg.m/s² | 2.2051590808893957 |
| Reference temperature | K | 298.15 |
| Reference pressure | Pa | 101325.0 |
| Temperature | K | 778.5219246537463 |
| Pressure | Pa | 101324.99999999999 |
| Density | kg/m³ | 0.453482022529039 |
| Specific enthalpy | J/(kg.K) | 500038.49937969144 |
| Specific heat capacity | J/(kg.K) | 1091.748083028256 |
| mass: AR | - | 0.013790127718329308 |
| mass: N2 | - | 0.7542602692440457 |
| mass: O2 | - | 0.23194960303762516 |source = HeatedGasEnergySource("airish.yaml", 500.0, mass_flow_rate=1.0,
cross_area=0.1, Y="N2: 0.79, O2: 0.21")
print(source.report())| Property | Unit | Value |
|------------------------|----------|-----------------------|
| Source kind | | HeatedGasEnergySource |
| Provided power | kW | 500.0 |
| Source | | airish.yaml |
| Phase | | air |
| Reference area | m² | 0.1 |
| Mass flow rate | kg/s | 1.0 |
| Volume flow rate | m³/s | 2.2090693490835482 |
| Momentum flux | kg.m/s² | 22.090693490835484 |
| Reference temperature | K | 298.15 |
| Reference pressure | Pa | 101325.0 |
| Temperature | K | 774.4142430033143 |
| Pressure | Pa | 101324.99999999999 |
| Density | kg/m³ | 0.45267931512193527 |
| Specific enthalpy | J/(kg.K) | 500040.3236104384 |
| Specific heat capacity | J/(kg.K) | 1100.3346580471573 |
| mass: N2 | - | 0.79 |
| mass: O2 | - | 0.21 |fuel_state = StateType("CH4: 1", 300, 101325)
oxid_state = StateType("N2: 0.79, O2: 0.21", 300, 101325)ops = CombustionPowerOp(500.0, 1.0, fuel_state, oxid_state)
source = CombustionEnergySource("ch4/bfer.yaml", operation=ops, cross_area=0.1)
# Recover values for other tests that follow
# TODO wrap these in properties for clean API:
mdot_fuel = source._qty_fuel.mass
mdot_oxid = source._qty_oxid.mass
nfr_fuel = NormalFlowRate.new_from_solution(source._qty_fuel)
nfr_oxid = NormalFlowRate.new_from_solution(source._qty_oxid)
qdot_fuel = 3600 * mdot_fuel / nfr_fuel.density
qdot_oxid = 3600 * mdot_oxid / nfr_oxid.density
print(source.report())| Property | Unit | Value |
|------------------------|----------|------------------------|
| Source kind | | CombustionEnergySource |
| Provided power | kW | 500.0 |
| Source | | ch4/bfer.yaml |
| Phase | | CH4_BFER |
| Reference area | m² | 0.1 |
| Mass flow rate | kg/s | 0.18117768922524702 |
| Volume flow rate | m³/s | 1.2210278725486656 |
| Momentum flux | kg.m/s² | 2.2122300842798666 |
| ******* FLUE: | | |
| Temperature | K | 2257.7969133307734 |
| Pressure | Pa | 101325.00000819027 |
| Density | kg/m³ | 0.14838128866548536 |
| Specific enthalpy | J/(kg.K) | -254492.89776989678 |
| Specific heat capacity | J/(kg.K) | 1516.597922218735 |
| mass: CH4 | - | 2.7662274427222582e-17 |
| mass: CO | - | 0.010559409120413024 |
| mass: CO2 | - | 0.13474113875913193 |
| mass: H2O | - | 0.12389490081229192 |
| mass: N2 | - | 0.7247731344386522 |
| mass: O2 | - | 0.006031416869510841 |
| ******* FUEL: | | |
| Temperature | K | 300.0 |
| Pressure | Pa | 101325.00000000001 |
| Density | kg/m³ | 0.6516985521312658 |
| Specific enthalpy | J/(kg.K) | -4645856.881890704 |
| Specific heat capacity | J/(kg.K) | 2229.0429122783953 |
| mass: CH4 | - | 1.0 |
| ******* OXIDIZER: | | |
| Temperature | K | 300.0 |
| Pressure | Pa | 101325.00000000001 |
| Density | kg/m³ | 1.171970349439655 |
| Specific enthalpy | J/(kg.K) | 1907.6015935135354 |
| Specific heat capacity | J/(kg.K) | 1010.0686132988887 |
| mass: N2 | - | 0.7670907820415769 |
| mass: O2 | - | 0.2329092179584231 |ops = CombustionFlowOp("mass", mdot_fuel, mdot_oxid, fuel_state, oxid_state)
source = CombustionEnergySource("ch4/bfer.yaml", operation=ops, cross_area=0.1)
print(source.report())| Property | Unit | Value |
|------------------------|----------|------------------------|
| Source kind | | CombustionEnergySource |
| Provided power | kW | 500.0 |
| Source | | ch4/bfer.yaml |
| Phase | | CH4_BFER |
| Reference area | m² | 0.1 |
| Mass flow rate | kg/s | 0.18117768922524702 |
| Volume flow rate | m³/s | 1.2210278725486656 |
| Momentum flux | kg.m/s² | 2.2122300842798666 |
| ******* FLUE: | | |
| Temperature | K | 2257.7969133307734 |
| Pressure | Pa | 101325.00000819027 |
| Density | kg/m³ | 0.14838128866548536 |
| Specific enthalpy | J/(kg.K) | -254492.89776989678 |
| Specific heat capacity | J/(kg.K) | 1516.597922218735 |
| mass: CH4 | - | 2.7662274427222582e-17 |
| mass: CO | - | 0.010559409120413024 |
| mass: CO2 | - | 0.13474113875913193 |
| mass: H2O | - | 0.12389490081229192 |
| mass: N2 | - | 0.7247731344386522 |
| mass: O2 | - | 0.006031416869510841 |
| ******* FUEL: | | |
| Temperature | K | 300.0 |
| Pressure | Pa | 101325.00000000001 |
| Density | kg/m³ | 0.6516985521312658 |
| Specific enthalpy | J/(kg.K) | -4645856.881890704 |
| Specific heat capacity | J/(kg.K) | 2229.0429122783953 |
| mass: CH4 | - | 1.0 |
| ******* OXIDIZER: | | |
| Temperature | K | 300.0 |
| Pressure | Pa | 101325.00000000001 |
| Density | kg/m³ | 1.171970349439655 |
| Specific enthalpy | J/(kg.K) | 1907.6015935135354 |
| Specific heat capacity | J/(kg.K) | 1010.0686132988887 |
| mass: N2 | - | 0.7670907820415769 |
| mass: O2 | - | 0.2329092179584231 |ops = CombustionFlowOp("volume", qdot_fuel, qdot_oxid, fuel_state, oxid_state)
source = CombustionEnergySource("ch4/bfer.yaml", operation=ops, cross_area=0.1)
print(source.report())| Property | Unit | Value |
|------------------------|----------|------------------------|
| Source kind | | CombustionEnergySource |
| Provided power | kW | 500.0 |
| Source | | ch4/bfer.yaml |
| Phase | | CH4_BFER |
| Reference area | m² | 0.1 |
| Mass flow rate | kg/s | 0.18117768922524705 |
| Volume flow rate | m³/s | 1.2210278725486658 |
| Momentum flux | kg.m/s² | 2.212230084279867 |
| ******* FLUE: | | |
| Temperature | K | 2257.7969133307734 |
| Pressure | Pa | 101325.00000819027 |
| Density | kg/m³ | 0.14838128866548536 |
| Specific enthalpy | J/(kg.K) | -254492.89776989678 |
| Specific heat capacity | J/(kg.K) | 1516.597922218735 |
| mass: CH4 | - | 2.7662274427222582e-17 |
| mass: CO | - | 0.010559409120413024 |
| mass: CO2 | - | 0.13474113875913193 |
| mass: H2O | - | 0.12389490081229192 |
| mass: N2 | - | 0.7247731344386522 |
| mass: O2 | - | 0.006031416869510841 |
| ******* FUEL: | | |
| Temperature | K | 300.0 |
| Pressure | Pa | 101325.00000000001 |
| Density | kg/m³ | 0.6516985521312658 |
| Specific enthalpy | J/(kg.K) | -4645856.881890704 |
| Specific heat capacity | J/(kg.K) | 2229.0429122783953 |
| mass: CH4 | - | 1.0 |
| ******* OXIDIZER: | | |
| Temperature | K | 300.0 |
| Pressure | Pa | 101325.00000000001 |
| Density | kg/m³ | 1.171970349439655 |
| Specific enthalpy | J/(kg.K) | 1907.6015935135354 |
| Specific heat capacity | J/(kg.K) | 1010.0686132988887 |
| mass: N2 | - | 0.7670907820415769 |
| mass: O2 | - | 0.2329092179584231 |Simplest of relaxation methods; assume you have a new updated solution \(A_{new}^{\star}\) for a problem whose past state was \(A_{old}\), then the manager will ensure the following relaxation will be applied to compute the next solution state \(A_{new}\) to be used in whatever you are computing:
\[ \begin{aligned} A_{new} &= \alpha{}A_{old} + (1-\alpha)A_{new}^{\star}\\ A_{old} &= A_{new} \end{aligned} \]
Notice that in this formulation, \(\alpha\) (or alpha in the API) represents the fraction of old solution to be used in the smearing process. Below we illustrate the effect of a step function \(H\) valued at 10 from the begining over consecutive updates (here we do not test for convergence, as that is problem specific and for this simple case the required number of steps could be evaluated by hand, take some time to try!).
alpha = 0.76
niter = 50
single = np.ones(1)
relaxer = RelaxUpdate(single, alpha)
opts = dict(n_vars=1, max_iter=niter, patience=3, rtol=0.001)
converged = StabilizeNvarsConvergenceCheck(**opts)
history = np.zeros(niter+1)
history[0] = single[0]
H = np.asarray([10])
for n in range(niter):
single[:] = relaxer(H)
history[n+1] = single[0]
if converged(single[0]):
history = history[:n+2]
break
_ = plot_history(history)
RelaxUpdate(
v_ini : numpy.ndarray[tuple[typing.Any, ...], numpy.dtype[numpy.float64]],
alpha : float = 0.5
) -> None:Relax solution for updating new iteration.
v_ini : numpy.ndarray[tuple[typing.Any, ...], numpy.dtype[numpy.float64]]alpha : float = 0.5update(
self : Any,
alpha : float
) -> None:Update relaxation coefficients.
alpha : float__call__(
self : Any,
v_new : numpy.ndarray[tuple[typing.Any, ...], numpy.dtype[numpy.float64]]
) -> numpy.ndarray[tuple[typing.Any, ...], numpy.dtype[numpy.float64]]:Evaluate new relaxed solution estimate.
v_new : numpy.ndarray[tuple[typing.Any, ...], numpy.dtype[numpy.float64]]StabilizeNvarsConvergenceCheck(
*,
n_vars : int,
min_iter : int = 1,
max_iter : int = 1000000,
patience : int = 10,
rtol : float = 1e-10,
atol : float = 1e-20,
equal_nan : bool = False,
log_iter : bool = False
) -> None:Check stabilization towards a constant value along iterations.
n_vars : intmin_iter : int = 1max_iter : int = 1000000patience : int = 10rtol : float = 1e-10numpy.isclose for details.
atol : float = 1e-20numpy.isclose for details.
equal_nan : bool = Falsenumpy.isclose for details.
log_iter : bool = False__call__(
self : Any,
state : numpy.ndarray[tuple[typing.Any, ...], numpy.dtype[numpy.float64]]
) -> bool:Check if all variables have stabilized at current iteration.
state : numpy.ndarray[tuple[typing.Any, ...], numpy.dtype[numpy.float64]]n_iterations(self : Any): -> int:Provides access to number of iterations performed.
ComposedStabilizedConvergence(
n_arrs : int,
**kwargs : Any
):Wrapper for checking stabilization of several arrays.
See StabilizeNvarsConvergenceCheck for keyword arguments; these are shared by all tested arrays. It is always possible to compose your own convergence checker using individual instances for more control over setup.
n_arrs : intkwargs : = NoneStabilizeNvarsConvergenceCheck for details.
__call__(
self : Any,
*arrs : tuple[numpy.ndarray[tuple[typing.Any, ...], numpy.dtype[numpy.float64]], ...]
) -> bool:Check if all arrays have stabilized at current iteration.
*arrs : tuple[numpy.ndarray[tuple[typing.Any, ...], numpy.dtype[numpy.float64]], ...]n_iterations(self : Any): -> int:Provides access to number of iterations performed.
This module contains utilities for working with symbolic functions, especially in the context of engineering problems treated by Majordome models. It is mostly a wrapper around CasADi, organizing symbolic operations and providing additional functionality for function manipulation in a physical context.
A major issue when working with symbolic expressions and real-world data is that they do not support flow-control constructs like if statements. This is a problem when we want to represent functions that are defined piecewise, especially when algorithmic differentiation is involved.
As an example, take the parameterization of thermodynamic properties of chemical species. For quantitative purposes, to ensure high accuracy, researchers have found that it is best to use different polynomial fits for different temperature ranges. This means that the thermodynamic properties of a species are defined as piecewise functions of temperature, which is a problem when we want to use them in a symbolic context.
To tackle this issue in certain contexts, PiecewiseSymbolicFunction provides an interpolation between different ranges of the domain of a function. Internally it makes use of Heaviside functions to ensure that the function is continuous. This is not a perfect solution, as it does not guarantee differentiability at the breakpoints, but it is a good compromise between accuracy and simplicity. It is also worth noting that this class is not meant to be used in a general context, but rather in specific contexts where the function is known to be well-behaved.
PiecewiseSymbolicFunction(
breakpoints : list[float],
functions : list[typing.Any]
) -> None:Compose a symbolic piecewise function with CasADi.
breakpoints : list[float]functions : list[typing.Any]Let’s illustrate its use with a simple example. We define two functions f and g, and we want to create a piecewise function that is equal to f for x < 0.5 and equal to g for x >= 0.5. Functions have being created such that they evaluate to the same value at the breakpoint, as the practical case that led to the class development, i.e. NASA-parameterization of thermochemical properties. Below we show how to create the piecewise function, evaluate it numerically and symbolically, and compute its derivative symbolically.
# Define two simple functions:
f = lambda x: x**2
g = lambda x: f(0.5) + (x - 0.5)
# Create the piecewise function:
F = PiecewiseSymbolicFunction([0, 0.5, 1], [f, g])
# Evaluate numerically:
x = np.linspace(-0.2, 1.2, 100)
y = F(x)
# Evaluate symbolically:
T = SX.sym("T")
Y = F(T)
# Create a function for evaluating the derivative:
ydot = Function("fdot", [T], [jacobian(Y, T)])Below we plot the piecewise function and its derivative. The function is continuous, but the derivative is not, as expected, although the transition is handled smoothly by the Heaviside functions. The red dashed line indicates the breakpoint at x = 0.5, where the function transitions from f to g.
symbolic_thermo_factory(
species : Species,
T : SX
) -> AbstractSymbolicThermo:Create an AbstractSymbolicThermo object.
species : SpeciesT : SXSuppose you need to create a symbolic representation of the NASA7 thermodynamic parameterization of molecular nitrogen, species of index 47 in the GRI-Mech 3.0 mechanism shipped with Cantera. Let’s start by loading and retrieving the thermodynamic data for this species.
gas = ct.Solution("gri30.yaml")
species = gas.species()[47]
thermo = species.thermoBelow we instantiate SymbolicThermo with the input data of the species. It is important to emphasize here that when building more complex systems, generally one should share the symbolic variables, e.g. temperature, across the different components of the system. For this reason, the constructor of symbolic thermodynamic classes expect the temperature symbol to be provided, instead of trying to create it internally; this is a design choice that will improve code robustness when handling mixtures and avoid many of the pitfalls of symbolic programming.
T = SX.sym("T")
nasa7 = symbolic_thermo_factory(species, T)You could instead directly instantiate Nasa7Thermo with the input data, but using the factory method SymbolicThermo.from_species is a better choice, as it will allow you to easily switch to a different thermodynamic model. Furthermore, it is the safe choice for loading entire databases.
nasa7 = Nasa7Thermo(T, species.thermo.input_data)
# This is just syntactic sugar for the above:
nasa7 = Nasa7Thermo.from_species(species, T)Here we illustrate the use of algorithmic differentiation to compute the derivative of the specific enthalpy with respect to temperature.
hdot = Function("hdot", [T], [jacobian(nasa7.h(T), T)])The above derivative is then confronted to the specific heat below:
Finally, we compare the deviations between symbolic and numerical (Cantera) evaluation of the thermodynamic functions. This is not for error analysis, but a proof of correctness of the symbolic implementation.
# XXX: Cantera evaluates quantities per kmol! Divide by 1000 to get
# per mol, as the symbolic implementation does.
c_ct = np.array([0.001 * species.thermo.cp(T) for T in T_num])
h_ct = np.array([0.001 * species.thermo.h(T) for T in T_num])
s_ct = np.array([0.001 * species.thermo.s(T) for T in T_num])
error_c = np.mean(np.abs(nasa7._cp(T_num) - c_ct))
error_h = np.mean(np.abs(nasa7._h(T_num) - h_ct))
error_s = np.mean(np.abs(nasa7._s(T_num) - s_ct))
error_c, error_h, error_s(np.float64(4.334310688136611e-15),
np.float64(6.311893230304122e-12),
np.float64(3.2400748750660566e-14))Nasa7Thermo(
T : SX,
input_data : dict[str, typing.Any]
) -> None:NASA7 thermodynamic parameterization.
This class does not implement Horner polynomial evaluation, as the main use case is to create symbolic expressions that are then evaluated by CasADi, which can handle polynomial evaluation efficiently. This is intentional and allows for easy verification.
It aims at providing a similar interface as SpeciesThermo from Cantera, from which it retrieves the data. That means, molar properties are provided through cp, h, and s properties.
T : SXinput_data : dict[str, typing.Any]SpeciesThermo property input_data.
specific_heat(
T : Any,
a : list[float],
symbolic : bool = False
) -> casadi.casadi.Function | casadi.casadi.SX:Compose NASA7 specific heat parameterization.
T : Anya : list[float]symbolic : bool = Falseenthalpy(
T : Any,
a : list[float],
symbolic : bool = False
) -> casadi.casadi.Function | casadi.casadi.SX:Compose NASA7 specific enthalpy parameterization.
T : Anya : list[float]symbolic : bool = Falseentropy(
T : Any,
a : list[float],
symbolic : bool = False
) -> casadi.casadi.Function | casadi.casadi.SX:Compose NASA7 specific entropy parameterization.
T : Anya : list[float]symbolic : bool = Falsecompose(
cls : Any,
name : str,
T : casadi.casadi.SX | numpy.ndarray[tuple[typing.Any, ...], numpy.dtype[numpy.float64]] | float,
data : list[list[float]],
symbolic : bool = False
) -> list[casadi.casadi.Function | casadi.casadi.SX | numpy.ndarray[tuple[typing.Any, ...], numpy.dtype[numpy.float64]] | float]:Compose a list of NASA7 functions for given data.
name : strT : casadi.casadi.SX | numpy.ndarray[tuple[typing.Any, ...], numpy.dtype[numpy.float64]] | floatdata : list[list[float]]symbolic : bool = FalseT is symbolic.
from_species(
cls : Any,
species : Species,
T : SX
) -> typing.Self:Create a Nasa7Thermo object from a Cantera species.
species : SpeciesT : SXsymbolic_transport_factory(
species : Species,
T : SX
) -> AbstractSymbolicTransport:Create an AbstractSymbolicTransport object.
species : SpeciesT : SXIn some fields of research, such as porous media or composite materials, the evaluation of effective thermal properties are key for simulation of macroscopic application cases. Class EffectiveThermalConductivity implements (static) thermal conductivity models for this sort of applications, including:
Maxwell-Garnett approximation (EffectiveThermalConductivity.maxwell_garnett) as exposed in (Hanein, Glasser, and Bannerman 2017); further discussion of its origin and applicability is provided by (Kiradjiev et al. 2019).
For porous media (packed bed in the context) with enhanced radiative effects the model by (Singh and Kaviany 1994) as discussed by (Kee et al. 2017) is implemented (EffectiveThermalConductivity.singh1994).
Its name is quite long, let’s start by getting an alias before evaluating the desired models:
etc = EffectiveThermalConductivity()Using Maxwell approximation, we could estimate the the effective thermal conductivity of a packed bed of particles in a matrix of air; assume particles loosly embeded in air and the following properties; the computed effective thermal conductivity is shown to approach the air limit:
phi = 0.30 # Loosely packed solids [30%v]
k_g = 0.025 # Air thermal conductivity [W/(m.K)]
k_s = 1.000 # Solids thermal conductivity [W/(m.K)]
etc.maxwell_garnett(phi, k_g, k_s)0.05396039603960396For the limit of high temperatures, it is usually important to account for particle-particle radiation heat transfer; this introduces a \(T^3\) dependence on temperature, as one should expect by linearizing Stefan-Boltzmann law.
Please notice that these models compute different things; while Maxwell approximation computes the medium properties (to approximate matrix-inclusion as a single domain), Singh’s model accounts only for solids properties. One might wish to combine them (warning: unverified validity!) to evaluate overall medium thermal conductivity. See the references in the class documentation for further discussion, specially the extension proposed by Kiradjiev (2019), which leads to a result similar to the assymptotic behavior displayed below.
d_p = 0.005
eps = 0.9
T = np.linspace(300, 1500, 50)
k_eff_s = etc.singh1994(T, phi, d_p, k_s, eps)
k_eff_m = etc.maxwell_garnett(phi, k_g, k_eff_s)
plot_etc((T, k_eff_s, k_eff_m)).resize(10, 5)
Due to kinetic theory implications, gas thermal conductivity tend to have a positive slope in temperature; the following illustrates how this can be accounted for and the roughly linear behaviour introduced when computing medium properties within air.
gas = ct.Solution("airish.yaml")
sol = ct.SolutionArray(gas, (T.shape[0],))
sol.TP = T, None
k_gas = sol.thermal_conductivity
k_eff_m = etc.maxwell_garnett(phi, k_gas, k_eff_s)
plot_etc((T, k_gas, k_eff_m)).resize(10, 5)
A dimensionless numbers calculator is provided for gas flows; it currently has a certain number of groups which are all evaluated by definition (which might change according to your field, please check the docs). The mechanics of using the class can be resumed to:
The meaning of the tuple of arguments provided to set_state is specified by tuple_name="TPX", which defaults to temperature, pressure, and molar proportions. Any triplet allowed by Cantera can be specified here.
Tw = 1000.0 # Wall temperature [K]
U = 10.0 # Characteristic velocity [m/s]
D = 0.05 # Pipe diameter [m]
L = 1.0 # Pipe length [m]
calculator = SolutionDimless("airish.yaml")
calculator.set_state(300.0, 101325.0, "N2: 1", tuple_name="TPX")
Re = calculator.reynolds(U, D)
Pr = calculator.prandtl()
Sc = calculator.schmidt()
print(calculator.report())-------- ------------ ---------------
Reynolds 31460.9 U=10.0, L=0.05
Prandtl 0.709327
Schmidt 0.761752 mix_diff_coeffs
-------- ------------ ---------------If you prefer to have direct access to the internal solution, you can set properties as usual in Cantera, but you need to keep in mind to call update() to refresh the internal state of the calculator. Every time the properties are updated, the internal buffer of computed dimensionless numbers is refreshed, as you migth notice in the following table.
calculator.solution.TPX = 300.0, 101325.0, "N2: 1"
calculator.update()
Pe_m = calculator.peclet_mass(U, L)
Pe_h = calculator.peclet_heat(U, L)
Gr = calculator.grashof(Tw, D)
Ra = calculator.rayleigh(Tw, D)
print(calculator.report())------------- ---------------- ------------------------------
Péclet (mass) 479308 U=10.0, L=1.0, mix_diff_coeffs
Péclet (heat) 446321 U=10.0, L=1.0
Grashof 1.13242e+07 Tw=1000.0, H=0.05, g=9.80665
Rayleigh 8.03259e+06 Tw=1000.0, H=0.05, g=9.80665
------------- ---------------- ------------------------------When performing CFD simulations of turbulent flows, it is important to ensure that the first layer of cells adjacent to the wall are properly graded to capture the boundary layer effects. The WallGradingCalculator class provides a simple interface for calculating the first layer thickness based on the desired dimensionless wall distance \(y^+\) and the flow properties.
The following example illustrates how to use the WallGradingCalculator to compute the first layer thickness for a turbulent flow of air at 700°C and atmospheric pressure, with a characteristic velocity of 10 m/s and a pipe diameter of 0.025 m. The skin friction factor is computed using the SkinFrictionFactor.smooth_wall method, which is appropriate for turbulent flows over smooth walls.
# First create a solution object with the desired state:
sol = SolutionDimless("airish.yaml")
sol.set_state(973.15, 101325, "N2: 0.79, O2: 0.21")
# Define a function to compute the skin friction factor for turbulent flow:
def f_tur(Re):
""" Turbulent skin friction factor for smooth walls."""
return SkinFrictionFactor.smooth_wall(Re, check=False) / 8
# Define the desired dimensionless wall distance y+:
y_plus = 1.0
# Instantiate the wall grading calculator from the solution and compute
# the first layer thickness with the provided skin friction factor:
calc = WallGradingCalculator.from_solution(sol, L=0.025, U=10.0)
y_first = calc.first_layer(y_plus, skin_factor=f_tur)
print(f"First layer thickness: {1e6*y_first:.0f} µm")First layer thickness: 153 µmWallGradingCalculator(
*,
L : float,
U : float,
rho : float,
mu : float,
skin_factor : typing.Optional[typing.Callable] = None
) -> None:Helper class for estimating first cell thickness given y+.
L : floatU : floatrho : floatmu : floatskin_factor : typing.Optional[typing.Callable] = Nonefirst_layer method will raise an error if called. If provided, it should be a callable that takes Reynolds number as input and returns the skin friction factor. Some examples of skin friction factors are provided in SkinFrictionFactor class.
set_skin_factor(
self : Any,
skin_factor : typing.Callable
) -> None:Set skin friction factor and compute related properties.
skin_factor : typing.Callable__init__ for details.
first_layer(
self : Any,
y_plus : Any,
skin_factor : typing.Optional[typing.Callable] = None
) -> float:Height of first cell for given y+ value [m].
y_plus : Anyskin_factor : typing.Optional[typing.Callable] = None__init__ for details.
from_solution(
cls : Any,
obj : SolutionDimless,
L : float,
U : float,
skin_factor : typing.Optional[typing.Callable] = None
) -> typing.Self:Alternative constructor from dimensionless solution.
obj : SolutionDimlessL : floatU : floatskin_factor : typing.Optional[typing.Callable] = None__init__ for details.
wall_shear_stress(
rho : Any,
U : Any,
Cf : Any
) -> float:Wall shear stress estimater from friction factor [Pa].
rho : AnyU : AnyCf : Anyfriction_velocity(
tw : Any,
rho : Any
) -> float:Dimensionless friction velocity.
tw : Anyrho : AnySkinFrictionFactor():Skin friction factors for y+ calculations.
laminar(Re : Any): -> float:Laminar limit theoretical value.
Re : Anysmooth_wall(
Re : Any,
check : bool = True
) -> float:Blasius smooth wall approximation.
As described here.
Re : Anycheck : bool = TrueTo the author’s knowledge, there is no standard tool to convert Cantera transport data to Sutherland parameters for use with OpenFOAM, what led to the motivation to develop SutherlandFitting. This simple class wraps a Cantera solution object, which it makes use for fitting Sutherland parameters to export as a table (to be used elswhere), and also allows for retrieving a converted database in OpenFOAM compatible format. The following example should be self-explanatory:
T = np.linspace(500, 2500, 100)
sutherland = SutherlandFitting("airish.yaml")
sutherland.fit(T, species_names=["O2", "N2"])
coef = sutherland.coefs_table
coef| species | As [uPa.s] | Ts [K] | RMSE [uPa.s] | |
|---|---|---|---|---|
| 0 | O2 | 1.873876 | 229.321211 | 0.527484 |
| 1 | N2 | 1.619350 | 226.141806 | 0.470152 |
If needed, it is also possible to have direct access to the viscosity data used in parameter fitting:
sutherland.viscosity.head()| T | O2 | N2 | |
|---|---|---|---|
| 0 | 500.000000 | 30.045456 | 26.123142 |
| 1 | 520.202020 | 30.886884 | 26.845049 |
| 2 | 540.404040 | 31.713820 | 27.554780 |
| 3 | 560.606061 | 32.527146 | 28.253075 |
| 4 | 580.808081 | 33.327662 | 28.940602 |
Because RMSE compresses all the error in a single value, you can check graphically where the deviations happen in the interval:
plot = sutherland.plot_species("N2")
Finally, it also responds to its initial goal of exporting values in OpenFOAM format:
print(sutherland.to_openfoam())
/* --- RMSE 0.5274842993902649 ---*/
"O2"
{
transport
{
As 1.8738763185e-06;
Ts 2.2932121082e+02;
}
}
/* --- RMSE 0.4701517988286822 ---*/
"N2"
{
transport
{
As 1.6193495479e-06;
Ts 2.2614180620e+02;
}
} Radiative properties of gases are often required when participating media characteristics imply so. Here we implement the WSGG model of Sadeghi et al. (2021) in class WSGGRadlibBordbar2020. The basic usage of the model to estimate total emissivity is done as follows:
model = WSGGRadlibBordbar2020()
model(L=1, T=1000, P=101325, x_h2o=0.18, x_co2=0.08, fvsoot=0.0)np.float64(0.18368577452587137)Evaluation of the model against original source of Bordbar (2014) is satisfactory, as follows:
Absorption coefficients for pure substances only reproduce approximately values reported by Bordbar (2020).
model(L=1, T=300, P=101325, x_h2o=0, x_co2=1)
model.absorption_coefs[1:]array([3.388079e-02, 4.544269e-01, 4.680226e+00, 1.038439e+02])model(L=1, T=300, P=101325, x_h2o=1, x_co2=0)
model.absorption_coefs[1:]array([ 0.07703541, 0.8242941 , 6.854761 , 65.93653 ])Below we illustrate the application of WSGGRadlibBordbar2020 to predict the emissivity of \(x\mathrm{H_2O}-(1-x)\mathrm{CO_2}\) mixtures over a broad temperature range applicable to the analysis of combustion processes. One observes an important participation of flue gases, especially at lower temperatures.
wsgg = WSGGRadlibBordbar2020()
@np.vectorize
def emissivity(T, X, L=1, P=101325):
return wsgg(L=L, T=T, P=P, x_h2o=X, x_co2=1-X, fvsoot=0.0)
T = np.linspace(800, 2400, 100)
X = np.linspace(0, 0.7, 100)
T, X = np.meshgrid(T, X)
eps = emissivity(T, X)
plot_emissivity(T, X, eps)
Upcoming…