diff --git a/src/cryo_tank/Index.md b/src/cryo_tank/Index.md
new file mode 100644
index 0000000..4795ff7
--- /dev/null
+++ b/src/cryo_tank/Index.md
@@ -0,0 +1,491 @@
+# cryo_tank 模型方程与 ODE 求解逻辑
+
+本文只保留当前 `src/cryo_tank/` 模型内部实际使用的方程，以及 ODE 方程组的求解逻辑。当前版本采用方案二：把气枕温度 `T_ull` 作为 ODE 状态量，氦气质量 `m_He` 由恒压状态方程派生。
+
+## 1. 模型假设
+
+当前模型是液氮贮箱的两区集总参数模型：
+
+- 液相区为液氮。
+- 气枕区为纯氦气。
+- 气枕压力固定为工作压力 `P_work`。
+- 不考虑氮气蒸发、冷凝和气液相间质量传递。
+- 氦气质量不是独立状态量，而是维持恒压所需的派生量。
+- 氦气只允许补入，不允许倒流；若计算得到负的氦气流量，则报告流量钳制为 0。
+
+对应代码：
+
+- 状态和派生量：`src/cryo_tank/tank_model.py:98-129`
+- ODE 右端项：`src/cryo_tank/tank_model.py:131-163`
+- 氦气补气温度方程：`src/cryo_tank/tank_model.py:165-214`
+
+## 2. 状态向量
+
+ODE 状态向量为：
+
+```text
+y = [m_liq, U_liq, T_ull]
+```
+
+其中：
+
+```text
+m_liq : 液氮质量 [kg]
+U_liq : 液相总内能 [J]
+T_ull : 气枕区温度 [K]
+```
+
+代码位置：`src/cryo_tank/tank_model.py:1-7`，`src/cryo_tank/tank_model.py:87-89`
+
+## 3. 初始状态方程
+
+初始液相体积：
+
+```text
+V_liq_0 = (1 - ullage_fraction) * V_total
+```
+
+初始气枕体积：
+
+```text
+V_ull_0 = ullage_fraction * V_total
+```
+
+初始液氮质量：
+
+```text
+rho_liq_0 = rho_LN2(T_init, P_work)
+m_liq_0 = rho_liq_0 * V_liq_0
+```
+
+初始液相总内能：
+
+```text
+U_liq_0 = m_liq_0 * u_LN2(T_init, P_work)
+```
+
+初始气枕温度：
+
+```text
+T_ull_0 = T_init
+```
+
+初始氦气质量是派生量：
+
+```text
+m_He_0 = rho_He(T_init, P_work) * V_ull_0
+```
+
+代码位置：`src/cryo_tank/tank_model.py:71-85`
+
+## 4. 几何方程
+
+横截面积：
+
+```text
+A_cross = V_total / H_tank
+```
+
+圆柱直径：
+
+```text
+D = sqrt(4 * A_cross / pi)
+```
+
+侧面积：
+
+```text
+A_side = pi * D * H_tank
+```
+
+端盖面积：
+
+```text
+A_cap = A_cross
+```
+
+总面积：
+
+```text
+A_total = A_side + 2 * A_cap
+```
+
+代码位置：`src/cryo_tank/tank_model.py:41-48`
+
+液位对应的湿壁和干壁面积：
+
+```text
+level = clamp(liquid_level, 0, H_tank)
+A_wet = A_cap + pi * D * level
+A_dry = A_cap + pi * D * (H_tank - level)
+```
+
+代码位置：`src/cryo_tank/tank_model.py:91-96`
+
+## 5. 派生量方程
+
+给定状态：
+
+```text
+y = [m_liq, U_liq, T_ull]
+```
+
+### 5.1 液相派生量
+
+```text
+u_liq = U_liq / m_liq
+T_liq = T_LN2_from_u(u_liq, P_work)
+rho_liq = rho_LN2(T_liq, P_work)
+V_liq = m_liq / rho_liq
+liquid_level = V_liq / A_cross
+fill_fraction = liquid_level / H_tank
+```
+
+代码位置：`src/cryo_tank/tank_model.py:107-113`
+
+### 5.2 气枕派生量
+
+```text
+V_ull = V_total - V_liq
+P_He = P_work
+rho_He = rho_He(T_ull, P_He)
+m_He = rho_He * V_ull
+U_ull = m_He * u_He(T_ull, P_He)
+```
+
+代码位置：`src/cryo_tank/tank_model.py:115-121`
+
+## 6. 换热方程
+
+### 6.1 液相与气枕界面换热
+
+```text
+Q_liq_to_ull = h_conv * A_cross * (T_liq - T_ull)
+```
+
+符号约定：
+
+```text
+Q_liq_to_ull > 0 : 热量从液相传给气枕
+Q_liq_to_ull < 0 : 热量从气枕传给液相
+```
+
+代码位置：`src/cryo_tank/tank_model.py:141-142`
+
+### 6.2 外界总漏热
+
+当前主程序默认使用 MLI 漏热模型：
+
+```text
+Q_leak = A_total * q_mli
+```
+
+代码位置：
+
+- 默认创建：`src/cryo_tank/main.py:32-33`
+- MLI 方程：`src/cryo_tank/heat_leak.py:21-40`
+
+代码中还提供 Foam 漏热模型：
+
+```text
+T_mean = (T_inner + T_env) / 2
+Q_leak = A_total * k_eff(T_mean) * (T_env - T_inner) / delta
+```
+
+如果 `k_eff` 是常数，则直接使用该常数。
+
+代码位置：`src/cryo_tank/heat_leak.py:43-70`
+
+### 6.3 漏热分配
+
+总漏热按湿壁面积和干壁面积分配：
+
+```text
+Q_leak_liq = Q_leak * A_wet / (A_wet + A_dry)
+Q_leak_ull = Q_leak * A_dry / (A_wet + A_dry)
+```
+
+代码位置：`src/cryo_tank/tank_model.py:144-148`
+
+## 7. 液氮质量方程
+
+液氮质量变化率：
+
+```text
+dm_liq/dt = mdot_in_ln2 - mdot_out_ln2
+```
+
+代码位置：
+
+- 净流量预计算：`src/cryo_tank/tank_model.py:64-65`
+- ODE 中使用：`src/cryo_tank/tank_model.py:150-151`
+
+## 8. 液相能量方程
+
+入口液氮焓：
+
+```text
+h_in_ln2 = h_LN2(T_in_ln2, P_work)
+```
+
+当前液相出口焓：
+
+```text
+h_liq = h_LN2(T_liq, P_work)
+```
+
+液相总内能变化率：
+
+```text
+dU_liq/dt =
+    mdot_in_ln2 * h_in_ln2
+  - mdot_out_ln2 * h_liq
+  - Q_liq_to_ull
+  + Q_leak_liq
+```
+
+代码位置：
+
+- 入口焓预计算：`src/cryo_tank/tank_model.py:60-62`
+- 能量方程：`src/cryo_tank/tank_model.py:152-156`
+
+## 9. 气枕温度和氦气质量导数方程
+
+当前模型把 `T_ull` 作为 ODE 状态。氦气质量由恒压关系派生：
+
+```text
+m_He = rho_He(T_ull, P_work) * V_ull
+```
+
+代码位置：`src/cryo_tank/tank_model.py:115-121`
+
+### 9.1 气枕体积变化率
+
+液氮质量变化导致气枕体积变化：
+
+```text
+dV_ull/dt = -dm_liq_dt / rho_liq
+```
+
+代码位置：`src/cryo_tank/tank_model.py:186-188`
+
+### 9.2 气枕区非氦气入口热量
+
+```text
+Q_ull_no_he = Q_liq_to_ull + Q_leak_ull
+```
+
+代码位置：`src/cryo_tank/tank_model.py:188`
+
+### 9.3 恒压氦气状态函数
+
+在固定压力 `P_work` 下，给定 `T_ull` 和 `V_ull`：
+
+```text
+rho = rho_He(T_ull, P_work)
+m(T, V) = rho * V
+U(T, V) = m(T, V) * u_He(T_ull, P_work)
+```
+
+代码位置：`src/cryo_tank/tank_model.py:115-121`
+
+### 9.4 解析偏导数
+
+代码不再对温度和体积做有限差分，而是使用恒压物性导数：
+
+```text
+drhodT_P = (partial rho / partial T)_P
+dudT_P = (partial u / partial T)_P
+```
+
+CoolProp 偏导接口位置：`src/cryo_tank/properties.py:100-109`
+
+由此得到：
+
+```text
+m_T = V_ull * drhodT_P
+m_V = rho
+U_T = V_ull * (u * drhodT_P + rho * dudT_P)
+U_V = rho * u
+```
+
+代码位置：`src/cryo_tank/tank_model.py:165-178`
+
+### 9.5 气枕温度变化率
+
+氦气入口焓：
+
+```text
+h_in_he = h_He(T_in_he, P_work)
+```
+
+气枕温度变化率：
+
+```text
+dT_ull/dt =
+    [Q_ull_no_he - (U_V + P_work - h_in_he * m_V) * dV_ull/dt]
+    / [U_T - h_in_he * m_T]
+```
+
+代码位置：`src/cryo_tank/tank_model.py:190-197`
+
+若分母过小：
+
+```text
+if abs(U_T - h_in_he * m_T) < 1e-30:
+    dT_ull/dt = 0
+    mdot_He = 0
+```
+
+代码位置：`src/cryo_tank/tank_model.py:190-193`
+
+### 9.6 氦气质量变化率
+
+```text
+dm_He/dt = mdot_He
+```
+
+其中：
+
+```text
+mdot_He = m_T * dT_ull/dt + m_V * dV_ull/dt
+```
+
+代码位置：`src/cryo_tank/tank_model.py:198`
+
+若 `mdot_He < 0`，报告流量钳制为 0：
+
+```text
+mdot_He = 0
+```
+
+代码位置：`src/cryo_tank/tank_model.py:200-205`
+
+## 10. 完整 ODE 方程组
+
+状态向量：
+
+```text
+y = [m_liq, U_liq, T_ull]
+```
+
+ODE 方程组：
+
+```text
+dm_liq/dt = mdot_in_ln2 - mdot_out_ln2
+```
+
+```text
+dU_liq/dt =
+    mdot_in_ln2 * h_in_ln2
+  - mdot_out_ln2 * h_liq
+  - Q_liq_to_ull
+  + Q_leak_liq
+```
+
+```text
+dT_ull/dt =
+    [Q_ull_no_he - (U_V + P_work - h_in_he * m_V) * dV_ull/dt]
+    / [U_T - h_in_he * m_T]
+```
+
+其中：
+
+```text
+h_liq = h_LN2(T_liq, P_work)
+Q_liq_to_ull = h_conv * A_cross * (T_liq - T_ull)
+Q_ull_no_he = Q_liq_to_ull + Q_leak_ull
+Q_leak_liq = Q_leak * A_wet / (A_wet + A_dry)
+Q_leak_ull = Q_leak * A_dry / (A_wet + A_dry)
+```
+
+氦气质量和质量流量是派生量：
+
+```text
+m_He = rho_He(T_ull, P_work) * V_ull
+mdot_He = m_T * dT_ull/dt + m_V * dV_ull/dt
+```
+
+代码位置：`src/cryo_tank/tank_model.py:131-214`
+
+## 11. ODE 求解逻辑
+
+求解器入口：
+
+```text
+run(tank, t_end, rtol=1e-8, atol=1e-10, max_step=10.0)
+```
+
+代码位置：`src/cryo_tank/solver.py:22-24`
+
+求解初值：
+
+```text
+y0 = tank.initial_state()
+```
+
+代码位置：`src/cryo_tank/solver.py:43`
+
+调用 `solve_ivp`：
+
+```text
+solve_ivp(
+    tank.rhs,
+    [0.0, t_end],
+    y0,
+    method="RK45",
+    rtol=rtol,
+    atol=atol,
+    max_step=max_step,
+    events=[_liquid_empty_event],
+    dense_output=True,
+)
+```
+
+代码位置：`src/cryo_tank/solver.py:45-55`
+
+液体排空事件：
+
+```text
+event(t, y) = m_liq
+```
+
+事件属性：
+
+```text
+terminal = True
+direction = -1
+```
+
+含义：
+
+- 当液氮质量下降到 0 时终止求解。
+- 只检测从正到零的方向。
+
+代码位置：`src/cryo_tank/solver.py:14-19`
+
+求解失败时抛出错误；若液体提前排空，则发出警告。
+
+代码位置：`src/cryo_tank/solver.py:57-64`
+
+## 12. 后处理逻辑
+
+求解完成后，`solver.run()` 对每个输出时刻执行：
+
+```text
+info = tank.derive(y_i)
+```
+
+并重新计算：
+
+```text
+Q_liq_to_ull
+Q_leak
+Q_leak_liq
+Q_leak_ull
+mdot_He
+```
+
+其中 `m_He`、`U_ull`、`mdot_He` 均为后处理派生结果，不是独立 ODE 状态。
+
+代码位置：`src/cryo_tank/solver.py:70-124`
diff --git a/src/cryo_tank/Index.pdf b/src/cryo_tank/Index.pdf
new file mode 100644
index 0000000..6f28992
Binary files /dev/null and b/src/cryo_tank/Index.pdf differ
diff --git a/src/cryo_tank/README.md b/src/cryo_tank/README.md
index ee16952..51f64b6 100644
--- a/src/cryo_tank/README.md
+++ b/src/cryo_tank/README.md
@@ -96,16 +96,16 @@ find results/cryo_tank -maxdepth 1 -type f -print
 状态量为：
 
 ```text
-y = [m_liq, U_liq, U_ull]
+y = [m_liq, U_liq, T_ull]
 ```
 
 含义分别是：
 
 - `m_liq`：液氮质量
 - `U_liq`：液相内能
