减少碳排放以遏制全球变暖是目前人类面临的巨大挑战。为解决传统化石燃料高碳排放的问题，有必要寻找更加清洁的新能源。氢能由于可持续、绿色清洁和能量密度高等优点，是未来低碳社会的理想能源之一。可再生能源由于间歇性和波动性的特点难以直接并网消纳。利用水电解将其转化为氢能，既便于可再生能源利用，又比化石燃料制氢更清洁，是未来制氢的理想途径^[1]。

静态供水式水电解最初因宇航员生命保障系统的需求而提出^{[2, 3]}。其实质是膜蒸馏与水电解的耦合。如图 1所示，疏水多孔膜将液态供水与电极隔开，在化学势差推动下水以蒸气的形式由供水传递至电解质以满足电解消耗。由于静态供水式水电解中供水不与电极和电解质直接接触，因此可避免常规水电解中可能存在的电解质污染与催化剂中毒等问题，大幅降低对水质的要求^[4]。与水蒸气电解相比，又无需额外的蒸汽发生过程。

由于供水方式的特殊性，水的供需变化会影响电解质的浓度与体积，进而影响电解性能。电解质层水含量的变化应处于适当范围内，防止系统因水淹或缺水而失稳。目前关于静态供水式水电解的研究较少，且多为实验研究^[2-6]，尤其缺少对其模拟工作的报道。数值建模有助于深入理解静态供水式水电解中发生的传递过程，是探究水的传递对电解性能影响的有效手段。关于膜蒸馏^[7-9]、常规水电解^[10-12]和水蒸气电解^{[13, 14]}的模拟研究已有诸多进展，但关于静态供水式水电解的模拟则仅见Lin等^[15]的报道。该工作建立了静态供水式质子交换膜水电解的一维等温模型，但并未对系统内膜蒸馏和水电解间的耦合关系进行建模，也没有考虑供水推动力与阻力变化对电解性能及系统稳定性的影响。

本研究首次建立了以KOH溶液为电解质的静态供水式碱性水电解(Static Feed Alkaline Water Electrolysis, SFAWE)二维数值模型，全面考虑了系统内部质量、动量、能量和电荷传递以及电化学反应之间的耦合关系。通过实验测定的极化曲线验证了模拟结果的准确性。计算得到了系统温度、压力、浓度和通量的分布，并对其合理性进行了讨论。着重考察了供水推动力和阻力的变化对电解性能的影响。分析了反应活化、浓度极化与欧姆阻抗导致的各部分过电位在总电压中的占比。探究了不同供水温度下电解质层KOH含量和阴极气体间隙层厚度对总电压与适宜运行工况范围的影响，为设计和过程优化提供参考。

1 模型与求解 1.1 基本概述

模型中SFAWE系统被划分为如下区域：液态供水流道(Feed Water Channel, FWC)、疏水多孔膜(Hydrophobic Porous Membrane, HPM)、阴极气体间隙层(Cathode Gas Gap Layer, CGGL)、阴极催化层(Cathode Catalyst Layer, CCL)、电解质层(Electrolyte Layer, EL)、阳极催化层(Anode Catalyst Layer, ACL)和阳极气体间隙层(Anode Gas Gap Layer, AGGL)，如图 1所示。图 1中粗黑线表示不可透壁面。气体间隙层填充多孔介质支撑。电解质层为多孔隔膜支撑的KOH溶液。电极反应为：

$\text{Cathode}:\ {{\text{H}}_{2}}\text{O}+{{\text{e}}^{-}}\to \frac{1}{2}{{\text{H}}_{2}}\uparrow +\text{O}{{\text{H}}^{-}}$

(1)

$\text{Anode:}\ \text{O}{{\text{H}}^{-}}\to \frac{1}{2}{{\text{H}}_{2}}\text{O}+\frac{1}{4}{{\text{O}}_{2}}\uparrow +{{\text{e}}^{-}}$

(2)

为简化计算，模型的基本假设如下：1)气相为理想气体；2)忽略重力影响；3)各区域流体流动为单相层流；4)仅液态供水流道与电解质层存在液相；5)在模型讨论的范围内，忽略电解质体积变化对气相质量与动量传递的影响；6)电极反应集中于电极催化层与电解质层的公共边上(CCL-EL边界和ACL-EL边界)；7)气态水在CCL-EL界面处冷凝进入电解质溶液后才会被电解反应消耗；8)忽略电解质层两侧气相压差对电解质溶液质量与动量传递的影响，不考虑电解质溶液的压力变化；9)忽略电解质层两侧H₂与O₂的交叉混合，阴极和阳极侧气相分别为湿氢气与湿氧气，透过疏水多孔膜溶解于液态供水的H₂的量也忽略不计。

1.2 方程与边界条件 1.2.1 质量与动量传递

对于流动相α(气相g，液相l)，总质量守恒与各组分i的质量守恒方程分别为：

$\nabla \cdot {\mathit{\boldsymbol{N}}_{\rm{ \mathsf{ α} }}} = \nabla \cdot \left( {{\rho _{\rm{ \mathsf{ α} }}}{\mathit{\boldsymbol{v}}_{\rm{ \mathsf{ α} }}}} \right) = 0$

(3)

$\nabla \cdot {\mathit{\boldsymbol{N}}_i} = \nabla \cdot \left[ {{\mathit{\boldsymbol{j}}_i} + {\rho _{\rm{ \mathsf{ α} }}}w(i){\mathit{\boldsymbol{v}}_{\rm{ \mathsf{ α} }}}} \right] = 0$

