读博那几年，我的研究方向是质子交换膜燃料电池冷启动的跨尺度建模，从催化层微结构一直做到单电池的温度场。毕业进了储能行业，这套手艺日常用在磷酸铁锂电池上，燃料电池的线偶尔也会想起，读博时的一些想法也没机会想完。现在纯粹出于兴趣，想把这些补上，做法是挑一篇当年给过我较多启发的文献，用 Python 从零把它的模型实现一遍，当年怎么算的，这回一行一行重新走。

就是 Huo、Jiao 和 Park 2019 年发表在 Applied Energy 上的那篇冷启动模型论文（原文）。读博搭自己模型的时候，这篇文章在机理组织和参数选择上给过我不少启发，这次换了身份重读，感受完全不同。当读者和当复现者要的东西不一样，读者要读懂主线，复现者要让每个符号都能落进代码，量纲、单位、参考态、分支条件，一处含糊后面都是坑。

先把总体判断放在这里。这篇文章的模型结构足以支撑一次严肃的复现，五种水状态、六类相变通道、电与热的双向耦合都摆在明面上；参数表和符号体系却不能直接抄进代码。我逐页核对了全部 36 个编号公式和 5 张参数表，量纲冲突、单位混用、数值排版问题都能找到，还有些闭合关系原文没有给。这篇札记就是这份对账单，公式编号与原文一致，方便对照。

按计划写六篇。本篇拆文献的水热机理和公式；第二篇先立符号与单位的规矩，再往后是离散骨架、物性与加载的坑、脚本框架与对照表；最后一篇记录调参和结果对标。

  1. 燃料电池冷启动复现札记（一）· 解构 Huo 2019 的水热机理与公式
  2. 燃料电池冷启动复现札记（二）· 符号转换表与 SI 摩尔基准

一、一块电池和五种水

质子交换膜燃料电池低温启动的核心矛盾很清楚。零下环境里电池发生氧化还原反应产水，水在离聚物、孔隙和冰之间反复分配，冰占住孔隙，氧气就到不了反应位点，电压跟着垮。启动成败取决于产水结冰和发热融冰谁跑得快。Huo 这篇的价值在于把过冷液水也纳入模型，用五种水状态把这条因果链完整摆开。

模型沿膜电极组件厚度方向建一维，完整层序是

$$ \mathrm{BP_a\,|\,Channel_a\,|\,GDL_a\,|\,MPL_a\,|\,CL_a\,|\,Membrane\,|\,CL_c\,|\,MPL_c\,|\,GDL_c\,|\,Channel_c\,|\,BP_c} $$

下标 a、c 表示阳极和阴极，BP 是双极板，GDL 是气体扩散层，MPL 是微孔层，CL 是催化层。层序明确不等于边界闭合，原文没有给出每层网格数，也没有讲清流道在一维传热传质里如何等效，复现时要自己决定流道是规定值边界还是独立控制体。

水状态符号主要位置迁移方式
未冻膜水$\lambda_{\rm nfm}$膜与两侧 CL 的离聚物扩散与电渗拖曳
冻结膜水$\lambda_{\rm fm}$膜内就地冻结与融化
水蒸气$C_{\rm vap}$CL、MPL、GDL 孔隙气相扩散
过冷液水$s_{\rm lq}$CL、MPL、GDL 孔隙压力驱动流动
孔隙冰$s_{\rm ice}$CL、MPL、GDL 孔隙就地累积

膜水这个名字容易误导。催化层里的离聚物同样装着膜水，名字说的是介质，不能据此只在膜层设变量。冻结膜水和孔隙冰是两个库存，前者影响还能导质子的未冻水量，后者直接占用气体通道，混在一起会把两条物理反馈搅乱。

我的做法是把五种状态统一折算成每单位总体积的摩尔库存

$$ q_{\rm nfm}=\omega\frac{\rho_m}{EW}\lambda_{\rm nfm},\quad q_{\rm lq}=\frac{\varepsilon\rho_{\rm lq}s_{\rm lq}}{M_w},\quad q_{\rm ice}=\frac{\varepsilon\rho_{\rm ice}s_{\rm ice}}{M_w},\quad q_{\rm vap}=\varepsilon g C_{\rm vap},\quad g=1-s_{\rm lq}-s_{\rm ice}. $$

