謝 龍，王玉璋，于建國

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

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

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

1 模型建立

1.1 物理模型

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

1.2 Lattice-Boltzmann模型

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

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

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

1.3 參數(shù)選取與多孔介質重構

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

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

最后將數(shù)值矩陣作為數(shù)值模擬的網(wǎng)格輸入。

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

2 模型計算與有效性分析

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

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

3 結果與分析

圖5 穩(wěn)態(tài)工況下多孔陽極內(nèi)參數(shù)分布

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

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

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

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

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

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

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

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

3.3 陽極結構對濃差極化的影響及電極結構優(yōu)化

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

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

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

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

4 結論

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

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

[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]王玉璋，于建國，惠宇，等.新型平板式固體氧化物燃料電池的開發(fā)和性能分析[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]王唯威，淮秀蘭.分形多孔介質導熱數(shù)值模擬分析[J].工程熱物理學報，2007，28(5)：835-837.