謝 龍，王玉璋，于建國

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

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

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

1 模型建立

1.1 物理模型

圖1以恒定工作溫度下的SOFC多孔陽極為研究對象。電極模型的上下側(cè)取為周期性邊界；左側(cè)為燃料氣體濃度邊界(concentration boundary,CB)，燃料氣體通過多孔介質(zhì)流道擴散至右側(cè)三相界面(triple phase boundary,TPB)發(fā)生電化學反應；右側(cè)為燃料氣體消耗與產(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的多孔介質(zhì)電極中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ù)選取與多孔介質(zhì)重構(gòu)

LB模擬中參數(shù)為格子單位，根據(jù)文獻[11]中物理單位與格子單位的轉(zhuǎn)換規(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 側(cè)J=0.0415 mol/(m2·s)，燃料側(cè) CT=11.4 mol/m3，燃料氣體摩爾比為 H2∶H2O∶N2=0.47∶0.03∶0.5，計算得J*=0.067。物理單位與格子單位轉(zhuǎn)換時需注意：格子單位中孔隙率Φ與物理單位一致；總體網(wǎng)格規(guī)模必須足夠大才能使多孔結(jié)構(gòu)具有足夠的“分辨率”，而計算機有限的計算能力又要求網(wǎng)格規(guī)模足夠小，本文選取120×120網(wǎng)格以滿足兩方面的需求，即格子單位中特征長度L=120；當互擴散系數(shù)過小時，互擴散松弛時間接近于0.5 s，將會導致數(shù)值穩(wěn)定性問題，因此將計算中格子單位的Dij擴大103數(shù)量級；J和CT的取值可調(diào)整以滿足J*的要求，如J選定為0.069，則此工況下CT值取為1140 mol/m3。由于采用相似原理進行數(shù)值模擬，故引入無量綱電流密度i*以進行模擬結(jié)果的無量綱分析，同時根據(jù)Nernst方程得到濃差極化表達式：

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

最后將數(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ù)與實際一致。本文采用連通情況較好二維多孔介質(zhì)為研究對象，其內(nèi)部擴散過程以分子擴散為主，主要關(guān)注LB方法模擬多孔電極內(nèi)微尺度流動的可行性。為了充分模擬多孔電極內(nèi)的微尺度流動過程，進一步的研究除了引入上述的無量綱擴散系數(shù)比之外，還需要以X射線CT掃描等技術(shù)[13]重構(gòu)出的三維多孔介質(zhì)為研究對象。

2 模型計算與有效性分析

孔隙尺度LB模型的每步迭代計算都需要對網(wǎng)格格點類型進行邏輯判斷，并對流體格點與固體邊界格點各分布函數(shù)進行計算，因此計算量巨大。一種提高計算效率的方法是采用稀疏矩陣，分理出固體格點，計算只在流體格點上進行，不需要進行格點的邏輯判斷和搜索，從而減少了存儲空間并能提高計算效率[14]；另一種方法是利用LB本身的并行特性，通過并行計算來提高計算效率[15]。本文采用MPI(message passing interface)并行方法計算LB模型，并行算法設(shè)計時按一維行分解將計算區(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(jié)果與文獻所報道結(jié)果較為一致，說明LB微尺度模型可有效模擬陽極支撐SOFC的濃差極化。三者之間的差異可能是由于采用的多孔介質(zhì)以及對多孔介質(zhì)的描述方式不同，以及計算中未充分考慮Knudsen擴散的影響。

3 結(jié)果與分析

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

圖5為穩(wěn)態(tài)工況下多孔電極內(nèi)燃料氣體摩爾分數(shù)的典型分布以及系統(tǒng)流線圖。H2從燃料側(cè)擴散至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)各截面上保持恒定，滿足質(zhì)量守恒，因此各組分分布云圖中摩爾分數(shù)偏小的截面上或狹窄流道中流速較大。

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

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

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

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

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

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

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

表1 不同陽極結(jié)構(gòu)的濃差極化

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

4 結(jié)論

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

(2)SOFC存在極限電流密度，該閥值隨著多孔陽極孔隙率的增大而增加；多孔陽極濃差極化隨無量綱電流密度i*或燃料側(cè)H2摩爾比的變化作近似線性變化，而在TPB處的H2被完全消耗時變化劇烈；通過優(yōu)化陽極多孔結(jié)構(gò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].電源技術(shù)，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]王唯威，淮秀蘭.分形多孔介質(zhì)導熱數(shù)值模擬分析[J].工程熱物理學報，2007，28(5)：835-837.