冻结膜水的记账与未冻膜水同形。这一步原文没有明说，是我为源项配平推出来的，但它决定了程序里源项会不会差一个孔隙率或摩尔质量因子。$\rho_m/EW$ 是每单位离聚物体积的酸性位点摩尔浓度，$\lambda$ 是每位点对应的水分子数，无量纲；$s$ 是相对孔隙体积的饱和度，$\varepsilon s$ 才是相对总体积的体积分数；$\omega$ 是离聚物占总体积的比例。顺带一个单位坑，表1的 $EW=1100\ \mathrm{kg/kmol}$ 在 SI 摩尔基准下是 $1.1\ \mathrm{kg/mol}$，直接把 1100 代进配着 $F=96487\ \mathrm{C/mol}$ 的体系会差一千倍，$\rho_m=1980\ \mathrm{kg/m^3}$ 时位点浓度是 $1800\ \mathrm{mol/m^3}$。

原文还有几处符号复用要当心。$r_i$ 是式(5)的浓度反应级数，式(28)的 $r$ 却是孔半径；$\sigma$ 一个下标是表面张力，另一个是电导率；$I$ 是面电流密度，$J_i$ 是体积反应电流密度。写代码时这些必须拆成不同的名字。

五条基本假设里，对复现影响最大的是两条。流道液水和水蒸气设为零，等价于极强的水移除边界，但它只约束流道，不等于 GDL 里也没有水。论文忽略体相对流，却给出了化学计量比，一维实现里若只规定流道浓度，化学计量比不会自动产生流量效应，这个缺口要自己补闭合说明。

二、电压从哪里来，损失在哪里

输出电压写成三项之和

$$ V_{\rm out}=V_{\rm nernst}+\eta_{\rm act}+V_{\rm ohm}, \tag{1} $$

其中 Nernst 电压为

$$ V_{\rm nernst}=1.23-0.9\times10^{-3}(T-T_0) +\frac{RT}{2F}\ln\left(\frac{p_{\rm H_2}p_{\rm O_2}^{0.5}}{p_{\rm H_2O}}\right). \tag{2} $$

这一处我逐字核对时撞上问题。对数必须作用在无量纲量上，原文把分压标成 Pa 直接放进比值，麻烦跟着来两条。其一，压力不归一化时数值依赖单位选择，在 253.15 K 下能引入约 0.063 V 的偏差，对冷启动这种电压动态本就不大的过程不算小数目。其二，初始流道是干的，$p_{\rm H_2O}$ 取零会让对数发散。我准备按活度形式实现，理想气体取 $a_i=p_i/p^0$，水的活度参考态单独声明，这属于登记在册的实现解释。

电子与质子两相的电荷守恒是

$$ \nabla\cdot(-\sigma_e^{\rm eff}\nabla V_e)=J_i, \tag{3} $$$$ \nabla\cdot(-\sigma_H^{\rm eff}\nabla V_H)=-J_i. \tag{4} $$

两式相加说明电荷只在两相之间转换。反应动力学用 Butler-Volmer 形式，叠上液水与冰对位点的阻塞

$$ J_i=(1-s_{\rm lq}-s_{\rm ice})J_{\rm ref}\left(\frac{C_i}{C_{i,\rm ref}}\right)^{r_i} \left[\exp\left(\frac{\alpha nF\eta_{\rm act}}{RT}\right) -\exp\left(-\frac{\alpha nF\eta_{\rm act}}{RT}\right)\right]. \tag{6} $$

两条堵塞通道同时存在，式(6)按 $g=1-s_{\rm lq}-s_{\rm ice}$ 线性压缩反应能力，表3的 $D^{\rm eff}\propto g^{1.5}$ 又压缩气体扩散，一个管位点一个管补给，可以并存，但不能再额外乘来源不明的冰覆盖因子。$s_{\rm lq}+s_{\rm ice}$ 超过 1 时反应因子变负，这是无效状态，程序里应该报错，静默截断会把 bug 藏进结果。

原文文字把 $J_i$ 称为交换电流密度，式(5)的两个对称指数表明它是净反应电流，零过电位时为零，这和交换电流密度的定义对不上。表1至表5也没有给两极的 $J_{\rm ref}$、$C_{i,\rm ref}$、$r_i$ 及其温度依赖，电压定量复现的关键缺口就在这里，只能回到原文引用的动力学文献里补。

