基于多尺度耦合的管式SOFC數(shù)值模擬

馬 旭，楊 晨

(1. 珠海城市職業(yè)技術(shù)學(xué)院，廣東 珠海 519090；2. 重慶大學(xué)能源與動力工程學(xué)院，重慶 400044)

1 引言

建模與仿真的發(fā)展趨勢是從計算機宏觀計算模擬到粒子團的模擬再到分子級的模擬，即從宏觀到介觀再到微觀。在宏觀層次上主要是依據(jù)三大守恒定律進行的CFD模擬；在介觀層次上近年來研究最多的是基于粒子團的格子玻爾茲曼方法，該方法通過一定的數(shù)學(xué)變換可以推導(dǎo)出三大守恒定律；在微觀層次上主要是依據(jù)分子動力學(xué)方法開展的模擬。在單尺度仿真中，盡管已經(jīng)獲得了不錯的成果，但是仍然存在一些局限性，例如模型無法準確恰當?shù)孛枋鲆呻y問題，現(xiàn)有算法難以適應(yīng)不斷增長的模型規(guī)模和精度要求以及計算機硬件的速度和效率的限制等問題。

多尺度耦合模擬是解決上述問題的有效途徑[1]。多尺度模擬在空間尺度上劃分出不同的層次，針對問題描述難易程度、模型規(guī)模精度、計算速度等問題，使用不同層次的模擬方法和手段[2]，以便能更恰當準確的描述疑難問題，提高模型計算的精度，實現(xiàn)并行計算，提高計算效率。多尺度模擬是仿真發(fā)展的趨勢，是大型復(fù)雜系統(tǒng)有力的描述手段[3]。

本文以管式固體氧化物燃料電池(SOFC)為例，從多尺度模擬的角度，在介觀層次上采取格子玻爾茲曼方法模擬電化學(xué)反應(yīng)擴散，在宏觀尺度上采取CFD方法模擬管式SOFC內(nèi)部的傳熱、流動過程，通過介觀層次和宏觀層次的耦合實現(xiàn)多尺度模擬。

2 格子玻爾茲曼方法的電化學(xué)反應(yīng)擴散模型

格子玻爾茲曼方法中模擬的是流體粒子團，流體粒子團可以理解為微觀充分大、宏觀充分小由若干流體分子組成的微團。其具體形式如下[4]

fi(x+eiΔt，t+Δt)=fi(x，t)+Ωi(f(x，t))

(1)

式中，fi是粒子速度分布函數(shù)；ei是當?shù)亓Ｗ铀俣龋沪竔(f(x，t))是碰撞算子，它表示粒子團發(fā)生碰撞后fi的變化率。

格子玻爾茲曼方法中關(guān)于電化學(xué)反應(yīng)擴散問題，其基本思想是對碰撞項進行修正，即在碰撞項中加入反應(yīng)項[5]，具體形式如下

(2)

其中

(3)

(4)

將(2)、(3)和(4)式代入(1)式，利用恰普曼—恩斯科格展開，可得到宏觀的反應(yīng)擴散方程：

(5)

式中，M是格子模型中的速度數(shù)；D′為擴散系數(shù)；R(ρ)是反應(yīng)項，形式由具體化學(xué)反應(yīng)決定。

3 管式SOFC的CFD模型

SOFC的基本工作原理如圖1。燃料和空氣需先壓縮并加熱到規(guī)定溫度。燃料被送到SOFC的陽極，空氣被送到SOFC的陰極。在陰極室中，隨著空氣的流動，O2逐漸溶入多孔電極中，在多孔電極中得到電子生成O2-，O2-通過電解質(zhì)內(nèi)的離子導(dǎo)體傳輸至陽極。在陽極室中，CH4與H2O反應(yīng)生成H2、CO、CO2的混合氣體，其中H2、CO與陰極傳送過來的O2-發(fā)生電化學(xué)反應(yīng)生成H2O、CO2，并產(chǎn)生電子，電子經(jīng)外部電路被再次送到陰極，如此循環(huán)往復(fù)，便產(chǎn)生了電能。

圖1 SOFC工作原理

過程中涉及到的主要化學(xué)反應(yīng)有[6]

CH4+H2O=3H2+CO

(重整反應(yīng))

CO+H2O=H2+CO2

(置換反應(yīng))

H2+O2-=H2O+2e

(陽極)

CO+O2-=CO2+2e

(陽極)

0.5O2+2e=O2-

(陰極)

總體反應(yīng)為

CH4+2H2O=4H2+CO2

CO+0.5O2=CO2

H2+0.5O2=H2O

3.1 質(zhì)量守恒方程

管式SOFC的質(zhì)量守恒方程可由下式表示[7]

?·(ερu)=0

(6)