(4)

式(3)和式(4)中：N_α为总质量通量；ρ_α为密度；v_α为质量平均速度；N_i为组分i的总质量通量；j_i为组分i的质量扩散通量；w(i)为组分i质量分数。

对于气相g，各组分的质量扩散通量通过Dusty Gas模型计算^[16]。对于孔径较小的疏水多孔膜，考虑Knudsen扩散的影响：

$\begin{array}{l} \;\;\;\;\;\;\;\;\sum\nolimits_{k \ne i} {\frac{{n(i)n(k)}}{{{D_{ik, {\rm{K}}, {\rm{eff}}}}}}} \left( {\frac{{{\mathit{\boldsymbol{j}}_i}}}{{{\rho _i}}} - \frac{{{\mathit{\boldsymbol{j}}_k}}}{{{\rho _k}}}} \right) = \\ - \nabla n(i) - \frac{{{\nabla _{\rm{g}}}}}{{{p_{\rm{g}}}}}n(i)\left\{ {1 - \frac{{1/{D_{{\rm{K}}, i, {\rm{eff}}}}}}{{\sum\nolimits_k {\left[ {n(k)/{D_{{\rm{K}}, k, {\rm{eff}}}}} \ \ \right]} }}} \right\} \end{array}$

(5)

${D_{ik, {\rm{K}}, {\rm{ eff }}}} = \frac{1}{{{D_{ik, {\rm{ eff }}}}}} + \frac{1}{{{D_{{\rm{K}}, i, {\rm{ eff }}}}{D_{{\rm{K}}, k, {\rm{ eff }}}}\sum\nolimits_j {\left[ {n(j)/{D_{{\rm{K}}, j, {\rm{ eff }}}}} \right]} }}$

(6)

式(5)和式(6)中：n(i)为组分i的摩尔分数；p_g为气相压力；D_{ik, K, eff}为包含Knudsen扩散影响组分i、k间有效二元扩散系数。对于孔径较大的气体间隙层与电极催化层，采用大孔条件简化方程^[16]：

$\begin{array}{c} \sum\nolimits_{k \ne i} {\frac{{n(i)n(k)}}{{{D_{ik, {\rm{eff}}}}}}} \left( {\frac{{{\mathit{\boldsymbol{j}}_i}}}{{{\rho _i}}} - \frac{{{\mathit{\boldsymbol{j}}_k}}}{{{\rho _k}}}} \right) = \\ - \nabla n(i) - \frac{{{\nabla _{{{\rm{g}}_{\rm{g}}}}}}}{{{p_{\rm{g}}}}}n(i)\left\{ {1 - \frac{{\sqrt {{M_i}} }}{{\sum\nolimits_k {\left[ {n(k)\sqrt {{M_k}} } \right]} }}} \right\} \end{array}$

(7)

有效二元扩散系数D_{ik, eff}与有效Knudsen扩散系数D_{K, i, eff}均按Bruggeman修正计算：

${D_{ik, {\rm{eff}}}} = \frac{{{\varepsilon _{\rm{g}}}}}{{{\tau _{\rm{g}}}}}{D_{ik}}, \quad {D_{{\rm{K}}, i, {\rm{eff}}}} = \frac{{{\varepsilon _{\rm{g}}}}}{{{\tau _{\rm{g}}}}}{D_{{\rm{K}}, i}}, \quad {\tau _{\rm{g}}} = \varepsilon _{\rm{g}}^{ - {\beta _{\rm{g}}}}$

(8)

式(8)中：ε_g为气相体积分数；τ_g为曲折因子；β_g为关联系数；D_ik为组分i、k间的二元扩散系数。Knudsen扩散系数D_{K, i}按式(9)计算。

${D_{{\rm{K}}, i}} = \frac{{{d_{\rm{s}}}}}{3}\sqrt {\frac{{8RT}}{{{\rm{ \mathsf{ π} }}{M_i}}}} $

(9)

式(9)中：d_s为多孔介质孔径；M_i为组分i的摩尔质量；R为气体常数；T为温度；π为圆周率。

气相的流动通过Brinkman方程描述：

$\begin{array}{c} \frac{{{\rho _{\rm{g}}}}}{{{\varepsilon _{\rm{g}}}}}\nabla \left( {\frac{{{\mathit{\boldsymbol{v}}_{\rm{g}}} \cdot {\mathit{\boldsymbol{v}}_{\rm{g}}}}}{{{\varepsilon _{\rm{g}}}}}} \right) = - \nabla {p_{\rm{g}}} + \nabla \cdot \\ \left\{ {\frac{1}{{{\varepsilon _{\rm{g}}}}}\left\{ {{\mu _{\rm{g}}}\left[ {\nabla {\mathit{\boldsymbol{v}}_{\rm{g}}} + {{\left( {\nabla {\mathit{\boldsymbol{v}}_{\rm{g}}}} \right)}^{\rm{T}}}} \right] - \frac{2}{3}{\mu _{\rm{g}}}\left( {\nabla \cdot {\mathit{\boldsymbol{v}}_{\rm{g}}}} \right)\mathit{\boldsymbol{I}}} \right\}} \right\} - \frac{{{\mu _{\rm{g}}}}}{{{\mathit{\boldsymbol{\kappa }}_{\rm{g}}}}}{\mathit{\boldsymbol{v}}_{\rm{g}}} \end{array}$

(10)

式(10)中：μ_g为气相黏度；κ_g为渗透率；I为单位张量。

对于电解质层KOH溶液，各组分质量扩散通量通过Newman等^{[17, 18]}的浓物质传递模型计算：

