基于棉线的微流体燃料电池阳极传质特性LB模拟

穆嫒萍，叶丁丁，陈蓉，朱恂，廖强

（1 重庆大学低品位能源利用技术及系统教育部重点实验室，重庆400030；2重庆大学能源与动力工程学院工程热物理研究所，重庆400030）

引 言

无膜微流体燃料电池因不存在由质子交换膜引起的与膜相关的一系列问题[1-3]，使其成为极具应用前景的新型微型电源之一。而传统的无膜微流体燃料电池常需外部驱动设备来维持燃料和氧化剂的层流流动，但该设备不仅会增加整个运行系统的体积，还会降低电池的净输出功率，限制其进一步向微型化发展。常见的可替代外部驱动设备的驱动方式有：流体的表面张力驱动[4]、重力驱动[5-6]、毛细力驱动[7-9]、化学反应驱动[10]等。和其他驱动方式相比，采用纸、纤维等多孔材料的毛细力驱动方式有着结构简单、制造成本低、占用空间小等优点[8]。但是当纸和纤维的孔隙被溶液充满后，其毛细驱动力减弱，难以形成持续稳定的流动，因此，在微流体装置中常采用与重力驱动相结合的方式来驱动流体以形成稳定流动[11-12]。因棉线有着天然三维的多孔结构、较好的柔韧性以及在浸润情况下具有较好的机械强度等优点，2018 年，Liu 等[13]构建了基于棉线并辅以重力驱动的微流体燃料电池，其最大功率密度和电流密度比文献[14]报道的基于织物类的微流体燃料电池的性能高3.7 和1.5 倍。2020 年，Wu等[15]直接将棉线平铺于电极表面作为流体流道，构建了基于棉线的微流体燃料电池，在碱性条件下电池的极限电流密度在同类电池中最高，而且该类电池的性能严重受到阳极侧燃料传质的限制，但目前对基于多孔纤维类的微流体燃料电池的相关理论研究工作还未见报道。因此，采用数值模拟的方法来研究该类燃料电池内部的传输机理将为电池性能的提升提供理论指导。

与其他数值模拟方法相比，格子Boltzmann方法（lattice Boltzmann method, LBM）能有效处理复杂的微观几何结构以及多孔介质中多相流体流动等方面的问题[16]。虽然近年来许多学者通过基于孔隙尺度的LBM 来研究燃料电池催化层或扩散层内的流动及传质问题[17-20]，但该方法对每个孔隙内的流动进行精确处理的要求限制了其模型计算域的尺度。与孔隙尺度不同，表征体元（representative elementary volume, REV）尺度是通过适当的作用力来描述流体与多孔介质之间的相互作用进而获取体积平均的宏观流动量，并不需要获取多孔介质的结构，因此常有学者使用该方法来研究复杂或无法直接获取的多孔介质结构内的流动特性[21-23]。而在采用多孔纤维的微流体燃料电池中，输运溶液的流道是由具有多尺度孔隙结构的棉线堆叠而成，具有复杂的多孔介质结构。因此本文采用LB 方法研究REV 尺度下基于棉线的微流体燃料电池阳极内流动及耦合电化学反应的传质特性。

本文首先介绍计算采用的物理模型；随后介绍LB 数值模拟方法并进行LBM 模型验证，证明该流动及传质模型的正确性；然后研究受重力驱动的基于棉线的微流体燃料电池阳极侧流动及传质特性，计算获得燃料浓度分布及阳极极化曲线，并讨论不同燃料进口浓度及流量对电池阳极性能及传质特性的影响。

1 物理模型

