气体扩散层亲疏水特性对流体分布的预测

摘要

采用伪势两相格子玻尔兹曼方法（LBM）研究了非增湿条件下碳纸和碳布气体扩散层（GDL）的流体流动状态。通过随机方法和X射线扫描法构建了气体扩散层样本。为确保模型中表面张力和接触角计算的准确性，采用玻璃微珠模型进行验证，随后通过调整气体扩散层的亲疏水特性，分析流体在气体扩散层中流动的实时状态，得到了亲疏水特性对孔隙结构内水饱和度的影响规律。结果表明：疏水性气体扩散层中的水分布与亲水性气体扩散层中的水分布明显不同，较大的疏水性更有利于气体扩散层内水的排出；疏水性更强的气体扩散层显著提高了液态水进入气体扩散层的入口压力，导致催化剂层（CL）更容易受到水渗透的影响。

关键词

质子交换膜燃料电池（PEMFC）; 碳纸; 碳布; 气体扩散层（GDL）; 流动传输; 疏水性和亲水性

质子交换膜燃料电池（PEMFC）作为一种新型清洁能源装置，在未来能源领域发挥着重要作用，而其中水管理对PEMFC的性能有很大影响。气体扩散层（GDL）及其亲疏水性对于确保反应气体均匀分布在催化剂层中以及保持气体通道畅通至关重要，疏水性强的气体扩散层具有良好的排水性能，可以保证反应产生的水及时排出燃料电池。然而，疏水性过强会导致膜水含量降低，脱水干燥使得质子传输受阻，内部电阻上升。因此，需要控制气体扩散层的亲疏水性以达到最好的性能，为了实现这一目标，研究人员试图通过改变基底材料的性能获取亲疏水性较为适中的气体扩散层。

碳纸和碳布是PEMFC中2种常用的气体扩散层基底材料，它们的结构和性能有很大不同。碳纸中的碳纤维与渗碳树脂结合，而碳布由于其本身就具有纤维编织结构，因此不需要此类黏合剂。

研究人员重点研究了碳纸和碳布材料在亲疏水性方面的性能^［1］。Kowal等^［2］测试了碳布和碳纸在0.2~1.0 A·cm^-2电流密度下的含水性能。结果表明，碳纸的单位体积含水率比碳布高174%。Wang等^［3］预测，碳布气体扩散层适合于较高湿度条件，因为其孔隙结构迂曲度较低，粗糙的纹理表面保证了有效的水分离，而碳纸气体扩散层由于其较高的迂曲度而保留了更多的水分，因此更适合干燥的工作条件。

为了进一步研究气体扩散层材料的特性对PEMFC性能的影响，研究人员从理论和实验两方面进行了分析。Jordan等^［4］、Lee等^［5］、Wang等^［6］都关注到了膜电极（MEA）和气体扩散层的特征参数，包括厚度、孔隙度和亲水性等，这些早期的理论研究通常基于疏水性气体扩散层中两相传输的连续模型。Park ^［7］在非增湿条件下对实验中的水凝结进行了可视化，实验结果表明气体扩散层的最佳含水饱和度为0.1~0.3。该研究还发现，液态水的分布依赖于气体扩散介质的亲水性，而不是气体流量，尤其是在非增湿条件下。

然而，由于上述模型均基于宏观尺度，在考虑疏水性或亲水性多孔介质中水传输过程的微尺度现象时存在一些局限性，因此研究人员倾向于使用其他方法来探究微孔中的两相流动过程。Lee等^［8］使用孔隙网络模型研究疏水性多孔传输层内的微观液态水传输，并利用该模型预测了气体扩散层内部空间中的水分布，此外还研究了气体扩散层的肋宽和厚度等几何形状参数的影响。Salah等^［9］使用格子玻尔兹曼方法（LBM）模拟了3种具有恒定气体流速的几何通道模型以及1个液滴模型，发现矩形通道具有最佳的除水特性。Chen等^［10］利用LBM研究了疏水性气体扩散层中液态水的两相动力学行为。模拟结果表明，气体扩散层的毛细特性和润湿性强烈影响液态水的传输，而且随着亲水性的增加，气体扩散层中液态水饱和度水平通常会降低。Molaeimanesh等^［11］研究了4种不同分布的聚四氟乙烯（PTFE）对气体扩散层的影响，结果表明碳纤维和水滴之间的毛细管力为主要作用力，通过适当的处理方法，可以使得气体扩散层内的碳纤维不被水覆盖。

