單顆粒多孔MgO吸附劑CO2吸附流動(dòng)傳輸?shù)母褡覤oltzmann方法模擬*

劉波，丁玉棟，廖強(qiáng)，朱恂，王宏

(重慶大學(xué)工程熱物理研究所 低品位能源利用及系統(tǒng)教育部重點(diǎn)實(shí)驗(yàn)室,重慶 400030)

近幾十年能源消耗量不斷增加，CO2排放量急劇上升，從而引發(fā)了溫室效應(yīng)等一系列的環(huán)境問(wèn)題。從排放源來(lái)看，中國(guó)火力發(fā)電CO2排放量約占工業(yè)CO2總排放量的近50%[1]，因此分離和捕集火力發(fā)電廠煙道氣中的CO2成為控制CO2排放最有效的手段之一?；痣姀S煙道氣往往排放溫度較高(100 ℃左右)、CO2體積比僅約10%～15%，且含有水蒸氣等雜質(zhì)性氣體成分。針對(duì)此特點(diǎn)，金屬氧化物及其鹽類通?？稍谳^高溫度下吸附CO2，且能夠耐受一定的水蒸氣，因此是一種非常理想的CO2吸附劑。在金屬氧化物及其鹽類材料中，Mg由于在自然界中儲(chǔ)量大，MgO易于制取，可在200 ℃以下吸附CO2，因此相對(duì)于CaO吸附劑和鋰鹽吸附劑，MgO吸附劑的CO2吸附溫度在200 ℃以下，是一種非常適于捕集電廠煙道氣中CO2的吸附劑。

現(xiàn)階段對(duì)于MgO吸附劑的研究多在于實(shí)驗(yàn)方面，包括提高比表面積[2]、優(yōu)化孔道結(jié)構(gòu)[3]、獲得顆粒成型方法[4]及其對(duì)材料吸附性能的影響等，如Bian等[5]采用熱分解法制備介孔MgO材料，和市面所售相比具有更大的比表面積和更高的吸附量；Song等[6]通過(guò)煅燒法合成介孔MgO材料，其比表面積大，吸附選擇性好；聞霞等[7]采用等體積浸漬法制取MgO/Al2O3吸附劑，獲得MgO負(fù)載量、吸附溫度等因素對(duì)吸附劑吸附CO2性能的影響；陳林等[4]制備碳參雜MgO顆粒，具有較好的機(jī)械強(qiáng)度與較高的比表面積。但現(xiàn)階段缺少數(shù)值模擬的手段分析吸附過(guò)程中多孔MgO顆粒內(nèi)煙道混合氣體的滲流及CO2傳輸反應(yīng)特性。

格子Boltzmann方法(lattice Boltzmann method,LBM)作為一種基于動(dòng)理學(xué)的介觀尺度數(shù)值理論，其粒子背景使得該方法在單組分流動(dòng)、傳熱傳質(zhì)、多組分多相流動(dòng)、多孔介質(zhì)內(nèi)流動(dòng)等領(lǐng)域有著廣泛的應(yīng)用[8-10]。國(guó)內(nèi)外已經(jīng)利用該方法開(kāi)展了大量氣固反應(yīng)的數(shù)值模擬研究工作，如：Sullivan等[11]利用格子Boltzmann方法開(kāi)展反應(yīng)床尺度流動(dòng)與傳質(zhì)的耦合建模，且該模型適用于低擴(kuò)散系數(shù)、高貝克萊數(shù)情況。張婷等[12]利用格子Boltzmann方法對(duì)孔隙尺度下堿性金屬氧化物CaO吸收CO2過(guò)程進(jìn)行模擬，分析摩爾體積比與達(dá)姆科勒數(shù)對(duì)CaO轉(zhuǎn)化率的影響。Bohn等[13]利用LBM法對(duì)CO與Fe2O3還原反應(yīng)進(jìn)行模擬，結(jié)果與實(shí)驗(yàn)數(shù)據(jù)吻合很好。Zhou等[14]利用格子Boltzmann方法研究孔隙尺度下氣體在固體表面吸附過(guò)程，并分析氣相速度分布及孔隙率對(duì)吸附平衡態(tài)的影響。已有研究針對(duì)MgO顆粒吸附法捕集電廠煙道氣中CO2的相關(guān)研究甚少，因此本文從表征體元(representative elementary volume,REV)尺度出發(fā)，采用格子Boltzmann方法對(duì)含化學(xué)反應(yīng)的單顆粒多孔MgO吸附劑吸附CO2的過(guò)程進(jìn)行數(shù)值模擬研究，其中假設(shè)顆粒為規(guī)則球體，混合氣體由N2與CO2組成，吸附過(guò)程絕熱，基于二級(jí)反應(yīng)動(dòng)力學(xué)模型研究單顆粒多孔MgO內(nèi)CO2的傳輸特性，分析孔隙率與顆粒粒徑變化對(duì)于顆粒內(nèi)滲流速度、CO2吸附速率及MgO顆粒轉(zhuǎn)化率的影響。