-- `U_ull`：气枕区总内能
+- `T_ull`：气枕区温度
 
-`derive(y)` 会根据状态量计算温度、体积、液位、充满率、氦气质量和氦气气枕压力等派生量。
+`derive(y)` 会根据状态量计算温度、体积、液位、充满率、氦气质量、气枕区总内能和氦气气枕压力等派生量。当前模型按恒压纯氦气枕处理，`m_He = rho_He(T_ull, P_work) * V_ull`，所以氦气质量是派生量，不是独立 ODE 状态。
 
 `rhs(t, y)` 是 ODE 右端函数，由 `scipy.integrate.solve_ivp` 调用。
 
@@ -132,9 +132,11 @@ run(tank, t_end, rtol=1e-8, atol=1e-10, max_step=10.0)
 主要功能：
 
 - 液氮密度、焓、内能计算。
-- 氮气饱和压力和饱和蒸气内能计算。
+- 根据液氮比内能和压力反推液相温度。
 - 氦气密度、焓、内能、比热计算。
-- 建立液氮内能到温度的查表插值，加速 ODE 求解。
+- 根据氦气压力和密度反推气枕温度。
+- 计算氦气恒压密度温度导数 `(partial rho / partial T)_P`。
+- 计算氦气恒压内能温度导数 `(partial u / partial T)_P`。
 
 该模块依赖 `CoolProp`。
 