在介观模型中，许多关于气体扩散层的研究都集中在疏水性对输运现象的影响上。Jinuntuya等^［12］利用基于LBM的X射线断层扫描研究了气体扩散层结构对水传输的影响，并测试了气体扩散层在不同接触角和不同压力下的传输过程。结果表明，稳定驱替和毛细管浸入不是一个渐进过程。当接触角为60°<θ<90°时，表现出稳定的驱替特性；当接触角为100°<θ<140°时，表现出毛细管指进特性。Kakaee等^［13］利用多相LBM模型研究了PTFE在气体扩散层上的分布对水滴去除的影响，采用3种不同的PTFE分布和3种不同的黏合剂含量进行了模拟。结果表明，当气体扩散层的黏合剂含量为30%或70%时，PTFE分布较为均匀，水管理效果也较好。Shakerinejad等^［14］利用LBM和气体扩散层计算机随机重构方法研究了气体扩散层中液态水的动力学。结果表明，气体扩散层界面附近结构的亲水性增加比气体扩散层纤维接触角的增加更加有效。

通过X射线扫描法构建气体扩散层结构体，随后利用Shan等^［15］提出的伪势两相格子玻尔兹曼模型构建了气体扩散层内的两相流动模型，并使用玻璃微珠结构进行了模型验证，随后利用该模型比较了不同亲疏水性碳布和碳纸的流动特性，以观察流体在气体扩散层中的传输现象，得出亲疏水性变化时不同材料气体扩散层的饱和度变化。

1 模型及算法

1.1　材料准备

1.1.1　随机生成碳纸气体扩散层

与其他方法相比，计算机随机生成方法具有成本低、易于实现等优点。随机生成碳纸的方法由Schladitz等^［16］提出，该方法也被其他格子玻尔兹曼模拟所采用^［17-21］。重构方法从随机生成直线开始，用这些直线表示具有固定和均匀半径的圆柱体。由于所绘制的线具有不同的分布，因此使用极坐标系中的方向密度函数$\rho (\theta ,\varphi )$、纬度$\theta$以及经度$\varphi$表示，如下所示：

$$\rho (\theta ,\varphi )=\frac{\mathrm{1}}{\mathrm{4}\mathrm{\pi }}\frac{\beta \mathrm{s}\mathrm{i}\mathrm{n}\theta }{(\mathrm{1}+({\beta }^{\mathrm{2}}-\mathrm{1})\mathrm{c}\mathrm{o}{\mathrm{s}}^{\mathrm{2}}{\theta )}^{\mathrm{3}/\mathrm{2}}}$$

（1）

式中：$\beta$为各向异性参数。随着$\beta$的增加，纤维将越来越平行于X‒Y平面。

1.1.2　X射线扫描碳纸和碳布

文献［22-23］中给出了三维气体扩散层图像的采集和重建细节。本研究中使用的玻璃微珠结构来源于文献［24］。碳纸和碳布气体扩散层具有较强的各向异性，其几何结构在不同方向上有较大的区别。在格子玻尔兹曼模拟中，孔隙空间用0表示，实体用1表示。后文将以上述文献中的碳纸和碳布气体扩散层结构作为基础构建模型。

1.2　伪势多相模型

伪势多相LBM^［15］计算需要2个步骤：第一步是流体粒子之间的碰撞，第二步是流体粒子的迁移。流体粒子分布函数的演变描述如下所示：

$$f_{a,k}^{}(x+c_a\mathrm{\Delta }t,t+\mathrm{\Delta }t)-f_{a,k}^{}(x,t)=-\mathit{M}\mathit{S}{\mathit{M}}^{-\mathrm{1}}(f_{a,k}^{}(x,t)-f_{a,k(\mathrm{e}\mathrm{q})}^{}(x,t))$$

（2）

式中：$f_{a,k}^{}$是流体k的粒子速度分布函数，沿$a$方向以速度$c_a$移动；$t$是量纲一时间；$x$是格点位置；$\mathit{M}$是变换矩阵；$\mathit{S}$是系数矩阵l；$f_{a,k(\mathrm{e}\mathrm{q})}^{}$是密度$\rho$和流体宏观速度$u$的函数，表示平衡状态下的粒子速度分布函数。$f_{a,k(\mathrm{e}\mathrm{q})}^{}$计算方法如下所示：