${{\mathit{\boldsymbol{j}}}_ + } = - \frac{{{M_ + }{\rho _1}{D_{{\rm{KOH, eff }}}}}}{{{M_{{\rm{KOH }}}}}}\nabla w({\rm{KOH}}) + \frac{{{\mathit{\boldsymbol{i}}_1}{t_ + }{M_ + }}}{F}$

(11)

${\mathit{\boldsymbol{j}}_ - } = - \frac{{{M_ - }{\rho _1}{D_{{\rm{KOH}}, {\rm{eff}}}}}}{{{M_{\rm{KOH}}}}}\nabla w({\rm{KOH}}) - \frac{{{\mathit{\boldsymbol{i}}_1}{t_ - }{M_ - }}}{F}$

(12)

${\mathit{\boldsymbol{j}}_{\rm{w}}} = - {\rho _1}{D_{{\rm{KOH}}, {\rm{eff}}}}\nabla w({\rm{w}}) - \left( {{t_ + }{M_ + } - {t_ - }{M_ - }} \right)\frac{{{\mathit{\boldsymbol{i}}_1}}}{F}$

(13)

式(11)至式(13)中：+、－、w分别代表K⁺、OH^－和水；i_l为液相电流密度；F为Faraday常数；t₊、t_－分别为K⁺和OH^－相对于质量平均速度的迁移数。

${t_ + } = w( - ) + w({\rm{w}})t_ + ^0, {t_ - } = w( + ) + w({\rm{w}})t_ - ^0$

(14)

式(14)中：t₊⁰、t_－⁰分别为K⁺和OH^－迁移数。有效整体扩散系数D_{KOH, eff}按Bruggeman修正计算：

${D_{{\rm{KOH}}, {\rm{eff}}}} = \frac{{{\varepsilon _1}}}{{{\tau _1}}}{D_{{\rm{KOH}}}}, \quad {\tau _1} = \varepsilon _1^{ - {\beta _1}}$

(15)

式(15)中：ε_l为液相体积分数；τ_l为曲折因子；β_l为关联系数；D_KOH为KOH的整体扩散系数。液相体积分数ε_l的大小反映了电解质层的水含量，由式(16)计算：

${\varepsilon _1} = \frac{{{V_1}}}{{{V_1} + {V_{\rm{s}}}}}$

(16)

隔膜骨架体积V_s和KOH溶液体积V_l分别为：

${V_{\rm{s}}} = \frac{{{S_{{\rm{EL}}}}{\rho _{{\rm{EL}}, {\rm{ area }}}}}}{{{\rho _{{\rm{EL}}}}}}, {V_1} = \frac{{{S_{{\rm{EL}}}}{\delta _{{\rm{EL}}}}{\varepsilon _{{\rm{EL}}}}{\rho _{1, {\rm{ fill }}}}w{{({\rm{KOH}})}_{{\rm{fill }}}}}}{{{\rho _1}w({\rm{KOH}})}}$

(17)

式(17)中：S_EL为隔膜面积；ρ_{EL, area}为隔膜面密度，ρ_EL为隔膜骨架密度；δ_EL为电解质层厚度；ε_EL为隔膜孔隙率；w(KOH)_fill为KOH溶液刚好充满隔膜孔隙时的质量分数，反映了电解质层的KOH含量；ρ_{l, fill}为w(KOH)_fill对应的KOH溶液的密度。

本模型中假定气液界面处无滑移，各组分均沿法向流入或流出。对于气相，FWC-HPM边界处H₂总通量为零，CCL-EL边界处的H₂和ACL-EL边界处的O₂的法向总通量大小均由Faraday定律计算：

${N_{{{\rm{H}}_2}}} = - \frac{{{i_{{\rm{loc}}, {\rm{c}}}}{M_{{{\rm{H}}_2}}}}}{{2F}}$

(18)

${N_{{{\rm{O}}_2}}} = \frac{{{i_{{\rm{loc}}, {\rm{a}}}}{M_{{{\rm{O}}_2}}}}}{{4F}}$

(19)

式(18)和式(19)中：i_loc为局部电流密度；下标a、c分别代表阳极和阴极。各气液界面处水蒸气的总通量大小N_vapor通过气液相变速率计算。FWC-HPM边界处：

${N_{{\rm{vapor}}, {\rm{wm}}}} = {k_{{\rm{g}} - 1}}\left[ {{p_{{\rm{sat}}}} - {p_{{\rm{g}}, {\rm{wm}}}}n{{({\rm{ vapor }})}_{{\rm{wm}}}}} \right]$

(20)

式(20)中：下标wm代表FWC-HPM边界；k_g-l为相变总传质系数。CCL-EL边界处：

${N_{{\rm{vapor, c}}}} = - {k_{{\rm{g}} - 1}}\left[ {{p_{{\rm{g}}, {\rm{c}}}}n{{({\rm{ vapor }})}_{\rm{c}}} - {p_{{\rm{sat}}, {\rm{KO}}{\rm{H}}, {\rm{c}}}}} \right]$

(21)

ACL-EL边界处：

${N_{{\rm{vapor}}, {\rm{a}}}} = {k_{{\rm{g}} - 1}}\left[ {{p_{{\rm{sat}}, {\rm{KOH}}, {\rm{a}}}} - {p_{{\rm{g}}, {\rm{a}}}}n{{({\rm{ vapor }})}_{\rm{a}}}} \right]$

(22)