@@ -373,6 +375,7 @@ pytest -q
 
 - 本模型当前将储箱分为液相区和气枕区两个区域，不是完整 CFD 模型。
 - 气枕区按纯氦气处理，物性由 CoolProp 计算。
+- 气枕温度 `T_ull` 是 ODE 状态量；氦气质量 `m_He` 和补气流量 `mdot_He` 是恒压条件下的派生结果。
 - 液氮物性依赖 `CoolProp`，缺少该包会导致程序无法启动。
 - `results/` 目录默认被 `.gitignore` 忽略，仿真输出不会自动上传到 Git。
 - 建议每次修改参数后保存对应工况说明，避免不同仿真结果混在同一个输出目录中。
diff --git a/src/cryo_tank/config.py b/src/cryo_tank/config.py
index 56be7bc..b634228 100644
--- a/src/cryo_tank/config.py
+++ b/src/cryo_tank/config.py
@@ -3,37 +3,38 @@
 Configuration for the cryogenic LN2 tank simulation.
 Pure data module -- no functions, no side effects.
 """
+
 import math
 
 # ---------- Tank geometry ----------
-V_TOTAL = 420.1e-3          # m^3  (420.1 L)
-H_TANK = 0.5                # m    (cylinder height)
+V_TOTAL = 420.1e-3  # m^3  (420.1 L)
+H_TANK = 0.5  # m    (cylinder height)
 A_CROSS = V_TOTAL / H_TANK  # m^2  (cross-section area)
 D_TANK = math.sqrt(4 * A_CROSS / math.pi)  # m (diameter)
-A_SIDE = math.pi * D_TANK * H_TANK         # m^2 (side wall)
-A_CAP = A_CROSS                             # m^2 (top or bottom cap)
-A_TOTAL = A_SIDE + 2 * A_CAP               # m^2 (total surface)
+A_SIDE = math.pi * D_TANK * H_TANK  # m^2 (side wall)
+A_CAP = A_CROSS  # m^2 (top or bottom cap)
+A_TOTAL = A_SIDE + 2 * A_CAP  # m^2 (total surface)
 
 # ---------- Tank limits ----------
-P_WORKING = 0.17e6          # Pa  (working pressure, absolute)
-P_MAX = 0.8e6               # Pa  (max bearing pressure)
+P_WORKING = 0.17e6  # Pa  (working pressure, absolute)
+P_MAX = 0.8e6  # Pa  (max bearing pressure)
 
 # ---------- Initial conditions ----------
-T_INIT = 78.0               # K
-ULLAGE_FRACTION = 0.30      # gas pocket = 30% of V_TOTAL
+T_INIT = 78.0  # K
+ULLAGE_FRACTION = 0.30  # gas pocket = 30% of V_TOTAL
 
 # ---------- Inlet / outlet ----------
-MDOT_IN_LN2 = 1.144         # kg/s
-T_IN_LN2 = 77.0             # K
-MDOT_OUT_LN2 = 1.1895       # kg/s
-T_IN_HE = 100.0             # K
+MDOT_IN_LN2 = 1.144  # kg/s
+T_IN_LN2 = 77.0  # K
+MDOT_OUT_LN2 = 1.1895  # kg/s
+T_IN_HE = 100.0  # K
 
 # ---------- Heat transfer ----------
-H_CONV_SURFACE = 50.0       # W/(m^2*K)  liquid-to-ullage surface convection
-T_ENV = 300.0                # K  ambient temperature
+H_CONV_SURFACE = 0.0  # W/(m^2*K)  liquid-to-ullage surface convection
+T_ENV = 300.0  # K  ambient temperature
 
 # ---------- Simulation control ----------
-T_END = 3600.0               # s  (1 hour)
+T_END = 3600.0  # s  (1 hour)
 RTOL = 1e-8
 ATOL = 1e-10
 
diff --git a/src/cryo_tank/properties.py b/src/cryo_tank/properties.py
index 66d1190..3fc02ab 100644
--- a/src/cryo_tank/properties.py
+++ b/src/cryo_tank/properties.py
@@ -96,3 +96,16 @@ 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()
+
+
+def he_drho_dT_const_p(T, P=_HE_P_REF):
+    """He density temperature derivative at constant pressure [kg/(m^3*K)]."""
+    _he_state.update(CP.PT_INPUTS, P, T)
+    return _he_state.first_partial_deriv(CP.iDmass, CP.iT, CP.iP)
+
+
+def he_du_dT_const_p(T, P=_HE_P_REF):
+    """He internal-energy temperature derivative at constant pressure [J/(kg*K)]."""
+    _he_state.update(CP.PT_INPUTS, P, T)
+    return _he_state.first_partial_deriv(CP.iUmass, CP.iT, CP.iP)
+
diff --git a/src/cryo_tank/solver.py b/src/cryo_tank/solver.py
index a0fe484..d69ef57 100644
--- a/src/cryo_tank/solver.py
+++ b/src/cryo_tank/solver.py
@@ -71,7 +71,7 @@ def run(tank, t_end, rtol=1e-8, atol=1e-10, max_step=10.0):
         't': t,
         'm_liq': sol.y[0],
         'U_liq': sol.y[1],
-        'U_ull': sol.y[2],
+        'U_ull': np.zeros(n),
         'T_liq': np.zeros(n),
         'T_ull': np.zeros(n),
         'm_He': np.zeros(n),
@@ -94,6 +94,7 @@ def run(tank, t_end, rtol=1e-8, atol=1e-10, max_step=10.0):
 
         history['T_liq'][i] = info['T_liq']
         history['T_ull'][i] = info['T_ull']
+        history['U_ull'][i] = info['U_ull']
         history['m_He'][i] = info['m_He']
         history['V_liq'][i] = info['V_liq']
         history['V_ull'][i] = info['V_ull']
diff --git a/src/cryo_tank/tank_model.py b/src/cryo_tank/tank_model.py
index d315d3c..352407a 100644
--- a/src/cryo_tank/tank_model.py
+++ b/src/cryo_tank/tank_model.py
@@ -2,7 +2,7 @@
 """
 CryoTank: two-zone (liquid + ullage) cryogenic tank model.
 
