质子交换膜水电解中气-液两相传输的模拟研究*

李 港，韩 波

质子交换膜水电解中气-液两相传输的模拟研究*

李 港，韩 波†

（浙江大学 航空航天学院，杭州 310027）

两相传输是影响质子交换膜水电解系统性能的关键因素。为掌握质子交换膜水电解单元中多孔扩散层内气液两相传输规律，基于数值重构的三维多孔扩散层结构，采用格子Boltzmann两相流动模型模拟研究了扩散层内两相传输过程，详细分析了扩散层孔隙率和表面接触角对气泡传输与分布的影响。数值模拟结果表明：孔隙率减小会明显降低气体渗透率，从而导致气泡难以在扩散层内找到有效传输通道。接触角的增大不仅增加了气泡在界面堆积的风险，也减缓了气泡在孔隙内的传输速度。从孔隙尺度水平初步掌握了质子交换膜水电解单元多孔扩散层内两相传输规律，可为高性能水电解系统设计和优化提供理论支撑。

质子交换膜水电解；多孔扩散层；两相传输；格子Boltzmann方法；多相流模型

0 引 言

高性能水电解技术可以实现电能与氢能之间的高效转换，为可再生能源存储提供了一种既高效又可持续的解决方案。常见的水电解制氢技术有碱性水电解、质子交换膜水电解和固体氧化物水电解三种。其中，质子交换膜水电解技术由于具有高效、气体制取纯度高、结构紧凑、工作电流密度高等优点[1-2]，被公认为是制氢和可再生能源存储领域极具发展前景的水电解技术。

气液两相传输过程是影响质子交换膜水电解系统性能的关键因素[3-4]。目前国内外针对该领域的研究主要集中在初步的实验观测和宏观建模方面。如LIM等[5]研究了碳纸扩散层中疏水聚合物的含量对质子交换膜水电解系统动力性能的影响。接触角测量结果表明，浸渍不同量氟化乙丙烯（fluorinated ethylene propylene, FEP）的多孔扩散层具有相似的疏水性。GRIGORIEV等[6]研究了烧结钛颗粒多孔扩散层，发现颗粒的平均孔径和钛板厚度对质子交换膜水电解性能有显著影响。SELAMET等[7]研究发现质子交换膜水电解中气泡的演化遵循一个周期性生长和脱落的规律；同时还发现水流速度、电流密度、温度对两相传输过程有显著影响。LAFMEJANI等[8]研究了单直流道内垂直向上的气液流动过程。结果表明，沿着流道的长度方向发生了气泡的聚合现象。另外，ZLOBINSKI等[9]使用中子成像技术观察了阳极扩散层内不同电流密度下气液分布情况，发现扩散层内两相流动都是由毛细管驱动的，而压力和电流密度对扩散层内水、气分布无显著影响。

在建模方面，一些宏观模型初步研究了质子交换膜水电解单元中的质量输运问题[10-11]。另外，AUBRAS等[12]设计了质子交换膜水电解电极的二维模型，用以研究电解单元内部传输过程。研究结果表明，气泡的聚合现象及传输性能与欧姆阻力密切相关。ARBABI等[13]使用流体体积（volume of fluid, VOF）方法模拟探究了多孔扩散层内气泡的传输过程。LAFMEJANI等[14]使用VOF方法预测了微流道内的泰勒流，并通过实验对其模拟结果进行了验证。

综上，目前针对质子交换膜水电解单元内部传输现象的研究还处于初步观察和预测阶段，而对复杂多孔扩散层内孔隙尺度下的气液两相传输规律还缺乏基本的认识。现有的实验手段还无法完全探测到扩散层尤其是催化层内部传输现象，而传统的宏观计算模型又很难有效模拟微介观尺度下的气液两相动态传输过程。本文采用基于介观理论的格子Boltzmann两相流动模型，数值模拟研究质子交换膜水电解单元扩散层内气液两相传输现象，详细分析不同参数对三维扩散层内气泡传输与分布的影响。

1 计算模型

1.1 格子Boltzmann两相模型

格子Boltzmann方法已被成功用于处理诸多两相传输问题[15-18]。本文采用Shan-Chen格子Boltzmann多组分多相流模型对质子交换膜水电解单元多孔扩散层内气液两相传输过程进行模拟研究。该模型的演化方程可以表示为：

式中：fαk(x,t)为t时刻位置x处组分k的分布函数，k = 1, 2 ... S表示组分，α = 1, 2 ... q表示离散速度方向；eα为离散速度；δt 为时间步长；τk为第k组分的松弛时间，与第k组分流体的黏度有关。fαk(eq)(x,t)为相应的平衡态分布函数，其离散速度模型采用QIAN等[19]提出的D3Q19模型，如图1所示。

D3Q19模型的离散速度矢量为：

第组分的平衡态分布函数为：

