SOFC多孔陽極內微尺度流動的Lattice-Boltzmann模擬

謝 龍，王玉璋，于建國

（上海交通大學動力機械與工程教育部重點實驗室，上海 200240）

固體氧化物燃料電池 (SOFC)一次發電效率高，易實現SOFC-GT聯合循環，電極材料成本低，燃料適用面廣，全固態結構適合進行模塊化設計和放大，是中長期最具吸引力的新型、綠色能源解決方案[1]。SOFC的一大發展瓶頸是較高的操作溫度要求電池的關鍵部件必須具有一定的化學和熱的相容性，局部熱不平衡引起的極大熱應力易導致微米級電解質層的密封性問題甚至斷裂，使得單電池難以做大[2]。目前主要采用具有低電池運行溫度的陽極支撐結構以改善SOFC工作溫度高這一技術難點。而在高電流密度或燃料利用率時，陽極支撐SOFC的性能主要受到陽極側濃度過電勢控制[3-5]。陽極中燃料氣體以擴散形式供應給電解質層，反應物氣體的擴散速率與氣體成分和陽極微結構有關[6-7]。已有的對電極微結構研究的局限性在于依據“容積平均”的概念將多孔介質看成一種在大尺度上均勻分布的虛擬連續介質，并采用難以測量（比表面積、迂曲度等）或難以承擔特性尺度角色（孔隙率、平均孔隙尺度等）的幾何參量來表征其幾何結構特性，難以模擬多孔電極內真實的微尺度流動。

本文建立了多孔電極內氣體輸運的Lattice-Boltzmann(LB)模型，結合多孔介質的數值構造方法，模擬陽極支撐SOFC的微尺度流動過程。在模型有效性驗證基礎上，進一步分析了不同孔隙率、電流密度、燃料組分及陽極結構對多孔電極內質量傳遞及濃差極化的影響。模型為多孔電極內氣體輸運過程的微尺度研究提供了思路，對于制備出性能優良、具有高輸出功率密度的電池具有重要的理論及實際工程意義。

1 模型建立

1.1 物理模型

圖1以恒定工作溫度下的SOFC多孔陽極為研究對象。電極模型的上下側取為周期性邊界；左側為燃料氣體濃度邊界(concentration boundary,CB)，燃料氣體通過多孔介質流道擴散至右側三相界面(triple phase boundary,TPB)發生電化學反應；右側為燃料氣體消耗與產物生成的流量邊界，固體顆粒表面處理為法向零流量、切向無滑移速度邊界。本文燃料選定為H2、H2O與N2的混合氣以模擬微尺度多組分流動系統。

1.2 Lattice-Boltzmann模型

LB方法是一種基于統計演化分布函數、研究粒子團碰撞和遷移的介觀方法，模擬適用范圍從Kn數小于10－3的宏觀連續流到101數量級的自由分子區內流動。SOFC的多孔介質電極中Kn≈1，發生分子擴散和Knudsen擴散[9]，同時LB方法具有獨特的邊界處理方式，使得LB方法對SOFC多孔電極的模擬具有極好的適應性。本文建立的二維等溫多組分LB模型由McCracken[10]的不同分子量氣體兩組分擴散DLS(different lattice speed)模型發展而來。模型演化方程為：

其中，從fαi(x,t)到為演化過程的遷移步，分布函數fαi遷移（或遷移后插值）到臨近格點；為演化過程的碰撞步，分為同組分氣體內的自我碰撞和異組分氣體間的交叉碰撞。

濃度邊界與流量邊界采用Zou-He格式，固體顆粒表面邊界采用非平衡反彈格式[11]。LB模型中的分子擴散由孔隙內分子間的自我碰撞和交叉碰撞體現，而Knudsen擴散則是由分子與固體顆粒的碰撞反映。

1.3 參數選取與多孔介質重構

