基于棉線的微流體燃料電池陽極傳質特性LB模擬

穆嬡萍，葉丁丁，陳蓉，朱恂，廖強

（1 重慶大學低品位能源利用技術及系統(tǒng)教育部重點實驗室，重慶400030；2重慶大學能源與動力工程學院工程熱物理研究所，重慶400030）

引 言

無膜微流體燃料電池因不存在由質子交換膜引起的與膜相關的一系列問題[1-3]，使其成為極具應用前景的新型微型電源之一。而傳統(tǒng)的無膜微流體燃料電池常需外部驅動設備來維持燃料和氧化劑的層流流動，但該設備不僅會增加整個運行系統(tǒng)的體積，還會降低電池的凈輸出功率，限制其進一步向微型化發(fā)展。常見的可替代外部驅動設備的驅動方式有：流體的表面張力驅動[4]、重力驅動[5-6]、毛細力驅動[7-9]、化學反應驅動[10]等。和其他驅動方式相比，采用紙、纖維等多孔材料的毛細力驅動方式有著結構簡單、制造成本低、占用空間小等優(yōu)點[8]。但是當紙和纖維的孔隙被溶液充滿后，其毛細驅動力減弱，難以形成持續(xù)穩(wěn)定的流動，因此，在微流體裝置中常采用與重力驅動相結合的方式來驅動流體以形成穩(wěn)定流動[11-12]。因棉線有著天然三維的多孔結構、較好的柔韌性以及在浸潤情況下具有較好的機械強度等優(yōu)點，2018 年，Liu 等[13]構建了基于棉線并輔以重力驅動的微流體燃料電池，其最大功率密度和電流密度比文獻[14]報道的基于織物類的微流體燃料電池的性能高3.7 和1.5 倍。2020 年，Wu等[15]直接將棉線平鋪于電極表面作為流體流道，構建了基于棉線的微流體燃料電池，在堿性條件下電池的極限電流密度在同類電池中最高，而且該類電池的性能嚴重受到陽極側燃料傳質的限制，但目前對基于多孔纖維類的微流體燃料電池的相關理論研究工作還未見報道。因此，采用數(shù)值模擬的方法來研究該類燃料電池內(nèi)部的傳輸機理將為電池性能的提升提供理論指導。

與其他數(shù)值模擬方法相比，格子Boltzmann方法（lattice Boltzmann method, LBM）能有效處理復雜的微觀幾何結構以及多孔介質中多相流體流動等方面的問題[16]。雖然近年來許多學者通過基于孔隙尺度的LBM 來研究燃料電池催化層或擴散層內(nèi)的流動及傳質問題[17-20]，但該方法對每個孔隙內(nèi)的流動進行精確處理的要求限制了其模型計算域的尺度。與孔隙尺度不同，表征體元（representative elementary volume, REV）尺度是通過適當?shù)淖饔昧砻枋隽黧w與多孔介質之間的相互作用進而獲取體積平均的宏觀流動量，并不需要獲取多孔介質的結構，因此常有學者使用該方法來研究復雜或無法直接獲取的多孔介質結構內(nèi)的流動特性[21-23]。而在采用多孔纖維的微流體燃料電池中，輸運溶液的流道是由具有多尺度孔隙結構的棉線堆疊而成，具有復雜的多孔介質結構。因此本文采用LB 方法研究REV 尺度下基于棉線的微流體燃料電池陽極內(nèi)流動及耦合電化學反應的傳質特性。

本文首先介紹計算采用的物理模型；隨后介紹LB 數(shù)值模擬方法并進行LBM 模型驗證，證明該流動及傳質模型的正確性；然后研究受重力驅動的基于棉線的微流體燃料電池陽極側流動及傳質特性，計算獲得燃料濃度分布及陽極極化曲線，并討論不同燃料進口濃度及流量對電池陽極性能及傳質特性的影響。

1 物理模型