假设不同流体粒子间存在非局部的相互作用，可以采用势函数来表征，其形式为：

根据相互作用势，可得第组分粒子与其他组分粒子间的相互作用力，其离散形式的表达式为：

在流固界面，考虑到壁面的润湿性，流体与固壁之间相互作用力的计算式可表示为：

式中：w表征了流体相与固体壁面间相互作用的强度。用()来标识固体格点，如果处是固体格点，则() = 0，否则() = 1。

重力的计算表达式为：

式中：为重力加速度。

作用力的影响通过下式引入，即：

而混合流体的宏观速度为：

为了满足系统的总动量守恒，需要对平衡态分布函数中的速度进行重新定义：

式中：为修正后的混合流体宏观速度，其表达式为：

混合流体的宏观密度和压力可以表示为：

1.2 两相模型验证

1.2.1 Laplace定律验证

以下验证采用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定律完全吻合。

1.2.2 接触角验证

如图3，计算区域网格设置为50 lu × 50 lu × 50 lu，初始在底部表面中心放置一个固定大小的气泡，通过调节强度系数进行数值模拟，当达到平衡状态时，气泡在固体表面上会形成特定的接触角。图中依次为w= −0.5、0、0.5时表面接触角的模拟结果，分别对应120°、90°、60°。因此目前的格子Boltzmann两相模型可以实现不同接触角的模拟。图4给出了流固间相互作用强度w与接触角之间的关系曲线。可以看到，接触角与强度w基本呈线性关系。

图3 接触角验证

图4 接触角θ与Gw的关系

2 结果与讨论

2.1 多孔扩散层孔隙率对气泡传输的影响

基于以上验证的格子Boltzmann两相流动模型，数值模拟研究了质子交换膜水电解单元扩散层内气泡传输现象，其中多孔扩散层采用非交织纤维随机生成方法[20-21]进行重构。

首先研究扩散层孔隙率对气泡传输与分布的影响。采用100 lu × 100 lu × 80 lu网格，在= 0处施加速度边界和反弹格式，其余边界条件设为非平衡态外推边界，扩散层厚度为30 lu，接触角为30°。图5给出了两种孔隙率扩散层中气泡传输过程的模拟结果对比。如图所示，孔隙率的变化会明显影响多孔扩散层内气泡传输与分布情况，孔隙率取0.78时，随着时间的推移，气泡能够顺利通过扩散层，以自由气泡的形态脱落传递；但是当孔隙率为0.56时，气泡在扩散层内的传输会受到明显的阻力，难以通过其多孔结构。这是由于孔隙率减小会缩小平均孔径，导致气体渗透率降低；另外，孔隙率减小会增加孔隙交叉复杂度，导致气泡难以找到有效通道。

图5 不同孔隙率的重构多孔扩散层中气泡传输与分布：（a、b）T = 500δt；（c、d）T = 1.5 × 103δt

进一步定量分析扩散层内气泡分布情况。图6给出了扩散层内不同截面上平均气体饱和度的变化情况。可以看到，当孔隙率为0.78时，在= 500δ时，位于= 10 lu、15 lu和20 lu截面处的气体饱和度分别为0.1804、0.1420和0.1395；在= 1.5 × 103δ时，对应截面处的气体饱和度变为0.180 3、0.214 7和0.259 9。而且在扩散层中上部乃至外流道内（> 20 lu）均有气体存在，如在= 1.5 × 103δ时，= 25 lu、35 lu、45 lu、55 lu截面处对应的气体饱和度分别为0.144 5、0.064 7、0.052 9和0.013 4。而当孔隙率为0.56时，气体则主要分布于扩散层的下半部分（= 10 ～ 20 lu），其相应的气体饱和度随着时间没有发生明显的变化。

2.2 多孔扩散层接触角对气泡传输的影响

图7给出了两种接触角条件下多孔扩散层内气泡传输过程的模拟结果，其中，扩散层厚度为30 lu，孔隙率为0.78。可以看到，接触角变化也会对扩散层内气泡传输产生影响。当接触角= 30°时，气泡受到的传输阻力相对较小，可以顺利通过扩散层。然而当接触角= 70°时，可以看到扩散层内部出现了气泡成团堆积的现象。这是由于接触角较小时，气泡与固体表面的相互作用力相对较弱，气泡较容易脱离表面，快速通过有效孔隙；反之，气泡会受到固体表面较大的阻力作用，很难在孔隙中进行传递，导致气泡堆积现象。