-State vector y = [m_liq, U_liq, U_ull] (3 components).
+State vector y = [m_liq, U_liq, T_ull] (3 components).
 Derived quantities (T, V, m_He, etc.) computed by derive(y).
 ODE right-hand side provided by rhs(t, y).
 """
@@ -10,7 +10,6 @@ import math
 import warnings
 
 import numpy as np
-from scipy.optimize import brentq
 
 from cryo_tank import properties as prop
 
@@ -60,7 +59,7 @@ class CryoTank:
 
         # Precompute constant inlet enthalpies
         self.h_in_ln2 = prop.ln2_h(T_in_ln2, P_work)
-        self.h_in_he = prop.he_h(T_in_he)
+        self.h_in_he = prop.he_h(T_in_he, P_work)
 
         # Net liquid flow (constant)
         self.dm_liq_dt = mdot_in_ln2 - mdot_out_ln2
@@ -81,14 +80,13 @@ class CryoTank:
         # 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._T_ull_0 = T_init
         self._m_He_0 = prop.he_rho(T_init, P_work) * V_ull_0
-
-        U_He_0 = self._m_He_0 * prop.he_u(T_init, P_work)
-        self._U_ull_0 = U_He_0
+        self._U_ull_0 = self._m_He_0 * prop.he_u(T_init, P_work)
 
     def initial_state(self):
-        """Return the ODE initial state vector y0 = [m_liq, U_liq, U_ull]."""
-        return np.array([self._m_liq_0, self._U_liq_0, self._U_ull_0])
+        """Return the ODE initial state vector y0 = [m_liq, U_liq, T_ull]."""
+        return np.array([self._m_liq_0, self._U_liq_0, self._T_ull_0])
 
     def wetted_areas(self, liquid_level):
         """Return (A_wet, A_dry) for the given liquid level [m]."""
@@ -101,86 +99,44 @@ 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_He, etc. Helium provides the full ullage
-        pressure in this no-evaporation model.
+        fill_fraction, m_He, U_ull, 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]
+        m_liq, U_liq, T_ull = y[0], y[1], y[2]
 
         # Liquid zone
-        u_liq = U_liq / m_liq  # specific internal energy
+        u_liq = U_liq / m_liq
         T_liq = prop.ln2_T_from_u(u_liq, self.P_work)
         rho_liq = prop.ln2_rho(T_liq, self.P_work)
         V_liq = m_liq / rho_liq
         liquid_level = V_liq / self.A_cross
         fill_fraction = liquid_level / self.H_tank
 
-        # Ullage zone
+        # Ullage zone: pure helium at tank pressure. T_ull is the ODE state;
+        # m_He is the amount required to maintain P_work in the current volume.
         V_ull = self.V_total - V_liq
-
-        # Solve T_ull from ullage energy
-        T_ull = self._solve_ullage_temperature(U_ull, V_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
+        rho_He = prop.he_rho(T_ull, P_He)
+        m_He = rho_He * V_ull
+        U_ull = m_He * prop.he_u(T_ull, P_He)
 
         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_He': P_He,
+            'rho_liq': rho_liq, 'rho_He': rho_He,
+            'm_He': m_He, 'U_ull': U_ull, 'P_He': P_He,
         }
 
-    def _solve_ullage_temperature(self, U_ull, V_ull):
-        """Solve pure-He ullage temperature using CoolProp EOS.
-
-        At fixed tank pressure, T is obtained from:
-            U_ull = rho_He(P_work, T) * V_ull * u_He(P_work, 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
-
-        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)
-
-        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
-
-        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].
+        """ODE right-hand side: dy/dt = [dm_liq/dt, dU_liq/dt, dT_ull/dt].
 
         This is called by scipy.integrate.solve_ivp.
         """
         info = self.derive(y)
