水淹被认为是制约质子交换膜燃料电池性能、耐久性以及紧凑性提升的重要问题，在高电流密度工况时水淹的影响更为严重^[1]。水淹问题主要出现在燃料电池气体流道(gas channels，GC)和气体扩散层(gas diffusion layer，GDL)中，通过优化GC和GDL的结构来增强其在高电流密度下的排水能力，从而避免水淹的出现，这对提升燃料电池的可靠性和耐久性有重要意义。

近年来，三维波浪形流道因其在提升除水性能方面的显著能力而被广泛关注。三维波浪形流道由Kuo等^[2-5]首次提出，并且实验、二维以及三维计算流体力学(computational fluid dynamics, CFD)模拟都表明该流道相比传统直流道可以提高气体流速以及改善流道内部的温度分布，从而具有更高的输出电压和功率密度。Kuo等^[6]对波浪形的障碍块的形状进行了优化，提出了圆弧形、梯形和楼梯形3种形状，结果表明：3种类型相比于传统直流道在对流传热和提高流速方面均有益，但是没有对3种类型进行相互对比。Byun等^[7]通过实验详细研究了单直流道中加入半圆形波浪的尺寸以及半圆的间距对燃料电池性能的影响，发现最大功率密度随着波浪凸起高度的增加而增加，随着相邻波浪间距的减小而减小。Li等^[8]设计了一种三维周期波浪形流道，并应用于蛇形流道，通过三维多相稳态的数值模拟发现可以在保持低压降的情况下提高当地流速从而增强除水能力。Kim等^[9]基于模拟分析丰田Mirai在阴极使用了类似于波浪形的三维复杂流道^[10]时发现，波浪形三维流道能在保持低压降的同时增强除水能力的主要原因是高电流密度条件下，波浪形流道引起的俯冲效应使得GDL中流体惯性力增大。而这种惯性力对于除水作用的机理却没有得到很好的解释。

同时，改造GDL的结构从而提升其排水能力也被众多学者研究。Gerteisen等^[11]首次在GDL上通过激光打孔提高排水性能。Wang等^[12]研究孔的直径和孔中心线之间的距离对GDL中水传输的影响。Fang等^[13]通过数值模拟获得GDL输运物性和孔直径的关系。然而目前针对在GDL上打孔的研究多为使用激光得到圆柱形的通孔。相比于圆柱形，圆锥形的孔具有更大侧面积，同时底部不易积水也能更好保护催化剂层。此外，目前研究多集中于流道或扩散层的独立优化，而流道和扩散层中的水输运则是相互紧密联系的。因此针对液态水传输，应将扩散层与波浪形流道进行协同优化，即将两者互为基础进行优化，这样更有利于实现多余液态水的加速排出。

介观尺度的格子Boltzmann方法(lattice Boltzmann method，LBM)近年来逐渐被发展完善，LBM可以用来模拟流体运动，尤其在复杂的多孔介质几何形状和多相流界面问题上有较大优势。LBM中可以用于多相流模拟的伪势模型(也被称为Shan-Chen模型)，由于其能直接模拟流体与流体、流体与壁面之间的作用力，算法简单，并行性能好，可实现相界面分离，因此被广泛应用^[14-15]，但传统的伪势模型难以实现较大密度比，而Li等^[16]通过在传统伪势模型中使用真实状态方程并修正了由此带来的热力学一致性问题，从而实现了高密度比。高密度比伪势模型更接近真实情况，也更能体现惯性力的作用，本文将采用此模型。

本文采用LBM-MRT高密度比伪势模型，模拟了高电流密度工况下质子膜燃料电池扩散层和流道内孔隙尺度下的水输运过程，得到GDL中Re、波浪形流道角度、GDL开孔形状和位置对燃料电池水管理的影响，并尝试得到最优化的GDL和波浪形流道的协同设计。

1 格子Boltzmann模型 1.1 多松弛时间模型

本文采用D2Q9多松弛时间格子Boltzmann模型来模拟气体扩散层和流道中的流动过程。该模型的演化方程表示如下^[17-18]：