同样，可以定量分析接触角对扩散层内气泡分布的影响。图8给出了两种接触角条件下扩散层不同截面处气体饱和度的变化曲线。可以看到，接触角变化对不同区域（= 0 ～ 5 lu、= 5 ～ 35 lu、> 35 lu）气体饱和度大小均有影响。当接触角= 70°，在= 0 ～ 5 lu和= 5 ～ 15 lu区域的气体饱和度都比接触角= 30°时更高。如= 70°时，在= 500δ，= 5 lu和10 lu截面处的气体饱和度分别为0.324 8和0.320 1；= 30°时对应截面的气体饱和度分别为0.183 1和0.180 4。从扩散层中部开始（≥ 15 lu），= 30°时的气体饱和度均高于= 70°时的对应数值。这是由于在接触角较大时，气泡与固体表面的相互作用力较强，气泡难以摆脱表面形成自由气泡，导致其在该区域发生堆积，难以进入扩散层内部；而接触角较小时，气泡可以顺利进入扩散层，因此在这个区域（≥ 15lu），小接触角扩散层内气体饱和度均高于大接触角扩散层内的气体饱和度。

图7 不同接触角的重构多孔扩散层中气泡传输与分布：（a、b）T = 500δt；（c、d）T = 1.5 × 103δt

3 结 论

对质子交换膜水电解单元扩散层结构进行数值重构，采用格子Boltzmann两相模型对其内部气泡传输现象进行了研究，结果表明：扩散层孔隙率会明显影响内部气泡传输过程，降低孔隙率会增加孔隙交叉复杂度，缩减孔隙尺度，导致气体渗透率降低，气泡难以顺利通过扩散层；接触角增大，会增加气泡的传输阻力，多个气泡在邻近孔隙不断合并，导致气泡发生堆积现象。

本研究从孔隙尺度水平初步掌握了质子交换膜水电解单元多孔扩散层内气泡传输规律。后续可以考虑催化层微观结构的影响，深入分析电解系统微介观两相传输规律。

[1] WANG X Y, ZHANG L S, LI G F, et al. The influence of Ferric ion contamination on the solid polymer electrolyte water electrolysis performance[J]. Electrochimica acta, 2015, 158: 253-257. DOI: 10.1016/j.electacta.2015.01.140.

[2] ARICÒ A S, SIRACUSANO S, BRIGUGLIO N, et al. Polymer electrolyte membrane water electrolysis: status of technologies and potential applications in combination with renewable power sources[J]. Journal of applied electrochemistry, 2013, 43(2): 107-118. DOI: 10.1007/s10800-012-0490-5.

[3] CARMO M, FRITZ D L, MERGEL J, et al. A comprehensive review on PEM water electrolysis[J]. International journal of hydrogen energy, 2013, 38(12): 4901-4934. DOI: 10.1016/j.ijhydene.2013.01.151.

[4] ZHANG X, YANG Y P, GUO L J, et al. Effects of carbon corrosion on mass transfer losses in proton exchange membrane fuel cells[J]. International journal of hydrogen energy, 2017, 42(7): 4699-4705. DOI: 10.1016/j.ijhydene.2016.08.223.

[5] LIM C, WANG C Y. Effects of hydrophobic polymer content in GDL on power performance of a PEM fuel cell[J]. Electrochimica acta, 2004, 49(24): 4149-4156. DOI: 10.1016/j.electacta.2004.04.009.

[6] GRIGORIEV S A, MILLET P, VOLOBUEV S A, et al. Optimization of porous current collectors for PEM water electrolysers[J]. International journal of hydrogen energy, 2009, 34(11): 4968-4973. DOI: 10.1016/j.ijhydene.2008.11.056.

[7] SELAMET O F, PASAOGULLARI U, SPERNJAK D, et al. Two-phase flow in a proton exchange membrane electrolyzer visualizedby simultaneous neutron radiography and optical imaging[J]. International journal of hydrogen energy, 2013, 38(14): 5823-5835. DOI: 10.1016/j.ijhydene.2013.02.087.

[8] LAFMEJANI S S, OLESEN A C, KÆR S K. Analysing gas-liquid flow in PEM electrolyser micro-channels[J]. ECS transactions, 2016, 75(14): 1121-1127. DOI: 10.1149/07514.1121ecst.

[9] ZLOBINSKI M, SCHULER T, BÜCHI F N, et al. Transient and steady state two-phase flow in anodic porous transport layer of proton exchange membrane waterelectrolyzer[J]. Journal of the electrochemical society, 2020, 167(8): 084509. DOI: 10.1149/1945-7111/ab8c89.

[10] HAN B, MO J K, KANG Z Y, et al. Effects of membrane electrode assembly properties on two-phase transport and performance in proton exchange membrane electrolyzer cells[J]. Electrochimica acta, 2016, 188: 317-326. DOI: 10.1016/j.electacta.2015.11.139.

[11] HAN B, MO J K, KANG Z Y, et al. Modeling of two-phase transport in proton exchange membrane electrolyzer cells for hydrogen energy[J]. International journal of hydrogen energy, 2017, 42(7): 4478-4489. DOI: 10.1016/j.ijhydene.2016.12.103.