电势部分还有一处必须在动手前想清楚，符号体系。式(7)、(8)把两极过电位写成相反符号的相电势差，却没有说明是否扣除了平衡电位参考；式(1)把损失当带符号加项，式(9)明确 $V_{\rm ohm}\le0$；后面热源又要求耗散非负。三套约定互相打架。我的打算是在程序内部用氧化为正的净法拉第电流 $j_F$ 和常规 $\eta=\phi_s-\phi_m-U_{\rm eq}$ 定义，保证 $j_F\eta\ge0$，再建一张转换表对接原式，Part2 里把这张表写完整。

欧姆损失按各层电阻串联汇总

$$ V_{\rm ohm}=-(\Omega_{\rm BP}+\Omega_{\rm GDL} +\Omega_{\rm MPL}+\Omega_{\rm CL}+\Omega_m)I. \tag{9} $$

$\Omega$ 是面积比电阻，单位 $\mathrm{\Omega\,m^2}$，乘上面电流密度才是伏特。催化层里电子电流逐渐转成质子电流，两条通路各自都小于总电流，把两种电阻直接串联再乘同一个 $I$ 是额外近似；双侧各层是逐层计数还是两侧合并，原文也没有写明。

三、水的储存与迁移

膜水守恒原文写成

$$ \frac{\rho_m}{EW}\frac{\partial(\omega\lambda_{\rm nfm})}{\partial t} =\frac{\rho_m}{EW} \frac{\partial^2(\omega^{1.5}D_{\rm nfm}\lambda_{\rm nfm})}{\partial x^2} +S_{\rm nfm}. \tag{10} $$

左端的 $\omega$ 把储量折到总体积，右端的 $\omega^{1.5}$ 是传输连通性修正，两者作用不同，不能约掉。原文把系数与状态的乘积整体放进了二阶导数，变系数时 $\partial_{xx}(a\lambda)$ 与 $\partial_x(a\,\partial_x\lambda)$ 相差一项 $\partial_x(\lambda\,\partial_xa)$。我打算按面通量守恒形式离散

$$ N_{\rm diff}=-\frac{\rho_m}{EW}\omega^{1.5}D_{\rm nfm} \frac{\partial\lambda_{\rm nfm}}{\partial x},\qquad \partial_t q_{\rm nfm}=-\partial_xN_{\rm diff}+S_{\rm nfm}, $$

并把这一修改登记为建模解释，不宣称逐字复制。

扩散系数本身是分段经验式，记 $f_T=\exp[2416(1/303-1/T)]$

$$ D_{\rm nfm}=\begin{cases} 2.69266\times10^{-10},&\lambda\le2,\\ 10^{-10}f_T[0.87(3-\lambda)+2.95(\lambda-2)],&2<\lambda\le3,\\ 10^{-10}f_T[2.9514(4-\lambda)+1.642454(\lambda-3)],&3<\lambda\le4,\\ 10^{-10}f_T(2.563-0.33\lambda+0.0264\lambda^2-0.000671\lambda^3),&\lambda>4. \end{cases} \tag{11} $$

单位 $\mathrm{m^2/s}$。冷启动关心低温，问题正好集中在第一段，它没有温度因子。在 253.15 K，$D(2)=2.693\times10^{-10}$，跨过 $\lambda=2$ 后 $D(2^+)\approx1.810\times10^{-11}$，跳变约 14.9 倍，且 $\lambda=3$ 两侧端点系数 2.95 与 2.9514 也不一致。我第一反应是想把这段平滑掉，忍住了。原始引用文献还没核到，擅自修好一段经验式，后面和作者结果对不上时就说清来路了。先保留原式，登记为敏感项，高含水量段的多项式同样有负值外推风险，也要挡住。

低温下未冻水的容量上限

$$ \lambda_{\rm sat}=\begin{cases} 4.837,&T\lt 223.15\ \mathrm K,\\ [-1.304+0.01479T-3.594\times10^{-5}T^2]^{-1},&223.15\le T\lt T_{\rm freeze}, \end{cases} \tag{13} $$