-        m_liq = y[0]
         T_liq = info['T_liq']
         T_ull = info['T_ull']
-        V_ull = info['V_ull']
-        rho_liq = info['rho_liq']
         liquid_level = info['liquid_level']
-        m_He = info['m_He']
 
         # --- Heat transfer ---
         Q_liq_to_ull = self.h_conv * self.A_cross * (T_liq - T_ull)
@@ -191,9 +147,6 @@ class CryoTank:
         Q_leak_liq = Q_leak * A_wet / A_total if A_total > 0 else 0.0
         Q_leak_ull = Q_leak * A_dry / A_total if A_total > 0 else 0.0
 
-        # --- He flow rate (analytical, per spec Section 2.6) ---
-        mdot_He = self._solve_he_flow_rate(info, Q_liq_to_ull, Q_leak_ull)
-
         # --- Liquid zone ---
         dm_liq_dt = self.dm_liq_dt  # constant: mdot_in - mdot_out
         h_liq = prop.ln2_h(T_liq, self.P_work)
@@ -203,59 +156,59 @@ class CryoTank:
                      + Q_leak_liq)
 
         # --- Ullage zone ---
-        dU_ull_dt = mdot_He * self.h_in_he + Q_liq_to_ull + Q_leak_ull
+        dT_ull_dt, _ = self._solve_ullage_temperature_rate(
+            info, Q_liq_to_ull, Q_leak_ull
+        )
 