LB模擬中參數為格子單位，根據文獻[11]中物理單位與格子單位的轉換規則，模擬1 mm厚的平板式SOFC多孔陽極（孔隙尺度為微米級）時，其厚度方向網格規模將達到103數量級。而根據相似原理，可對計算參數做無量綱處理，通過保證無量綱準則數的一致性來模擬真實流動，從而減小計算規模。選取無量綱準則數為無量綱流量J*=JL/(CTD13)，互擴散系數比D12/D23、D13/D23和孔隙率Φ。研究對象選取為1073 K、0.1 MPa環境下SOFC多孔陽極中H2、H2O與N2多組分理想氣體流動系統，物理參數為D12=3.37×10－4m2/s，D13=0.692×10－4m2/s，D23=1.085×10－4m2/s[12]，Φ=0.5，L=0.002 m，TPB 側J=0.0415 mol/(m2·s)，燃料側 CT=11.4 mol/m3，燃料氣體摩爾比為 H2∶H2O∶N2=0.47∶0.03∶0.5，計算得J*=0.067。物理單位與格子單位轉換時需注意：格子單位中孔隙率Φ與物理單位一致；總體網格規模必須足夠大才能使多孔結構具有足夠的“分辨率”，而計算機有限的計算能力又要求網格規模足夠小，本文選取120×120網格以滿足兩方面的需求，即格子單位中特征長度L=120；當互擴散系數過小時，互擴散松弛時間接近于0.5 s，將會導致數值穩定性問題，因此將計算中格子單位的Dij擴大103數量級；J和CT的取值可調整以滿足J*的要求，如J選定為0.069，則此工況下CT值取為1140 mol/m3。由于采用相似原理進行數值模擬，故引入無量綱電流密度i*以進行模擬結果的無量綱分析，同時根據Nernst方程得到濃差極化表達式：

模擬時將多孔介質信息以數字化網格形式讀入。首先以二值圖形式建立多孔介質，然后將孔隙結構二值圖轉換為數值0、1形式的數字化網格，即通過標志函數來表示多孔介質的空間分布：

最后將數值矩陣作為數值模擬的網格輸入。

需要注意的是，為了實現對多孔電極微尺度流動的完全模擬，還需要引入表征Knudsen擴散與分子擴散比值的無量綱參數D1,k/D1,3（其中D1,k為組分1的Knudsen擴散系數）來保證數值模擬與真實流動的完全相似。一種方法是保持現有網格形式不變（即固定D1,k），通過增大D1,3（即弱化分子間的碰撞，減小連續區域效應）來間接模擬實際多孔電極內的非連續區域效應；另一種方法是修正數值網格（通常是減少孔隙間格子數）以保證網格的Knudsen數與實際一致。本文采用連通情況較好二維多孔介質為研究對象，其內部擴散過程以分子擴散為主，主要關注LB方法模擬多孔電極內微尺度流動的可行性。為了充分模擬多孔電極內的微尺度流動過程，進一步的研究除了引入上述的無量綱擴散系數比之外，還需要以X射線CT掃描等技術[13]重構出的三維多孔介質為研究對象。

2 模型計算與有效性分析

孔隙尺度LB模型的每步迭代計算都需要對網格格點類型進行邏輯判斷，并對流體格點與固體邊界格點各分布函數進行計算，因此計算量巨大。一種提高計算效率的方法是采用稀疏矩陣，分理出固體格點，計算只在流體格點上進行，不需要進行格點的邏輯判斷和搜索，從而減少了存儲空間并能提高計算效率[14]；另一種方法是利用LB本身的并行特性，通過并行計算來提高計算效率[15]。本文采用MPI(message passing interface)并行方法計算LB模型，并行算法設計時按一維行分解將計算區域分配給各個處理器，如圖2。值得注意的是，若按一維列分解或二維塊分解劃分計算區域，臨近分區互相決定了對方的邊界條件，不利于數值穩定性，影響計算精度。采用8核、單核2.5 GHz工作站進行并行計算，圖3給出了處理器數量與加速比（串行運算時間與并行運算時間之比），可見隨著處理器數的增加，加速效果趨于穩定。

為驗證LB多組分流動模型的正確性，采用Yakabe[5]的陽極支撐SOFC性能試驗數據進行模擬計算，同時與Suwanwarangkul[16]采用的Dusty-Gas Model(DGM)模型預測進行比較，如圖4?？梢姳灸Ｐ陀嬎憬Y果與文獻所報道結果較為一致，說明LB微尺度模型可有效模擬陽極支撐SOFC的濃差極化。三者之間的差異可能是由于采用的多孔介質以及對多孔介質的描述方式不同，以及計算中未充分考慮Knudsen擴散的影響。

3 結果與分析

圖5 穩態工況下多孔陽極內參數分布