p_sat和p_{sat, KOH}分别为纯水和KOH溶液的饱和蒸汽压。本模型中通过调整k_g-l取值足够大(10^-5 s·m^-1)，即水的气液相变传质阻力很小，来近似水的相平衡。气相出口压力为常压，即101.325 kPa。

对于电解质层KOH溶液，结合式(27)、式(28)以及Faraday定律可知，液态水在CCL-EL边界和ACL-EL边界处的法向总通量大小分别为：

${N_{{\rm{w}}, {\rm{c}}}} = \frac{{{i_{{\rm{loc}}, {\rm{c}}}}{M_{{{\rm{H}}_2}{\rm{O}}}}}}{F} + {k_{{\rm{g}} - 1}}\left[ {{p_{{\rm{g}}, {\rm{c}}}}n{{({\rm{ vapor }})}_{\rm{c}}} - {p_{{\rm{sat}}, {\rm{KOH}}, {\rm{c}}}}} \right]$

(23)

${N_{{\rm{w}}, {\rm{a}}}} = \frac{{{i_{{\rm{loc }}, {\rm{a}}}}{M_{{{\rm{H}}_2}0}}}}{{2F}} - {k_{{\rm{g}} - 1}}\left[ {{p_{{\rm{sat}}, {\rm{K}}0{\rm{H}}, {\rm{a}}}} - {p_{{\rm{g}}, {\rm{a}}}}n{{({\rm{ vapor }})}_{\rm{a}}}} \right]$

(24)

K⁺仅存于溶液中，在CCL-EL边界和ACL-EL边界处总通量均为零。CCL-EL边界和ACL-EL边界处OH^－法向总通量大小均由Faraday定律计算：

${N_{ - , {\rm{c}}}} = - \frac{{{i_{{\rm{loc}}, {\rm{c}}}}{M_ - }}}{F}$

(25)

${N_{ - , {\rm{a}}}} = - \frac{{{i_{{\rm{loc}}, {\rm{a}}}}{M_ - }}}{F}$

(26)

气相和液相所在区域各壁面均设为无滑移。

1.2.2 电荷传递

多孔电极电导率很高，视为等势体，各处电势均为φ_s。电解质层KOH溶液电荷守恒方程为：

$\nabla \cdot {\mathit{\boldsymbol{i}}_1} = 0$

(27)

离子电流密度i_l为^[18]：

$\begin{array}{c} {\mathit{\boldsymbol{i}}_1} = - {\sigma _{1, {\rm{ eff }}}}\nabla {\varphi _1} - \frac{{\left( {{\nu _ + } + {\nu _ - }} \right)RT{\sigma _{1, {\rm{ eff }}}}}}{{w({\rm{W}})w({\rm{KOH}})F}}\\ \left[ {\frac{{t_ + ^0}}{{{z_ + }{\nu _ + }}} + \frac{{w({\rm{KOH}}){M_{{{\rm{H}}_2}{\rm{O}}}}}}{{2w({\rm{W}}){M_{{\rm{KOH}}}}}}} \right]\left( {1 + \frac{{\partial \ln {\gamma _ \pm }}}{{\partial \ln {m_{{\rm{KOH}}}}}}} \right)\nabla w({\rm{KOH}}) \end{array}$

(28)

式(28)中：φ_l为KOH溶液的电势；γ_±为KOH的平均活度系数；m_KOH为KOH的质量摩尔浓度。KOH溶液的有效电导率σ_{l, eff}按Bruggeman修正计算：

${\sigma _{1, {\rm{eff}}}} = \frac{{{\varepsilon _1}}}{{{\tau _1}}}{{{\sigma }}_1}$

(29)

其中：σ_l为KOH溶液本征电导率。电解质层欧姆过电位η_ohm通过阳极和阴极的液相电势差计算：

${\eta _{{\rm{ohm}}}} = {\varphi _{1, {\rm{a}}}} - {\varphi _{1, {\rm{c}}}}$

(30)

阴极电位设为0 V，阳极电位为总电压U_cell。

1.2.3 电化学反应

本模型以参照浓度c_ref和实际温度T下的平衡电极电位E_{eq, T}作为过电位的计算基准：

${E_{{\rm{eq }}, T}} = {E_{{\rm{eq }}, {\rm{ ref }}}} + \frac{{{\rm{d}}{E_{{\rm{eq }}}}}}{{{\rm{d}}T}}\left( {T - {T_{{\rm{ref }}}}} \right)$

(31)

式(31)中：E_{eq, ref}为参照状态(参照浓度c_ref、参照温度T_ref)下的平衡电极电位；dE_eq/dT为平衡电极电位的温度系数。阴阳极过电位均由式(32)计算：

$\eta = {\varphi _{\rm{s}}} - {\varphi _1} - {E_{{\rm{eq}}, T}}$

(32)

不同于电极完全浸没于溶液中的常规水电解，SFAWE中电极有效面积会受到电解质层水含量影响。本模型对有效交换电流密度i_{0, T, eff}进行了修正：

${i_{0, T, {\rm{ eff }}}} = a\varepsilon _1^b \cdot {i_{0, T}}$

(33)

${i_{0, T}} = {i_{0, {\rm{ ref }}}}\exp \left[ { - \frac{{{E_{{\rm{act }}}}}}{R}\left( {\frac{1}{T} - \frac{1}{{{T_{{\rm{ref }}}}}}} \right)} \right]$

(34)

式(34)中：i_{0, T}和i_{0, ref}分别为基准状态和参照状态下的本征交换电流密度；E_act为电极反应的活化能；关联系数a、b通过拟合实验数据获取。交换电流密度按电极的表观面积计算。