$ \begin{array} [c]{c} f_{a}\left(x+e_{a} \delta_{t}, t+\delta_{t}\right)-f_{a}(x, t)=\\-\overline{\boldsymbol{\varLambda}}_{\alpha \beta}\left.\left(f_{\beta}-f_{\beta}^{e q}\right)\right|_{(x, t)}+\delta_{t}\left.\left(\boldsymbol{S}_{a}-0.5 \overline{\boldsymbol{\varLambda}}_{a \beta} \boldsymbol{S}_{\beta}\right)\right|_{(x, t)}, \end{array} $

(1)

$ \boldsymbol{\varLambda}=\operatorname{diag}\left(\tau_{\rho}^{-1}, \tau_{\mathrm{e}}^{-1}, \tau_{\xi}^{-1}, \tau_{\mathrm{j}}^{-1}, \tau_{\mathrm{q}}^{-1}, \tau_{\mathrm{j}}^{-1}, \tau_{\mathrm{q}}^{-1}, \tau_{v}^{-1}, \tau_{v}^{-1}\right). $

(2)

其中：$\overline{\boldsymbol{\varLambda}}=\boldsymbol{M}^{-1} \boldsymbol{\varLambda} \boldsymbol{M}$是碰撞矩阵，Λ是松弛时间矩阵。

式(2)中松弛时间表示如下^[19]：

$ \begin{array}{c}{\tau_{\rho}^{-1}=\tau_{\mathrm{j}}^{-1}=1} , \\ {\tau_{\mathrm{e}}^{-1}=\tau_{\xi}^{-1}=\tau_{v}^{-1}=0.65^{-1}}, \\ {\tau_{\mathrm{q}}^{-1}=\frac{8\left(2-\tau_{v}^{-1}\right)}{8-\tau_{v}^{-1}}}.\end{array} $

(3)

利用m=Mf可将密度分布函数从速度空间映射到动量空间。因此演化方程中右侧的碰撞项转化为

$ \boldsymbol{m}^{\prime}=\boldsymbol{m}-\boldsymbol{\varLambda}\left[\boldsymbol{m}-\boldsymbol{m}^{eq}\right]+\delta_{t}\left(\boldsymbol{I}-\frac{\boldsymbol{\varLambda}}{2}\right) \overline{\boldsymbol{S}}. $

(4)

其中：m′是碰撞后的动量矢量，I是单位矩阵，$\overline{\pmb{S}}$是动量空间的外力。

1.2 高密度比伪势多相模型

本文采用Shan-Chen单组份两相伪势模型来实现两相流模拟。该模型中流体之间作用力和流体与壁面间的相互作用力分别表示分别表示为^[20-22]：

$ \boldsymbol{F}=-G \psi(\boldsymbol{x}) \sum\limits_{a=1}^{N-1} w\left(\left|\boldsymbol{e}_{a}\right|^{2}\right) \psi\left(\boldsymbol{x}+\boldsymbol{e}_{a} \delta_{t}\right) \boldsymbol{e}_{a}, $

(5)

$ \boldsymbol{F}_{\mathrm{ads}}=-G_{\mathrm{w}} \psi(\boldsymbol{x}) \sum\limits_{\alpha=1}^{N-1} w\left(\left|\boldsymbol{e}_{\alpha}\right|^{2}\right) s\left(\boldsymbol{x}+\boldsymbol{e}_{\alpha} \delta_{t}\right) \boldsymbol{e}_{\alpha}. $

(6)

其中：G_w为可实现接触角的调节；w(|e_α|²)为权重系数，定义为w(1)=1/3和w(2)=1/12。ψ(x)为有效质量函数，为了实现高密度比模型，可定义为^[23]

$ \psi=\sqrt{\frac{2\left(P_{\mathrm{EOS}}-\rho c _{\mathrm{s}}^{2}\right)}{Gc_{\mathrm{s}}^{2}}}. $

(7)

其中P_EOS为状态方程。在本文选用分段线性状态方程为^[24]