在 $-20\ ^\circ\mathrm{C}$ 得 $\lambda_{\rm sat}\approx7.306$。原文还有一行 $T\ge T_{\rm freeze}$ 的分支，只写了 $\lambda_{\rm sat}>\lambda_{\rm nfm}$ 这个不等式，没有返回值，翻译成数值函数时要当成停止冻结的分支处理，不能随手挑一个极大含水量。这个阈值管的是膜内冻结驱动力，表3由水活度算出的 $\lambda_{\rm equil}$ 管的是膜气交换，两套经验式约束两个过程，名字看着像，不能互相覆盖。

液水走 Darcy 流动

$$ \frac{\partial(\varepsilon\rho_{\rm lq}s_{\rm lq})}{\partial t} =\frac{\partial^2\left(\rho_{\rm lq}\frac{K_{\rm lq}}{\mu_{\rm lq}}p_{\rm lq}\right)}{\partial x^2}+S_{\rm lq}, \tag{15} $$

原文用压力做主变量，好处是材料界面上压力协调、饱和度允许跳变，代价是各层 $K_{\rm lq}/\mu_{\rm lq}$ 不同时，乘积二阶导数和 Darcy 通量散度不等价。表1给的是本征渗透率 $K_0$，没有给 $K_{\rm lq}=K_0k_r(s)$ 的相对渗透率函数，干孔隙里液水怎么停流也缺闭合，这条留到 Part2 单列。

蒸气方程的麻烦在单位

$$ \frac{\partial[\varepsilon(1-s_{\rm lq}-s_{\rm ice})\rho_gC_{\rm vap}]}{\partial t} =\frac{\partial^2(\rho_gD_{\rm vap}^{\rm eff}C_{\rm vap})}{\partial x^2}+S_{\rm vap}. \tag{16} $$

$C_{\rm vap}$ 定义成 $\mathrm{mol/m^3}$，再乘气体密度 $\rho_g$，存储项的量纲就和质量源对不上。我的理解是直接采用气相体积摩尔浓度做状态量

$$ \frac{\partial(\varepsilon gC_{\rm vap})}{\partial t} =\frac{\partial}{\partial x}\left(D_{\rm vap}^{\rm eff}\frac{\partial C_{\rm vap}}{\partial x}\right)+S_{\rm vap}^{\rm mol}, $$

注意 $D^{\rm eff}$ 已经含了表3的孔隙与堵塞修正，面通量不能重复乘同样的 $\varepsilon^{1.5}g^{1.5}$。存储项里的 $g$ 本身随时间变，这个记账细节容易被忽略，液水冻结时密度又从 990 变到 $920\ \mathrm{kg/m^3}$，等质量的冰占的孔体积更大，只积分 $C_{\rm vap}$ 而不管 $g$ 的变化，库存迟早对不上。把 $\dot s_{\rm ice}$ 直接取成 $-\dot s_{\rm lq}$ 当质量守恒也是错的，等质量冻结时两者幅度差一个密度比，配平要在质量或摩尔库存上做。

四、产水、相变与拖曳的源项网络

电化学产水写作

$$ S_{\rm Prod}=\frac{J_i}{nF}. \tag{25} $$

这个 $n$ 在原文里承担了太多身份。产水计量应该用 $J_c^+/(2F)$，氢气消耗用 $J_a^+/(2F)$，氧气消耗用 $J_c^+/(4F)$，把阴极动力学的 $n=4$ 原封不动传给产水，水会少生成一半。表4用 kmol、正文用 mol，$EW$、$M_w$、$F$ 必须整套统一到同一物质量基准，式(20)里乘 $M_w$ 的换算才成立。在 $0.1\ \mathrm{A/cm^2}$ 下面积产水速率是 $I/(2F)\approx5.182\times10^{-3}\ \mathrm{mol/(m^2\,s)}$，这是不用解方程就该满足的整电池约束，我把它列为复现验收的第一条。

论文比较了三种产水归宿，生成于离聚物、孔隙液相或孔隙气相，作者的结论是膜水产水更合理，液相或气相产水会让冻结更早、运行时间更短。这个结论在它的方程和工况下有明确含义，拿来比较建模假设可以，推广成所有催化层都不会产气相液相水就过头了。实现上我打算用分配系数 $(f_m,f_v,f_l)$，$f_m+f_v+f_l=1$，三种机制就是三个纯分配端点，混合分配留给以后扩展。

