nicolas
/
gaspype
mirror of https://github.com/DLR-Institute-of-Future-Fuels/gaspype.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

Solver updated with start value estimation and log-error for elemental species balance

This commit is contained in:
Nicolas Kruse 2025-06-24 09:45:36 +02:00 committed by Nicolas Kruse
parent a5f806aa31
commit 63b1808341
2 changed files with 22 additions and 12 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

32
src/gaspype/_solver.py
View File

@ -69,6 +69,7 @@ def equilibrium_eq(fs: fluid_system, element_composition: FloatArray, t: float,
"""Calculate the equilibrium composition of a fluid based on equilibrium equations""" """Calculate the equilibrium composition of a fluid based on equilibrium equations"""
el_max = np.max(element_composition) el_max = np.max(element_composition)
element_norm = element_composition / el_max element_norm = element_composition / el_max
element_norm_log = np.log(element_norm + epsy)
a = fs.array_stoichiometric_coefficients a = fs.array_stoichiometric_coefficients
a_sum = np.sum(a) a_sum = np.sum(a)
@ -81,35 +82,44 @@ def equilibrium_eq(fs: fluid_system, element_composition: FloatArray, t: float,
bp = b - np.sum(a * np.log(p / p0), axis=1) bp = b - np.sum(a * np.log(p / p0), axis=1)
# Calculating the maximum possible amount for each species based on the elements # Calculating the maximum possible amount for each species based on the elements
logn_start = np.min(element_norm / (fs.array_species_elements + epsy), axis=1) species_max = np.min(element_norm / (fs.array_species_elements + epsy), axis=1)
logn_start = np.log(species_max + epsy)
def residuals(logn: FloatArray): # type: ignore # global count
# count = 0
def residuals(logn: FloatArray) -> tuple[FloatArray, FloatArray]:
# global count
# count += 1
# print('------', count)
# assert count < 100
n = np.exp(logn) n = np.exp(logn)
n_sum = np.sum(n) n_sum = np.sum(n)
# Residuals from equilibrium equations: # Residuals from equilibrium equations:
eq_resid = np.dot(a, logn - np.log(n_sum)) - bp resid_eq = np.dot(a, logn - np.log(n_sum)) - bp
# Derivative: # Jacobian:
j_eq = a - a_sum * n / n_sum j_eq = a - a_sum * n / n_sum
# Residuals from elemental balance: # Residuals from elemental balance:
el_error = np.dot(el_matrix, n) - element_norm el_error = np.dot(el_matrix, n)
ab_resid = np.log1p(el_error) resid_ab = np.log(el_error) - element_norm_log
# Derivative: # Jacobian
j_ab = el_matrix * n / np.expand_dims(el_error + 1, axis=1) j_ab = el_matrix * n / el_error[:, np.newaxis]
return (np.hstack([eq_resid, ab_resid]), np.concatenate([j_eq, j_ab], axis=0)) return (np.hstack([resid_eq, resid_ab]), np.concatenate([j_eq, j_ab], axis=0))
ret = root(residuals, logn_start, jac=True, tol=1e-30) ret = root(residuals, logn_start, jac=True, tol=1e-10)
n = np.exp(np.array(ret['x'], dtype=NDFloat)) n = np.exp(np.array(ret['x'], dtype=NDFloat))
# print(ret)
return n * el_max return n * el_max
def equilibrium(f: fluid | elements, t: float | FloatArray, p: float = 1e5) -> fluid: def equilibrium(f: fluid | elements, t: float | FloatArray, p: float = 1e5) -> fluid:
"""Calculate the equilibrium composition of a fluid at a given temperature and pressure" """Calculate the isobaric equilibrium composition of a fluid at a given temperature and pressure"
Args: Args:
f: Fluid or elements object f: Fluid or elements object

2
src/gaspype/constants.py
View File

@ -21,4 +21,4 @@ p0 = 1e5 # Pa
t0 = 273.15 + 25 # K t0 = 273.15 + 25 # K
p_atm = 101325 # Pa p_atm = 101325 # Pa
epsy = 1e-18 epsy = 1e-30