圖5為穩態工況下多孔電極內燃料氣體摩爾分數的典型分布以及系統流線圖。H2從燃料側擴散至TPB參與反應，TPB上H2不斷消耗生成H2O，H2O則從TPB擴散至燃料流道；N2在TPB上不參與電化學反應，在多孔電極內表現為梯度較小的濃度差分布；電極的總速度分布則是從TPB指向燃料流道，電極內的電化學反應熱與歐姆熱不斷被帶出電極，以達到良好的冷卻效果，避免局部過熱區的產生。經計算發現各組分流量在電極內各截面上保持恒定，滿足質量守恒，因此各組分分布云圖中摩爾分數偏小的截面上或狹窄流道中流速較大。

3.1 孔隙率、電流密度對質量傳遞和濃差極化的影響

圖6給出了不同孔隙率與電流密度時TPB上H2的平均摩爾分數分布，可見固定J*時，隨著孔隙率的增大，電極內氣體擴散通道增加，維持一定摩爾流量所需的濃度差減小；每個孔隙率對應著一定的J*閥值（或極限電流密度i*），該流量閥值表示擴散到TPB處的H2被完全消耗，該閥值隨著孔隙率的增大而增加；固定孔隙率時，濃度差隨著J*的增加而增大以滿足大摩爾流量的需求。

圖7給出了不同孔隙率與電流密度情況下多孔陽極的濃差極化特性曲線，可見固定孔隙率時，濃差極化隨著i*的增大而增大，固定i*時，濃差極化隨著孔隙率的減小而增大。原因是隨著孔隙率的減小，H2難以擴散到TPB參與反應且電化學反應生成的H2O難以往外擴散；或者隨著i*的增大，TPB上電化學反應速度加快從而H2大量消耗且H2O大量生成，因此TPB上H2的摩爾濃度減小且H2O濃度增大，從而引起較大的極化損失。同時還觀察到濃差極化在很大i*范圍間呈近似線性分布，而當i*接近該孔隙率對應的極限電流密度時濃差極化出現較大的階躍，i*恒定時濃差極化隨孔隙率減小而增大的趨勢也越發明顯。因此高i*或低孔隙率情況下濃差極化的變化更加敏感。

3.2 氣體組分對濃差極化的影響

圖8給出了孔隙率為0.45時，固定燃料氣體中N2摩爾分數為50%，變化H2與H2O摩爾比時各i*對應的濃差極化曲線，可見濃差極化隨燃料側H2摩爾比的增大而減小，且除了大的i*或小的H2摩爾體積時表現出的濃差極化陡增以外，其余區域為線性變化。原因是當i*、孔隙率以及N2摩爾分數一定時，多孔電極兩側H2濃度差保持恒定；電極內N2純擴散不參與反應，在TPB上H2與H2O摩爾濃度之和保持恒定。因此隨著燃料側H2摩爾分數的降低，TPB處H2摩爾分數線性減小，而H2O則線性增加，從而引起濃差極化的近似線性變化。

圖9給出了孔隙率為0.45時，固定燃料氣體中H2與H2O摩爾比為4∶1，變化H2濃度時各i*對應的濃差極化曲線。i*、孔隙率以及H2/H2O摩爾比一定時，多孔電極兩側H2濃度差隨H2摩爾分數的增加而線性減?。欢S著N2摩爾分數的減小，TPB上H2與H2O摩爾濃度之和則線性增加。因此燃料側H2摩爾分數降低會引起TPB處H2摩爾分數的線性減小以及H2O摩爾分數的線性增加，從而引起濃差極化在大H2摩爾體積時的近似線性變化以及小H2摩爾體積時的陡增。

圖9 變H2濃度時的濃差極化曲線

3.3 陽極結構對濃差極化的影響及電極結構優化

表1給出了i*=0.064，孔隙率為0.45時不同陽極結構對濃差極化的影響。其中多孔電極結構有二級規則分形Sierpinski地毯結構多孔介質[17]，以等半徑圓為基本構造單位的順排、錯排和隨機多孔介質，以及真實多孔介質。

順排、錯排與隨機三種多孔陽極結構的濃差極化值比較接近說明不同電極結構的主要特性尺度為孔隙率，而濃差極化的細微差異則是由不同迂曲度引起，說明迂曲度為多孔電極的次要特性尺度。分形與真實多孔介質得到了相同的濃差極化值，則說明真實多孔介質具有分形特性。各種形式陽極結構在極化特性上的不同表現說明對多孔結構的合理優化可以減弱電極的濃差極化特性。

表1 不同陽極結構的濃差極化