电渗拖曳把水从阳极搬向阴极

$$ S_{\rm EOD}=n_d\frac{J_i}{nF},\qquad n_d=\frac{2.5}{22}\lambda_{\rm nfm}. \tag{23,24} $$

$\lambda=3.4$ 时 $n_d\approx0.386$。同样的分母问题，拖曳该跟质子通量 $i_H/F$ 走，用氧气的 $n=4$ 计算会小四倍。原文还只在两侧 CL 里放配对源项，膜区没有拖曳散度项，两侧 $J_i(x)$ 与 $n_d(\lambda)$ 各不相同，逐点代入并不自动等量搬运。我准备用贯穿离聚物通路的面通量 $N_{\rm EOD}=n_di_H/F$，让边界和电荷守恒自然配平两侧，这也登记为改动。

表4的六类相变是这篇文章的心脏，我逐条核下来，问题也集中在这里，挑三条说。

蒸气与液水之间的蒸发支

$$ S_{\rm vap\to lq}=\kappa_{\rm vap-lq}(s_{\rm lq}+s_{\rm ice})\frac{C_{\rm vap}-C_{\rm sat}}{RT}M_w, \qquad C_{\rm vap}\lt C_{\rm sat}. $$

两个问题。浓度差再除 $RT$，量纲对不上，若 $C$ 实为压力，$(p-p_{\rm sat})/RT$ 才是摩尔浓度差，我按后者理解并让 $C_{\rm sat}=p_{\rm sat}/RT$。蒸发供体乘的是 $s_{\rm lq}+s_{\rm ice}$，孔隙里只有冰没有液水时，这条支路还在通过式(20)消耗液水，能把液水库存推成负值。若意图是模拟升华，该从冰库存扣并用相应潜热；首轮实现我会把蒸发供体限制为液水，升华单独登记。

凝华支只允许蒸气变冰，没有反向升华，也没有 $T\lt T_{\rm freeze}$ 的关闭条件，模型算到冰点以上时这条通道何时停要自己规定。蒸气对过冷液水和对冰的饱和压力一般不同，原文用一个 $C_{\rm sat}$ 同时驱动凝结与凝华，是个额外简化，两条正源还竞争同一份蒸气库存，联立时要做供体耗尽检查。

膜水冻结那条最典型

$$ S_{\rm nfm\to fm}=\begin{cases} \kappa_{\rm nfm-fm}\dfrac{\rho_m}{EW}(\lambda_{\rm nfm}-\lambda_{\rm sat}), &\lambda_{\rm nfm}\ge\lambda_{\rm sat},\\[6pt] \kappa_{\rm nfm-fm}\dfrac{\rho_m}{EW}\lambda_{\rm fm}, &\lambda_{\rm nfm}<\lambda_{\rm sat}. \end{cases} $$

PDF 里下支没有负号。$\lambda_{\rm nfm}<\lambda_{\rm sat}$ 且 $\lambda_{\rm fm}>0$ 时，它仍让冻结水增长、未冻水减少，方向反了，按式(18)、(19)的守恒配对，融化支应为负。表2也没有单列 $\kappa_{\rm nfm-fm}$ 的值，另外几个速率都是 $1\ \mathrm{s^{-1}}$，不能顺推这个也是。

$\kappa=1\ \mathrm{s^{-1}}$ 的物理含义值得写一句。液冰转换两支都正比于供体库存，孤立冻结过程是一阶衰减，$1\ \mathrm{s^{-1}}$ 表示秒级松弛，并非一秒内全部冻完。有限冻结速率意味着低温下液水可以一直存在，只要产水或凝结在补。

所有内部交换加起来应该满足

$$ \sum_{k\in\{\rm nfm,fm,vap,lq,ice\}}S_k^{\rm internal}=0, $$

每条转换给供体减 $\Gamma$、给受体加 $\Gamma$。全电池积分后电渗拖曳也只是内部搬运，总水量变化只能由产水和外排解释。这条检查不依赖任何拟合参数，比曲线长得像更基础，我会把它做成程序里的运行时断言。

五、热量收支与潜热物性

能量方程

$$ \frac{\partial[(\rho C_p)_{\rm fl,sl}^{\rm eff}T]}{\partial t} =\frac{\partial^2(k_{\rm fl,sl}^{\rm eff}T)}{\partial x^2}+S_T. \tag{31} $$