-        return np.array([dm_liq_dt, dU_liq_dt, dU_ull_dt])
+        return np.array([dm_liq_dt, dU_liq_dt, dT_ull_dt])
 
-    def _solve_he_flow_rate(self, info, Q_liq_to_ull, Q_leak_ull):
-        """Solve He inlet mass flow from constant-pressure CoolProp EOS.
-
-        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.
-        """
+    def _he_mass_energy_partials(self, info):
+        """Return local partials for m(T,V) and U(T,V) at constant pressure."""
         T_ull = info['T_ull']
         V_ull = info['V_ull']
-        rho_liq = info['rho_liq']
+        rho = info['rho_He']
+        u = prop.he_u(T_ull, self.P_work)
+        drho_dT = prop.he_drho_dT_const_p(T_ull, self.P_work)
+        du_dT = prop.he_du_dT_const_p(T_ull, self.P_work)
 
+        m_T = V_ull * drho_dT
+        m_V = rho
+        U_T = V_ull * (u * drho_dT + rho * du_dT)
+        U_V = rho * u
+        return m_T, U_T, m_V, U_V
+
+    def _solve_ullage_temperature_rate(self, info, Q_liq_to_ull, Q_leak_ull):
+        """Solve dT_ull/dt and He inlet flow from constant-pressure EOS.
+
+        T_ull is the ODE state. The pure-He ullage is maintained at P_work,
+        so m_He = rho(P_work, T_ull) * V_ull is a derived quantity.
+        """
+        rho_liq = info['rho_liq']
         dV_ull_dt = -self.dm_liq_dt / rho_liq
         Q_ull_no_he = Q_liq_to_ull + Q_leak_ull
 