阴极和阳极的过电位η与局部电流密度i_loc的关系均通过Butler-Volmer方程描述：

$\begin{array}{c} {i_{1{\rm{oc}}}} = a\varepsilon _1^b \cdot \\ {i_{0, T}}\left[ {{C_{\rm{R}}}\exp \left( {\frac{{{\alpha _{\rm{f}}}F\eta }}{{RT}}} \right) - {C_0}\exp \left( {\frac{{ - \left( {1 - {\alpha _{\rm{f}}}} \right)F\eta }}{{RT}}} \right)} \right] \end{array}$

(35)

其中α_f为电极反应的正向传递系数，阴极和阳极反应中还原态和氧化态组分的浓度商C_R、C_O为：

${C_{{\rm{R}}, {\rm{c}}}} = \left( {\frac{{{c_{ - , {\rm{c}}}}}}{{{c_{ - , {\rm{ref}}}}}}} \right){\left( {\frac{{{p_{{{\rm{H}}_2}}}}}{{{p_{{{\rm{H}}_2}, {\rm{ref}}}}}}} \right)^{\frac{1}{2}}}, {C_{{\rm{R}}, {\rm{a}}}} = {\left( {\frac{{{c_{ - , {\rm{a}}}}}}{{{c_{ - , {\rm{ref}}}}}}} \right)}$

(36)

${C_{0, {\rm{c}}}} = \left( {\frac{{{c_{{\rm{w}}, {\rm{c}}}}}}{{{c_{{\rm{w}}, {\rm{ref}}}}}}} \right), \quad {C_{0, {\rm{a}}}} = {\left( {\frac{{{c_{{\rm{w}}, {\rm{a}}}}}}{{{c_{{\rm{w}}, {\rm{ref}}}}}}} \right)^{\frac{1}{2}}}{\left( {\frac{{{p_{{0_2}}}}}{{{p_{{{\rm{O}}_2}, {\rm{ref}}}}}}} \right)^{\frac{1}{4}}}$

(37)

阴极和阳极的过电位η均可分解为活化过电位η_act和浓度极化过电位η_conc 2部分。

$\eta = {\eta _{{\rm{conc}}}} + {\eta _{{\rm{act}}}}$

(38)

若以电极表面附近参与反应各组分的实际浓度c和实际温度T对应的平衡电极电位E_eq作为过电位的计算基准，此时的过电位即为考虑电解质层水含量影响的实际活化过电位η_act：

${\eta _{{\rm{act}}}} = {\varphi _{\rm{s}}} - {\varphi _1} - {E_{{\rm{eq}}}}$

(39)

根据Nernst方程：

${E_{{\rm{eq}}}} = {E_{{\rm{eq}}, T}} - \frac{{RT}}{F}\ln \frac{{{C_{\rm{R}}}}}{{{C_{\rm{O}}}}}$

(40)

由式(32)、(38)~(40)，浓度极化过电位η_conc为：

${\eta _{{\rm{conc }}}} = \eta - {\eta _{{\rm{act }}}} = - \frac{{RT}}{F}\ln \frac{{{C_{\rm{R}}}}}{{{C_{\rm{O}}}}}$

(41)

忽略浓度极化影响后B-V方程转化为：

$\begin{array}{c} {i_{{\rm{loc }}}} = a\varepsilon _1^b \cdot \\ {i_0}\left[ {\exp \left( {\frac{{{\alpha _{\rm{f}}}F{\eta _{{\rm{act}}}}}}{{RT}}} \right) - \exp \left( {\frac{{ - \left( {1 - {\alpha _{\rm{f}}}} \right)F{\eta _{{\rm{act}}}}}}{{RT}}} \right)} \right] \end{array}$

(42)

${i_0} = {i_{0, T}}C_{\rm{R}}^{\left( {1 - {\alpha _{\rm{f}}}} \right)}C_{\rm{O}}^{{\alpha _{\rm{f}}}}$

(43)

式(42)中：i₀为实际浓度c和实际温度T对应的本征交换电流密度。若再忽略电解质层内水含量对活化过电位的影响，此时的过电位即为本征活化过电位η_{act, 0}，B-V方程的形式转化为：

${i_{{\rm{loc }}}} = {i_0}\left[ {\exp \left( {\frac{{{\alpha _{\rm{f}}}F{\eta _{{\rm{act}}, 0}}}}{{RT}}} \right) - \exp \left( {\frac{{ - \left( {1 - {\alpha _{\rm{f}}}} \right)F{\eta _{{\rm{act}}, 0}}}}{{RT}}} \right)} \right]$

(44)

η_act与η_{act, 0}的差值，即活化过电位中与电解质层水含量相关联的部分记为：

$\mathit{\Delta }{\eta _{{\rm{act}}}} = {\eta _{{\rm{act}}}} - {\eta _{{\rm{act}}, 0}}$

(45)

阴极析氢反应使用高活性的Pt/C催化剂，其本征活化过电位η_{act, 0, c}可忽略。阳极析氧反应的本征活化过电位η_{act, 0, a}较高，B-V方程可简化为Tafel形式：

${\eta _{{\rm{act}}, 0, {\rm{a}}}} = \frac{{RT}}{{n{\alpha _{{\rm{f}}, {\rm{a}}}}F}}\ln \left( {\frac{{{i_{{\rm{loc}}, {\rm{a}}}}}}{{{i_{0, {\rm{a}}}}}}} \right)$

(46)

由式(46)，阴极和阳极的过电位均可分解为：