根据前期构建的基于棉线的微流体燃料电池[15]建立物理模型。实验中采用了医用亲水纱布上的棉线作为可输运阳极和阴极溶液的流道，通过控制溶液进出口高度差形成的重力维持溶液的连续流动。在该电池结构中，阳极和阴极流道分别用10根棉线平铺堆叠于负载有催化剂（有效面积为1 mm×10 mm）的碳纸电极上，为防止阴、阳极溶液发生扩散混合，在二者中间加入由40 根棉线平铺堆叠而成的隔离液通道。在电池运行过程中，阳极液（甲酸钠与氢氧化钾的混合溶液）和阴极液（含饱和溶解氧的氢氧化钾溶液）在重力驱动下输运到负载有催化剂的碳纸电极上发生电化学反应，阳极甲酸根与氢氧根反应生成碳酸根、水和电子，阴极在碱性条件下反应生成氢氧根，氢氧根经隔离液到达阳极，电子经外电路到达阴极，形成回路产生电流。

基于棉线的微流体燃料电池运行初期，由棉线多孔介质的毛细力和重力驱动输运溶液，当多孔介质被溶液充满后，电池靠重力的作用持续供液，使其稳定运行。溶液在棉线流道内输运的过程中，其多孔流道主要有三个尺度的孔隙结构：一是棉线与棉线之间的间隙，二是单根棉线中棉线纤维之间的间隙，三是单根棉线纤维内部的间隙[15]。因此，本文在理论计算模型中采用的基本假设如下：

（1）微流体燃料电池在常温常压下运行，燃料溶液为不可压缩流体，且在多孔流道内为等温稳态流动。

（2）不考虑棉线未被溶液充满的不稳定阶段，即基于棉线的微流体燃料电池在重力驱动的作用下稳定运行。

（3）单根棉线内部结构复杂且难以进行三维重构，因此，忽略单根棉线内部复杂的多孔结构，假设单根棉线为具有一定孔隙率ɛ的圆柱。

（4）催化层的厚度仅为几十微米，本文主要考虑棉线流道内的流动和传质，不考虑催化层内部的传输，即将催化层考虑为反应界面。

模型计算区域如图1（a）所示，选取棉线组成的基本单元进行计算，考虑阳极电极布置在计算区域底部位置，棉线与棉线之间为等边三角形排列，如图1（b）所示，图中绿色表示棉线，其余为棉线之间的间隙。

图1 模型计算区域示意图（a）及棉线堆叠方式示意图（b）Fig.1 Schematics of computational domain(a)and cotton thread array(b)

甲酸钠溶液在阳极电极处发生电化学反应

其物质守恒方程为

式中，u为速度；C和Deff分别为物质浓度和有效扩散系数。在该计算模型中，溶液在棉线内部的传质阻力要比棉线与棉线间间隙内的传质阻力大，采用该方程描述溶液在棉线与棉线之间的传质时，其扩散系数为D，而溶液在单根棉线内部的扩散系数则修正为Deff=ɛ1.5D。其反应源项Rs=-Ilocal/2F，局部电流Ilocal采用Butler-Volmer方程描述为

其中，Iref为参考电流密度；C为反应界面上的当地浓度；Cref为参考浓度；β为反应级数；αa为阳极传输系数；F为法拉第常数；R为通用气体常数；T为温度；η为过电势。

2 格子Boltzmann方法

2.1 流动模型

不可压缩流体在多孔介质内流动的控制方程为[24]：

式中，u、P为表观速度和压力；ρ0为参考密度；υe为有效运动黏性系数；F为包括了介质阻力及其他外力的总体积力，即

式中，K、υ分别为渗透率和流体运动黏性系数，与平均孔隙率ɛ相关的结构函数Fɛ可用式（6）表示[25]

低速流动时，总体积力F[式（5）]中包含结构函数Fɛ的非线性介质阻力项可以忽略。

对于式（4）不可压缩流体在多孔介质中的等温流动模型，Guo 和Zhao[26]于2002 年提出求解该问题的广义LBE模型（GLBE），演化方程为

其中,

在本文采用的D3Q19 的流动模型参数中的权重因子ωi及各方向粒子离散速度ei为

因速度定义项的外力F包含速度，根据Guo 和Zhao[26]提出的临界速度求解得出

式中c0与c1为

