-
Notifications
You must be signed in to change notification settings - Fork 1
/
solution.py
122 lines (108 loc) · 3.9 KB
/
solution.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
from typing import Tuple, Union
import cantera as ct
import jax.numpy as np
# local imports
from .dataclass import (
GasInfo,
KineticsCoeffs,
KineticsData,
NASAPolynomials,
ProductionRates,
)
from .kinetics import (
get_equilibirum_constants,
get_forward_rate_constants,
get_reverse_rate_constants,
)
from .thermo import calculate_cp, calculate_enthalpy, calculate_entropy
def get_initial_mole_fractions_and_state(
path: str,
T: Union[float, np.ndarray],
P: Union[float, np.ndarray],
X: Union[str, dict],
) -> Tuple[np.ndarray, GasInfo]:
"""
Helper function to get Cantera's Solution object
"""
gas = ct.Solution(path)
gas.TPX = T, P, X
gas_info = GasInfo(
reactant_stioch_coeffs=gas.reactant_stoich_coeffs(),
product_stioch_coeffs=gas.product_stoich_coeffs(),
molecular_weights=gas.molecular_weights,
)
return np.array(gas.Y, dtype=np.float64), gas_info
def Y2X(Y: np.ndarray, mol_wts: np.ndarray, mean_mol_wt: np.ndarray):
"""
helper function to calculate mass fractions from mole fractions
"""
return np.divide(np.multiply(Y, mean_mol_wt), mol_wts)
def Y2C(Y: np.ndarray, mol_wts: np.ndarray, density_mass: np.ndarray):
"""
helper function to calculate concentrations from mole fractions
"""
return np.divide(np.multiply(Y, density_mass), mol_wts)
def get_mean_molecular_weight(Y: np.ndarray, mol_wts: np.ndarray) -> np.ndarray:
"""
helper function to calculate concentrations from mole fractions
"""
return 1.0 / (np.matmul(Y, 1.0 / mol_wts))
def get_production_rates(
kf: np.ndarray, kr: np.ndarray, C: np.ndarray, gas_info: GasInfo
) -> np.ndarray:
"""
calculate net production rates from rate constants and concentration
"""
vk = -gas_info.reactant_stioch_coeffs + gas_info.product_stioch_coeffs
forward_rates_of_progress = kf * np.exp(
np.matmul(np.log(C + 1e-300), gas_info.reactant_stioch_coeffs)
)
reverse_rates_of_progress = kr * np.exp(
np.matmul(np.log(C + 1e-300), gas_info.product_stioch_coeffs)
)
qdot = np.subtract(forward_rates_of_progress, reverse_rates_of_progress)
wdot = np.matmul(qdot, vk.T)
return ProductionRates(
forward_rates_of_progress=forward_rates_of_progress,
reverse_rates_of_progress=reverse_rates_of_progress,
qdot=qdot,
wdot=wdot,
)
def get_dYdt(
P: Union[float, np.ndarray],
R: float,
gas_info: GasInfo,
nasa_poly: NASAPolynomials,
kinetics_coeffs: KineticsCoeffs,
kinetics_data: KineticsData,
t: Union[float, np.ndarray],
state_vec: np.ndarray,
) -> np.ndarray:
"""
get dYdt where dYdt[0] is dTdt and rest is dYdt
"""
T = state_vec[0]
Y = state_vec[1:]
T = np.clip(T, a_min=200.0, a_max=1e5)
Y = np.clip(Y, a_min=0.0, a_max=1.0)
mean_molecular_weight = get_mean_molecular_weight(Y, gas_info.molecular_weights)
density_mass = P / R / T * mean_molecular_weight
X = Y2X(Y, gas_info.molecular_weights, mean_molecular_weight)
C = Y2C(Y, gas_info.molecular_weights, density_mass)
cp_data = calculate_cp(T, X, R, nasa_poly)
enthalpy_data = calculate_enthalpy(T, X, R, nasa_poly)
entropy_data = calculate_entropy(T, P, X, R, nasa_poly)
cp_mass = cp_data.cp_mole / mean_molecular_weight
Kc = get_equilibirum_constants(
T, P, R, entropy_data.sdivR, enthalpy_data.hdivRT, gas_info
)
kf = get_forward_rate_constants(T, R, C, kinetics_coeffs, kinetics_data)
kr = get_reverse_rate_constants(kf, Kc, kinetics_data.is_reversible)
production_rates = get_production_rates(kf, kr, C, gas_info)
Ydot = (production_rates.wdot * gas_info.molecular_weights) / density_mass
Tdot = -(
np.matmul(enthalpy_data.partial_molar_enthalpies, production_rates.wdot)
) / (density_mass * cp_mass)
dYdt = np.hstack((Tdot, Ydot))
# dYdt = delete_small_numbers(dYdt)
return dYdt