$$f_{a,k(\mathrm{e}\mathrm{q})}^{}={\omega }_a{\rho }_k\left(\mathrm{1}-\frac{\mathrm{3}}{\mathrm{2}}\frac{u_{\mathrm{e}\mathrm{q}}{u}_{\mathrm{e}\mathrm{q}}}{c_{\mathrm{s}}^{\mathrm{2}}}\right),a=\mathrm{0}$$

$$f_{a,k(\mathrm{e}\mathrm{q})}^{}={\omega }_a{\rho }_k\left(\mathrm{1}+\mathrm{3}\frac{c_a{u}_{\mathrm{e}\mathrm{q}}}{c_{\mathrm{s}}^{\mathrm{2}}}+\frac{\mathrm{9}}{\mathrm{2}}\frac{(c_a{u}_{\mathrm{e}\mathrm{q}}{)}^{\mathrm{2}}}{c_{\mathrm{s}}^{\mathrm{4}}}-\frac{\mathrm{3}}{\mathrm{2}}\frac{u_{\mathrm{e}\mathrm{q}}{u}_{\mathrm{e}\mathrm{q}}}{c_{\mathrm{s}}^{\mathrm{2}}}\right),a=\mathrm{1},\cdots ,\mathrm{18}$$

（3）

式中：${\rho }_k$为流体$k$的密度，${\rho }_k={\displaystyle \sum _{a=\mathrm{0}}^{\mathrm{18}}}{f}_{a,k}^{}$；$c_{\mathrm{s}}$为格子声速，采用D3Q19模型，速度分布如图1所示；$u_{\mathrm{e}\mathrm{q}}$为宏观平衡速度；${\omega }_a$为基于D3Q19模型的权重系数。$c_a$计算式如下所示：

图1  D3Q19格子玻尔兹曼模型中19个格子速度

Fig.1  Nineteen lattice velocities used in the D3Q19 lattice Boltzmann model