$\eta = {\eta _{{\rm{conc }}}} + {\eta _{{\rm{act }}}} = {\eta _{{\rm{conc }}}} + {\eta _{{\rm{act }}, 0}} + \Delta {\eta _{{\rm{act }}}}$

(47)

从而总电压U_cell可分解为：

${U_{{\rm{cell }}}} = {E_{{\rm{eq }}, T}} + {\eta _{{\rm{conc }}}} + {\eta _{{\rm{act }}, 0}} + \Delta {\eta _{{\rm{act }}}} + {\eta _{{\rm{ohm }}}}$

(48)

1.2.4 能量传递

各区域的总能量守恒方程为：

${\rho _{\rm{ \mathsf{ α} }}}{C_{{\rm{p}}, {\rm{ \mathsf{ α} }}}}{\mathit{\boldsymbol{v}}_{\rm{ \mathsf{ α} }}} \cdot \nabla T - \nabla \cdot \left( {{k_{{\rm{eff}}}}\nabla T} \right) = {Q_{{\rm{total}}}}$

(49)

${k_{{\rm{eff}}}} = {\varepsilon _{\rm{s}}}{k_{\rm{s}}} + {\varepsilon _{\rm{ \mathsf{ α} }}}{k_{\rm{ \mathsf{ α} }}}$

(50)

式(49)中：Q_total为区域总热源；C_{p, α}为流动相α的热容；k_eff为有效平均热导率；ε_s为固相体积分数；k_s为固相热导率；k_α为流动相α热导率。本模型中假定液态供水流速足够快，可及时交换系统内部传递出的热量，故FWC被视作等温域，温度等于供水温度T_water，不作为模型实际计算区域。

电极焦耳热可忽略。电解质层焦耳热为：

${Q_{\rm{J}}} = - {\mathit{\boldsymbol{i}}_1} \cdot \nabla {\varphi _1}$

(51)

CCL-EL边界和ACL-EL边界处由于水的相变以及电极反应和过电位存在，边界热源大小分别为：

$\begin{array}{c} {q_{\rm{c}}} = {k_{{\rm{g}} - 1}}\left[ {{p_{{\rm{g}}, {\rm{c}}}}n{\rm{ (vapor) }} - {p_{{\rm{sat}}, {\rm{KOH}}, {\rm{c}}}}} \right] \cdot \Delta {H_{{\rm{g}} - 1}} + \\ \left[ {{\eta _{\rm{c}}} + {T_{\rm{c}}}{{\left( {\frac{{{\rm{d}}{E_{{\rm{eq}}}}}}{{{\rm{d}}T}}} \right)}_{\rm{c}}}} \right]{i_{{\rm{loc}}, {\rm{c}}}} \end{array}$

(52)

$\begin{array}{c} {q_{\rm{a}}} = - {k_{{\rm{g}} - 1}}\left[ {{p_{{\rm{sat}}, {\rm{KOH}}, {\rm{a}}}} - {p_{{\rm{g}}, {\rm{a}}}}n{{({\rm{ vapor }})}_{\rm{a}}}} \right] \cdot \Delta {H_{{\rm{g}} - 1}} + \\ \left[ {{\eta _{\rm{a}}} + {T_{\rm{a}}}{{\left( {\frac{{{\rm{d}}{E_{{\rm{eq}}}}}}{{{\rm{d}}T}}} \right)}_{\rm{a}}}} \right]{i_{{\rm{loc}}, {\rm{a}}}} \end{array}$

(53)

式(52)和式(53)中：ΔH_g-l为水的气液相变焓。

本模型针对串联多池中的典型单池进行模拟。在AGGL竖直外壁处，相邻单池液态供水会以对流和热传导方式，经金属外框(传热阻力可忽略)与系统交换热量。SFAWE系统内部产生的焦耳热、电极反应热和相变热，一部分通过对流被出口气体带出，另一部分则通过电极催化层、气体间隙层和疏水多孔膜内的对流与热传导，传递至FWC-HPM边界与AGGL竖直外壁处被液态供水交换。供水流速较快时主要以对流方式在FWC-HPM边界与AGGL竖直外壁处与系统交换热量。因此FWC-HPM边界处的边界热源大小为：

$\begin{array}{c} {q_{{\rm{wm}}}} = - {k_{{\rm{g}} - 1}}\left[ {{p_{{\rm{sat}}}} - {p_{{\rm{g}}, {\rm{wm}}}}n{{({\rm{ vapor }})}_{{\rm{wm}}}}} \right] \cdot \Delta {H_{{\rm{g}} - 1}} + \\ {h_{{\rm{wm}}}} \cdot \left( {{T_{{\rm{water}}}} - {T_{{\rm{wm}}}}} \right) \end{array}$

(54)

AGGL竖直外壁处的边界热源大小为：

${q_{{\rm{AGGL}}, {\rm{ wall }}}} = {h_{{\rm{AGGL}}, {\rm{ wall }}}} \cdot \left( {{T_{{\rm{water }}}} - {T_{{\rm{AGGL}}, {\rm{ wall }}}}} \right)$

(55)

式(54)和式(55)中：h_wm与h_{AGGL, wall}分别为FWC-HPM边界与AGGL外壁处的对流传热系数。由于假定液态供水流速足够快，即对流传热系数h_wm与h_{AGGL, wall}足够大，FWC-HPM边界、AGGL竖直外壁处温度均等于供水温度T_water。除AGGL竖直外壁的其余壁面均设为热绝缘，热通量为0。

1.3 模型参数与求解