1 格子-Boltzmann模型

1.1 流動(dòng)格子-Boltzmann模型

本文采用Guo和Zhao[15]提出的表征體元尺度下滲流廣義LBE模型模擬單顆粒多孔MgO吸附劑內(nèi)主體氣體(N2)的流動(dòng)過(guò)程，并采用D2Q9格子模型，其控制方程為

(1)

式中：fi(r,t)為網(wǎng)格點(diǎn)r處t時(shí)刻粒子的分布函數(shù)；δt為時(shí)間步長(zhǎng)；i=0～8；ei為粒子離散速度，e0=0，e1=-e3=(1,0)c，e2=-e4=(0,1)c，e5=-e7=(1,1)c，e6=-e8=(-1,1)c；其中c為格子速度,c=δx/δt，δx為格子步長(zhǎng),τv為流動(dòng)格子方程的無(wú)量綱松弛時(shí)間。速度的平衡態(tài)分函數(shù)為

(2)

外力項(xiàng)Fi定義為

(3)

式中:當(dāng)?shù)芈曀贋閏s；ε為孔隙率；ωi為權(quán)系數(shù)，ω0=4/9，ωi=1/9(i=1～4)，ωi=1/36(i=5～8)；F是包含有外力和介質(zhì)阻力的合力，其表達(dá)式為

(4)

式中：υ為流體運(yùn)動(dòng)黏度;G為外部體積力;K為滲透率;Fε為結(jié)構(gòu)參數(shù)，利用Ergun方程得到：

(5)

滲流的密度ρ與宏觀滲流速度U定義為

(6)

式中速度v為一個(gè)臨時(shí)速度，公式為

(7)

a0，a1為參數(shù)，其取值為：

(8)

對(duì)式(1)做Chanpman-Enskog多尺度展開(kāi)得到其對(duì)應(yīng)的宏觀動(dòng)量方程如下：

(9)

(10)

1.2 傳質(zhì)格子-Boltzmann模型

1.2.1 傳質(zhì)模型

利用格子Boltzmann傳質(zhì)模型模擬多孔MgO顆粒內(nèi)的CO2流動(dòng)傳輸與吸附過(guò)程。CO2溶質(zhì)視為被動(dòng)標(biāo)量，認(rèn)為其與主體氣體(N2)之間是單向作用，溶質(zhì)的對(duì)流擴(kuò)散過(guò)程不影響主體氣體N2的流動(dòng)。根據(jù)Zeiser等[16]的研究，在不影響計(jì)算精度的情況下，可將速度離散方向由9個(gè)降為5個(gè)(D2Q9模型變?yōu)镈2Q5模型，i取0～4)，其控制方程為

gi(r+ciδt)-gi(r,t)=

(11)

式中：gi為溶質(zhì)濃度的分布函數(shù)，τc為傳質(zhì)格子方程的無(wú)量綱松弛時(shí)間，R為反應(yīng)源項(xiàng)，Ji為常數(shù)。二維情況下Ji=J0(i=0)，Ji=(1-J0)/4(i=1～4)，J0為可調(diào)參數(shù)，0≤J0≤1，geq為濃度的平衡態(tài)分布函數(shù),公式為

(12)

對(duì)濃度分布函數(shù)gi求和可以獲得其濃度

C=∑gi.

(13)

與流動(dòng)方程一樣，對(duì)式(11)做Chanpman-Enskog多尺度展開(kāi)得到其宏觀的傳質(zhì)對(duì)流擴(kuò)散方程

(14)