基于此，對多孔電極結構的優化路線應該從減少濃差極化出發，充分考慮歐姆極化以及活化極化的影響，使極化損失之和達到最小值。i*一定時，孔隙率越大濃差極化越小，而較大的孔隙率會導致電極的結構性能下降，結構整體的歐姆極化損失增加；選擇合適孔隙率的同時，多孔結構應滿足燃料附近孔隙率較大而TPB附近孔隙率較小，既使得燃料氣體可以順利的往TPB擴散，又保證TPB附近足夠的有效反應面積從而減小活化極化損失。

4 結論

(1)LB方法適用于燃料電池多孔電極內多組分氣體輸運過程的微尺度模擬，通過LB模型可以建立多孔陽極微結構及燃料組分與質量傳遞和濃差極化的關系。進一步耦合電極內傳熱過程以及TPB上電化學反應的LB模型將是設計與分析SOFC多孔陽極更為準確有效的工具。

(2)SOFC存在極限電流密度，該閥值隨著多孔陽極孔隙率的增大而增加；多孔陽極濃差極化隨無量綱電流密度i*或燃料側H2摩爾比的變化作近似線性變化，而在TPB處的H2被完全消耗時變化劇烈；通過優化陽極多孔結構可以減小濃差極化，以提升單電池性能。

[1]WANG Y Z,YOSHIBA F,WATANABE T,et al.Numerical analysis of electrochemical characteristics and heat/species transport for planar porous-electrode-supported SOFC[J]. Journal of Power Sources,2007,170(1):101-110.

[2]毛宗強，黃建兵，王誠，等.低溫固體氧化物燃料電池研究進展[J].電源技術，2008，32(2)：75-79.

[3]王玉璋，于建國，惠宇，等.新型平板式固體氧化物燃料電池的開發和性能分析[J].動力工程，2008，28(6)：949-954.

[4]YI J,ANIL V V.Fuel composition and diluent effect on gas transport and performance of anode-supported SOFCs[J].Journal of theElectrochemical Society,2003,150(7):942-951.

[5]YAKABE H,HISHINUMA M,URATANI M,et al.Evaluation and modeling of performance of anode-supported solid oxide fuel cell[J].Journal of Power Sources,2000,86(1/2):423-431.

[6]VIRKAR A V,CHEN J,TANNER C W,et al.The role of electrode microstructure on activation and concentration polarizations in solid oxide fuel cells[J].Solid State Ionics,2000,131(1/2):189-198.

[7]ZHAO F,VIRKAR A V.Dependence of polarization in anode-supported solid oxide fuel cells on various cell parameters[J].Journal of Power Sources,2005,141(1):79-95.

[8]BIRD R B,STEWART W E,LIGHTFOOT E N.Transport phenomena 2nd EDN[M].New York:Wiley,2002.

[9]于建國，王玉璋，翁史烈.以煤氣化合成氣為燃料的平板式固體氧化物燃料電池性能[J].中國電機工程學報，2011，30：88-93.

[10]MCCRACKEN M E,ABRAHAM J.Lattice Boltzmann methods for binary mixtures with different molecular weights[J].Phys Rev E,2005,71(4):8.

[11]何雅玲，王勇，李慶，等.格子Boltzmann方法的理論及應用[M].北京：科學出版社，2009.

[12]JOSHI A S,PERACCHIO A A,GREW K N,et al.Lattice Boltzmann method for multi-component,non-continuum mass diffusion[J].Journal of Physics D:Applied Physics,2007,40(23):7593-7600.

[13]LEFEBVRE,PHILIPPE L,BANHART,et al.Porous metals and metallic foams:current status and recent developments[J].Advanced Engineering Materials,2008,10(9):775-787.

[14]郭照立，鄭楚光.格子Boltzmann方法的原理及應用[M].北京：科學出版社，2009.

[15]VIDALA D,ROYA R,BERTRAND F.On improving the performance of large parallel Lattice Boltzmann flow simulations in heterogeneous porous media[J].Computers&Fluids,2010,39(2):324-337.

[16]SUWANWARANGKUL R,CROISET E,FOWLER M W,et al.Performance comparison of Fick's,dusty-gas and Stefan-Maxwell models to predict the concentration overpotential of a SOFC anode[J].Journal of Power Sources,2003,122(1):9-18.

[17]王唯威，淮秀蘭.分形多孔介質導熱數值模擬分析[J].工程熱物理學報，2007，28(5)：835-837.