-        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
-
-        dT = max(1e-3, abs(T_ull) * 1e-5)
-        dV = max(1e-8, abs(V_ull) * 1e-5)
-
-        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
-
-        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
-
+        m_T, U_T, m_V, U_V = self._he_mass_energy_partials(info)
         denominator = U_T - self.h_in_he * m_T
         if abs(denominator) < 1e-30:
-            return 0.0
+            return 0.0, 0.0
 
-        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
+        dT_ull_dt = (Q_ull_no_he
+                     - (U_V + self.P_work - self.h_in_he * m_V)
+                     * dV_ull_dt) / denominator
+        mdot_He = m_T * dT_ull_dt + m_V * dV_ull_dt
 
         if mdot_He < 0:
             warnings.warn(
                 f"He backflow requested (m_dot_He={mdot_He:.4e} kg/s); "
-                "clamping to 0. Pressure may drift above target."
+                "clamping reported flow to 0. Pressure control may be invalid."
             )
             mdot_He = 0.0
 
+        return dT_ull_dt, mdot_He
+
+    def _solve_he_flow_rate(self, info, Q_liq_to_ull, Q_leak_ull):
+        """Return derived He inlet mass flow for the current state [kg/s]."""
+        _, mdot_He = self._solve_ullage_temperature_rate(
+            info, Q_liq_to_ull, Q_leak_ull
+        )
         return mdot_He