式中，擴(kuò)散系數(shù)為

D=δtc2CQ(1-J0)(τm-0.5).

(15)

1.2.2 傳質(zhì)源項(xiàng)

固體MgO顆粒吸附CO2過(guò)程是在其表面進(jìn)行，酸性CO2氣體與MgO上的堿性位點(diǎn)O2-發(fā)生反應(yīng)，生成了單齒碳酸根、雙齒碳酸根[17]?？紤]到反應(yīng)本身的復(fù)雜性，將該反應(yīng)視為基元反應(yīng)，顆粒內(nèi)每個(gè)網(wǎng)格節(jié)點(diǎn)為反應(yīng)位點(diǎn)，并參考Bohn等[13]采用二級(jí)反應(yīng)動(dòng)力學(xué)模型對(duì)傳質(zhì)源項(xiàng)進(jìn)行處理：

(16)

(17)

2 模型驗(yàn)證及邊界條件

2.1 流動(dòng)耦合傳質(zhì)模型的驗(yàn)證

對(duì)二維矩形通道內(nèi)充滿多孔介質(zhì)的Poiseuille流動(dòng)及傳質(zhì)過(guò)程進(jìn)行模擬，其沿Y軸方向的高度為H，沿X軸方向最大流速為U0，已知其流動(dòng)宏觀控制方程為

(18)

式(18)的解析解為

U=GK/υ×(1-cosh(γ(Y-H/2))/cosh(γH/2)),

(19)

U0=GK/υ×(1-cosh-1(γH/2)).

(20)

式中：γ=(υε/Kυe)。假設(shè)通道的上、下壁面濃度為定值，分別為Cu=1.0、Cb=0.0，則對(duì)流擴(kuò)散方程解析解為

(21)

模擬中，矩形計(jì)算域采用的網(wǎng)格數(shù)為76×76，參數(shù)設(shè)定為：ε=0.1，Da=0.000 1，Re=0.1，Sc=100。圖1(a)和圖1(b)分別為出口處，沿Y方向速度和濃度的分布曲線，結(jié)果表明采用格子Boltzmann方法得到的模擬果與數(shù)值解吻合得很好。

圖1 出口處沿Y軸方向速度和濃度的分布曲線Fig.1 The velocity(a) and concentration(b) profiles along the Y-direction at the outlet

2.2 數(shù)值模擬條件

為節(jié)省計(jì)算時(shí)間，本文采用多塊網(wǎng)格技術(shù)，將計(jì)算域劃分為上中下3個(gè)部分，其上、下部分均為粗網(wǎng)格，中間部分為細(xì)網(wǎng)格，計(jì)算域?yàn)槎S矩形計(jì)算域，顆粒的無(wú)量綱直徑為d，X軸方向無(wú)量綱長(zhǎng)度Lx為顆粒直徑d的7倍，Y軸方向無(wú)量綱長(zhǎng)度Ly為顆粒直徑d的10倍，混合氣體(N2與CO2組成)沿X軸方向流入，圖2為計(jì)算域示意圖。

圖2 計(jì)算域示意圖Fig.2 Computational domain diagram

初始條件：參照課題組前期實(shí)驗(yàn)結(jié)果[18]，N2分壓0.85 bar，CO2分壓0.15 bar,氣體宏觀溫度取50 ℃，N2密度0.895 1 kg/m3，CO2濃度5.64 mol/m3，MgO顆粒密度3.58 g/cm3(濃度為89 500 mol/m3)，入口氣體速度0.5 m/s，入口N2密度、濃度及速度均保持恒定。粒徑dp1=100 μm，參考長(zhǎng)度Xc1=5×10-6m，參考速度Uc=0.1 m/s，參考濃度Cc=5.64 mol/m3,參考密度ρc=1.053 kg/m3。當(dāng)粒徑發(fā)生變化時(shí)，改變參考長(zhǎng)度(計(jì)算域大小相應(yīng)變化，但其無(wú)量綱值不變，其余物理量值不變)。對(duì)于參考量取值與無(wú)量綱換算見(jiàn)表1。通過(guò)計(jì)算可知粒徑dp1=100 μm時(shí)，Re=3.2，粒徑dp2=200 μm時(shí)，Re=6.4，粒徑dp3=400 μm時(shí)，Re=12.8，施密特?cái)?shù)均為Sc=υp/Dp=0.96。