对式（7）做Chapman-Enskog分析可知流体运动黏性系数与松弛时间的关系为

当ɛ=1 时，该GLBE 模型演变为用于求解自由流动区域速度分布情况的LBE模型。

2.2 传质模型

在LBM 传质模型中，当物质的扩散系数D与流动模型中的运动黏性系数υ相差较大时，常用的格子Bhatnagar-Gross-Krook（LBGK）模型会出现数值不稳定现象。Inamuro[27]提出在平衡分布函数中添加修正项的格子动理格式（lattice kinetic scheme，LKS）模型中，通过调节该参数可使松弛因子稳定在合理范围内，改善数值不稳定现象。基于该模型，Yang 等[28]提出用于求解对流扩散方程的修正LBGK模型，其演化方程为[28-29]

其时间导数可由向后差分获得

在该模型中宏观参数C为

对式（15）做Chapman-Enskog 分析可知扩散系数与松弛时间的关系及浓度梯度∇C为

为保证数值稳定性，本文中τυ及τD均设为1，则可调参数A可根据τD及扩散系数D来确定。

边界条件设定：进口（X=0）处已知压力Pin和浓度Cin；出口（X=L）为已知出口压力Pout，在该边界处浓度为充分发展边界；左右（Y=0，Y=W）为周期性边界；上边界（Z=H）处速度u= 0，浓度梯度为0；下边界（Z=0）处速度u=0 且在该界面发生电化学反应时，采用Kang等[30]报道的反应边界方法来求解。

2.3 LB模型验证

采用两平板间充满多孔介质的流体流动来验证REV-LBM 模型的正确性。流体在流道内的速度分布表述为[26]

在本文中υeυ=1，ɛ、K分别为孔隙率和渗透率，G为驱动流体流动的外力，其值为G= ΔP/(ρL)。解析解为

图2 沿Y方向的速度分布Fig.2 Velocity profile along with Y direction

所采用的边界条件为：进出口为已知压力边界，上下边界速度u=0。图2 为孔隙率ɛ为1.0、0.5、0.2 时沿Y方向的速度分布。从图中可以看出LB 的模拟结果与解析解吻合较好，证明了REV-LBM 模型的正确性。

采用D3Q7修正后的LBGK传质模型验证对流-扩散-反应模型。其边界条件为：进出口为已知浓度边界（Cin=1，Cout=0），上下边界的浓度梯度为0，左右为周期性边界。反应源项Rs=-KC（K=0、2、10）。图3为沿中心线上的浓度分布，可以看出当可调参数A=0 和A≠0 时的数值解与其解析解[31]吻合较好，可证明该模型的正确性。

图3 可调参数A=0、A≠0时沿X方向的平均浓度分布Fig.3 Average concentration distribution along with X direction when additional parameter A=0 and A≠0

3 计算结果与讨论

3.1 模型验证

通过LB 方法计算得到的阳极极化曲线与实验结果对比如图4所示，其计算参数取值见表1。从图4 可知，在电位较低时，理论结果与实验结果有一定差异，这主要因为电流密度较小时，Butler-Volmer公式不能准确描述此时的活化极化损失。但整体上计算结果与实验结果吻合较好，证明该模型及参数的适用性。

表1 计算参数Table 1 Calculation parameters used in the model

图4 LB计算与实验[15]得到的阳极极化曲线对比Fig.4 Comparison of anode polarization curve of LB method and experimental data[15]

3.2 流道内的速度与浓度分布

通过计算获得流道内的速度分布如图5 所示。由于不考虑单根棉线内部的多孔结构，仅求解REV尺度下体积平均的宏观流量，所以棉线内部流速分布均匀，且由于棉线内部多孔介质阻力使得棉线内部与棉线间隙之间的速度存在数量级上的差异。

图5 X=L/2截面处的速度分布Fig.5 Velocity profile at X=L/2

图6 反应界面（Z=0）上的燃料浓度分布Fig.6 Fuel concentration distribution on reaction interface(Z=0)