根據(jù)前期構建的基于棉線的微流體燃料電池[15]建立物理模型。實驗中采用了醫(yī)用親水紗布上的棉線作為可輸運陽極和陰極溶液的流道，通過控制溶液進出口高度差形成的重力維持溶液的連續(xù)流動。在該電池結構中，陽極和陰極流道分別用10根棉線平鋪堆疊于負載有催化劑（有效面積為1 mm×10 mm）的碳紙電極上，為防止陰、陽極溶液發(fā)生擴散混合，在二者中間加入由40 根棉線平鋪堆疊而成的隔離液通道。在電池運行過程中，陽極液（甲酸鈉與氫氧化鉀的混合溶液）和陰極液（含飽和溶解氧的氫氧化鉀溶液）在重力驅動下輸運到負載有催化劑的碳紙電極上發(fā)生電化學反應，陽極甲酸根與氫氧根反應生成碳酸根、水和電子，陰極在堿性條件下反應生成氫氧根，氫氧根經(jīng)隔離液到達陽極，電子經(jīng)外電路到達陰極，形成回路產(chǎn)生電流。

基于棉線的微流體燃料電池運行初期，由棉線多孔介質的毛細力和重力驅動輸運溶液，當多孔介質被溶液充滿后，電池靠重力的作用持續(xù)供液，使其穩(wěn)定運行。溶液在棉線流道內(nèi)輸運的過程中，其多孔流道主要有三個尺度的孔隙結構：一是棉線與棉線之間的間隙，二是單根棉線中棉線纖維之間的間隙，三是單根棉線纖維內(nèi)部的間隙[15]。因此，本文在理論計算模型中采用的基本假設如下：

（1）微流體燃料電池在常溫常壓下運行，燃料溶液為不可壓縮流體，且在多孔流道內(nèi)為等溫穩(wěn)態(tài)流動。

（2）不考慮棉線未被溶液充滿的不穩(wěn)定階段，即基于棉線的微流體燃料電池在重力驅動的作用下穩(wěn)定運行。

（3）單根棉線內(nèi)部結構復雜且難以進行三維重構，因此，忽略單根棉線內(nèi)部復雜的多孔結構，假設單根棉線為具有一定孔隙率?的圓柱。

（4）催化層的厚度僅為幾十微米，本文主要考慮棉線流道內(nèi)的流動和傳質，不考慮催化層內(nèi)部的傳輸，即將催化層考慮為反應界面。

模型計算區(qū)域如圖1（a）所示，選取棉線組成的基本單元進行計算，考慮陽極電極布置在計算區(qū)域底部位置，棉線與棉線之間為等邊三角形排列，如圖1（b）所示，圖中綠色表示棉線，其余為棉線之間的間隙。

圖1 模型計算區(qū)域示意圖（a）及棉線堆疊方式示意圖（b）Fig.1 Schematics of computational domain(a)and cotton thread array(b)

甲酸鈉溶液在陽極電極處發(fā)生電化學反應

其物質守恒方程為

式中，u為速度；C和Deff分別為物質濃度和有效擴散系數(shù)。在該計算模型中，溶液在棉線內(nèi)部的傳質阻力要比棉線與棉線間間隙內(nèi)的傳質阻力大，采用該方程描述溶液在棉線與棉線之間的傳質時，其擴散系數(shù)為D，而溶液在單根棉線內(nèi)部的擴散系數(shù)則修正為Deff=?1.5D。其反應源項Rs=-Ilocal/2F，局部電流Ilocal采用Butler-Volmer方程描述為

其中，Iref為參考電流密度；C為反應界面上的當?shù)貪舛龋籆ref為參考濃度；β為反應級數(shù)；αa為陽極傳輸系數(shù)；F為法拉第常數(shù)；R為通用氣體常數(shù)；T為溫度；η為過電勢。

2 格子Boltzmann方法

2.1 流動模型

不可壓縮流體在多孔介質內(nèi)流動的控制方程為[24]：

式中，u、P為表觀速度和壓力；ρ0為參考密度；υe為有效運動黏性系數(shù)；F為包括了介質阻力及其他外力的總體積力，即

式中，K、υ分別為滲透率和流體運動黏性系數(shù)，與平均孔隙率?相關的結構函數(shù)F?可用式（6）表示[25]