式中，ε是管式SOFC電極的孔隙率。

燃料電池內(nèi)部是由CH4、H2O、H2、CO、CO2等氣體形成的混合氣體，其密度為：

(7)

式中，mi是混合氣體中某氣體的質(zhì)量百分比，Mi是混合氣體中某氣體的摩爾質(zhì)量。

3.2 動量守恒方程

管式SOFC的動量守恒方程可由下式表示

?·(ερuu)=-ε?P+?·(εμeff?u)+Su

(8)

式中，Su為源項，在流道中Su數(shù)值取零。在多孔電極中發(fā)生電化學(xué)反應(yīng)使得氣體濃度發(fā)生改變，進而引起傳質(zhì)現(xiàn)象的產(chǎn)生。在傳熱傳質(zhì)過程中，由于流道內(nèi)的速度分布差異很大，在電極界面將產(chǎn)生很大的速度梯度，因此不能忽視粘性力和慣性力的影響，此時應(yīng)用布林克曼方程描述動量，源項為

(9)

式中，μeff是氣體的有效粘性系數(shù)，k是多孔電極的滲透率。

3.3 能量守恒方程

管式SOFC的能力守恒方程可由下式表示[8-9]

?·(ερCpuT)=?·(keff?T)+ST

(10)

式中，ST為熱源項，keff為有效導(dǎo)熱系數(shù)(W/m*K)。

3.4 邊界條件

在流道入口處初始流速和溫度參數(shù)采用定值，即

uin=u0；Tin=T0

在流道出口處假設(shè)氣體流速和溫度已達到穩(wěn)定狀態(tài)，即

4 尺度耦合的技術(shù)路線

對管式SOFC采用多尺度耦合模擬的技術(shù)路線[8]如圖2。在宏觀尺度上，采用三大守恒定律建立CFD模型模擬流動和傳熱問題。在介觀尺度上，采用格子玻爾茲曼方法模擬反應(yīng)擴散問題，其原因是SOFC陽極為多孔介質(zhì)，分布不連續(xù)，采用宏觀模型難以進行準確描述。

圖2 管式SOFC多尺度耦合模擬的技術(shù)路線圖

本文的多尺度耦合計算流程如圖3。假設(shè)任意初始溫度T，根據(jù)式(1)～式(4)計算反應(yīng)放熱量，并以此作為邊界條件，輸入流動和傳熱的CFD模型中計算溫度場，取其平均值Tc， 將T和Tc進行比較，當相對誤差關(guān)系式err=(T-Tc)/Tc成立，則退出計算，計算流程如圖3。

圖3 計算流程圖

5 數(shù)值模擬結(jié)果及分析

5.1 格子玻爾茲曼電化學(xué)反應(yīng)模擬結(jié)果及分析

本文采用格子玻爾茲曼D2Q9模型模擬管式SOFC的電化學(xué)反應(yīng)擴散現(xiàn)象，二維網(wǎng)格劃分為256×256，H2、O2的初始濃度分別為15.2250mol和9.9167mol，網(wǎng)格初始分布為正態(tài)隨機分布。管式SOFC的電流密度可由下式表示[10]

(11)

式中，z為反應(yīng)的H2的量。

圖4所示為1000℃下反應(yīng)達到平衡后H2的濃度分布。圖中x軸表示管式SOFC陽極的厚度，y軸表示管式SOFC陽極的長度，z軸表示反應(yīng)達到平衡后H2的濃度分布。圖4結(jié)果表明反應(yīng)達到平衡后陽極內(nèi)H2濃度呈多個正態(tài)隨機分布，這種正態(tài)分布與所采取的初始分布相關(guān)，也與反應(yīng)進行的程度相關(guān)，反應(yīng)進行的越徹底，各網(wǎng)格節(jié)點上濃度的偏差越小，因此隨著反應(yīng)溫度的升高，反應(yīng)進行的程度越高，陽極內(nèi)各節(jié)點達到平衡后，其濃度分布越均勻。

圖4 1000℃時反應(yīng)平衡后H2濃度分布

圖6 電流密度與反應(yīng)溫度關(guān)系

圖7 反應(yīng)時間與反應(yīng)溫度關(guān)系

圖5、6、7所示為H2的濃度、電流密度、反應(yīng)時間與反應(yīng)溫度的關(guān)系。結(jié)果表明，H2的濃度隨著反應(yīng)溫度的增加而逐漸減少；電流密度隨著反應(yīng)溫度的增加而逐漸變大，電流密度的大小和H2量有關(guān)，H2量越多電流密度就越大；反應(yīng)達到平衡所需要的時間隨著反應(yīng)溫度的增加變得越來越短；本文在模擬時，模擬所需要的相關(guān)參數(shù)，如指前因子、活化能等參數(shù)都是依據(jù)有關(guān)手冊選取的，這可能與實際情況存在偏差，導(dǎo)致800℃時反應(yīng)所需時間較長。