底部发生电化学反应时，由于棉线内部的速度比棉线间隙内的速度低，且燃料在棉线内部的扩散阻力大于棉线间隙处的，使得棉线与底部接触处的燃料传输比在棉线间隙内传输得慢。因此当底部发生电化学反应时，燃料传输快的地方强化了燃料传质，而传输慢的地方燃料无法及时得到补充，因此，在底部反应界面上燃料浓度呈现曲线式分布（图6）。

阳极过电位对燃料浓度分布的影响如图7 所示。从图中可以看出，在不同阳极过电位下，燃料浓度沿X方向均呈下降趋势。因为阳极过电位的大小反映了电化学极化发生的程度，所以阳极过电位较低时，电化学反应速率低，导致沿X方向上的燃料浓度降低的趋势较慢；阳极过电位较高时，电化学反应速率较高，燃料消耗快，沿X方向上燃料浓度下降的趋势更加明显。

图7 不同阳极过电位下沿X方向的平均燃料浓度分布Fig.7 Average fuel concentration distribution along with X direction at different anodic overpotentials

3.3 进口燃料浓度的影响

不同进口燃料浓度影响电池阳极燃料的传质，进而影响其阳极性能。在探究进口燃料浓度对电池阳极传质及性能的影响时，入口燃料流量保持为29.24 μl/min，分别选取进口浓度为2、1 和0.5 mol/L进行计算分析。计算得到在同一过电位、不同进口燃料浓度下在反应界面上沿流动方向的平均浓度分布，结果如图8所示。可以看出，在各个进口燃料浓度下，沿X方向上燃料浓度均下降，而且燃料浓度越高，浓度降低的趋势越快。在X方向的后半段，燃料浓度降低的速率均在减缓，且进口浓度越低时，其出现燃料浓度速率减缓时的位置越靠近入口。这主要是由于进口浓度越高，电化学反应速率越快，燃料消耗也就越快。

图8 反应界面（Z=0）上沿X方向的平均燃料浓度分布（η=0.49 V）Fig.8 Average fuel concentration distribution on reaction interface along with X direction（η=0.49 V）

图9 反应界面（Z=0）上沿X方向的平均电流密度分布（η=0.49 V）Fig.9 Average concentration distribution on reaction interface along with X direction（η=0.49 V）

图10 不同进口燃料浓度下阳极极化曲线Fig.10 Anode polarization curves under different inlet fuel concentrations

在同一过电位下沿流道X方向的平均电流密度分布如图9 所示。从图中可以看出，平均电流密度的变化趋势和图8 中燃料浓度分布一致，由Butler-Volmer 公式可知，一定过电位下，电化学反应速率与当地燃料浓度呈正比，因此，不同进口燃料浓度下平均电流密度沿X方向不断降低。此外，从图9还可以看出，当进口燃料浓度为0.5 mol/L时，沿X方向的平均电流密度变化不大，而进口燃料浓度越高，变化越明显，浓度为2 mol/L 时进出口的平均电流密度差值约为0.5 mol/L 时的4倍左右。不同进口燃料浓度下电池阳极极化曲线如图10 所示。从图10 可知，对不同进口燃料浓度，随着阳极电势升高，电流密度均在升高，符合阳极极化曲线的特征。在同一电位下，燃料传质强化使得阳极电流密度随着进口燃料浓度的增加而升高。

3.4 进口燃料流量的影响

图11 在反应界面（Z=0）上不同进口燃料流量下的燃料浓度分布Fig.11 Fuel concentration distribution on reaction interface(Z=0)under different inlet fuel flow rates