低速流動時，總體積力F[式（5）]中包含結構函數(shù)F?的非線性介質阻力項可以忽略。

對于式（4）不可壓縮流體在多孔介質中的等溫流動模型，Guo 和Zhao[26]于2002 年提出求解該問題的廣義LBE模型（GLBE），演化方程為

其中,

在本文采用的D3Q19 的流動模型參數(shù)中的權重因子ωi及各方向粒子離散速度ei為

因速度定義項的外力F包含速度，根據(jù)Guo 和Zhao[26]提出的臨界速度求解得出

式中c0與c1為

對式（7）做Chapman-Enskog分析可知流體運動黏性系數(shù)與松弛時間的關系為

當?=1 時，該GLBE 模型演變?yōu)橛糜谇蠼庾杂闪鲃訁^(qū)域速度分布情況的LBE模型。

2.2 傳質模型

在LBM 傳質模型中，當物質的擴散系數(shù)D與流動模型中的運動黏性系數(shù)υ相差較大時，常用的格子Bhatnagar-Gross-Krook（LBGK）模型會出現(xiàn)數(shù)值不穩(wěn)定現(xiàn)象。Inamuro[27]提出在平衡分布函數(shù)中添加修正項的格子動理格式（lattice kinetic scheme，LKS）模型中，通過調(diào)節(jié)該參數(shù)可使松弛因子穩(wěn)定在合理范圍內(nèi)，改善數(shù)值不穩(wěn)定現(xiàn)象。基于該模型，Yang 等[28]提出用于求解對流擴散方程的修正LBGK模型，其演化方程為[28-29]

其時間導數(shù)可由向后差分獲得

在該模型中宏觀參數(shù)C為

對式（15）做Chapman-Enskog 分析可知擴散系數(shù)與松弛時間的關系及濃度梯度?C為

為保證數(shù)值穩(wěn)定性，本文中τυ及τD均設為1，則可調(diào)參數(shù)A可根據(jù)τD及擴散系數(shù)D來確定。

邊界條件設定：進口（X=0）處已知壓力Pin和濃度Cin；出口（X=L）為已知出口壓力Pout，在該邊界處濃度為充分發(fā)展邊界；左右（Y=0，Y=W）為周期性邊界；上邊界（Z=H）處速度u= 0，濃度梯度為0；下邊界（Z=0）處速度u=0 且在該界面發(fā)生電化學反應時，采用Kang等[30]報道的反應邊界方法來求解。

2.3 LB模型驗證

采用兩平板間充滿多孔介質的流體流動來驗證REV-LBM 模型的正確性。流體在流道內(nèi)的速度分布表述為[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為沿中心線上的濃度分布，可以看出當可調(diào)參數(shù)A=0 和A≠0 時的數(shù)值解與其解析解[31]吻合較好，可證明該模型的正確性。

圖3 可調(diào)參數(shù)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所示，其計算參數(shù)取值見表1。從圖4 可知，在電位較低時，理論結果與實驗結果有一定差異，這主要因為電流密度較小時，Butler-Volmer公式不能準確描述此時的活化極化損失。但整體上計算結果與實驗結果吻合較好，證明該模型及參數(shù)的適用性。

表1 計算參數(shù)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 流道內(nèi)的速度與濃度分布

通過計算獲得流道內(nèi)的速度分布如圖5 所示。由于不考慮單根棉線內(nèi)部的多孔結構，僅求解REV尺度下體積平均的宏觀流量，所以棉線內(nèi)部流速分布均勻，且由于棉線內(nèi)部多孔介質阻力使得棉線內(nèi)部與棉線間隙之間的速度存在數(shù)量級上的差異。

圖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)

底部發(fā)生電化學反應時，由于棉線內(nèi)部的速度比棉線間隙內(nèi)的速度低，且燃料在棉線內(nèi)部的擴散阻力大于棉線間隙處的，使得棉線與底部接觸處的燃料傳輸比在棉線間隙內(nèi)傳輸?shù)寐Ｒ虼水數(shù)撞堪l(fā)生電化學反應時，燃料傳輸快的地方強化了燃料傳質，而傳輸慢的地方燃料無法及時得到補充，因此，在底部反應界面上燃料濃度呈現(xiàn)曲線式分布（圖6）。