$$c_a=\left\{\begin{array}{l}(\mathrm{0,0},\mathrm{0})c/\mathrm{\Delta }t,a=\mathrm{0}\\ (\pm \mathrm{1,0},\mathrm{0})c/\mathrm{\Delta }t,(\mathrm{0},\pm \mathrm{1,0})c/\mathrm{\Delta }t,(\mathrm{0,0},\pm \mathrm{1})c/\mathrm{\Delta }t,a=\mathrm{1},\cdots ,\mathrm{6}\\ (\pm \mathrm{1},\pm \mathrm{1,0})c/\mathrm{\Delta }t,(\pm \mathrm{1,0},\pm \mathrm{1})c/\mathrm{\Delta }t,(\mathrm{0},\pm \mathrm{1},\pm \mathrm{1})c/\mathrm{\Delta }t,a=\mathrm{7},\cdots ,\mathrm{18}\end{array}\right.$$

（4）

式中：$c$是单位格子长度。${\omega }_a$计算式如下所示：

$${\omega }_a=\left\{\begin{array}{l}\mathrm{1}/\mathrm{3},a=\mathrm{0}\\ \mathrm{1}/\mathrm{18},a=\mathrm{1,2},\cdots ,\mathrm{6}\\ \mathrm{1}/\mathrm{36},a=\mathrm{7,8},\cdots ,\mathrm{18}\end{array}\right.$$

（5）

在式（2）中，$\mathit{S}={\mathit{M}}^{-\mathrm{1}}{\mathit{S}}_{\mathit{k}}\mathit{M}$是描述流体粒子之间碰撞的矩阵，其中$\mathit{M}$是将粒子分布函数转换为D’Humières等^［25］给出的力矩变换矩阵，${\mathit{S}}_{\mathit{k}}$是流体k在力矩空间中执行碰撞的对角矩阵。${\mathit{S}}_{\mathit{k}}$的具体形式如下所示：

$${\mathit{S}}_{\mathit{k}}=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}(s_{k\mathrm{0}},s_{k\mathrm{1}},s_{k\mathrm{2}},s_{k\mathrm{3}},s_{k\mathrm{4}},s_{k\mathrm{5}},s_{k\mathrm{6}},s_{k\mathrm{7}},s_{k\mathrm{8}},s_{k\mathrm{9}},s_{k\mathrm{10}},s_{k\mathrm{11}},s_{k\mathrm{12}},s_{k\mathrm{13}},s_{k\mathrm{14}},s_{k\mathrm{15}},s_{k\mathrm{16}},s_{k\mathrm{17}},s_{k\mathrm{18}})$$

（6）

$$s_{k\mathrm{0}}=s_{k\mathrm{2}}=s_{k\mathrm{5}}=s_{k\mathrm{7}}=\mathrm{0},s_{k\mathrm{1}}=s_{k\mathrm{3}}=s_{k\mathrm{9}}=s_{k\mathrm{10}}=s_{k\mathrm{11}}=s_{k\mathrm{12}}=s_{k\mathrm{13}}=s_{k\mathrm{14}}=s_{k\mathrm{15}}=\frac{\mathrm{1}}{{\tau }_k},s_{k\mathrm{4}}=s_{k\mathrm{6}}=s_{k\mathrm{8}}=s_{k\mathrm{16}}=s_{k\mathrm{17}}=s_{k\mathrm{18}}=\mathrm{8}\left(\left(\mathrm{2}-\frac{\mathrm{1}}{{\tau }_k}\right)\left(\mathrm{8}-\frac{\mathrm{1}}{{\tau }_k}\right)\right)$$

（7）

式中：${\tau }_k$是与流体$k$黏度相关的参数，也被称为松弛时间。在上述格子玻尔兹曼模型中，碰撞在力矩中计算如下所示：

$$f_{a,k}(x_a+c_a\mathrm{\Delta }t,t+\mathrm{\Delta }t)=f_{a,k}(x_a,t)-{\mathit{M}}^{-\mathrm{1}}\mathit{S}(m_{a,k}(x_a,t)-m_{a,k(\mathrm{e}\mathrm{q})}(x_a,t))$$

（8）

式中：$m_{a,k}=\mathit{M}{f}_{a,k}$为动量矩阵，且$m_{a,k(\mathrm{e}\mathrm{q})}=\mathit{M}{f}_{a,k(\mathrm{e}\mathrm{q})}$。变换矩阵$\mathit{M}$的计算式如下所示：

$$M_{ij}={\phi }_{ij}$$

（9）

式中：${\phi }_{ij}$代表矩阵第$i$行第$j$列的元素。$\mathit{M}$由正交基构成，如下所示：

$${\phi }_{\mathrm{0}a}={\left|e_a\right|}^{\mathrm{0}},{\phi }_{\mathrm{1}a}=\mathrm{19}{\left|e_a\right|}^{\mathrm{2}}-\mathrm{30},{\phi }_{\mathrm{2}a}=\frac{\mathrm{21}{\left|e_a\right|}^{\mathrm{4}}-\mathrm{53}{\left|e_a\right|}^{\mathrm{2}}+\mathrm{24}}{\mathrm{2}},{\phi }_{\mathrm{3}a}=e_{ax},{\phi }_{\mathrm{4}a}=(\mathrm{5}{\left|e_a\right|}^{\mathrm{2}}-\mathrm{9})e_{ax},{\phi }_{\mathrm{5}a}=e_{ay},{\phi }_{\mathrm{6}a}=(\mathrm{5}{\left|e_a\right|}^{\mathrm{2}}-\mathrm{9})e_{ay},{\phi }_{\mathrm{7}a}=e_{az},{\phi }_{\mathrm{8}a}=(\mathrm{5}{\left|e_a\right|}^{\mathrm{2}}-\mathrm{9})e_{az},{\phi }_{\mathrm{9}a}=\mathrm{3}{e}_{ax}^{\mathrm{2}}-{\left|e_a\right|}^{\mathrm{2}},{\phi }_{\mathrm{10}a}=\mathrm{3}{e}_{ax}^{\mathrm{2}}-{\left|e_a\right|}^{\mathrm{2}},{\phi }_{\mathrm{11}a}=\mathrm{3}{e}_{ax}^{\mathrm{2}}-{\left|e_a\right|}^{\mathrm{2}},{\phi }_{\mathrm{12}a}=-\mathrm{5}(e_{ay}^{\mathrm{2}}-e_{az}^{\mathrm{2}}),{\phi }_{\mathrm{13}a}=e_{ax}{e}_{ay},{\phi }_{\mathrm{14}a}=e_{ay}{e}_{az},{\phi }_{\mathrm{15}a}=e_{ax}{e}_{az},{\phi }_{\mathrm{16}a}=(e_{ay}^{\mathrm{2}}-e_{az}^{\mathrm{2}})e_{ax},{\phi }_{\mathrm{17}a}=(e_{az}^{\mathrm{2}}-e_{ax}^{\mathrm{2}})e_{ay},{\phi }_{\mathrm{18}a}=(e_{ax}^{\mathrm{2}}-e_{ay}^{\mathrm{2}})e_{az}$$

（10）

式中：$e_{ax}\mathrm{、}{e}_{ay}\mathrm{、}{e}_{az}$分别为x、y、z 3个方向上的速度分量；$e_a$为合速度。式（3）中流体$k$的平衡速度$u_{\mathrm{e}\mathrm{q}}$通过下式计算：

$${\rho }_k{u}_{\mathrm{e}\mathrm{q}}={\rho }_{k}u\mathrm{\text{'}}+{\tau }_k{F}_k$$

（11）

式中：$u\mathrm{\text{'}}$为体积速度；$F_k$为流体作用力。$u\mathrm{\text{'}}$可通过下式计算：

$$u\mathrm{\text{'}}=\frac{\sum _k{\rho }_k{u}_k/{\tau }_k}{\sum _k{\rho }_k/{\tau }_k}$$

（12）

流体作用力$F_k$包括流体‒流体相互作用力$F_{\mathrm{1},k}$和流体‒固体相互作用力$F_{\mathrm{2},k}$。$F_k=F_{\mathrm{1},k}+F_{\mathrm{2},k}$是每个粒子上的总作用力。流体‒流体相互作用力定义如下所示：

$$F_{\mathrm{1},k}(x)=-{\rho }_k(x)\sum _{x\mathrm{\text{'}}}\sum _{\overline{k}}{G}_{k\overline{k}}(x,x\mathrm{\text{'}}){\rho }_{\overline{k}}(x^{\mathrm{\text{'}}})(x\mathrm{\text{'}}-x)$$

（13）

$G_{k\overline{k}}(x,x\mathrm{\text{'}})=G_{\overline{k}k}(x,x\mathrm{\text{'}})$，其中$x$和$x\mathrm{\text{'}}$分别为流体粒子原位置和下一时刻的位置。本研究中，$G_{k\overline{k}}(x,x\mathrm{\text{'}})$由以下分段函数定义：

$$G_{k\overline{k}}(x,x\mathrm{\text{'}})=\left\{\begin{array}{l}{g}_{k\overline{k}},\left|x-x\mathrm{\text{'}}\right|=\mathrm{1}\\ \frac{g_{k\overline{k}}}{\mathrm{2}},\left|x-x\mathrm{\text{'}}\right|=\sqrt[]{\mathrm{2}}\\ \mathrm{0},\mathrm{其他}\end{array}\right.$$

（14）

式中：$g_{k\overline{k}}$为流体之间的作用力常数。流体‒固体作用力定义如下所示：

$$F_{\mathrm{2},k}(x)=-{\rho }_k(x)G_{k\mathrm{s}}(x,x\mathrm{\text{'}})n_{\mathrm{s}}(x\mathrm{\text{'}})(x\mathrm{\text{'}}-x)$$

（15）

式中：$n_{\mathrm{s}}(x\mathrm{\text{'}})$在流固界面处为常数，若不在流固界面处则为零；$G_{k\mathrm{s}}(x,x\mathrm{\text{'}})$为流体‒固体相互作用强度，定义固体壁面的亲疏水性。

为确保与流体‒流体相互作用的一致性，流体‒固体相互作用强度的描述$G_{ks}(x,x\mathrm{\text{'}})$如下所示：

$$G_{k\mathrm{s}}(x,x\mathrm{\text{'}})=\left\{\begin{array}{l}{g}_{ks},\left|x-x\mathrm{\text{'}}\right|=\mathrm{1}\\ \frac{g_{ks}}{\mathrm{2}},\left|x-x\mathrm{\text{'}}\right|=\sqrt[]{\mathrm{2}}\\ \mathrm{0},\mathrm{其他}\end{array}\right.$$

（16）

式中：$g_{k\mathrm{s}}$为流体和固体之间的作用力常数。

1.3　模型假设

在建模过程中，作以下假设：

（1） 模型中的流体包括气体和液体，是不可压缩的，具有恒定的密度。

（2） 模型处于恒温状态。

（3） 气体分子或液体分子与固体表面之间没有摩擦。

2 模型验证

2.1　表面张力验证

为了确保该模型中表面张力计算的准确性，使用数值方法验证两相平衡后液滴半径和液滴内外压差是否符合拉普拉斯定律。网格的每个边界均采用周期边界条件。在远离界面的节点处进行压力测量，测量数据以点阵单位表示。在模拟中，2个松弛时间设置为${\tau }_{\mathrm{1}}=\mathrm{1.0}$和${\tau }_{\mathrm{2}}=\mathrm{1.2}$。量纲一流体‒流体相互作用系数设置为0.000 1。从图2可以看出，所有点都与直线吻合良好，可以保证模型计算的准确性。

图2  气泡表面压力随气泡半径R的变化

Fig.2  Variation of bubble surface pressure with bubble radius R

2.2　接触角验证

随着表面流体润湿性的增加，非润湿相开始从表面分离，并与表面形成接触角。接触角可定义为液滴与所在固相平面之间的角度，如图3所示。液滴的接触角$\theta$可通过下式计算：

$$\mathrm{t}\mathrm{a}\mathrm{n}\theta =\frac{L}{\mathrm{2}(R-H)}$$

（17）

式中：$H$是液滴的高度；$L$是液滴与壁面的接触长度；$R$是液滴的半径。

图3  接触角示意图

Fig.3  Schematic diagram of contact angle

水‒固体接触角取决于固液接触常数$g_{k\mathrm{s}}$的值。为了确定接触角与$g_{k\mathrm{s}}$的关系，将半径为10个格子单元的水滴放置在固体壁上。然后，分别选择$g_{k\mathrm{s}}$为-0.02、0、0.02进行模拟，最终达到平衡时计算的接触角为120°、90°、60°。.

3 结果和讨论

3.1　疏水性和亲水性多孔介质的影响

使用X射线断层扫描技术获取碳纸和碳布气体扩散层结构体，结合伪势两相LBM计算模型计算非增湿条件下不同材料表面亲疏水性对水传输的影响。

3.1.1　亲疏水性不同的玻璃微珠结构

由于玻璃微珠的结构均匀，3个维度都具有一定的相似性，各向同性较高，因此首先使用玻璃微珠结构结合格子玻尔兹曼模型计算了表面疏水性和亲水性对水传输的影响。图4为液态水在结构体中的浸入过程。

图4  接触角为120°和60^o的玻璃微珠结构中的水分布

Fig.4  Water distribution in the microstructure of glass bead with contact angles of 120 °and 60 °

从图4可以看出，玻璃微珠的微孔结构呈明显的球形堆积结构，孔隙主要存在于球之间。当注入水时，水沿着球体边缘进入孔隙。比较120°和60°接触角可以发现，当接触角为120°时，球体表面的水分布更均匀，水覆盖的区域更宽；当接触角为60°时，水更多地集中于固体表面，并且总量较少。总的来说，在疏水性结构中流体表现出更多的碎片化形态，而亲水性结构中流体表现为连续状态。

3.1.2　碳纸亲疏水性对于两相流动的影响

碳纸的疏水性较高，其中的液态水更容易呈现破碎状分布，疏水性结构中的液体更倾向于停留在结构中的较大孔隙中，如图5所示。

图5  碳纸结构体表面水分布

Fig.5  Surface water distribution of carbon paper structure

图6a、b展示了碳纸结构体水分布等值线图及截面图。疏水性气体扩散层中的水分布与亲水性气体扩散层中的水分布明显不同。亲水性气体扩散层中的水浸入类似活塞流，几乎占据了所有可用的孔隙空间。在疏水性气体扩散层中，流体呈现碎片分布，留下较多的孔隙空间供气体流动。

图6  碳纸结构体水分布等值线图及截面图

Fig.6  Surface water distribution contour map and cross-sectional view of carbon paper structure

通过观察图6a中等值线图的放大图，以及图6b中的截面图，可以清楚地看到疏水结构表面的水状态为碎片状。当固液作用力常数为正时，接触角大于90°，表现出疏水性，液态水需要克服较大的毛细力才能在面内方向膨胀，因此只能占据相对较小的孔隙。当固液作用力常数为负时，接触角小于90°，表现出亲水性，液态水进入孔隙的毛细管阻力降低。因此，对于亲水性气体扩散层结构，液态水不仅会占据气体扩散层中较大的孔，还会占据气体扩散层中更多较小孔径的孔，液态水更大程度地存留在孔隙中。

3.1.3　碳布亲疏水性对于两相流动的影响

虽然碳布和碳纸均为基于碳纤维的多孔材料，但碳纸内纤维处于无序分布状态，而碳布中的纤维分布更加有序，因此碳布比碳纸迂曲度更小，孔隙结构复杂度相对较低，如图7、8所示。

图7  碳布结构示意图

Fig.7  Schematic diagram of carbon cloth structure

图8  碳布结构体水分布等值线图及截面图

Fig.8  Surface water distribution contour map and cross-sectional view of carbon cloth structure

与碳纸相比，疏水性碳布中的水分布更均匀，表面的分布面积更大。从图9和图10可以发现，压缩条件下疏水性碳布的水分散度仍高于亲水性碳布，与未压缩情况相比效果较为明显。观察相同位置可以看出，纤维结构发生了巨大变化，整体压缩变形^［26］。同时，孔隙面积也变得更小、更复杂，使流体传输阻力上升，传输速率降低。由于碳布结构在压缩后本身孔隙空间更小，因此疏水性对液体输送的影响减小，液体的传输速率主要取决于压缩程度。总的来说，不管是碳布还是碳纸，疏水性更强的材料更有利于液体的传输，高效的液体传输能够保证气体不被堵塞，但过高的疏水性也会导致气体扩散层干燥，使得传输阻力增加。

图9  120° 接触角碳布不同压力下浸入情况

Fig.9  Immersion of carbon cloth with 120° contact angle under different pressures

图10  60° 接触角碳布不同压力下浸入情况

Fig.10  Immersion of carbon cloth with 60° contact angle under different pressures

3.2　疏水性和亲水性气体扩散层的液体饱和曲线

饱和度定义为水所占孔隙面积与总孔隙面积之比。对于疏水性/亲水性碳纸和碳布气体扩散层，通过格子玻尔兹曼模型中水的格点数与总孔隙格点数的比例计算饱和度。图11给出了稳态时沿气体扩散层厚度方向的液体饱和度分布。沿着气体扩散层厚度方向，亲水性/疏水性气体扩散层的液体饱和度都会逐渐降低。亲水性气体扩散层中液态水进入孔隙的毛细管阻力较小，因此亲水性气体扩散层中除了较大的孔，更多较小孔径的孔也被液态水占据，液态水进入气体扩散层后，会储存在更多的孔隙中，从而提高气体扩散层中的液体饱和度。对于疏水性气体扩散层，液态水需要克服较大的毛细管阻力，只能进入较少的半径较大的孔径，大部分的孔受限于较大的毛细管阻力难以浸入，导致整体饱和度较低，含水量较少。

图11  沿气体扩散层厚度方向的液体饱和度分布

Fig.11  Liquid saturation profile along thickness direction of GDL

4 结论

（1）疏水性气体扩散层中的水分布与亲水性气体扩散层中的水分布明显不同。亲水性气体扩散层中的水侵入类似活塞流，几乎占据了所有可用的孔隙空间。在疏水性气体扩散层中，流体则呈现碎片分布，留下较多孔隙空间供气体流动。此外，对比碳纸和碳布气体扩散层表面的液态水状态可以发现，液态水覆盖在碳纸气体扩散层表面，碳布比碳纸更粗糙，因此液态水覆盖率也更低。

（2）使用格子波尔兹曼模型计算了亲水性/疏水性气体扩散层中的液体饱和度，液态水进入亲水性气体扩散层后，会在更多的孔隙中停留，从而提高气体扩散层中的液体饱和度。对于疏水性气体扩散层，液态水需要克服较大的毛细管力，且只能占据较少的半径较大的孔隙，因此疏水性更强的气体扩散层显著提高了液态水进入气体扩散层的入口阻力，使得气体扩散层保持较低的饱和度。

作者贡献声明

高 源：前期材料整理，催化层结构扫描，指导论文写作与修改。

丁兆丰：论文模型构建，论文写作与修改。