[12] AUBRAS F, DESEURE J, KADJO J J A, et al. Two-dimensional model of low-pressure PEM electrolyser: two-phase flow regime, electrochemical modelling and experimental validation[J]. International journal of hydrogen energy, 2017, 42(42): 26203-26216. DOI: 10.1016/j.ijhydene.2017.08.211.

[13] ARBABI F, MONTAZERI H, ABOUATALLAH R, et al. Three-dimensional computational fluid dynamics modelling of oxygen bubble transport in polymer electrolyte membrane electrolyzer porous transport layers[J]. Journal of the electrochemical society, 2016, 163(11): F3062-F3069. DOI: 10.1149/2.0091611jes.

[14] LAFMEJANI S S, OLESEN A C, KÆR S K. VOF modelling of gas–liquid flow in PEM water electrolysis cell micro-channels[J]. International journal of hydrogen energy, 2017, 42(26): 16333-16344. DOI: 10.1016/j.ijhydene.2017.05.079.

[15] HAN B, MENG H. Numerical studies of interfacial phenomena in liquid water transport in polymer electrolytemembrane fuel cells using the lattice Boltzmann method[J]. International journal of hydrogen energy, 2013, 38(12): 5053-5059. DOI: 10.1016/j.ijhydene.2013.02.055.

[16] HAO L, CHENG P. Lattice Boltzmann simulations of water transport in gas diffusion layer of a polymer electrolyte membrane fuel cell[J]. Journal of power sources, 2010, 195(12): 3870-3881. DOI: 10.1016/j.jpowsour.2009.11.125.

[17] PARK J, MATSUBARA M, LI X. Application of lattice Boltzmann method to a micro-scale flow simulation in the porous electrode of a PEM fuel cell [J]. Journal of Power Sources, 2007, 173(1): 404-414. DOI: 10.1016/j. jpowsour.2007.04.021.

[18] HAO L, CHENG P. Pore-scale simulations on relative permeabilities of porous media by lattice Boltzmann method[J]. International journal of heat and mass transfer, 2010, 53(9/10): 1908-1913. DOI: 10.1016/j.ijheatmasstransfer.2009.12.066.

[19] QIAN Y H, D'HUMIÈRES D, LALLEMAND P. Lattice BGK models for navier-stokes equation[J]. Europhysics letters, 1992, 17(6): 479-484. DOI: 10.1209/0295-5075/17/6/001.

[20] SCHULZ V P, MUKHERJEE P P, BECKER J, et al. Numerical evaluation of effective gas diffusivity - Saturation dependence of uncompressed and compressed gas diffusion media in PEFCs[J]. ECS transactions, 2006, 3(1): 1069-1075. DOI: 10.1149/1.2356226.

[21] SCHULZ V P, BECKER J, WIEGMANN A, et al. Modeling of two-phase behavior in the gas diffusion medium of PEFCs via full morphology approach[J]. Journal of the electrochemical society, 2007, 154(4): B419-B426. DOI: 10.1149/1.2472547.

Simulation of Gas/Liquid Two-Phase Transport in Proton Exchange Membrane Water Electrolysis

LI Gang, HAN Bo

(School of Aeronautics and Astronautics, Zhejiang University, Hangzhou 310027, China)

Two-phase transport is the key factor that affects the performance of the proton exchange membrane (PEM) water electrolysis system. The lattice Boltzmann method (LBM) was employed to simulate the two-phase transport process in a three-dimensional reconstructed porous diffusion layer for the mechanism of gas/liquid two-phase transport inside the porous diffusion layer of a PEM water electrolysis cell. The relevant simulations were conducted to study the effects of the porosity and surface contact angle of the diffusion layer on the gas bubble transport and distribution. Numerical results indicated that the gas permeability of the diffusion layer significantly reduced as the porosity decreased, which made it difficult for gas bubbles to get the effective transport pathway inside the diffusion layer. Additionally, as the contact angle increased, the risk of gas bubble accumulation on the interface increased, and the velocity of gas bubble transport in the pores also slowed down. A fundamental understanding of the pore-scale two-phase transport mechanism in the porous diffusion layer of the PEM water electrolysis cell was obtained, which might provide theoretical support for the design and optimization of a high-performance water electrolysis system.

proton exchange membrane water electrolysis; porous diffusion layer; two-phase transport; lattice Boltzmann method; multiphase flow model

2095-560X（2022）05-0393-07

TK91；O359+.1

A

10.3969/j.issn.2095-560X.2022.05.001

2022-04-22

2022-07-20

国家自然科学基金项目（11702245）

韩 波，E-mail：bohan@zju.edu.cn

李 港（1997-），男，硕士研究生，主要从事微介观气液两相传输过程的数值模拟研究。

韩 波（1984-），男，博士，副教授，主要从事气液两相流动机理方面的研究。