fix(cryo_tank): use helium-only CoolProp pressure model

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。
 - 建议每次修改参数后保存对应工况说明，避免不同仿真结果混在同一个输出目录中。

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)

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()
 

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]')

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()

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

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); "