0 引言
气液两相传输过程是影响质子交换膜水电解系统性能的关键因素[3,4]。目前国内外针对该领域的研究主要集中在初步的实验观测和宏观建模方面。如LIM等[5]研究了碳纸扩散层中疏水聚合物的含量对质子交换膜水电解系统动力性能的影响。接触角测量结果表明,浸渍不同量氟化乙丙烯（fluorinated ethylene propylene, FEP）的多孔扩散层具有相似的疏水性。GRIGORIEV等[6]研究了烧结钛颗粒多孔扩散层,发现颗粒的平均孔径和钛板厚度对质子交换膜水电解性能有显著影响。SELAMET等[7]研究发现质子交换膜水电解中气泡的演化遵循一个周期性生长和脱落的规律;同时还发现水流速度、电流密度、温度对两相传输过程有显著影响。LAFMEJANI等[8]研究了单直流道内垂直向上的气液流动过程。结果表明,沿着流道的长度方向发生了气泡的聚合现象。另外,ZLOBINSKI等[9]使用中子成像技术观察了阳极扩散层内不同电流密度下气液分布情况,发现扩散层内两相流动都是由毛细管驱动的,而压力和电流密度对扩散层内水、气分布无显著影响。
综上,目前针对质子交换膜水电解单元内部传输现象的研究还处于初步观察和预测阶段,而对复杂多孔扩散层内孔隙尺度下的气液两相传输规律还缺乏基本的认识。现有的实验手段还无法完全探测到扩散层尤其是催化层内部传输现象,而传统的宏观计算模型又很难有效模拟微介观尺度下的气液两相动态传输过程。本文采用基于介观理论的格子Boltzmann两相流动模型,数值模拟研究质子交换膜水电解单元扩散层内气液两相传输现象,详细分析不同参数对三维扩散层内气泡传输与分布的影响。
1 计算模型
1.1 格子Boltzmann两相模型
$\begin{matrix} f_{\alpha }^{k}\left( \mathbf{x}+{{\mathbf{e}}_{\alpha }}{{\delta }_{t}},t+{{\delta }_{t}} \right)-f_{\alpha }^{k}\left( \mathbf{x},t \right)= \\ -\frac{1}{{{\tau }_{k}}}\left[ f_{\alpha }^{k}\left( \mathbf{x},t \right)-f_{\alpha }^{k(eq)}\left( \mathbf{x},t \right) \right] \\ \end{matrix}$ (1)
Fig. 1 D3Q19 discrete velocity model图1 D3Q19离散速度模型 |
D3Q19模型的离散速度矢量E为：
$E=\left[ \begin{matrix} 0 & 1 & -1 & 0 & 0 & 0 & 0 & 1 & -1 & 1 & -1 & 1 & -1 & 1 & -1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 & -1 & 0 & 0 & 1 & -1 & -1 & 1 & 0 & 0 & 0 & 0 & 1 & -1 & 1 & -1 \\ 0 & 0 & 0 & 0 & 0 & 1 & -1 & 0 & 0 & 0 & 0 & 1 & -1 & -1 & 1 & 1 & -1 & -1 & 1 \\\end{matrix} \right]$ (2)
第k组分的平衡态分布函数为：
$f_{\alpha }^{k(eq)}={{w}_{\alpha }}{{\rho }_{k}}\left[ 1+\frac{{{\mathbf{e}}_{\alpha }}\cdot {{u}_{k}}}{{{c}_{s}}^{2}}+\frac{{{\left( {{\mathbf{e}}_{\alpha }}\cdot {{u}_{k}} \right)}^{2}}}{2{{c}_{s}}^{4}}-\frac{{{u}_{k}}^{2}}{2{{c}_{s}}^{2}} \right]$ (3)
式中：${{w}_{\alpha }}$为权系数,取值为1/3（α = 0）、1/18（α = 1 ~ 6）、1/36（α = 7~18）;cs为格子声速[19],取值为1/$\sqrt{3}$;而模型中第k组分的密度和速度表达式如下：
${{\rho }_{k}}=\sum{f_{\alpha }^{k}}$ (4)
${{u}_{k}}=\frac{1}{{{\rho }_{k}}}\sum{{{\mathbf{e}}_{\alpha }}f_{\alpha }^{k}}$ (5)
假设不同流体粒子间存在非局部的相互作用,可以采用势函数来表征,其形式为：
${{V}_{k{k}'}}\left( x,{x}' \right)={{G}_{k{k}'}}{{\varphi }_{k}}\left( x \right){{\varphi }_{{{k}'}}}\left( {{x}'} \right)\left( \left| x- \right.\left. {{x}'} \right| \right)$ (6)
式中：$x\prime =x+{{\mathbf{e}}_{\alpha }}{{\delta }_{t}}$,表示格子位置x的临近格点位置;φ(x)为第k组分粒子的伪势函数;${{G}_{k{k}'}}\left( x,{x}' \right)$为格林函数,表示流体粒子之间相互作用的强度,具体表达式如下：
${{G}_{k{k}'}}\left( \left| x- \right.\left. {{x}'} \right| \right)=\left\{ \begin{matrix} G\text{ }\left| x- \right.\left. {{x}'} \right|={{\delta }_{x}} \\ 0\text{ }\left| x- \right.\left. {{x}'} \right|>{{\delta }_{x}} \\\end{matrix} \right\}$ (7)
根据相互作用势,可得第k组分粒子与其他组分粒子间的相互作用力,其离散形式的表达式为：
${{F}_{f}}\left( x \right)=-{{G}_{k{k}'}}\left( x,{x}' \right){{\varphi }_{k}}\left( x \right)\sum\limits_{\alpha }{{{w}_{\alpha }}}{{\varphi }_{{{k}'}}}\left( x+{{e}_{\alpha }}{{\delta }_{t}} \right){{e}_{\alpha }}$ (8)
在流固界面,考虑到壁面的润湿性,流体与固壁之间相互作用力的计算式可表示为：
${{F}_{w}}\left( x \right)=-{{\varphi }_{k}}\left( x \right)\sum\limits_{\alpha }{{{G}_{w}}{{w}_{\alpha }}}s\left( x+{{e}_{\alpha }}{{\delta }_{t}} \right){{e}_{\alpha }}$ (9)
式中：Gw表征了流体相与固体壁面间相互作用的强度。用s(x)来标识固体格点,如果x处是固体格点,则s(x) = 0,否则s(x) = 1。
重力的计算表达式为：
${{F}_{g}}(x)=-{{\rho }_{k}}\cdot g$ (10)
式中：g为重力加速度。
作用力Fk的影响通过下式引入,即：
${{u}_{k}}^{eq}={{u}_{k}}+\frac{{{\tau }_{k}}}{{{\rho }_{k}}}{{F}_{k}}$ (11)
而混合流体的宏观速度u为：
$u=\frac{\sum{{{\rho }_{k}}{{u}_{k}}}}{\rho }$ (12)
为了满足系统的总动量守恒,需要对平衡态分布函数中的速度进行重新定义：
${{u}_{k}}^{eq}={u}'+\frac{{{\tau }_{k}}}{{{\rho }_{k}}}{{F}_{k}}$ (13)
式中：u' 为修正后的混合流体宏观速度,其表达式为：
${u}'=\frac{\sum{\frac{{{\rho }_{k}}{{u}_{k}}}{{{\tau }_{k}}}}}{\sum{\frac{{{\rho }_{k}}}{{{\tau }_{k}}}}}$ (14)
混合流体的宏观密度ρ和压力p可以表示为：
$\rho =\sum{{{\rho }_{k}}}$ (15)
$p(x)=c_{s}^{2}\left[ \sum{{{\rho }_{k}}\left( x \right)+{{G}_{k{k}'}}\left( x,{x}' \right)}{{\varphi }_{k}}\left( x \right){{\varphi }_{{{k}'}}}\left( x \right) \right]$ (16)
1.2 两相模型验证
1.2.1 Laplace定律验证
根据Laplace定律,气泡与周围液体处于平衡时的内外压力差（${{p}_{in}}-{{p}_{out}}$）应与作用于气液界面上的表面张力σ成正比,与气泡的半径R成反比,即满足下式：
${{p}_{in}}-{{p}_{out}}=\frac{\sigma }{R}$ (17)
以下验证采用D3Q19离散速度模型,有效密度为φ(x) = ρ,计算网格为100 lu × 100 lu × 100 lu（lu表示格子单位）,四周采用非平衡态外推边界,1相为气态,2相为液态。初始时刻将 (50, 50, 50) 设置为球心,半径为R的区域内密度为ρ1 = 1 × 10-3,ρ2 = 1.00;区域外密度为ρ1 = 1.00,ρ2 = 1 × 10-3;松弛时间取1.0,G12和G21取值为 -1.8,G11和G22取值为0。图2给出了稳定状态下气泡内外压力差与曲率之间的关系,可以看到格子Boltzmann模拟结果与Laplace定律完全吻合。
Fig. 2 Validations of Laplace’s law图2 Laplace定律验证 |
1.2.2 接触角验证
Fig. 3 Validation of contact angle图3 接触角验证 |
Fig. 4 Relationship between contact angle θ and Gw图4 接触角θ与Gw的关系 |
2 结果与讨论
2.1 多孔扩散层孔隙率对气泡传输的影响
首先研究扩散层孔隙率对气泡传输与分布的影响。采用100 lu × 100 lu × 80 lu网格,在Z = 0处施加速度边界和反弹格式,其余边界条件设为非平衡态外推边界,扩散层厚度为30 lu,接触角为30°。图5给出了两种孔隙率扩散层中气泡传输过程的模拟结果对比。如图所示,孔隙率的变化会明显影响多孔扩散层内气泡传输与分布情况,孔隙率取0.78时,随着时间的推移,气泡能够顺利通过扩散层,以自由气泡的形态脱落传递;但是当孔隙率为0.56时,气泡在扩散层内的传输会受到明显的阻力,难以通过其多孔结构。这是由于孔隙率减小会缩小平均孔径,导致气体渗透率降低;另外,孔隙率减小会增加孔隙交叉复杂度,导致气泡难以找到有效通道。
Fig. 5 Gas bubble transport and distribution in a reconstructed diffusion layer with different porosities: (a, b) T = 500δt; (c, d) T = 1.5 × 103δt图5 不同孔隙率的重构多孔扩散层中气泡传输与分布：（a、b）T = 500δt;（c、d）T = 1.5 × 103δt |
进一步定量分析扩散层内气泡分布情况。图6给出了扩散层内不同截面上平均气体饱和度的变化情况。可以看到,当孔隙率为0.78时,在T = 500δt时,位于Z = 10 lu、15 lu和20 lu截面处的气体饱和度分别为0.180 4、0.142 0和0.139 5;在T = 1.5 × 103δt时,对应截面处的气体饱和度变为0.180 3、0.214 7和0.259 9。而且在扩散层中上部乃至外流道内（Z > 20 lu）均有气体存在,如在T = 1.5 × 103δt时,Z = 25 lu、35 lu、45 lu、55 lu截面处对应的气体饱和度分别为0.144 5、0.064 7、0.052 9和0.013 4。而当孔隙率为0.56时,气体则主要分布于扩散层的下半部分（Z = 10 ~ 20 lu）,其相应的气体饱和度随着时间没有发生明显的变化。
Fig. 6 Gas distribution in each section of diffusion layer with different porosities图6 不同孔隙率的扩散层各截面内气相分布情况 |
2.2 多孔扩散层接触角对气泡传输的影响
图7给出了两种接触角条件下多孔扩散层内气泡传输过程的模拟结果,其中,扩散层厚度为30 lu,孔隙率为0.78。可以看到,接触角变化也会对扩散层内气泡传输产生影响。当接触角θ = 30°时,气泡受到的传输阻力相对较小,可以顺利通过扩散层。然而当接触角θ = 70°时,可以看到扩散层内部出现了气泡成团堆积的现象。这是由于接触角较小时,气泡与固体表面的相互作用力相对较弱,气泡较容易脱离表面,快速通过有效孔隙;反之,气泡会受到固体表面较大的阻力作用,很难在孔隙中进行传递,导致气泡堆积现象。
Fig. 7 Gas bubble transport and distribution in a reconstructed diffusion layer with different contact angles: (a, b) T = 500δt; (c, d) T = 1.5 × 103δt图7 不同接触角的重构多孔扩散层中气泡传输与分布：（a、b）T = 500δt;（c、d）T = 1.5 × 103δt |
同样,可以定量分析接触角对扩散层内气泡分布的影响。图8给出了两种接触角条件下扩散层不同截面处气体饱和度的变化曲线。可以看到,接触角变化对不同区域（Z = 0 ~ 5 lu、Z = 5 ~ 35 lu、Z > 35 lu）气体饱和度大小均有影响。当接触角θ = 70°,在Z = 0 ~ 5 lu和Z = 5 ~ 15 lu区域的气体饱和度都比接触角θ = 30°时更高。如θ = 70°时,在T = 500δt,Z = 5 lu和10 lu截面处的气体饱和度分别为0.324 8和0.320 1;θ = 30°时对应截面的气体饱和度分别为0.183 1和0.180 4。从扩散层中部开始（Z ≥ 15 lu）,θ = 30°时的气体饱和度均高于θ = 70°时的对应数值。这是由于在接触角较大时,气泡与固体表面的相互作用力较强,气泡难以摆脱表面形成自由气泡,导致其在该区域发生堆积,难以进入扩散层内部;而接触角较小时,气泡可以顺利进入扩散层,因此在这个区域（Z ≥ 15lu）,小接触角扩散层内气体饱和度均高于大接触角扩散层内的气体饱和度。
Fig. 8 Gas distribution in each section of diffusion layer with different contact angles图8 不同接触角的扩散层各截面内气相分布情况 |
3 结论
对质子交换膜水电解单元扩散层结构进行数值重构,采用格子Boltzmann两相模型对其内部气泡传输现象进行了研究,结果表明：扩散层孔隙率会明显影响内部气泡传输过程,降低孔隙率会增加孔隙交叉复杂度,缩减孔隙尺度,导致气体渗透率降低,气泡难以顺利通过扩散层;接触角增大,会增加气泡的传输阻力,多个气泡在邻近孔隙不断合并,导致气泡发生堆积现象。
本研究从孔隙尺度水平初步掌握了质子交换膜水电解单元多孔扩散层内气泡传输规律。后续可以考虑催化层微观结构的影响,深入分析电解系统微介观两相传输规律。