陽極過電位對燃料濃度分布的影響如圖7 所示。從圖中可以看出，在不同陽極過電位下，燃料濃度沿X方向均呈下降趨勢。因為陽極過電位的大小反映了電化學極化發(fā)生的程度，所以陽極過電位較低時，電化學反應速率低，導致沿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方向的后半段，燃料濃度降低的速率均在減緩，且進口濃度越低時，其出現(xiàn)燃料濃度速率減緩時的位置越靠近入口。這主要是由于進口濃度越高，電化學反應速率越快，燃料消耗也就越快。

圖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 公式可知，一定過電位下，電化學反應速率與當?shù)厝剂蠞舛瘸收龋虼耍煌M口燃料濃度下平均電流密度沿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所示。從圖中可以看出，流量越大時，棉線與反應界面接觸部位的燃料濃度與其他區(qū)域濃度之間的差異越大；而流量越小，該濃度差異越小且流道后段的濃度越低。這主要是因為進口流量越大時，棉線內(nèi)部多孔介質的阻力使得棉線內(nèi)部與外部速度之間的差異較大，因此，發(fā)生電化學反應時，由于反應界面上流速低的地方其燃料補充不及時，進而導致此處燃料濃度與棉線間隙內(nèi)濃度的差異較大[如圖11（a）所示]。當進口流量較小時，由于棉線內(nèi)外部速度差異較小，所以從圖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 結 論

本文針對基于棉線的微流體燃料電池建立了陽極棉線流道內(nèi)含電化學反應的理論計算模型，采用格子Boltzmann 方法研究了流道內(nèi)燃料的流動及耦合電化學反應的傳質特性，討論了燃料的進口濃度及流量對電池陽極性能及傳質特性的影響，得到的主要結論如下：

（1）通過格子Boltzmann方法計算獲得的陽極極化曲線與實驗結果吻合較好。

（2）進口燃料濃度相同時，不同陽極過電位下沿流動方向上的燃料濃度均呈下降趨勢，且過電位越大，燃料濃度下降得越多。而進口燃料濃度越高時，反應界面上燃料濃度越高，電流密度及陽極性能越高。

（3）隨著燃料進口流量的增加，陽極性能增強。當流量較高時，棉線與反應界面接觸部位的濃度與其他區(qū)域之間的差異較大，而流量較低時，該濃度差異減小且流道后段的濃度較低。

符 號 說 明

A——可調(diào)參數(shù)

C——濃度，mol/L

Cref——參考濃度，mol/L

c——LBM中的格子速度

cs——LBM中的格子聲速

D——擴散系數(shù)，m2/s

Deff——有效擴散系數(shù)，m2/s

ei——LBM中的離散速度

F——法拉第常數(shù)，C/mol

F?——結構函數(shù)

fi——i方向的粒子速度分布函數(shù)

——i方向的粒子速度平衡分布函數(shù)

gi——i方向的粒子濃度分布函數(shù)

——i方向的粒子濃度平衡分布函數(shù)

H——模型高度，m

Iref——參考電流密度，mA/cm2

K——滲透率，m2

L——模型長度，m

P——壓力，Pa

Q——流量，μl/min

R——通用氣體常數(shù)，J/(mol·K)

T——溫度,K

u——速度，m/s

W——模型寬度，m

αa——陽極傳輸系數(shù)

β——反應級數(shù)

δt——時間步

δx——網(wǎng)格寬度

?——孔隙率

η——過電位，V

υ——運動黏性系數(shù),m2/s

υe——有效運動黏性系數(shù),m2/s

ρ——密度，kg/m3

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

τD——與擴散系數(shù)相關的松弛時間

τυ——與運動黏性系數(shù)相關的松弛時間

ωi——平衡分布函數(shù)中的權系數(shù)

下角標

in——進口

out——出口