模型通过COMSOL Multiphysics 5.6软件求解。采用细化自由三角形网格以确保精度。求解器相对误差设为10^-5。求解分3步进行，首先初始化电流分布，然后求解气相质量与动量传递耦合解，最后以这些计算结果为初值，对所有传递过程进行耦合求解。各步骤均使用全耦合求解器。

各计算区域的材料与几何结构参数列于表 1。

模型相关的物性参数如下：对于气相，密度ρ_g按理想气体计算。黏度μ_g、热导率k_g和质量热容C_{p, g}以及二元扩散系数D_ik均通过COMSOL中的热力学系统模块获取。μ_g选用Wilke公式计算，k_g采用动力学理论计算，D_ik采用Fuller公式计算。气相所在多孔介质的曲折因子关联系数β_g取值为1.5。

KOH溶液的密度ρ_l和本征电导率σ_l均采用Gilliam等^[19]的关联式计算。KOH的平均活度系数γ_±采用Akerlof等^[20]的关联式计算。K⁺、OH^－的离子迁移数t₊⁰、t_－⁰均由离子极限摩尔电导率计算：

$t_ + ^0 = \frac{{\lambda _ + ^0}}{{\lambda _ + ^0 + \lambda _ - ^0}}, t_ - ^0 = \frac{{\lambda _ - ^0}}{{\lambda _ + ^0 + \lambda _ - ^0}}$

(56)

极限摩尔电导率λ₊⁰、λ_－⁰采用Robinson等^[21]的关联式计算。KOH整体扩散系数D_KOH采用De Vidts等^[22]的关联式计算。溶液热导率k_l和质量热容C_{p, l}采用Zaytsev等^[19]的关联式计算。纯水和KOH溶液的饱和蒸汽压p_sat和p_{sat, KOH}采用Balej等^[23]的关联式计算。气液相变焓通过COMSOL材料数据库获取。电解质层曲折因子关联系数β_l取值为0.5。

对于固相，PVDF膜的热导率取值为0.262 2 W·m^-1·K^-1，PPS膜的热导率取值为0.28 W·m^-1·K^-1，镍的热导率由COMSOL材料数据库获取。

电极反应动力学参数列于表 2。

模型求解的基准条件列于表 3。

2 结果与讨论 2.1 模拟结果的验证与分析

模拟结果与实验数据^[25]的对比如图 2所示。除供水温度外，实验数据的测定条件均与表 1~表 3中所给出的参数相同。可见模拟结果与实验数据较为吻合，说明模型具有较好的准确性。

各部分过电位在总电压中的占比如图 3所示。由于阳极反应活性较低，本征活化过电位η_{act, 0}在总电压中占比很高。活化过电位中与电解质层水含量相关联的部分Δη_act仅次于η_{act, 0}，意味着水的传递对电解性能具有重要影响。KOH溶液的高电导率使欧姆过电位η_ohm较低。浓度极化过电位η_conc占比最低，这是由于阴极和阳极反应中OH^－均为还原态组分，根据式(36)、(37)和(41)，阴极和阳极侧的浓度极化对过电位的影响相互抵消，因而总的浓度极化过电位很小。

如图 4所示，液态供水中水的化学势高于KOH溶液中的水，水在化学势差的驱动下供入电解质。反应生成的H₂和O₂会夹带部分水蒸气，但阴极侧H₂夹带出的水量(1.73×10^-4 g·s^-1)和阳极侧O₂夹带出的水量(8.28×10^-5 g·s^-1)与供入电解质层的水量(1.48×10^-3 g·s^-1)相比很少，意味着供水效率较高。阳极侧气相主要以对流的方式从系统内流出，因此其水蒸气浓度基本无变化，由阳极侧KOH溶液饱和蒸汽压决定。

如图 5所示，气体间隙层压降较小，电解质层两侧气相压差不大。由于疏水多孔膜的渗透率远小于气体间隙层与多孔电极(见表 1)，因此阴极侧气相压降主要集中于疏水多孔膜。

SFAWE系统的温度分布如图 6所示。水电解反应吸收的总热量来自可逆反应热之和(阴极4.06 W，阳极-8.15 W)。但由于过电位的存在产生大量不可逆反应热(阴极1.20 W，阳极6.72 W)，因此水电解反应总体释放热量(3.83 W)。阳极侧不可逆反应热与可逆反应热相互中和，故反应热主要集中于阴极侧。CCL-EL边界处水蒸气冷凝放热(3.56 W)，ACL-EL边界处O₂带走的水量较少，蒸发吸热量很小(-0.20 W)。此外，电解质层还会产生焦耳热(0.57 W)，将所有热量进行汇总，可知SFAWE系统总体上释放热量(7.76 W)。这些热量大多通过热传导传递至FWC-HPM边界，被液态供水以蒸发吸热(-3.98 W)和对流传热(-0.88 W)的方式带走，还有一部分传导至AGGL竖直外壁后被相邻单池的液态供水交换(-2.37 W)，但由于电解质层热阻较大，这部分热量相对于传递至FWC-HPM边界的热量较少。气相以对流的方式从气体出口带出的热量则相对较少(-0.53 W)。

2.2 适宜运行工况的判定

供水不足会导致电解质层缺水，可能出现电极与电解质局部接触不良、气体穿透电解质层发生交叉返混等问题，使系统性能降低甚至失稳。因此操作电流不宜过高，从系统实际运行的角度考虑，设定电解质层KOH溶液体积V_l减小至电解质层多孔隔膜孔隙体积的70%时对应的电流值为上限。