热源按区域组合，式(32)的分区里 GDL、MPL、CL 那一行与阳极 CL、阴极 CL 两行重叠，合理解释是那一行只管 GDL 和 MPL，这条要登记修订，程序里用多个独立条件分支会重复加热。可逆热全部放在阴极也是空间分配近似，阳极反应快并不能从热力学上推出阳极熵热严格为零。导热项同样建议写成守恒面热流形式，理由和第三节相同。

欧姆热

$$ S_{T,\rm ohm}=\sigma_e^{\rm eff}\|\nabla V_e\|^2 +\sigma_H^{\rm eff}\|\nabla V_H\|^2 \tag{35} $$

在催化层外的区域只列了电子相，那一栏却包含不导电子的膜，膜内质子欧姆热被漏掉了。含水量决定膜电导率，这里正是水与热的一条直接连接，漏掉这项会同时破坏热场和能量核算。

可逆热的熵变，表5给 $-163110\ \mathrm{J/(kmol\,K)}$，以每摩尔氢气为反应基准应换算成 $-163.110\ \mathrm{J/(mol\,K)}$ 并用 $n=2$ 配反应速率。直接拿 kmol 数值配 $F=96487\ \mathrm{C/mol}$ 是一千倍误差，拿 $n=4$ 又少一半。

表5的潜热数值存在排版问题，字面值不能用于热计算，复现时应换成有来源的物性，并核对 $L_{\rm sub}=L_{\rm vap}+L_{\rm fus}$ 的相焓关系。我拿作者团队早年公开博士论文中的 $h_{\rm fusn}=3.336\times10^5\ \mathrm{J/kg}$ 做过一次上限检查。$0.1\ \mathrm{A/cm^2}$、阴极 CL 厚 $10\ \mathrm{\mu m}$ 时平均体积电流约 $10^8\ \mathrm{A/m^3}$，可逆热源约 $2.14\times10^7\ \mathrm{W/m^3}$；体积产水质量速率约 $9.3\ \mathrm{kg/(m^3\,s)}$，乘上潜热后，产水全部冻结的潜热热源约 $3.11\times10^6\ \mathrm{W/m^3}$，与可逆热源同量级，复现时相变热这一项值得保留。

六、两个不用求解器的数量级检查

动手写代码前我算了两个手算例，用途是查计量数、孔隙率和微米换算有没有出错。

阴极 CL 离聚物从初始含水量装到 $\lambda_{\rm sat}$ 的储水容量

$$ \Delta n_{m,\rm CL}=\delta_{\rm CL}\,\omega\frac{\rho_m}{EW}(\lambda_{\rm sat}-\lambda_0), $$

在 253.15 K、$\lambda_0=3.4$ 时约 $0.02812\ \mathrm{mol/m^2}$，按基准产水速率填满约需 5.4 秒，$\lambda_0=2.9$ 时约 6.1 秒。若阴极 CL 孔隙全部由冰占据

$$ n_{\rm ice,max}=\frac{\delta_{\rm CL}\varepsilon\rho_{\rm ice}}{M_w} \approx0.153\ \mathrm{mol/m^2}, $$

全部产水留在这里并成冰，填满约需 29.6 秒。两个数都忽略了跨层输运、相变、拖曳和升温，不能拿去和图4的运行总时长直接比。它们的价值在另一处，如果程序算出的填充时间比这两个数差出一个量级，先查自己的电子数和单位，再怀疑模型。

七、从方程里读出的物理图景

把方程连起来，冷启动的四段演化自然浮出来。产水先填未饱和的离聚物，膜水沿膜扩散、被电渗拖曳重新分配；$\lambda_{\rm sat}$ 随温度、$\lambda_{\rm nfm}$ 随位置都在变，各点到达阈值的时间不同；有限脱附速率移不走新增水，出现 $\lambda_{\rm nfm}>\lambda_{\rm sat}$ 的过饱和，膜内离子通道还湿着，孔隙冰还没完全建立；最后冰占孔隙，$g$ 下降，氧气补给和反应位点一起受限，恒流要求更大的损失，直到停机。四个阶段在局部时间上重叠，实现成四段互斥、硬切换的程序逻辑就错了。

