diff --git a/src/cryo_tank/Index.md b/src/cryo_tank/Index.md index 4795ff7..0a1a638 100644 --- a/src/cryo_tank/Index.md +++ b/src/cryo_tank/Index.md @@ -1,7 +1,3 @@ -# cryo_tank 模型方程与 ODE 求解逻辑 - -本文只保留当前 `src/cryo_tank/` 模型内部实际使用的方程，以及 ODE 方程组的求解逻辑。当前版本采用方案二：把气枕温度 `T_ull` 作为 ODE 状态量，氦气质量 `m_He` 由恒压状态方程派生。 - ## 1. 模型假设 当前模型是液氮贮箱的两区集总参数模型： @@ -15,25 +11,29 @@ 对应代码： -- 状态和派生量：`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` +- 状态和派生量：`src/cryo_tank/tank_model.py:98-130` +- ODE 右端项：`src/cryo_tank/tank_model.py:132-162` +- 液相温度、气枕体积和氦气补气方程：`src/cryo_tank/tank_model.py:172-257` ## 2. 状态向量 ODE 状态向量为： -```text -y = [m_liq, U_liq, T_ull] -``` +$$ +\mathbf{y}=\left[m_{\mathrm{liq}},\ U_{\mathrm{liq}},\ T_{\mathrm{ull}}\right] +$$ 其中： -```text -m_liq : 液氮质量 [kg] -U_liq : 液相总内能 [J] -T_ull : 气枕区温度 [K] -``` +$$ +\begin{aligned} +m_{\mathrm{liq}} &: \text{液氮质量}\ [\mathrm{kg}] \\ +U_{\mathrm{liq}} &: \text{液相总内能}\ [\mathrm{J}] \\ +T_{\mathrm{ull}} &: \text{气枕区温度}\ [\mathrm{K}] +\end{aligned} +$$ + +式中，`m_liq` 为液氮质量，`U_liq` 为液相总内能，`T_ull` 为氦气气枕温度。`m_liq` 和 `U_liq` 是液相守恒方程的状态量，`T_ull` 是气枕能量方程的状态量。 代码位置：`src/cryo_tank/tank_model.py:1-7`，`src/cryo_tank/tank_model.py:87-89` @@ -41,40 +41,54 @@ T_ull : 气枕区温度 [K] 初始液相体积： -```text -V_liq_0 = (1 - ullage_fraction) * V_total -``` +$$ +V_{\mathrm{liq},0}=\left(1-f_{\mathrm{ull},0}\right)V_{\mathrm{total}} +$$ + +式中，`V_liq_0` 为初始液相体积，`ullage_fraction` 为初始气枕体积分数，`V_total` 为贮箱总容积。 初始气枕体积： -```text -V_ull_0 = ullage_fraction * V_total -``` +$$ +V_{\mathrm{ull},0}=f_{\mathrm{ull},0}V_{\mathrm{total}} +$$ + +式中，`V_ull_0` 为初始气枕体积。 初始液氮质量： -```text -rho_liq_0 = rho_LN2(T_init, P_work) -m_liq_0 = rho_liq_0 * V_liq_0 -``` +$$ +\begin{aligned} +\rho_{\mathrm{liq},0} &= \rho_{\mathrm{LN2}}\left(T_{\mathrm{init}},P_{\mathrm{work}}\right) \\ +m_{\mathrm{liq},0} &= \rho_{\mathrm{liq},0}V_{\mathrm{liq},0} +\end{aligned} +$$ + +式中，`rho_liq_0` 为初始液氮密度，`rho_LN2(T, P)` 为液氮在温度 `T`、压力 `P` 下的密度物性函数，`T_init` 为初始温度，`P_work` 为工作压力，`m_liq_0` 为初始液氮质量。 初始液相总内能： -```text -U_liq_0 = m_liq_0 * u_LN2(T_init, P_work) -``` +$$ +U_{\mathrm{liq},0}=m_{\mathrm{liq},0}u_{\mathrm{LN2}}\left(T_{\mathrm{init}},P_{\mathrm{work}}\right) +$$ + +式中，`U_liq_0` 为初始液相总内能，`u_LN2(T, P)` 为液氮比内能。 初始气枕温度： -```text -T_ull_0 = T_init -``` +$$ +T_{\mathrm{ull},0}=T_{\mathrm{init}} +$$ + +式中，`T_ull_0` 为初始气枕温度。 初始氦气质量是派生量： -```text -m_He_0 = rho_He(T_init, P_work) * V_ull_0 -``` +$$ +m_{\mathrm{He},0}=\rho_{\mathrm{He}}\left(T_{\mathrm{init}},P_{\mathrm{work}}\right)V_{\mathrm{ull},0} +$$ + +式中，`m_He_0` 为初始氦气质量，`rho_He(T, P)` 为氦气密度物性函数。 代码位置：`src/cryo_tank/tank_model.py:71-85` @@ -82,43 +96,57 @@ m_He_0 = rho_He(T_init, P_work) * V_ull_0 横截面积： -```text -A_cross = V_total / H_tank -``` +$$ +A_{\mathrm{cross}}=\frac{V_{\mathrm{total}}}{H_{\mathrm{tank}}} +$$ + +式中，`A_cross` 为圆柱贮箱横截面积，`H_tank` 为贮箱高度。 圆柱直径： -```text -D = sqrt(4 * A_cross / pi) -``` +$$ +D=\sqrt{\frac{4A_{\mathrm{cross}}}{\pi}} +$$ + +式中，`D` 为等效圆柱直径，`pi` 为圆周率。 侧面积： -```text -A_side = pi * D * H_tank -``` +$$ +A_{\mathrm{side}}=\pi D H_{\mathrm{tank}} +$$ + +式中，`A_side` 为圆柱侧壁面积。 端盖面积： -```text -A_cap = A_cross -``` +$$ +A_{\mathrm{cap}}=A_{\mathrm{cross}} +$$ + +式中，`A_cap` 为单个端盖面积。 总面积： -```text -A_total = A_side + 2 * A_cap -``` +$$ +A_{\mathrm{total}}=A_{\mathrm{side}}+2A_{\mathrm{cap}} +$$ + +式中，`A_total` 为贮箱外表面总面积，包含侧壁和两个端盖。 代码位置：`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) -``` +$$ +\begin{aligned} +L &= \operatorname{clamp}\left(H_{\mathrm{liq}},0,H_{\mathrm{tank}}\right) \\ +A_{\mathrm{wet}} &= A_{\mathrm{cap}}+\pi D L \\ +A_{\mathrm{dry}} &= A_{\mathrm{cap}}+\pi D\left(H_{\mathrm{tank}}-L\right) +\end{aligned} +$$ + +式中，`liquid_level` 为由液相体积换算得到的液位，`level` 为限制在 `[0, H_tank]` 范围内的液位，`A_wet` 为与液相接触的湿壁面积，`A_dry` 为与气枕接触的干壁面积。 代码位置：`src/cryo_tank/tank_model.py:91-96` @@ -126,32 +154,40 @@ A_dry = A_cap + pi * D * (H_tank - level) 给定状态： -```text -y = [m_liq, U_liq, T_ull] -``` +$$ +\mathbf{y}=\left[m_{\mathrm{liq}},\ U_{\mathrm{liq}},\ T_{\mathrm{ull}}\right] +$$ ### 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 -``` +$$ +\begin{aligned} +u_{\mathrm{liq}} &= \frac{U_{\mathrm{liq}}}{m_{\mathrm{liq}}} \\ +T_{\mathrm{liq}} &= T_{\mathrm{LN2}}\left(u_{\mathrm{liq}},P_{\mathrm{work}}\right) \\ +\rho_{\mathrm{liq}} &= \rho_{\mathrm{LN2}}\left(T_{\mathrm{liq}},P_{\mathrm{work}}\right) \\ +V_{\mathrm{liq}} &= \frac{m_{\mathrm{liq}}}{\rho_{\mathrm{liq}}} \\ +H_{\mathrm{liq}} &= \frac{V_{\mathrm{liq}}}{A_{\mathrm{cross}}} \\ +f_{\mathrm{fill}} &= \frac{H_{\mathrm{liq}}}{H_{\mathrm{tank}}} +\end{aligned} +$$ + +式中，`u_liq` 为液氮比内能，`T_liq` 为由比内能和压力反算得到的液相温度，`rho_liq` 为当前液氮密度，`V_liq` 为液相体积，`liquid_level` 为液位高度，`fill_fraction` 为液位高度占贮箱高度的比例。 代码位置：`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) -``` +$$ +\begin{aligned} +V_{\mathrm{ull}} &= V_{\mathrm{total}}-V_{\mathrm{liq}} \\ +P_{\mathrm{He}} &= P_{\mathrm{work}} \\ +\rho_{\mathrm{He}} &= \rho_{\mathrm{He}}\left(T_{\mathrm{ull}},P_{\mathrm{He}}\right) \\ +m_{\mathrm{He}} &= \rho_{\mathrm{He}}V_{\mathrm{ull}} \\ +U_{\mathrm{ull}} &= m_{\mathrm{He}}u_{\mathrm{He}}\left(T_{\mathrm{ull}},P_{\mathrm{He}}\right) +\end{aligned} +$$ + +式中，`V_ull` 为气枕体积，`P_He` 为氦气压力，`rho_He` 为当前氦气密度，`m_He` 为维持恒压所需的氦气质量，`U_ull` 为气枕氦气总内能，`u_He(T, P)` 为氦气比内能。 代码位置：`src/cryo_tank/tank_model.py:115-121` @@ -159,26 +195,32 @@ U_ull = m_He * u_He(T_ull, P_He) ### 6.1 液相与气枕界面换热 -```text -Q_liq_to_ull = h_conv * A_cross * (T_liq - T_ull) -``` +$$ +Q_{\mathrm{liq}\to\mathrm{ull}}=h_{\mathrm{conv}}A_{\mathrm{cross}}\left(T_{\mathrm{liq}}-T_{\mathrm{ull}}\right) +$$ + +式中，`Q_liq_to_ull` 为液相到气枕的界面换热率，`h_conv` 为界面对流换热系数，`A_cross` 为气液界面面积，`T_liq` 和 `T_ull` 分别为液相温度和气枕温度。 符号约定： -```text -Q_liq_to_ull > 0 : 热量从液相传给气枕 -Q_liq_to_ull < 0 : 热量从气枕传给液相 -``` +$$ +\begin{aligned} +Q_{\mathrm{liq}\to\mathrm{ull}} &> 0 && \text{热量从液相传给气枕} \\ +Q_{\mathrm{liq}\to\mathrm{ull}} &< 0 && \text{热量从气枕传给液相} +\end{aligned} +$$ -代码位置：`src/cryo_tank/tank_model.py:141-142` +代码位置：`src/cryo_tank/tank_model.py:142-143` ### 6.2 外界总漏热 当前主程序默认使用 MLI 漏热模型： -```text -Q_leak = A_total * q_mli -``` +$$ +Q_{\mathrm{leak}}=A_{\mathrm{total}}q_{\mathrm{MLI}} +$$ + +式中，`Q_leak` 为外界向贮箱的总漏热率，`q_mli` 为 MLI 单位面积漏热热流密度。 代码位置： @@ -187,10 +229,14 @@ Q_leak = A_total * q_mli 代码中还提供 Foam 漏热模型： -```text -T_mean = (T_inner + T_env) / 2 -Q_leak = A_total * k_eff(T_mean) * (T_env - T_inner) / delta -``` +$$ +\begin{aligned} +T_{\mathrm{mean}} &= \frac{T_{\mathrm{inner}}+T_{\mathrm{env}}}{2} \\ +Q_{\mathrm{leak}} &= A_{\mathrm{total}}k_{\mathrm{eff}}\left(T_{\mathrm{mean}}\right)\frac{T_{\mathrm{env}}-T_{\mathrm{inner}}}{\delta} +\end{aligned} +$$ + +式中，`T_mean` 为泡沫层平均温度，`T_inner` 为贮箱内侧温度，`T_env` 为环境温度，`k_eff` 为有效导热系数，`delta` 为泡沫保温层厚度。 如果 `k_eff` 是常数，则直接使用该常数。 @@ -200,92 +246,141 @@ Q_leak = A_total * k_eff(T_mean) * (T_env - T_inner) / delta 总漏热按湿壁面积和干壁面积分配： -```text -Q_leak_liq = Q_leak * A_wet / (A_wet + A_dry) -Q_leak_ull = Q_leak * A_dry / (A_wet + A_dry) -``` +$$ +\begin{aligned} +Q_{\mathrm{leak,liq}} &= Q_{\mathrm{leak}}\frac{A_{\mathrm{wet}}}{A_{\mathrm{wet}}+A_{\mathrm{dry}}} \\ +Q_{\mathrm{leak,ull}} &= Q_{\mathrm{leak}}\frac{A_{\mathrm{dry}}}{A_{\mathrm{wet}}+A_{\mathrm{dry}}} +\end{aligned} +$$ -代码位置：`src/cryo_tank/tank_model.py:144-148` +式中，`Q_leak_liq` 为分配到液相的漏热，`Q_leak_ull` 为分配到气枕的漏热。 + +代码位置：`src/cryo_tank/tank_model.py:145-149` ## 7. 液氮质量方程 液氮质量变化率： -```text -dm_liq/dt = mdot_in_ln2 - mdot_out_ln2 -``` +$$ +\frac{dm_{\mathrm{liq}}}{dt}=\dot{m}_{\mathrm{in,LN2}}-\dot{m}_{\mathrm{out,LN2}} +$$ + +式中，`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` +- ODE 中使用：`src/cryo_tank/tank_model.py:151-152` ## 8. 液相能量方程 入口液氮焓： -```text -h_in_ln2 = h_LN2(T_in_ln2, P_work) -``` +$$ +h_{\mathrm{in,LN2}}=h_{\mathrm{LN2}}\left(T_{\mathrm{in,LN2}},P_{\mathrm{work}}\right) +$$ + +式中，`h_in_ln2` 为入口液氮比焓，`h_LN2(T, P)` 为液氮比焓物性函数，`T_in_ln2` 为入口液氮温度。 当前液相出口焓： -```text -h_liq = h_LN2(T_liq, P_work) -``` +$$ +h_{\mathrm{liq}}=h_{\mathrm{LN2}}\left(T_{\mathrm{liq}},P_{\mathrm{work}}\right) +$$ + +式中，`h_liq` 为当前液相出口比焓，模型假定出口液氮状态等于贮箱液相主体状态。 液相总内能变化率： -```text -dU_liq/dt = - mdot_in_ln2 * h_in_ln2 - - mdot_out_ln2 * h_liq - - Q_liq_to_ull - + Q_leak_liq -``` +$$ +\frac{dU_{\mathrm{liq}}}{dt}=\dot{m}_{\mathrm{in,LN2}}h_{\mathrm{in,LN2}}-\dot{m}_{\mathrm{out,LN2}}h_{\mathrm{liq}}-Q_{\mathrm{liq}\to\mathrm{ull}}+Q_{\mathrm{leak,liq}} +$$ + +式中，`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` +- 能量方程调用：`src/cryo_tank/tank_model.py:153-155` +- 能量方程实现：`src/cryo_tank/tank_model.py:164-170` ## 9. 气枕温度和氦气质量导数方程 当前模型把 `T_ull` 作为 ODE 状态。氦气质量由恒压关系派生： -```text -m_He = rho_He(T_ull, P_work) * V_ull -``` +$$ +m_{\mathrm{He}}=\rho_{\mathrm{He}}\left(T_{\mathrm{ull}},P_{\mathrm{work}}\right)V_{\mathrm{ull}} +$$ + +式中，`m_He` 是派生氦气质量，不是 ODE 独立状态。 代码位置：`src/cryo_tank/tank_model.py:115-121` ### 9.1 气枕体积变化率 -液氮质量变化导致气枕体积变化： +液相体积由液氮质量和恒压液氮密度决定： -```text -dV_ull/dt = -dm_liq_dt / rho_liq -``` +$$ +V_{\mathrm{liq}}=\frac{m_{\mathrm{liq}}}{\rho_{\mathrm{liq}}\left(T_{\mathrm{liq}},P_{\mathrm{work}}\right)} +$$ -代码位置：`src/cryo_tank/tank_model.py:186-188` +气枕体积为： + +$$ +V_{\mathrm{ull}}=V_{\mathrm{total}}-V_{\mathrm{liq}} +$$ + +对液相体积求全导数： + +$$ +\frac{dV_{\mathrm{liq}}}{dt}=\frac{1}{\rho_{\mathrm{liq}}}\frac{dm_{\mathrm{liq}}}{dt}-\frac{m_{\mathrm{liq}}}{\rho_{\mathrm{liq}}^2} +\left(\frac{\partial\rho_{\mathrm{liq}}}{\partial T}\right)_P\frac{dT_{\mathrm{liq}}}{dt} +$$ + +其中液相温度变化率由 `u_liq = U_liq / m_liq` 得到： + +$$ +\frac{dT_{\mathrm{liq}}}{dt}=\frac{\frac{dU_{\mathrm{liq}}}{dt}-u_{\mathrm{liq}}\frac{dm_{\mathrm{liq}}}{dt}}{m_{\mathrm{liq}}\left(\frac{\partial u_{\mathrm{liq}}}{\partial T}\right)_P} +$$ + +因此气枕体积变化率为： + +$$ +\frac{dV_{\mathrm{ull}}}{dt}=-\frac{1}{\rho_{\mathrm{liq}}}\frac{dm_{\mathrm{liq}}}{dt}+\frac{m_{\mathrm{liq}}}{\rho_{\mathrm{liq}}^2} +\left(\frac{\partial\rho_{\mathrm{liq}}}{\partial T}\right)_P\frac{dT_{\mathrm{liq}}}{dt} +$$ + +式中，`dV_ull/dt` 为气枕体积变化率，第一项为液氮质量变化导致的体积变化，第二项为液氮温度变化导致密度变化后产生的体积修正项。液氮的恒压密度温度导数和恒压内能温度导数由 CoolProp 提供。 + +代码位置： + +- LN2 恒压偏导：`src/cryo_tank/properties.py:47-56` +- 液相温度变化率：`src/cryo_tank/tank_model.py:172-180` +- 气枕体积变化率：`src/cryo_tank/tank_model.py:182-194` ### 9.2 气枕区非氦气入口热量 -```text -Q_ull_no_he = Q_liq_to_ull + Q_leak_ull -``` +$$ +Q_{\mathrm{ull,noHe}}=Q_{\mathrm{liq}\to\mathrm{ull}}+Q_{\mathrm{leak,ull}} +$$ -代码位置：`src/cryo_tank/tank_model.py:188` +式中，`Q_ull_no_he` 为不含氦气入口焓流的气枕净受热率，包括液相传入气枕的界面换热和分配到气枕干壁的外界漏热。 + +代码位置：`src/cryo_tank/tank_model.py:221` ### 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) -``` +$$ +\begin{aligned} +\rho &= \rho_{\mathrm{He}}\left(T_{\mathrm{ull}},P_{\mathrm{work}}\right) \\ +m(T,V) &= \rho V \\ +U(T,V) &= m(T,V)u_{\mathrm{He}}\left(T_{\mathrm{ull}},P_{\mathrm{work}}\right) +\end{aligned} +$$ + +式中，`rho` 为氦气密度，`T` 表示 `T_ull`，`V` 表示 `V_ull`，`m(T, V)` 为恒压下由温度和体积确定的氦气质量，`U(T, V)` 为恒压下由温度和体积确定的气枕总内能。 代码位置：`src/cryo_tank/tank_model.py:115-121` @@ -293,120 +388,226 @@ U(T, V) = m(T, V) * u_He(T_ull, P_work) 代码不再对温度和体积做有限差分，而是使用恒压物性导数： -```text -drhodT_P = (partial rho / partial T)_P -dudT_P = (partial u / partial T)_P -``` +$$ +\begin{aligned} +\left(\frac{\partial\rho}{\partial T}\right)_P &= \rho_T\big|_P \\ +\left(\frac{\partial u}{\partial T}\right)_P &= u_T\big|_P +\end{aligned} +$$ -CoolProp 偏导接口位置：`src/cryo_tank/properties.py:100-109` +式中，`drhodT_P` 为恒压氦气密度对温度的偏导数，`dudT_P` 为恒压氦气比内能对温度的偏导数。 + +CoolProp 偏导接口位置：`src/cryo_tank/properties.py:113-122` 由此得到： -```text -m_T = V_ull * drhodT_P -m_V = rho -U_T = V_ull * (u * drhodT_P + rho * dudT_P) -U_V = rho * u -``` +$$ +\begin{aligned} +m_T &= V_{\mathrm{ull}}\left(\frac{\partial\rho}{\partial T}\right)_P \\ +m_V &= \rho \\ +U_T &= V_{\mathrm{ull}}\left[u\left(\frac{\partial\rho}{\partial T}\right)_P+\rho\left(\frac{\partial u}{\partial T}\right)_P\right] \\ +U_V &= \rho u +\end{aligned} +$$ -代码位置：`src/cryo_tank/tank_model.py:165-178` +式中，`m_T = (partial m / partial T)_V`，`m_V = (partial m / partial V)_T`，`U_T = (partial U / partial T)_V`，`U_V = (partial U / partial V)_T`。这些偏导用于把 `dm_He/dt` 和 `dU_ull/dt` 展开为 `dT_ull/dt` 与 `dV_ull/dt` 的线性组合。 + +推导过程如下。恒压下 `rho = rho(T, P_work)`，`u = u(T, P_work)`，所以： + +$$ +m(T,V)=\rho\left(T,P_{\mathrm{work}}\right)V +$$ + +对温度和体积求偏导： + +$$ +\begin{aligned} +m_T &= \left(\frac{\partial m}{\partial T}\right)_V + = V\left(\frac{\partial\rho}{\partial T}\right)_P + = V_{\mathrm{ull}}\rho_T\big|_P \\ +m_V &= \left(\frac{\partial m}{\partial V}\right)_T + = \rho +\end{aligned} +$$ + +气枕总内能为： + +$$ +U(T,V)=m(T,V)u\left(T,P_{\mathrm{work}}\right)=\rho\left(T,P_{\mathrm{work}}\right)Vu\left(T,P_{\mathrm{work}}\right) +$$ + +对温度和体积求偏导： + +$$ +\begin{aligned} +U_T &= \left(\frac{\partial U}{\partial T}\right)_V + = V\left[u\left(\frac{\partial\rho}{\partial T}\right)_P+\rho\left(\frac{\partial u}{\partial T}\right)_P\right] \\ + &= V_{\mathrm{ull}}\left(u\rho_T\big|_P+\rho u_T\big|_P\right) \\ +U_V &= \left(\frac{\partial U}{\partial V}\right)_T + = \rho u +\end{aligned} +$$ + +代码位置：`src/cryo_tank/tank_model.py:196-209` ### 9.5 气枕温度变化率 氦气入口焓： -```text -h_in_he = h_He(T_in_he, P_work) -``` +$$ +h_{\mathrm{in,He}}=h_{\mathrm{He}}\left(T_{\mathrm{in,He}},P_{\mathrm{work}}\right) +$$ + +式中，`h_in_he` 为入口氦气比焓，`h_He(T, P)` 为氦气比焓物性函数，`T_in_he` 为入口氦气温度。 + +氦气气枕能量方程采用开口、变体积控制体的一阶热力学形式： + +$$ +\frac{dU_{\mathrm{ull}}}{dt}=\dot{m}_{\mathrm{He}}h_{\mathrm{in,He}}+Q_{\mathrm{ull,noHe}}-P_{\mathrm{work}}\frac{dV_{\mathrm{ull}}}{dt} +$$ + +式中，`dU_ull/dt` 为气枕氦气总内能变化率，`mdot_He * h_in_he` 为补入氦气带入的焓流，`Q_ull_no_he` 为界面换热与外界漏热给气枕的热输入，`P_work * dV_ull/dt` 为气枕边界膨胀功项。当前模型无气枕出口流量，因此没有出口焓流项。 + +由恒压状态函数的全微分： + +$$ +\begin{aligned} +\frac{dm_{\mathrm{He}}}{dt} &= m_T\frac{dT_{\mathrm{ull}}}{dt}+m_V\frac{dV_{\mathrm{ull}}}{dt} \\ +\frac{dU_{\mathrm{ull}}}{dt} &= U_T\frac{dT_{\mathrm{ull}}}{dt}+U_V\frac{dV_{\mathrm{ull}}}{dt} +\end{aligned} +$$ + +把质量导数代入气枕能量方程： + +$$ +U_T\frac{dT_{\mathrm{ull}}}{dt}+U_V\frac{dV_{\mathrm{ull}}}{dt}=h_{\mathrm{in,He}}\left(m_T\frac{dT_{\mathrm{ull}}}{dt}+m_V\frac{dV_{\mathrm{ull}}}{dt}\right)+Q_{\mathrm{ull,noHe}}-P_{\mathrm{work}}\frac{dV_{\mathrm{ull}}}{dt} +$$ + +整理 `dT_ull/dt` 项： + +$$ +\left(U_T-h_{\mathrm{in,He}}m_T\right)\frac{dT_{\mathrm{ull}}}{dt}=Q_{\mathrm{ull,noHe}}-\left(U_V+P_{\mathrm{work}}-h_{\mathrm{in,He}}m_V\right)\frac{dV_{\mathrm{ull}}}{dt} +$$ 气枕温度变化率： -```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] -``` +$$ +\frac{dT_{\mathrm{ull}}}{dt}=\frac{Q_{\mathrm{ull,noHe}}-\left(U_V+P_{\mathrm{work}}-h_{\mathrm{in,He}}m_V\right)\frac{dV_{\mathrm{ull}}}{dt}}{U_T-h_{\mathrm{in,He}}m_T} +$$ -代码位置：`src/cryo_tank/tank_model.py:190-197` +式中，分子表示扣除体积变化、边界功和入口质量变化耦合项后的气枕有效热输入，分母表示在恒压约束和入口补气耦合下气枕对温度变化的等效热容项。 + +代码位置：`src/cryo_tank/tank_model.py:223-230` 若分母过小： -```text -if abs(U_T - h_in_he * m_T) < 1e-30: - dT_ull/dt = 0 - mdot_He = 0 -``` +$$ +\left|U_T-h_{\mathrm{in,He}}m_T\right|<10^{-30} +\quad\Longrightarrow\quad +\begin{cases} +\dfrac{dT_{\mathrm{ull}}}{dt}=0 \\ +\dot{m}_{\mathrm{He}}=0 +\end{cases} +$$ -代码位置：`src/cryo_tank/tank_model.py:190-193` +式中，`1e-30` 是数值保护阈值，用于避免除以过小分母。 + +代码位置：`src/cryo_tank/tank_model.py:223-226` ### 9.6 氦气质量变化率 -```text -dm_He/dt = mdot_He -``` +$$ +\frac{dm_{\mathrm{He}}}{dt}=\dot{m}_{\mathrm{He}} +$$ + +式中，`dm_He/dt` 为恒压约束下氦气质量变化率，`mdot_He` 为需要补入的氦气质量流量。 其中： -```text -mdot_He = m_T * dT_ull/dt + m_V * dV_ull/dt -``` +$$ +\dot{m}_{\mathrm{He}}=m_T\frac{dT_{\mathrm{ull}}}{dt}+m_V\frac{dV_{\mathrm{ull}}}{dt} +$$ -代码位置：`src/cryo_tank/tank_model.py:198` +式中，第一项为气枕温度变化导致的氦气质量变化，第二项为气枕体积变化导致的氦气质量变化。 + +代码位置：`src/cryo_tank/tank_model.py:231` 若 `mdot_He < 0`，报告流量钳制为 0： -```text -mdot_He = 0 -``` +$$ +\dot{m}_{\mathrm{He}}=0 +$$ -代码位置：`src/cryo_tank/tank_model.py:200-205` +式中，钳制只作用于报告的补气流量；当前 ODE 状态 `T_ull` 的导数仍由未钳制前的恒压能量方程求得。 + +代码位置：`src/cryo_tank/tank_model.py:233-238` ## 10. 完整 ODE 方程组 状态向量： -```text -y = [m_liq, U_liq, T_ull] -``` +$$ +\mathbf{y}=\left[m_{\mathrm{liq}},\ U_{\mathrm{liq}},\ T_{\mathrm{ull}}\right] +$$ + +式中，`y` 为 ODE 状态向量，三个分量分别为液氮质量、液相总内能和气枕温度。 ODE 方程组： -```text -dm_liq/dt = mdot_in_ln2 - mdot_out_ln2 -``` +$$ +\frac{dm_{\mathrm{liq}}}{dt}=\dot{m}_{\mathrm{in,LN2}}-\dot{m}_{\mathrm{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] -``` +$$ +\frac{dU_{\mathrm{liq}}}{dt}=\dot{m}_{\mathrm{in,LN2}}h_{\mathrm{in,LN2}}-\dot{m}_{\mathrm{out,LN2}}h_{\mathrm{liq}}-Q_{\mathrm{liq}\to\mathrm{ull}}+Q_{\mathrm{leak,liq}} +$$ + +式中，液相总内能变化率由入口焓流、出口焓流、液相向气枕换热和液相漏热共同决定。 + +$$ +\frac{dT_{\mathrm{ull}}}{dt}=\frac{Q_{\mathrm{ull,noHe}}-\left(U_V+P_{\mathrm{work}}-h_{\mathrm{in,He}}m_V\right)\frac{dV_{\mathrm{ull}}}{dt}}{U_T-h_{\mathrm{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) -``` +$$ +\begin{aligned} +h_{\mathrm{liq}} &= h_{\mathrm{LN2}}\left(T_{\mathrm{liq}},P_{\mathrm{work}}\right) \\ +Q_{\mathrm{liq}\to\mathrm{ull}} &= h_{\mathrm{conv}}A_{\mathrm{cross}}\left(T_{\mathrm{liq}}-T_{\mathrm{ull}}\right) \\ +Q_{\mathrm{ull,noHe}} &= Q_{\mathrm{liq}\to\mathrm{ull}}+Q_{\mathrm{leak,ull}} \\ +Q_{\mathrm{leak,liq}} &= Q_{\mathrm{leak}}\frac{A_{\mathrm{wet}}}{A_{\mathrm{wet}}+A_{\mathrm{dry}}} \\ +Q_{\mathrm{leak,ull}} &= Q_{\mathrm{leak}}\frac{A_{\mathrm{dry}}}{A_{\mathrm{wet}}+A_{\mathrm{dry}}} +\end{aligned} +$$ + +气枕体积变化率使用完整链式法则： + +$$ +\begin{aligned} +\frac{dT_{\mathrm{liq}}}{dt} &= \frac{\frac{dU_{\mathrm{liq}}}{dt}-u_{\mathrm{liq}}\frac{dm_{\mathrm{liq}}}{dt}}{m_{\mathrm{liq}}\left(\frac{\partial u_{\mathrm{liq}}}{\partial T}\right)_P} \\ +\frac{dV_{\mathrm{ull}}}{dt} &= -\frac{1}{\rho_{\mathrm{liq}}}\frac{dm_{\mathrm{liq}}}{dt}+\frac{m_{\mathrm{liq}}}{\rho_{\mathrm{liq}}^2} +\left(\frac{\partial\rho_{\mathrm{liq}}}{\partial T}\right)_P\frac{dT_{\mathrm{liq}}}{dt} +\end{aligned} +$$ + +式中，`h_liq` 为当前液相比焓，`Q_liq_to_ull` 为液相到气枕界面换热，`Q_ull_no_he` 为不含入口氦气焓流的气枕热输入，`Q_leak_liq` 和 `Q_leak_ull` 分别为分配到液相和气枕的漏热；`dV_ull/dt` 同时包含液氮质量变化和液氮密度随温度变化的贡献。 氦气质量和质量流量是派生量： -```text -m_He = rho_He(T_ull, P_work) * V_ull -mdot_He = m_T * dT_ull/dt + m_V * dV_ull/dt -``` +$$ +\begin{aligned} +m_{\mathrm{He}} &= \rho_{\mathrm{He}}\left(T_{\mathrm{ull}},P_{\mathrm{work}}\right)V_{\mathrm{ull}} \\ +\dot{m}_{\mathrm{He}} &= m_T\frac{dT_{\mathrm{ull}}}{dt}+m_V\frac{dV_{\mathrm{ull}}}{dt} +\end{aligned} +$$ -代码位置：`src/cryo_tank/tank_model.py:131-214` +式中，`m_He` 为恒压状态方程给出的当前氦气质量，`mdot_He` 为由质量全微分得到的所需补气流量。 + +代码位置：`src/cryo_tank/tank_model.py:132-257` ## 11. ODE 求解逻辑 @@ -416,13 +617,17 @@ mdot_He = m_T * dT_ull/dt + m_V * dV_ull/dt run(tank, t_end, rtol=1e-8, atol=1e-10, max_step=10.0) ``` +式中，`tank` 为 `CryoTank` 模型实例，`t_end` 为仿真终止时间，`rtol` 和 `atol` 为相对、绝对误差容限，`max_step` 为求解器最大时间步长。 + 代码位置：`src/cryo_tank/solver.py:22-24` 求解初值： -```text -y0 = tank.initial_state() -``` +$$ +\mathbf{y}_0=\operatorname{initial\_state}\left(\mathrm{tank}\right) +$$ + +式中，`y0` 为初始 ODE 状态向量。 代码位置：`src/cryo_tank/solver.py:43` @@ -442,20 +647,28 @@ solve_ivp( ) ``` +式中，`tank.rhs` 为 ODE 右端函数，`[0.0, t_end]` 为积分时间区间，`method="RK45"` 指四、五阶 Runge-Kutta 自适应算法，`events` 用于注册液体排空终止事件，`dense_output=True` 表示生成连续插值解。 + 代码位置：`src/cryo_tank/solver.py:45-55` 液体排空事件： -```text -event(t, y) = m_liq -``` +$$ +g(t,\mathbf{y})=m_{\mathrm{liq}} +$$ + +式中，`event(t, y)` 为事件函数；当状态中的 `m_liq` 到达零时，事件函数到达零。 事件属性： -```text -terminal = True -direction = -1 -``` +$$ +\begin{aligned} +\mathrm{terminal} &= \mathrm{True} \\ +\mathrm{direction} &= -1 +\end{aligned} +$$ + +式中，`terminal=True` 表示事件触发后终止积分，`direction=-1` 表示只检测事件函数从正值下降到零的穿越。 含义： @@ -472,19 +685,25 @@ direction = -1 求解完成后，`solver.run()` 对每个输出时刻执行： -```text -info = tank.derive(y_i) -``` +$$ +\mathrm{info}_i=\operatorname{derive}\left(\mathbf{y}_i\right) +$$ + +式中，`y_i` 为某一输出时刻的状态向量，`info` 为根据该状态计算出的温度、体积、质量和换热等派生量字典。 并重新计算： -```text -Q_liq_to_ull -Q_leak -Q_leak_liq -Q_leak_ull -mdot_He -``` +$$ +\left\{ +Q_{\mathrm{liq}\to\mathrm{ull}},\ +Q_{\mathrm{leak}},\ +Q_{\mathrm{leak,liq}},\ +Q_{\mathrm{leak,ull}},\ +\dot{m}_{\mathrm{He}} +\right\} +$$ + +式中，这些量分别为界面换热、总漏热、液相漏热、气枕漏热和派生氦气补气质量流量，用于输出和结果分析。 其中 `m_He`、`U_ull`、`mdot_He` 均为后处理派生结果，不是独立 ODE 状态。 diff --git a/src/cryo_tank/Index.pdf b/src/cryo_tank/Index.pdf deleted file mode 100644 index 6f28992..0000000 Binary files a/src/cryo_tank/Index.pdf and /dev/null differ diff --git a/src/cryo_tank/config.py b/src/cryo_tank/config.py index b634228..d07a669 100644 --- a/src/cryo_tank/config.py +++ b/src/cryo_tank/config.py @@ -34,7 +34,7 @@ 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 = 4818.884042 # s (liquid level reaches 5%) RTOL = 1e-8 ATOL = 1e-10 diff --git a/src/cryo_tank/properties.py b/src/cryo_tank/properties.py index 3fc02ab..18688e3 100644 --- a/src/cryo_tank/properties.py +++ b/src/cryo_tank/properties.py @@ -44,6 +44,18 @@ def ln2_T_from_u(u, P): return _n2_state.T() +def ln2_drho_dT_const_p(T, P): + """LN2 density temperature derivative at constant pressure [kg/(m^3*K)].""" + _n2_state.update(CP.PT_INPUTS, P, T) + return _n2_state.first_partial_deriv(CP.iDmass, CP.iT, CP.iP) + + +def ln2_du_dT_const_p(T, P): + """LN2 internal-energy temperature derivative at constant pressure [J/(kg*K)].""" + _n2_state.update(CP.PT_INPUTS, P, T) + return _n2_state.first_partial_deriv(CP.iUmass, CP.iT, CP.iP) + + # --------------------------------------------------------------------------- # Helium properties via CoolProp # --------------------------------------------------------------------------- diff --git a/src/cryo_tank/solver.py b/src/cryo_tank/solver.py index d69ef57..6093268 100644 --- a/src/cryo_tank/solver.py +++ b/src/cryo_tank/solver.py @@ -120,7 +120,7 @@ def run(tank, t_end, rtol=1e-8, atol=1e-10, max_step=10.0): 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 + info, Q_liq_to_ull, Q_leak_ull, Q_leak_liq ) diff --git a/src/cryo_tank/tank_model.py b/src/cryo_tank/tank_model.py index 352407a..8b44ee8 100644 --- a/src/cryo_tank/tank_model.py +++ b/src/cryo_tank/tank_model.py @@ -121,6 +121,7 @@ class CryoTank: U_ull = m_He * prop.he_u(T_ull, P_He) return { + 'm_liq': m_liq, 'U_liq': U_liq, 'u_liq': u_liq, 'T_liq': T_liq, 'T_ull': T_ull, 'V_liq': V_liq, 'V_ull': V_ull, 'liquid_level': liquid_level, 'fill_fraction': fill_fraction, @@ -149,19 +150,49 @@ class CryoTank: # --- Liquid zone --- dm_liq_dt = self.dm_liq_dt # constant: mdot_in - mdot_out - h_liq = prop.ln2_h(T_liq, self.P_work) - dU_liq_dt = (self.mdot_in_ln2 * self.h_in_ln2 - - self.mdot_out_ln2 * h_liq - - Q_liq_to_ull - + Q_leak_liq) + dU_liq_dt = self._liquid_energy_rate( + info, Q_liq_to_ull, Q_leak_liq + ) # --- Ullage zone --- dT_ull_dt, _ = self._solve_ullage_temperature_rate( - info, Q_liq_to_ull, Q_leak_ull + info, Q_liq_to_ull, Q_leak_ull, dU_liq_dt ) return np.array([dm_liq_dt, dU_liq_dt, dT_ull_dt]) + def _liquid_energy_rate(self, info, Q_liq_to_ull, Q_leak_liq): + """Return liquid-zone total internal-energy rate [W].""" + h_liq = prop.ln2_h(info['T_liq'], self.P_work) + return (self.mdot_in_ln2 * self.h_in_ln2 + - self.mdot_out_ln2 * h_liq + - Q_liq_to_ull + + Q_leak_liq) + + def _liquid_temperature_rate(self, info, dm_liq_dt, dU_liq_dt): + """Return dT_liq/dt from u_liq=U_liq/m_liq at constant tank pressure.""" + m_liq = info['m_liq'] + u_liq = info['u_liq'] + du_dT = prop.ln2_du_dT_const_p(info['T_liq'], self.P_work) + denominator = m_liq * du_dT + if abs(denominator) < 1e-30: + return 0.0 + return (dU_liq_dt - u_liq * dm_liq_dt) / denominator + + def _ullage_volume_rate(self, info, dm_liq_dt, dU_liq_dt): + """Return dV_ull/dt including liquid density variation with temperature.""" + rho_liq = info['rho_liq'] + m_liq = info['m_liq'] + T_liq = info['T_liq'] + dT_liq_dt = self._liquid_temperature_rate( + info, dm_liq_dt, dU_liq_dt + ) + drho_dT = prop.ln2_drho_dT_const_p(T_liq, self.P_work) + + dV_liq_dt = (dm_liq_dt / rho_liq + - m_liq * drho_dT * dT_liq_dt / rho_liq ** 2) + return -dV_liq_dt + 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'] @@ -177,14 +208,16 @@ class CryoTank: 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): + def _solve_ullage_temperature_rate(self, info, Q_liq_to_ull, Q_leak_ull, + dU_liq_dt): """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 + dV_ull_dt = self._ullage_volume_rate( + info, self.dm_liq_dt, dU_liq_dt + ) Q_ull_no_he = Q_liq_to_ull + Q_leak_ull m_T, U_T, m_V, U_V = self._he_mass_energy_partials(info) @@ -206,9 +239,19 @@ class CryoTank: return dT_ull_dt, mdot_He - def _solve_he_flow_rate(self, info, Q_liq_to_ull, Q_leak_ull): + def _solve_he_flow_rate(self, info, Q_liq_to_ull, Q_leak_ull, + Q_leak_liq=None): """Return derived He inlet mass flow for the current state [kg/s].""" + if Q_leak_liq is None: + Q_leak = self.heat_leak_model.compute(info['T_liq'], self.T_env) + A_wet, A_dry = self.wetted_areas(info['liquid_level']) + A_total = A_wet + A_dry + Q_leak_liq = Q_leak * A_wet / A_total if A_total > 0 else 0.0 + + dU_liq_dt = self._liquid_energy_rate( + info, Q_liq_to_ull, Q_leak_liq + ) _, mdot_He = self._solve_ullage_temperature_rate( - info, Q_liq_to_ull, Q_leak_ull + info, Q_liq_to_ull, Q_leak_ull, dU_liq_dt ) return mdot_He diff --git a/tests/cryo_tank/test_properties.py b/tests/cryo_tank/test_properties.py index dadf977..371340a 100644 --- a/tests/cryo_tank/test_properties.py +++ b/tests/cryo_tank/test_properties.py @@ -6,6 +6,7 @@ sys.path.insert(0, "src") from cryo_tank.properties import ( ln2_rho, ln2_h, ln2_u, ln2_T_from_u, + ln2_drho_dT_const_p, ln2_du_dT_const_p, he_u, he_h, he_cp, he_cv, ) from cryo_tank.config import P_WORKING @@ -32,6 +33,21 @@ class TestLN2Properties: T_recovered = ln2_T_from_u(u, P_WORKING) assert abs(T_recovered - T_orig) < 0.01 + def test_ln2_density_derivative_matches_finite_difference(self): + T = 78.0 + step = 1e-3 + actual = ln2_drho_dT_const_p(T, P_WORKING) + expected = ( + ln2_rho(T + step, P_WORKING) + - ln2_rho(T - step, P_WORKING) + ) / (2.0 * step) + assert actual < 0.0 + assert abs(actual - expected) / abs(expected) < 1e-5 + + def test_ln2_internal_energy_derivative_is_positive(self): + du_dT = ln2_du_dT_const_p(78.0, P_WORKING) + assert du_dT > 0.0 + class TestHeliumProperties: """Helium properties.""" diff --git a/tests/cryo_tank/test_tank_model.py b/tests/cryo_tank/test_tank_model.py index fef21e9..0fa354d 100644 --- a/tests/cryo_tank/test_tank_model.py +++ b/tests/cryo_tank/test_tank_model.py @@ -12,10 +12,9 @@ from cryo_tank.config import ( ) -def _make_tank(): +def _make_tank(**overrides): """Create a CryoTank with default config and MLI heat leak.""" - heat_leak = MLIHeatLeak(A_total=A_TOTAL, q_mli=1.0) - return CryoTank( + kw = dict( V_total=V_TOTAL, H_tank=H_TANK, P_work=P_WORKING, T_init=T_INIT, ullage_fraction=ULLAGE_FRACTION, @@ -23,8 +22,10 @@ def _make_tank(): mdot_out_ln2=MDOT_OUT_LN2, T_in_he=T_IN_HE, h_conv=H_CONV_SURFACE, T_env=T_ENV, - heat_leak_model=heat_leak, + heat_leak_model=MLIHeatLeak(A_total=A_TOTAL, q_mli=1.0), ) + kw.update(overrides) + return CryoTank(**kw) class TestGeometry: @@ -75,3 +76,37 @@ class TestInitialState: info = tank.derive(y0) P_He = info['P_He'] assert abs(P_He - P_WORKING) / P_WORKING < 1e-12 + + +class TestVolumeRate: + + def test_ullage_volume_rate_includes_liquid_thermal_expansion(self): + tank = _make_tank( + mdot_in_ln2=0.0, + mdot_out_ln2=0.0, + h_conv=0.0, + ) + y0 = tank.initial_state() + info = tank.derive(y0) + + Q_liq_to_ull = 0.0 + Q_leak = tank.heat_leak_model.compute(info['T_liq'], tank.T_env) + A_wet, A_dry = tank.wetted_areas(info['liquid_level']) + Q_leak_liq = Q_leak * A_wet / (A_wet + A_dry) + dU_liq_dt = tank._liquid_energy_rate( + info, Q_liq_to_ull, Q_leak_liq + ) + + dV_ull_dt = tank._ullage_volume_rate( + info, tank.dm_liq_dt, dU_liq_dt + ) + + y_next = y0.copy() + dt = 1.0 + y_next[1] += dU_liq_dt * dt + finite_difference = ( + tank.derive(y_next)['V_ull'] - info['V_ull'] + ) / dt + + assert dV_ull_dt < 0.0 + assert abs(dV_ull_dt - finite_difference) / abs(finite_difference) < 1e-5