为研究燃料进口流量对电池阳极性能的影响，针对进口燃料浓度2 mol/L 时，进口燃料流量分别为29.24、15.90 和2.92 μl/min 进行计算，得到不同进口燃料流量下在底部反应界面上的浓度分布，如图11所示。从图中可以看出，流量越大时，棉线与反应界面接触部位的燃料浓度与其他区域浓度之间的差异越大；而流量越小，该浓度差异越小且流道后段的浓度越低。这主要是因为进口流量越大时，棉线内部多孔介质的阻力使得棉线内部与外部速度之间的差异较大，因此，发生电化学反应时，由于反应界面上流速低的地方其燃料补充不及时，进而导致此处燃料浓度与棉线间隙内浓度的差异较大[如图11（a）所示]。当进口流量较小时，由于棉线内外部速度差异较小，所以从图11（c）中可以看出在Y方向上的浓度差异并不明显。

图12 不同进口燃料流量下沿X方向的平均燃料浓度分布Fig.12 Average concentration distribution along with X direction at different inlet fuel flow rates(η=0.49 V)

图13 不同进口燃料流量下的阳极极化曲线Fig.13 Anode polarization curves at different inlet fuel flow rates

在反应界面上沿X方向的平均燃料浓度分布如图12 所示。可以看出，不同进口燃料流量下，燃料的平均浓度均沿X方向减小，在同一位置处，进口燃料流量越大，燃料浓度越高。这主要是由于进口燃料流量较小时，流道中的流速较慢导致反应界面上的燃料被消耗之后无法及时得到补充，而进口燃料流量较大时，对流传质强化使得反应界面上的燃料浓度增加。不同进口燃料流量下的阳极极化曲线如图13所示，从图中可以看出，当阳极电势较低时，由于反应速率较慢，所以不同进口燃料流量下电流密度相差不大，但当阳极电势高于-0.7 V时，阳极性能随着燃料流量的增加而升高，这是由于流量升高强化了对流传质，引起反应界面上燃料浓度升高，最终使得电流密度升高。

4 结 论

本文针对基于棉线的微流体燃料电池建立了阳极棉线流道内含电化学反应的理论计算模型，采用格子Boltzmann 方法研究了流道内燃料的流动及耦合电化学反应的传质特性，讨论了燃料的进口浓度及流量对电池阳极性能及传质特性的影响，得到的主要结论如下：

（1）通过格子Boltzmann方法计算获得的阳极极化曲线与实验结果吻合较好。

（2）进口燃料浓度相同时，不同阳极过电位下沿流动方向上的燃料浓度均呈下降趋势，且过电位越大，燃料浓度下降得越多。而进口燃料浓度越高时，反应界面上燃料浓度越高，电流密度及阳极性能越高。

（3）随着燃料进口流量的增加，阳极性能增强。当流量较高时，棉线与反应界面接触部位的浓度与其他区域之间的差异较大，而流量较低时，该浓度差异减小且流道后段的浓度较低。

符 号 说 明

A——可调参数

C——浓度，mol/L

Cref——参考浓度，mol/L

c——LBM中的格子速度

cs——LBM中的格子声速

D——扩散系数，m2/s

Deff——有效扩散系数，m2/s

ei——LBM中的离散速度

F——法拉第常数，C/mol

Fɛ——结构函数

fi——i方向的粒子速度分布函数

——i方向的粒子速度平衡分布函数

gi——i方向的粒子浓度分布函数

——i方向的粒子浓度平衡分布函数

H——模型高度，m

Iref——参考电流密度，mA/cm2

K——渗透率，m2

L——模型长度，m

P——压力，Pa

Q——流量，μl/min

R——通用气体常数，J/(mol·K)

T——温度,K

u——速度，m/s

W——模型宽度，m

αa——阳极传输系数

β——反应级数

δt——时间步

δx——网格宽度

ɛ——孔隙率

η——过电位，V

υ——运动黏性系数,m2/s

υe——有效运动黏性系数,m2/s

ρ——密度，kg/m3

ρ0——参考密度，kg/m3

τD——与扩散系数相关的松弛时间

τυ——与运动黏性系数相关的松弛时间

ωi——平衡分布函数中的权系数

下角标

in——进口

out——出口