网络里有几组竞争关系，读文献时容易看漏。膜水增多改善离子导电，同时让离聚物更快逼近未冻阈值，电压的早期改善与后期恶化可以来自同一条产水路径的不同阶段。提高电流同时提高产水与发热，更快升温帮助融冰，更快结冰提前堵塞，输赢取决于热惯量、散热、传水距离和相变速率，从发热更多推不出更容易启动成功。初始膜更干储水空间更大，电导率也更低，图4、5的比较发生在特定工况，不能扩展成越干越好的普遍结论。

作者提出的平衡模型版本也值得记一笔。它按膜水、蒸气、液水的顺序安置新增水，容量满了转入下一态。实现时如果只把 $\lambda$ 截到 $\lambda_{\rm sat}$，却不把超额水挪进其他库存，平衡处理就变成凭空删水；阈值随温度和气相体积变化时，已有库存也要参与重分配。还有一个直接的阈值检查，$-20\ ^\circ\mathrm{C}$ 时 $\lambda_{\rm sat}\approx7.306$，而气相活度接近 1 时 $\lambda_{\rm equil}\approx14$，把 7.306 代回平衡关系反解得到水活度约 0.804。膜水到未冻阈值时气相并未饱和，两套经验关系的共同作用推不出统一化学势意义上的全相平衡，作者的顺序分配是一个明确的近似选择，复现时就按它做，标注清楚。

$\kappa_{\rm nfm-lq}$ 是图8主要扫描的参数，取 1、2、3 $\mathrm{s^{-1}}$，调大它膜水过饱和减小、孔隙液水供应增加、冻结堵塞更早。但 $g\to0$ 时原式又把脱附通道关掉，接近堵塞的局部会保留过多膜水。停机前电压偏高到底归因于这个参数，还是归因于传输闭合或初态，要等 Part3 的结果说话。

八、对标口径与下一篇要做的事

原文的四张靶图，对标时口径要先钉死。图4是 $-20\ ^\circ\mathrm{C}$、$0.1\ \mathrm{A/cm^2}$、$\lambda_0=3.4$，输出电压和 CL 五个位置的含水量，比较时要区分局部值与平均值；图5换成 $\lambda_0=2.9$，输出平均含水量和冰占比，冰占比的空间统计口径原文没有写明，不能默认当成局部峰值；图8比较平衡与非平衡及三个脱附速率，其中 CL 最大含水量那条不能换成平均量；图9是三种产水机制下的空间分布，CL 横坐标左侧接膜、右侧接 MPL。图9里液水积累位置会随时间向 CL 内部迁移，这正说明平均量方程承载不了空间过程，也是我认定这篇文章适合做一维复现模板的原因之一。

对标实验还要记得 Tajiri 等 2007 年那组数据的条件。实验用大热容夹具，启动期间近似等温；干吹扫也会留下孔隙残水和不均匀分布。模型的理想初值和均匀初值都带近似，停机时间的偏差不能全算在脱附系数头上。

留给 Part2 的清单我按优先级记在这里。建符号转换表统一三套正负号约定；全模型换 SI 摩尔基准；传输方程全部采用面通量守恒形式离散；补相对渗透率与干湿界面闭合；饱和蒸气压二次项指数按 $10^{-5}$ 修正，作者团队的博士论文可以佐证；潜热换成有来源物性；定加载协议，解决零反应物初值与恒流加载的冲突；准备一份原文表达与实现解释的对照表，能修的记依据，多解的留分支。

这篇文章读到这里，我对它的评价没有变。模型骨架是好的，坑集中在书面上。把书面的账理清，比写代码本身更接近复现的一半。动 Python 之前，下一篇先把符号和单位钉死。

来源

  1. Huo S, Jiao K, Park J W. On the water transport behavior and phase transition mechanisms in cold start operation of PEM fuel cell. Applied Energy, 2019, 233-234, 776-788. doi.org/10.1016/j.apenergy.2018.10.068
  2. Jiao K. Experimental and Modelling Studies of Cold Start Processes in Proton Exchange Membrane Fuel Cells. PhD thesis, University of Waterloo, 2011. 大学资料库
  3. Tajiri K, Wang Y, Tabuchi Y, et al. Effects of operating and design parameters on PEFC cold start. Journal of Power Sources, 2007, 165, 279-286. 作者机构全文