SOFC多孔陽(yáng)極內(nèi)微尺度流動(dòng)的Lattice-Boltzmann模擬

謝 龍，王玉璋，于建國(guó)

（上海交通大學(xué)動(dòng)力機(jī)械與工程教育部重點(diǎn)實(shí)驗(yàn)室，上海 200240）

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

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

1 模型建立

1.1 物理模型

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

1.2 Lattice-Boltzmann模型

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

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

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

1.3 參數(shù)選取與多孔介質(zhì)重構(gòu)

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

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

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

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

2 模型計(jì)算與有效性分析

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

為驗(yàn)證LB多組分流動(dòng)模型的正確性，采用Yakabe[5]的陽(yáng)極支撐SOFC性能試驗(yàn)數(shù)據(jù)進(jìn)行模擬計(jì)算，同時(shí)與Suwanwarangkul[16]采用的Dusty-Gas Model(DGM)模型預(yù)測(cè)進(jìn)行比較，如圖4?？梢?jiàn)本模型計(jì)算結(jié)果與文獻(xiàn)所報(bào)道結(jié)果較為一致，說(shuō)明LB微尺度模型可有效模擬陽(yáng)極支撐SOFC的濃差極化。三者之間的差異可能是由于采用的多孔介質(zhì)以及對(duì)多孔介質(zhì)的描述方式不同，以及計(jì)算中未充分考慮Knudsen擴(kuò)散的影響。

3 結(jié)果與分析

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

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

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

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

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

3.2 氣體組分對(duì)濃差極化的影響

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

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

圖9 變H2濃度時(shí)的濃差極化曲線

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

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

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

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

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

4 結(jié)論

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

(2)SOFC存在極限電流密度，該閥值隨著多孔陽(yáng)極孔隙率的增大而增加；多孔陽(yáng)極濃差極化隨無(wú)量綱電流密度i*或燃料側(cè)H2摩爾比的變化作近似線性變化，而在TPB處的H2被完全消耗時(shí)變化劇烈；通過(guò)優(yōu)化陽(yáng)極多孔結(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]毛宗強(qiáng)，黃建兵，王誠(chéng)，等.低溫固體氧化物燃料電池研究進(jìn)展[J].電源技術(shù)，2008，32(2)：75-79.

[3]王玉璋，于建國(guó)，惠宇，等.新型平板式固體氧化物燃料電池的開(kāi)發(fā)和性能分析[J].動(dòng)力工程，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]于建國(guó)，王玉璋，翁史烈.以煤氣化合成氣為燃料的平板式固體氧化物燃料電池性能[J].中國(guó)電機(jī)工程學(xué)報(bào)，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方法的理論及應(yīng)用[M].北京：科學(xué)出版社，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方法的原理及應(yīng)用[M].北京：科學(xué)出版社，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ì)導(dǎo)熱數(shù)值模擬分析[J].工程熱物理學(xué)報(bào)，2007，28(5)：835-837.