$ p(\rho)=\left\{\begin{array}{ll}{\rho \theta_{\mathrm{V}}, } & {\rho \leqslant \rho_{1}} ;\\ {\rho_{1} \theta_{\mathrm{V}}+\left(\rho-\rho_{1}\right) \theta_{\mathrm{M}}, } & {\rho_{1}<\rho \leqslant \rho_{2}} ;\\ {\rho_{1} \theta_{\mathrm{V}}+\left(\rho_{2}-\rho_{1}\right) \theta_{\mathrm{M}}+\left(\rho-\rho_{2}\right) \theta_{\mathrm{L}}, } & {\rho>\rho_{2} .}\end{array}\right. $

(8)

其中：θ_V=0.49c_s²，θ_L=c_s²，θ_M=-0.04c_s²。本文取ρ₁=0.012和ρ₂=0.972来实现密度比100(对应两相的密度分别为ρ_V=0.01和ρ_L=1.0)。

为了保持热力学一致性，式(4)中的作用力项$\overline{\boldsymbol{S}}$修正为^[25]：

$ \begin{aligned} \overline{\boldsymbol{S}}_{1} &=\overline{\boldsymbol{S}}_{1}^{\prime}+\frac{12 \sigma\left(\boldsymbol{F}_{x}^{2}+\boldsymbol{F}_{y}^{2}\right)}{\psi^{2} \delta_{t}\left(\tau_{\mathrm{e}}-0.5\right)} ,\\ \overline{\boldsymbol{S}}_{2} &=\overline{\boldsymbol{S}}_{2}^{\prime}-\frac{12 \sigma\left(\boldsymbol{F}_{x}^{2}+\boldsymbol{F}_{y}^{2}\right)}{\psi^{2} \delta_{t}\left(\tau_{\xi}-0.5\right)}. \end{aligned} $

(9)

1.3 计算区域和边界条件

本文主要模拟GC和GDL中水的输运过程，系统截取400×250 lu的计算区域如图 1中虚线框所示, 这里的lu指无量纲格子单位(lattice unit，lu)。计算区域y方向小于200 lu处为GDL区域，大于200 lu处为GC区域。本文二维多孔介质采用随机方法生成(如图 1中灰色曲线所围成的区域)，包含了碳纤维及其涂覆物，其中平均纤维直径10 lu，孔隙率为0.88。初始时刻，一个密度为1.0 lu的液滴(如图 1中的黑色圆形)置于GDL区域内，其圆心坐标为(100, 100) lu，直径为30 lu，该液滴被周围密度0.01 lu的空气包裹(图 1中除去液滴和多孔介质的区域)。系统的上下边界为壁面无滑移边界，左右边界为周期边界条件。

为了简化波浪形流道引起的GDL中惯性力的影响，本文将惯性力以体积力的形式施加到流体上，使得计算区域内的流体产生波浪形运动，且此处设置惯性力(体积力)的方向与波浪形流道方向一致，即在计算区域左侧(x≤200 lu)流体的体积力指向右下方，计算区域右侧(x>200 lu)体积力则指向右上方。

该系统由两组无量纲数来描述:密度比ρ_liquid/ρ_gas和多孔介质的$R e_{p}=U \overline{D} / \nu$，其中：$\overline{D}$为碳纤维直径，ν为流体运动黏度。Re描述了流体惯性力和黏性力的关系，相关文献表明，当Re_p/(1-ϕ)≈1~10(其中ϕ表示孔隙率)时，惯性力的影响将不可忽略^[26]。故本文多孔介质中水的运动不仅需考虑表面张力，还应考虑惯性力的影响^[9]。

1.4 模型验证

本文通过在200×200 lu区域内改变液滴的初始半径r，来验证Young-Laplace定律。计算结果如图 2a所示，可以看出数值结果的液滴内外压力差与半径的倒数成线性关系，与Young-Laplace定律基本吻合。液滴与壁面的接触角大小通过改变G_w的值来实现，如图 2b所示，本文选取接触角为120°。

2 工况设置

本文主要模拟分析Re、波浪形流道角度、GDL上的开孔形状以及开孔位置对除水性能的影响，具体设置如表 1给出。由于需要考虑多孔介质中惯性力影响的临界Re_p约为1，因此本文所采用的Re范围为1至5。波浪形流道的角度与惯性力的方向保持一致，因此可以通过改变惯性力的方向从而实现波浪形角度的改变。GDL中开孔形状和位置如图 3所示。锥形孔的形状使用三角形替代，在保证三角形面积不变的情况下改变长宽比。开孔位置也包含了在波浪形流道上游、中间、下游3种。

3 结果与讨论 3.1 Re对除水性能的影响

Re的大小对应着所需施加惯性力的大小。图 4为不同Re工况下GDL中水的质量与初始总水量之比随时间的变化。Re=1时由于惯性力太小，液滴在GDL中缓慢运动，无法排出；随着Re的增加，惯性力增大，水逐渐从GDL中被排出，且排出速度逐渐增加。当Re=4时，GDL中水的质量出现剧烈波动。之后随着Re增加(Re=5)，惯性力持续增大，一方面使得液滴很快被吹出GDL，另一方面使得在流道中的液滴能停留在流道中，因此除水性能增强。

图 5表示水开始排出GDL的无量纲时刻(开始排出时间/总计算时间)随Re的关系。可以看出水开始排出GDL的时刻随着Re的增加而减小，这是因为Re增大，惯性力增大后，水的运动速度更快，从而更早地开始排出。从图 6的液滴动态行为中发现Re=4时GDL中水的质量出现剧烈波动是因为当水被吹到流道中后，由于液滴破碎程度大，导致不断有小液滴又重新进入GDL中引起其中的水量产生剧烈波动。

3.2 波浪形流道角度对除水性能的影响

本文分别研究对称波浪形流道角度和非对称波浪形流道角度大小对GDL除水性能的影响。图 7表示不同两侧对称夹角时GDL中水的质量与初始总水量之比随时间的变化关系。从图中可以看出不同角度时，水被吹出GDL的起始时刻几乎相同。两侧均为60°时，水更快速地排出GDL。两侧均为30°时，GDL中水量变化波动较大，且排出速度较慢。这是由于夹角越小，当水排出GDL后，在流道中向下的分力越大，因此越不稳定，导致GDL中水量的波动。

图 8和9表示两侧为非对称夹角时GDL中水的质量与初始总水量之比随时间的变化关系。图 8中，当左区域角度相同时，右区域夹角越小，排水性能越好。图 9中，当右区域角度相同时，左区域夹角越大，排水性能越好。从图 8和9的对比可以看出，左侧角度对排水性能的影响大于右侧角度的影响。因此应增大波浪形流道的入口角度(左侧)，减小波浪形流道的出口角度(右侧)，这样会更加利于排水。

3.3 GDL的开孔形状和位置对除水性能的影响

图 10定量统计了不同开孔形状对GDL排水性能的影响。从图中结果可以看出，开孔面积相同时，通孔、锥孔3以及不开孔的出水性能差别不大，而当开孔形状为锥孔2时，GDL的排水速率显著增加，锥孔1虽然排水速度和锥孔2差别不大，但波动较为剧烈。这是因为随着锥形孔的底边增长、高度减小，GDL中水的运动阻力逐渐减小，从而加快排水速率，但随着底边的继续增大、高度减小，锥孔难以深入GDL内部，水运动阻力减少程度变弱，因此排水速率几乎和无孔一样。同时注意到，孔形状对水开始从GDL中被排出的时刻几乎无影响。图 11表示锥孔位置对GDL排水性能的影响，结果表明，当孔越靠近GDL右侧时，排水速率越快。这是因为右侧即水的出口位置，右侧开孔使得GDL中水的运动阻力减小，从而增快排水速率。

4 结论

本文采用LBM-MRT高密度比伪势模型，分析了水的雷诺数Re、波浪形流道角度、GDL上的开孔形状、开孔位置对除水性能的影响，可以得到如下结论：1)随着Re的增加，GDL中除水性能增加；继续增加Re，会引起GDL中水量的剧烈波动；之后Re继续增加，除水性能增强；2)增大波浪形流道的入口角度，减小出口角度，会更加利于排水；3) GDL中开锥形孔2时，排水速率最快；在右侧开孔会增快GDL排水速率；4) GDL中水开始排出的时刻随着Re的增加而减小，而与夹角、开孔形状以及位置几乎没有关系。在GC和GDL的协同设计中，不仅要考虑水从GDL的排出，还应考虑GC中的水被重新吹入GDL中。

综上，当使用大出口角度小入口角度的波浪形流道和具有放置于波浪右侧的锥形孔2的GDL时燃料电池的排水性能最佳。本文所提出的锥孔型GDL以及针对锥孔型GDL和波浪形流道的优化可为未来的燃料电池在高电流密度下的水管理优化设计提供指导。