邊界條件：入口處N2密度與速度恒定，采用Zou、He入口邊界；上、下邊界處認(rèn)為速度邊界層已充分發(fā)展，沿Y方向速度不再變化，采用非平衡外推格式；出口認(rèn)為速度不再變化，仍采用非平衡外推格式。對(duì)于濃度邊界，入口CO2濃度不變，采用恒定濃度邊界；上、下及出口邊界處認(rèn)為濃度梯度為零，同采用非平衡外推格式。

表1 物理參數(shù)的無(wú)量綱轉(zhuǎn)化表Table 1 Non-dimensional transformation table for physical parameters

3 結(jié)果分析

3.1 粒徑大小影響

圖3為粒徑dp2=200 μm時(shí)不同孔隙率下的局部流場(chǎng)圖，從圖中可以看到孔隙率不同時(shí)流線穿過(guò)MgO顆粒的偏折程度也不同，孔隙率越小流線偏折越厲害，迂曲度越大；由于流動(dòng)本身雷諾數(shù)較小，所以不存在尾渦與氣體回流。流場(chǎng)結(jié)構(gòu)的差異與流速的不同會(huì)直影響到對(duì)流傳質(zhì)過(guò)程，進(jìn)而影響顆粒內(nèi)氣體濃度，并最終影響氣體與固體的反應(yīng)吸附。圖4為不同粒徑下通過(guò)圓心截面速度的分布情況，可以看出當(dāng)粒徑較小(dp1=100 μm)時(shí)，顆粒內(nèi)的主流體速度越大且其最小速度越接近球心；顆粒粒徑越大時(shí)，顆粒內(nèi)速度越小，且其最小速度越靠近沿流動(dòng)方向顆粒的后部。

圖3 不同孔隙率下局部流場(chǎng)分布情況Fig.3 Local streamline distribution at different porosities

圖4 不同粒徑下通過(guò)圓心截面處速度分布Fig.4 Velocity distribution along the X-direction through the center of the particle with different sizes

圖5為不同粒徑下，t=1.0 s，ε=0.3時(shí)細(xì)網(wǎng)格計(jì)算域內(nèi)濃度分布云圖；圖6為不同粒徑下通過(guò)圓心截面處沿X軸方向的濃度分布。通過(guò)圖5和圖6可以看到，粒徑越小顆粒內(nèi)平均濃度值越大；并且隨著顆粒直徑的增大，CO2最低濃度位置逐漸沿流動(dòng)方向向顆粒后部移動(dòng)。這是由于入口初始速度相同，相同孔隙率下粒徑越大，流動(dòng)阻力越大，顆粒內(nèi)流速越小，對(duì)流傳質(zhì)越弱，傳質(zhì)阻力增加，自然顆粒內(nèi)濃度值更低。

圖5 相同時(shí)間不同粒徑下細(xì)網(wǎng)格內(nèi)CO2濃度分布Fig.5 CO2 concentration distribution in fine grids at the same time with different particle sizes

圖6 不同粒徑下通過(guò)圓心截面CO2濃度分布Fig.6 CO2 concentration distribution along the X-direction through the center of the particle with different sizes

圖7為粒徑dp2=200 μm孔隙率為0.5時(shí)，不同時(shí)步MgO顆粒內(nèi)固體轉(zhuǎn)化率X(r,t)的分布情況，可以看出由于CO2氣體濃度分布的不均勻，固體顆粒的轉(zhuǎn)化率也不均勻，顆粒內(nèi)沿流動(dòng)方向顆粒前部濃度較高，其局部固體轉(zhuǎn)化率也更高，由于其本身雷諾數(shù)較小，不存在尾渦與氣體回流，而顆?？亢筇幍牧鲃?dòng)及傳質(zhì)阻力最大所以濃度最低點(diǎn)及轉(zhuǎn)化率最低點(diǎn)均出現(xiàn)在沿流動(dòng)方向顆粒后部。

圖7 不同時(shí)步下顆粒內(nèi)固體轉(zhuǎn)化率分布情況Fig.7 Adsorbent conversion rate distributions in particle at different calculation time steps

3.2 孔隙率影響