5.2 耦合數(shù)值模擬結(jié)果及分析

圖8為空氣入口流速為1.634m/s時管式SOFC多尺度耦合模擬的溫度場分布。結(jié)果表明，空氣溫度和管壁溫度在入口處變化較大，隨著空氣的流動，溫度很快達到穩(wěn)定。管式SOFC內(nèi)部的溫度可以通過控制空氣的流動狀態(tài)來控制，因為空氣起著冷卻的作用，隨著空氣流速的增加，被空氣帶走的電化學(xué)反應(yīng)的放熱量越多，溫度降低的越多，達到平衡后便處于穩(wěn)定狀態(tài)，且溫度分布趨于均勻。

圖8 管式SOFC多尺度耦合模擬的溫度場分布

在流動方向上，隨著電化學(xué)反應(yīng)的不斷進行，放出的熱量不斷增多，流道中的氣體被不斷加熱，氣體溫度逐漸升高，這加速了電化學(xué)反應(yīng)的進行，使得放熱量快速增加。另外，在這過程中管壁也被不斷加熱，溫度逐漸升高。管式SOFC整體溫度分布情況為入口溫度最低，出口溫度達到最高值。因此，管式SOFC內(nèi)部的溫度分布與反應(yīng)放熱量相關(guān)。

綜上所述，可以得出影響管式SOFC內(nèi)部溫度分布的兩個重要因素，一個為反應(yīng)放熱量，另一個為空氣的流動狀態(tài)。

通過多尺度耦合模擬，發(fā)現(xiàn)隨著燃氣的流動，反應(yīng)發(fā)生的位置不斷向后移動，反應(yīng)產(chǎn)生的熱量不斷增加，對流道和管壁加熱效果越明顯。因此沿著管式SOFC流動方向，流道和管壁的溫度越來越高。

圖9所示為計算穩(wěn)定后管式SOFC的電解質(zhì)壁面的溫度分布，圖中以1.634m/s和1.0m/s為例，分別以熱量和對流換熱作為邊界條件計算溫度分布。結(jié)果表明，采用熱量作為邊界條件的溫度整體上比采用對流換熱作為邊界條件的溫度高。當采用熱量作為邊界條件時，流速高的溫度高，其原因在于熱量邊界對溫度分布影響大，流速不足以將熱量全部帶走。當采用對流換熱作為邊界條件時，流速低的溫度高，其原因在于對流換熱邊界中對流效應(yīng)影響大，空氣流速對溫度分布起主要作用，流速越高，溫度越低，因此穩(wěn)定后流速低的溫度高。

圖9 管式SOFC電解質(zhì)壁面溫度分布

圖10、11所示為空氣流道入口流速u=1.634m/s時的速度場分布和管式SOFC中截面的速度分布。結(jié)果表明，管式SOFC的中部流速最大，壁面處流速最小。

圖10 管式SOFC流場分布

圖11 管式SOFC中截面速度分布

圖12是管式SOFC管壁中軸線上的速度分布，結(jié)果表明多孔介質(zhì)中的流速分布與多孔介質(zhì)的準周期性分布一致。

圖12 管式SOFC管壁中軸線速度分布

6 結(jié)論

本文從介觀層次和宏觀層次出發(fā)，介觀層次上采用格子玻爾茲曼方法建立了管式SOFC的電化學(xué)反應(yīng)擴散模型，宏觀層次上采用計算流體力學(xué)方法建立了管式SOFC的CFD模型，并將格子玻爾茲曼方法計算的反應(yīng)放熱量作為邊界條件求解CFD模型，實現(xiàn)了多尺度耦合模擬，得到了管式SOFC的穩(wěn)態(tài)數(shù)值模擬結(jié)果。從結(jié)果中可知，H2的濃度和反應(yīng)時間隨反應(yīng)溫度逐漸降低，電流密度隨反應(yīng)溫度逐漸增加。空氣流道中的氣體溫度和管壁溫度沿著流動方向逐漸升高。影響溫度分布的兩個重要因素，一個是反應(yīng)放熱量，另一個是空氣的流動狀態(tài)。反應(yīng)放熱量越多，管式SOFC內(nèi)部的總體溫度越高。另外，空氣流速越高，管式SOFC內(nèi)部溫度降低的越多，達到平衡后便處于穩(wěn)定狀態(tài)，且溫度分布趨于均勻。

結(jié)果表明采用多尺度耦合方法模擬管式SOFC是行之有效的，在一定程度上解決了管式SOFC介觀尺度和宏觀尺度多物理場耦合的問題。