供水过多会使气体间隙层水淹，过量的水携带KOH进入气体间隙，增大气相传质阻力的同时，还会使电解质浓度降低，导致系统性能下降。若过量的碱液同疏水多孔膜接触，还可能使膜被腐蚀，这些因素都可能导致系统失稳。因此操作电流不宜过低，从系统实际运行的角度考虑，设定电解质层KOH溶液体积V_l增加至多孔电极与电解质层多孔隔膜的总孔隙体积时对应的电流值为下限。

2.3 不同供水温度下，电解质层KOH含量对电解性能的影响

提高电流密度会增加耗水量，使平衡时电解质层水含量降低。根据电解质层液相体积变化可确定系统适宜运行工况范围，如图 7所示。可见适当升高供水温度或提高电解质层KOH含量均利于提高电解性能，扩大适宜运行工况范围，提高系统操作弹性。但为确保系统稳定运行，电流密度与供水温度和电解质层KOH含量应尽量匹配。

电解质层KOH含量变化对电解性能的影响可通过平衡时电解质层KOH浓度与水含量变化进行分析。如图 8所示，一定供水温度与电流密度下，电解质层KOH含量变化几乎不会影响平衡时阴极侧KOH浓度w(KOH)_c。电流密度的大小决定了达到平衡所需的供水量，进而决定平衡时供水推动力与阻力的匹配关系。供水温度一定时，若忽略电解质层KOH含量变化对系统温度分布的影响，则供水推动力大小由阴极侧KOH浓度决定。多孔介质参数给定时，供水阻力的大小由阴极侧气相组成决定。而阴极侧的气相组成受供水温度、阴极侧KOH浓度以及氢气产量共同影响，氢气产量又由电流密度确定下来，因此一定供水温度与电流密度下，阴极侧气相组成与供水阻力也由阴极侧KOH浓度决定。故其他条件一定时，供水温度与电流密度决定了平衡时的供水推动力与阻力，进而决定了与之匹配的阴极侧KOH浓度。

另一方面，电解质层KOH含量的提高会增大KOH溶液的平均浓度w(KOH)_avg，使其中水的化学势降低，供水推动力变大，供水量增加。但由于一定电流密度下所需的水量已经确定，因此电解质层的水含量便会提高，以抵消电解质层KOH含量变化对供水量的影响，从而电解质层液相体积分数ε_l增大。于是KOH溶液的有效电导率修正系数项增大，有效电导率提高，一定电流密度下，为满足传导电流的需求，KOH溶液的浓度梯度便会减小，而平衡时阴极侧KOH浓度几乎不发生变化，因此阳极侧KOH溶液浓度便会提高。

电解质层KOH浓度与水含量的变化使过电位发生改变，进而影响了电解性能，如图 9所示。KOH溶液的有效电导率增大使欧姆过电位η_ohm降低。水含量的提高使有效交换电流密度的修正系数项增大，活化过电位中与水含量相关联的部分Δη_act变小。由于过电位降低，产热变少，平衡时系统整体温度会下降，为保持供水推动力，阴极侧KOH溶液浓度便会小幅下降，使阴极侧KOH溶液的饱和蒸汽压保持不变，这也是图 8中w(KOH)_c会有极小幅度降低的原因。

2.4 不同供水温度下，阴极气体间隙层厚度对电解性能的影响

如图 10所示，适当地减小阴极气体间隙厚度有利于提高电解性能，扩大适宜运行工况范围，但减小阴极气体间隙层厚度会使适宜运行工况对应电流密度区间上下限数值同时增大。为确保系统稳定运行，电流密度与供水温度和阴极气体间隙层厚度应尽量匹配。阴极气体间隙层厚度变化对电解性能的影响亦可通过平衡时电解质层KOH浓度与水含量的变化来分析。在适宜运行工况范围内，降低阴极气体间隙层厚度使电解性能提高主要原因在于水蒸气向阴极传递距离减小，供水阻力降低，平衡时电解质层水含量更高，如图 11(a)所示。由于水含量提高，液相体积分数ε_l增大，有效交换电流密度修正系数项增大，因此活化过电位中与水含量相关联部分Δη_act减小，如图 11(b)所示。虽然电解质层水含量提高使平衡时电解质层KOH溶液的平均浓度降低，但KOH溶液有效电导率修正系数会随电解质层液相体积分数提高而增大，二者总效应使KOH溶液有效电导率先增大后降低，因此欧姆过电位η_ohm呈现先降低后增大的趋势。但由于欧姆过电位变化相较于Δη_act较为平缓，因此在适宜的运行工况范围内，总电压随阴极气体间隙层厚度的降低而下降。

3 结论

利用COMSOL软件首次建立了以KOH溶液为电解质的静态供水式碱性水电解二维数值模型，结合质量、动量、能量和电荷传递以及电化学反应，对系统内的膜蒸馏过程和水电解反应间的复杂耦合关系进行了完整描述，并结合极化曲线实测数据验证了模型可靠性。水的传递对静态供水式水电解性能具有重要影响。供水推动力和阻力变化会显著地改变电解质层水含量，在适宜运行工况范围内，更加充足的水供应有利于改善电解性能。供水温度适当升高、电解质层KOH质量分数略高于30%的工业标准以及水蒸气向阴极传递距离的降低均有利于降低总电压，增大系统操作弹性。但为防止水的供需失衡使系统缺水或水淹而失稳，应尽可能地使供水温度、电流密度、电解质层KOH含量和阴极气体间隙层厚度相互匹配，确保系统始终运行于适宜工况条件下。