圖8為粒徑dp1=100 μm，相同計(jì)算時(shí)步，不同孔隙率下通過(guò)顆粒球心截面處CO2氣體濃度的分布情況。從圖8可以發(fā)現(xiàn)，同一粒徑、相同計(jì)算時(shí)步下，隨著孔隙率的增加，顆粒內(nèi)的CO2濃度越高。主要原因?yàn)榭紫堵试酱髮?dǎo)致顆粒內(nèi)的滲透率越大(Ergun方程算得)，氣體傳輸?shù)淖枇Ω?，并且相同粒徑時(shí)固體的量隨孔隙率增加而減小，使相同時(shí)步下吸附的CO2氣體量更少，CO2在顆粒內(nèi)更容易聚集，因而MgO顆粒內(nèi)的CO2濃度越高。

圖8 不同孔隙率下通過(guò)球心截面的CO2濃度分布Fig.8 CO2 concentration distribution along the X-direction through the center of the particle with different porosities

圖9展示粒徑分別為100、200和400 μm時(shí)孔隙率變化對(duì)顆粒內(nèi)吸附劑平均轉(zhuǎn)化率Xa的影響。由圖可知，吸附過(guò)程初始階段，固體轉(zhuǎn)化率較小，轉(zhuǎn)化率隨時(shí)間基本呈現(xiàn)線性變化；隨著時(shí)間推移，固體轉(zhuǎn)化率曲線斜率越來(lái)越小，這是因?yàn)殡S著反應(yīng)的進(jìn)行，MgO顆粒吸附位點(diǎn)不斷消耗，雖固體轉(zhuǎn)化率越來(lái)越高，但剩余固體的量越來(lái)越少，根據(jù)反應(yīng)動(dòng)力學(xué)模型可知吸附速率受到顆粒固體轉(zhuǎn)化率的限制，顆粒的固體轉(zhuǎn)化率變大，反應(yīng)速率降低，從而導(dǎo)致固體轉(zhuǎn)化率曲線斜率減小，CO2吸附逐漸達(dá)到飽和。此外，由前文分析結(jié)果知，相同粒徑、相同計(jì)算時(shí)步下，孔隙率的增加使MgO顆粒內(nèi)CO2濃度升高，根據(jù)反應(yīng)動(dòng)力學(xué)模型可知，CO2濃度越高會(huì)導(dǎo)致顆粒內(nèi)反應(yīng)位點(diǎn)的吸附反應(yīng)速率更大，固體的量消耗得更快。以上因素共同作用，導(dǎo)致圖9中相同粒徑下，孔隙率越大，反應(yīng)速率越快，固體顆粒轉(zhuǎn)化率的變化速率越快，CO2達(dá)到飽和的過(guò)程也越快。

圖9 不同孔隙率下平均轉(zhuǎn)化率隨時(shí)間的變化Fig.9 Variation in average adsorbent conversion rate with time at different porosities

4 結(jié)論

本文從表征體元尺度出發(fā)，基于二級(jí)反應(yīng)動(dòng)力學(xué)模型，利用格子Boltzmann方法對(duì)單顆粒多孔MgO吸附劑化學(xué)吸附CO2過(guò)程流動(dòng)傳輸特性進(jìn)行數(shù)值模擬，分析孔隙率與顆粒粒徑變化對(duì)于顆粒內(nèi)滲流速度、CO2吸附速率及MgO顆粒內(nèi)固體轉(zhuǎn)化率的影響。主要結(jié)論如下：

1)固體多孔MgO顆粒吸附CO2過(guò)程可將吸附位點(diǎn)處理為網(wǎng)格節(jié)點(diǎn)，節(jié)點(diǎn)處的反應(yīng)視為基元反應(yīng)，符合二級(jí)反應(yīng)動(dòng)力學(xué)模型。

2)沿氣體流動(dòng)方向顆粒內(nèi)CO2濃度與顆粒轉(zhuǎn)化率逐漸遞減；在相同時(shí)間下顆粒粒徑越大，顆粒內(nèi)滲流速度越小，CO2濃度越低，轉(zhuǎn)化率越小，吸附速率越慢。

3)相同粒徑和相同時(shí)間下，MgO顆粒孔隙率越大，顆粒內(nèi)CO2濃度越高，MgO固體轉(zhuǎn)化率越高，平均吸附速率越快。