diff --git a/src/cryo_tank/README.md b/src/cryo_tank/README.md index 7c2b017..ee16952 100644 --- a/src/cryo_tank/README.md +++ b/src/cryo_tank/README.md @@ -46,7 +46,7 @@ OUTPUT_DIR = "results/cryo_tank" - `cryo_tank_level.png`：液位和充满率曲线。 - `cryo_tank_he_flow.png`：氦气质量和氦气流量曲线。 - `cryo_tank_heat.png`：漏热和气液换热曲线。 -- `cryo_tank_pressure.png`：总压、氮气分压和氦气分压曲线。 +- `cryo_tank_pressure.png`：总压和氦气气枕压力曲线。 如果没有看到 CSV 文件，先确认本地代码里是否有： @@ -105,7 +105,7 @@ y = [m_liq, U_liq, U_ull] - `U_liq`：液相内能 - `U_ull`：气枕区总内能 -`derive(y)` 会根据状态量计算温度、体积、液位、充满率、氦气质量、氮气分压和氦气分压等派生量。 +`derive(y)` 会根据状态量计算温度、体积、液位、充满率、氦气质量和氦气气枕压力等派生量。 `rhs(t, y)` 是 ODE 右端函数，由 `scipy.integrate.solve_ivp` 调用。 @@ -133,7 +133,7 @@ run(tank, t_end, rtol=1e-8, atol=1e-10, max_step=10.0) - 液氮密度、焓、内能计算。 - 氮气饱和压力和饱和蒸气内能计算。 -- 氦气理想气体焓、内能、比热计算。 +- 氦气密度、焓、内能、比热计算。 - 建立液氮内能到温度的查表插值，加速 ODE 求解。 该模块依赖 `CoolProp`。 @@ -348,7 +348,6 @@ CSV 第一行是列名，列名来自 `history` 字典，包括： - `V_ull` - `liquid_level` - `fill_fraction` -- `P_N2` - `P_He` - `P_total` - `Q_leak` @@ -373,7 +372,7 @@ pytest -q ## 注意事项 - 本模型当前将储箱分为液相区和气枕区两个区域，不是完整 CFD 模型。 -- 气枕区氮气和氦气采用简化处理，氦气按理想气体热力学关系计算。 +- 气枕区按纯氦气处理，物性由 CoolProp 计算。 - 液氮物性依赖 `CoolProp`，缺少该包会导致程序无法启动。 - `results/` 目录默认被 `.gitignore` 忽略，仿真输出不会自动上传到 Git。 - 建议每次修改参数后保存对应工况说明，避免不同仿真结果混在同一个输出目录中。 diff --git a/src/cryo_tank/config.py b/src/cryo_tank/config.py index 01f0b71..56be7bc 100644 --- a/src/cryo_tank/config.py +++ b/src/cryo_tank/config.py @@ -5,13 +5,6 @@ Pure data module -- no functions, no side effects. """ import math -# ---------- Gas constants ---------- -R_UNIVERSAL = 8314.46 # J/(kmol*K) -M_N2 = 28.014 # kg/kmol -M_HE = 4.0026 # kg/kmol -R_HE = R_UNIVERSAL / M_HE # 2077.1 J/(kg*K) -R_N2 = R_UNIVERSAL / M_N2 # 296.8 J/(kg*K) - # ---------- Tank geometry ---------- V_TOTAL = 420.1e-3 # m^3 (420.1 L) H_TANK = 0.5 # m (cylinder height) diff --git a/src/cryo_tank/main.py b/src/cryo_tank/main.py index 0a353fb..d486d34 100644 --- a/src/cryo_tank/main.py +++ b/src/cryo_tank/main.py @@ -49,7 +49,7 @@ def main(): print(f" T_liq = {info0['T_liq']:.2f} K, T_ull = {info0['T_ull']:.2f} K") print(f" fill_fraction = {info0['fill_fraction']:.1%}") print(f" m_He = {info0['m_He']*1000:.2f} g") - print(f" P_N2 = {info0['P_N2']/1e6:.4f} MPa, P_He = {info0['P_He']/1e6:.4f} MPa") + print(f" P_He = {info0['P_He']/1e6:.4f} MPa") print(f"Running to t_end = {T_END:.0f} s ...") print() diff --git a/src/cryo_tank/output.py b/src/cryo_tank/output.py index c22b948..11498e4 100644 --- a/src/cryo_tank/output.py +++ b/src/cryo_tank/output.py @@ -116,12 +116,11 @@ def plot_heat_fluxes(history, path): def plot_pressure(history, path): - """Plot tank pressure (P_total, P_N2, P_He) vs time.""" + """Plot tank pressure (P_total and P_He) vs time.""" fig, ax = plt.subplots(figsize=(10, 5)) t = history['t'] ax.plot(t, history['P_total'] / 1e6, label='P_total', linewidth=2) - ax.plot(t, history['P_N2'] / 1e6, label='P_N2', linestyle='--') ax.plot(t, history['P_He'] / 1e6, label='P_He', linestyle='--') ax.set_xlabel('Time [s]') ax.set_ylabel('Pressure [MPa]') diff --git a/src/cryo_tank/properties.py b/src/cryo_tank/properties.py index 8d23d0d..66d1190 100644 --- a/src/cryo_tank/properties.py +++ b/src/cryo_tank/properties.py @@ -1,13 +1,11 @@ # src/cryo_tank/properties.py """ -Fluid property wrappers for liquid nitrogen, N2 vapor, and helium. +Fluid property wrappers for liquid nitrogen and helium. -Performance strategy (per spec Section 8.1): - - N2 (liquid & vapor): CoolProp with persistent AbstractState objects - - He: analytical ideal gas (cp=5196.2 J/(kg*K), cv=3117.1 J/(kg*K)) - - Lookup tables for the ODE hot path (built at import time) +Performance strategy: + - N2 liquid: CoolProp with a persistent AbstractState object + - He: CoolProp with a persistent AbstractState object """ -import numpy as np import CoolProp.CoolProp as CP from CoolProp import AbstractState @@ -41,81 +39,60 @@ def ln2_u(T, P): def ln2_T_from_u(u, P): - """Recover LN2 temperature from specific internal energy [K]. - - Uses lookup table interpolation for speed; falls back to CoolProp - if outside the table range. - """ - return float(np.interp(u, _ln2_u_table, _ln2_T_table)) + """Recover LN2 temperature from pressure and specific internal energy [K].""" + _n2_state.update(CP.PUmass_INPUTS, P, u) + return _n2_state.T() # --------------------------------------------------------------------------- -# N2 vapor properties (at saturation or specified conditions) +# Helium properties via CoolProp # --------------------------------------------------------------------------- -def n2_sat_pressure(T): - """N2 saturation pressure [Pa] at temperature T.""" - _n2_state.update(CP.QT_INPUTS, 1.0, T) - return _n2_state.p() +_HE_P_REF = 170000.0 -def n2_vapor_u(T): - """N2 saturated vapor specific internal energy [J/kg] at temperature T.""" - _n2_state.update(CP.QT_INPUTS, 1.0, T) - return _n2_state.umass() +def he_rho(T, P=_HE_P_REF): + """He density [kg/m^3].""" + _he_state.update(CP.PT_INPUTS, P, T) + return _he_state.rhomass() -def n2_vapor_rho(T): - """N2 saturated vapor density [kg/m^3] at temperature T.""" - _n2_state.update(CP.QT_INPUTS, 1.0, T) - return _n2_state.rhomass() - - -# --------------------------------------------------------------------------- -# Helium properties (ideal gas: cp=5/2 R, cv=3/2 R, monatomic) -# --------------------------------------------------------------------------- -_HE_CP = 5196.2 # J/(kg*K), = 5/2 * R_He -_HE_CV = 3117.1 # J/(kg*K), = 3/2 * R_He - -# Reference state: CoolProp He at T_ref=0K gives u_ref, h_ref -# We match CoolProp's reference by computing offset at a known point. -_he_state.update(CP.PT_INPUTS, 170000.0, 100.0) -_HE_H_REF = _he_state.hmass() - _HE_CP * 100.0 # h = cp*T + h_ref -_HE_U_REF = _he_state.umass() - _HE_CV * 100.0 # u = cv*T + u_ref - - -def he_cp(): +def he_cp(T=100.0, P=_HE_P_REF): """He specific heat at constant pressure [J/(kg*K)].""" - return _HE_CP + _he_state.update(CP.PT_INPUTS, P, T) + return _he_state.cpmass() -def he_cv(): +def he_cv(T=100.0, P=_HE_P_REF): """He specific heat at constant volume [J/(kg*K)].""" - return _HE_CV + _he_state.update(CP.PT_INPUTS, P, T) + return _he_state.cvmass() -def he_h(T): - """He specific enthalpy [J/kg] (ideal gas).""" - return _HE_CP * T + _HE_H_REF +def he_h(T, P=_HE_P_REF): + """He specific enthalpy [J/kg].""" + _he_state.update(CP.PT_INPUTS, P, T) + return _he_state.hmass() -def he_u(T): - """He specific internal energy [J/kg] (ideal gas).""" - return _HE_CV * T + _HE_U_REF +def he_u(T, P=_HE_P_REF): + """He specific internal energy [J/kg].""" + _he_state.update(CP.PT_INPUTS, P, T) + return _he_state.umass() -def he_T_from_u(u): - """Recover He temperature from specific internal energy [K].""" - return (u - _HE_U_REF) / _HE_CV +def he_T_from_u(u, P=_HE_P_REF): + """Recover He temperature from pressure and specific internal energy [K].""" + _he_state.update(CP.PUmass_INPUTS, P, u) + return _he_state.T() -# --------------------------------------------------------------------------- -# Lookup table for LN2: u(T) -> T at P = 0.17 MPa (built at import time) -# --------------------------------------------------------------------------- -_LN2_T_MIN = 65.0 -_LN2_T_MAX = 82.0 # stay below saturation at 0.17 MPa (~82.03 K) -_LN2_TABLE_N = 200 -_P_WORK = 170000.0 +def he_T_from_rho(P, rho): + """Recover He temperature from pressure and density [K].""" + _he_state.update(CP.DmassP_INPUTS, rho, P) + return _he_state.T() -_ln2_T_table = np.linspace(_LN2_T_MIN, _LN2_T_MAX, _LN2_TABLE_N) -_ln2_u_table = np.array([ln2_u(T, _P_WORK) for T in _ln2_T_table]) -# _ln2_u_table is monotonically increasing, so np.interp works for inverse lookup + +def he_u_from_rho(P, rho): + """Recover He specific internal energy from pressure and density [J/kg].""" + _he_state.update(CP.DmassP_INPUTS, rho, P) + return _he_state.umass() diff --git a/src/cryo_tank/solver.py b/src/cryo_tank/solver.py index 90fa3a3..a0fe484 100644 --- a/src/cryo_tank/solver.py +++ b/src/cryo_tank/solver.py @@ -80,7 +80,6 @@ def run(tank, t_end, rtol=1e-8, atol=1e-10, max_step=10.0): 'V_ull': np.zeros(n), 'liquid_level': np.zeros(n), 'fill_fraction': np.zeros(n), - 'P_N2': np.zeros(n), 'P_He': np.zeros(n), 'P_total': np.zeros(n), 'Q_leak': np.zeros(n), @@ -100,9 +99,8 @@ def run(tank, t_end, rtol=1e-8, atol=1e-10, max_step=10.0): history['V_ull'][i] = info['V_ull'] history['liquid_level'][i] = info['liquid_level'] history['fill_fraction'][i] = info['fill_fraction'] - history['P_N2'][i] = info['P_N2'] history['P_He'][i] = info['P_He'] - history['P_total'][i] = info['P_N2'] + info['P_He'] + history['P_total'][i] = info['P_He'] # Recompute heat terms for recording T_liq = info['T_liq'] @@ -120,11 +118,9 @@ def run(tank, t_end, rtol=1e-8, atol=1e-10, max_step=10.0): history['Q_leak'][i] = Q_leak history['Q_leak_liq'][i] = Q_leak_liq history['Q_leak_ull'][i] = Q_leak_ull + history['mdot_He'][i] = tank._solve_he_flow_rate( + info, Q_liq_to_ull, Q_leak_ull + ) - # He flow rate via finite difference on m_He (post-processing only, not in RHS) - dt = np.diff(t) - dm_He = np.diff(history['m_He']) - history['mdot_He'][0] = dm_He[0] / dt[0] if len(dt) > 0 else 0.0 - history['mdot_He'][1:] = dm_He / dt return history diff --git a/src/cryo_tank/tank_model.py b/src/cryo_tank/tank_model.py index 95c311b..d315d3c 100644 --- a/src/cryo_tank/tank_model.py +++ b/src/cryo_tank/tank_model.py @@ -10,8 +10,8 @@ import math import warnings import numpy as np +from scipy.optimize import brentq -from cryo_tank.config import R_HE, R_N2 from cryo_tank import properties as prop @@ -77,21 +77,14 @@ class CryoTank: self._m_liq_0 = rho_liq_0 * V_liq_0 self._U_liq_0 = self._m_liq_0 * prop.ln2_u(T_init, P_work) - # Ullage initial state: N2 vapor at saturation + He to fill pressure - # Use ideal gas law for N2 mass (consistent with derive/rhs which - # treat ullage N2 as ideal gas: P_N2 = m_N2 * R_N2 * T / V) - P_N2_0 = prop.n2_sat_pressure(T_init) - self.m_N2_ull = P_N2_0 * V_ull_0 / (R_N2 * T_init) # FIXED for all time + # Ullage initial state: helium pressurization only. + # Nitrogen evaporation is intentionally not modeled, so the ullage + # pressure is provided entirely by helium and the liquid sees the same + # tank pressure for property lookup. + self._m_He_0 = prop.he_rho(T_init, P_work) * V_ull_0 - P_He_0 = P_work - P_N2_0 - self._m_He_0 = P_He_0 * V_ull_0 / (R_HE * T_init) - - U_N2_ull_0 = self.m_N2_ull * prop.n2_vapor_u(T_init) - U_He_0 = self._m_He_0 * prop.he_u(T_init) - self._U_ull_0 = U_N2_ull_0 + U_He_0 - - # Saturation temperature warning threshold - self._T_sat = 82.0 # approximate, from CoolProp: ~82.03 K at 0.17 MPa + U_He_0 = self._m_He_0 * prop.he_u(T_init, P_work) + self._U_ull_0 = U_He_0 def initial_state(self): """Return the ODE initial state vector y0 = [m_liq, U_liq, U_ull].""" @@ -108,7 +101,8 @@ class CryoTank: """Compute all derived quantities from state vector y. Returns a dict with T_liq, T_ull, V_liq, V_ull, liquid_level, - fill_fraction, m_He, P_N2, P_He, etc. + fill_fraction, m_He, P_He, etc. Helium provides the full ullage + pressure in this no-evaporation model. """ m_liq, U_liq, U_ull = y[0], y[1], y[2] @@ -126,62 +120,53 @@ class CryoTank: # Solve T_ull from ullage energy T_ull = self._solve_ullage_temperature(U_ull, V_ull) - # Recompute P_N2 and m_He with correct T_ull - P_N2 = self.m_N2_ull * R_N2 * T_ull / V_ull - P_He = self.P_work - P_N2 - m_He = P_He * V_ull / (R_HE * T_ull) + # No nitrogen evaporation: helium provides all gas pressure. + P_He = self.P_work + m_He = prop.he_rho(T_ull, P_He) * V_ull return { 'T_liq': T_liq, 'T_ull': T_ull, 'V_liq': V_liq, 'V_ull': V_ull, 'liquid_level': liquid_level, 'fill_fraction': fill_fraction, 'rho_liq': rho_liq, - 'm_He': m_He, 'P_N2': P_N2, 'P_He': P_He, + 'm_He': m_He, 'P_He': P_He, } def _solve_ullage_temperature(self, U_ull, V_ull): - """Solve for T_ull given total ullage internal energy and volume. + """Solve pure-He ullage temperature using CoolProp EOS. - U_ull = m_N2_ull * u_N2_vap(T) + m_He(T) * u_He(T) - where m_He(T) = (P_work - m_N2_ull * R_N2 * T / V_ull) * V_ull / (R_He * T) - - Solved by Newton iteration. + At fixed tank pressure, T is obtained from: + U_ull = rho_He(P_work, T) * V_ull * u_He(P_work, T) """ - T = self._T_init # initial guess - for _ in range(50): - P_N2 = self.m_N2_ull * R_N2 * T / V_ull - P_He = self.P_work - P_N2 - if P_He < 0: - P_He = 0.0 - m_He = P_He * V_ull / (R_HE * T) + def residual(T): + rho = prop.he_rho(T, self.P_work) + u = prop.he_u(T, self.P_work) + return rho * V_ull * u - U_ull - # N2 vapor internal energy (ideal gas approx): u = cv_N2 * T + u_ref - cv_N2 = 743.0 - u_N2_ref = prop.n2_vapor_u(78.0) - cv_N2 * 78.0 - u_N2 = cv_N2 * T + u_N2_ref + T_grid = np.geomspace(20.0, 1000.0, 160) + T_prev = T_grid[0] + f_prev = residual(T_prev) + if abs(f_prev) < 1e-8: + return float(T_prev) - U_calc = self.m_N2_ull * u_N2 + m_He * prop.he_u(T) - residual = U_calc - U_ull + best_T = T_prev + best_abs_f = abs(f_prev) + for T in T_grid[1:]: + f = residual(T) + if abs(f) < best_abs_f: + best_T = T + best_abs_f = abs(f) + if f_prev * f <= 0.0: + return brentq(residual, T_prev, T, xtol=1e-8, rtol=1e-10) + T_prev = T + f_prev = f - if abs(residual) < 1e-3: # converged (< 1 mJ) - return T - - # Numerical derivative - dT = 0.01 - P_N2_p = self.m_N2_ull * R_N2 * (T + dT) / V_ull - P_He_p = max(0.0, self.P_work - P_N2_p) - m_He_p = P_He_p * V_ull / (R_HE * (T + dT)) - u_N2_p = cv_N2 * (T + dT) + u_N2_ref - U_calc_p = self.m_N2_ull * u_N2_p + m_He_p * prop.he_u(T + dT) - dU_dT = (U_calc_p - U_calc) / dT - - if abs(dU_dT) < 1e-20: - break - T = T - residual / dU_dT - T = max(50.0, min(T, 400.0)) # clamp - - warnings.warn(f"_solve_ullage_temperature did not converge: T={T:.2f}, residual={residual:.2e}") - return T + warnings.warn( + "Unable to bracket ullage temperature from CoolProp He state: " + f"U_ull={U_ull:.6g}, V_ull={V_ull:.6g}; " + f"using nearest T={best_T:.6g} K" + ) + return float(best_T) def rhs(self, t, y): """ODE right-hand side: dy/dt = [dm_liq/dt, dU_liq/dt, dU_ull/dt]. @@ -220,58 +205,52 @@ class CryoTank: # --- Ullage zone --- dU_ull_dt = mdot_He * self.h_in_he + Q_liq_to_ull + Q_leak_ull - # --- Warnings --- - if T_liq > self._T_sat - 1.0: - warnings.warn( - f"T_liq={T_liq:.2f}K approaching saturation ({self._T_sat:.1f}K); " - "evaporation effects may be significant." - ) - return np.array([dm_liq_dt, dU_liq_dt, dU_ull_dt]) def _solve_he_flow_rate(self, info, Q_liq_to_ull, Q_leak_ull): - """Analytically solve for m_dot_He from dP/dt = 0 constraint. + """Solve He inlet mass flow from constant-pressure CoolProp EOS. - The key equation: P_total = P_N2 + P_He = const. - - P_N2 = m_N2_ull * R_N2 * T_ull / V_ull (N2 ideal gas, m_N2_ull = const) - P_He = m_He * R_HE * T_ull / V_ull - - dP/dt = 0 implies dP_He/dt = -dP_N2/dt. - - Expanding and solving for m_dot_He gives a linear equation. - See spec Section 2.6 for full derivation. + The ullage is pure He at P_work. With V_ull changing due to liquid + volume change, m_dot_He is found by differentiating: + m(T, V) = rho(P_work, T) * V + U(T, V) = m(T, V) * u(P_work, T) + and combining it with dU/dt = m_dot_He*h_in + Q_ull. """ T_ull = info['T_ull'] V_ull = info['V_ull'] - m_He = info['m_He'] rho_liq = info['rho_liq'] - # dV_ull/dt = -dV_liq/dt = -dm_liq/dt / rho_liq = -(mdot_in - mdot_out) / rho_liq - dV_ull_dt = -self.dm_liq_dt / rho_liq # positive when liquid drains - - # Total ullage heat input (excluding He inlet, which we're solving for) + dV_ull_dt = -self.dm_liq_dt / rho_liq Q_ull_no_he = Q_liq_to_ull + Q_leak_ull - # Ullage total cv*mass (for dT_ull/dt estimation) - cv_N2 = 743.0 # J/(kg*K), N2 vapor - cv_He = prop.he_cv() - C_ull = self.m_N2_ull * cv_N2 + m_He * cv_He # total heat capacity [J/K] + def state_at(T, V): + rho = prop.he_rho(T, self.P_work) + m = rho * V + U = m * prop.he_u(T, self.P_work) + return m, U - R_mix = self.m_N2_ull * R_N2 + m_He * R_HE # effective "mR" [J/K] + dT = max(1e-3, abs(T_ull) * 1e-5) + dV = max(1e-8, abs(V_ull) * 1e-5) - # Coefficient of m_dot_He in dP/dt = 0: - # A * m_dot_He + B = 0 - # m_dot_He = -B / A - A = R_HE * T_ull / V_ull + R_mix / (V_ull * C_ull) * self.h_in_he - B = R_mix / (V_ull * C_ull) * Q_ull_no_he - R_mix * T_ull / (V_ull ** 2) * dV_ull_dt + m_plus, U_plus = state_at(T_ull + dT, V_ull) + m_minus, U_minus = state_at(max(20.0, T_ull - dT), V_ull) + actual_dT = (T_ull + dT) - max(20.0, T_ull - dT) + m_T = (m_plus - m_minus) / actual_dT + U_T = (U_plus - U_minus) / actual_dT - if abs(A) < 1e-30: + m_plus, U_plus = state_at(T_ull, V_ull + dV) + m_minus, U_minus = state_at(T_ull, max(1e-12, V_ull - dV)) + actual_dV = (V_ull + dV) - max(1e-12, V_ull - dV) + m_V = (m_plus - m_minus) / actual_dV + U_V = (U_plus - U_minus) / actual_dV + + denominator = U_T - self.h_in_he * m_T + if abs(denominator) < 1e-30: return 0.0 - mdot_He = -B / A + dT_dt = (Q_ull_no_he - (U_V - self.h_in_he * m_V) * dV_ull_dt) / denominator + mdot_He = m_T * dT_dt + m_V * dV_ull_dt - # Clamp: He can only flow in (strict mode per spec Section 10.3) if mdot_He < 0: warnings.warn( f"He backflow requested (m_dot_He={mdot_He:.4e} kg/s); " diff --git a/tests/cryo_tank/test_properties.py b/tests/cryo_tank/test_properties.py index 984c4f6..dadf977 100644 --- a/tests/cryo_tank/test_properties.py +++ b/tests/cryo_tank/test_properties.py @@ -6,7 +6,6 @@ sys.path.insert(0, "src") from cryo_tank.properties import ( ln2_rho, ln2_h, ln2_u, ln2_T_from_u, - n2_vapor_u, n2_sat_pressure, he_u, he_h, he_cp, he_cv, ) from cryo_tank.config import P_WORKING @@ -34,24 +33,12 @@ class TestLN2Properties: assert abs(T_recovered - T_orig) < 0.01 -class TestN2VaporProperties: - """N2 vapor properties.""" - - def test_n2_sat_pressure_at_78K(self): - P_sat = n2_sat_pressure(78.0) - assert 0.10e6 < P_sat < 0.12e6 # ~0.1093 MPa - - def test_n2_vapor_internal_energy_at_78K(self): - u = n2_vapor_u(78.0) - assert 50000 < u < 60000 # ~55547 J/kg - - class TestHeliumProperties: - """Helium (ideal gas) properties.""" + """Helium properties.""" def test_he_cp_near_5196(self): cp = he_cp() - assert abs(cp - 5196.2) < 10 # monatomic ideal gas + assert abs(cp - 5196.2) < 10 def test_he_cv_near_3117(self): cv = he_cv() diff --git a/tests/cryo_tank/test_tank_model.py b/tests/cryo_tank/test_tank_model.py index 3a821fa..fef21e9 100644 --- a/tests/cryo_tank/test_tank_model.py +++ b/tests/cryo_tank/test_tank_model.py @@ -69,10 +69,9 @@ class TestInitialState: assert abs(info['T_liq'] - T_INIT) < 0.1 assert abs(info['T_ull'] - T_INIT) < 1.0 - def test_initial_pressure_components_sum_to_P_working(self): + def test_initial_pressure_is_provided_entirely_by_helium(self): tank = _make_tank() y0 = tank.initial_state() info = tank.derive(y0) - P_N2 = info['P_N2'] P_He = info['P_He'] - assert abs(P_N2 + P_He - P_WORKING) / P_WORKING < 1e-6 + assert abs(P_He - P_WORKING) / P_WORKING < 1e-12