氣固流化床啟動階段擋板內(nèi)構(gòu)件受力特性的CFD-DEM 模擬

李鐵男，趙碧丹，趙鵬，張永民，王軍武，4

（1 中國石油大學(xué)（北京）重質(zhì)油國家重點實驗室，北京 102249； 2 中國科學(xué)院過程工程研究所多相復(fù)雜系統(tǒng)國家重點實驗室，北京 100190； 3 中國科學(xué)院大學(xué)化學(xué)工程學(xué)院，北京 100049；4中國科學(xué)院綠色過程制造創(chuàng)新研究院，北京 100190）

引 言

流化床反應(yīng)器具有良好的傳質(zhì)傳熱性能，且能夠較好輸送固體顆粒物料，因此廣泛用于工業(yè)生產(chǎn)過程[1]。反應(yīng)器內(nèi)設(shè)置合適的內(nèi)構(gòu)件可強化傳質(zhì)與傳熱并提高反應(yīng)器性能[2-3]。基于以往研究經(jīng)驗，采用傾斜葉片的流化床內(nèi)構(gòu)件破碎氣泡和改善氣固接觸的效果更好，而且還可以通過調(diào)整結(jié)構(gòu)參數(shù)有效調(diào)節(jié)氣固相的停留時間分布，具有比垂直構(gòu)件和具有豎直葉片的網(wǎng)狀格柵更好的反應(yīng)強化效果，因此此類內(nèi)構(gòu)件在工業(yè)中應(yīng)用最廣泛[3]，典型的例子如百葉窗格柵、脊型和塔型內(nèi)構(gòu)件等。由于工業(yè)流化床內(nèi)存在復(fù)雜的動態(tài)兩相流動以及不同操作狀態(tài)之間的切換，內(nèi)構(gòu)件會受到不同形式和強度的作用力，如果內(nèi)構(gòu)件結(jié)構(gòu)設(shè)計不合理，很可能會致使內(nèi)構(gòu)件發(fā)生損傷甚至被破壞，進(jìn)而導(dǎo)致工業(yè)裝置出現(xiàn)故障和經(jīng)濟損失。因此，為了保障流化床內(nèi)構(gòu)件的長周期可靠性，系統(tǒng)掌握內(nèi)構(gòu)件在流化床內(nèi)的受力特性尤為重要[4]。

早期針對流化床中內(nèi)構(gòu)件受力特性的研究主要采用實驗手段，大多數(shù)關(guān)注的是正常流化狀態(tài)下B 類顆粒流化床內(nèi)圓形水平換熱管的受力特性[5-9]，其中Grace等[7-8]的研究為揭示正常流化狀態(tài)下圓管受力的機理奠定了良好基礎(chǔ)。近年來，Zhang等[4,10-13]系統(tǒng)研究了化工領(lǐng)域A 類顆粒流化床反應(yīng)器中經(jīng)常使用的斜片擋板內(nèi)構(gòu)件的受力特性，發(fā)現(xiàn)啟動階段流化床層內(nèi)部構(gòu)件會受到一個很大的向上載荷脈沖，其峰值可達(dá)正常流化狀態(tài)下內(nèi)構(gòu)件受到的平均載荷值的數(shù)倍甚至一個數(shù)量級以上[11]，這對工業(yè)裝置內(nèi)構(gòu)件可靠性是一個巨大的潛在威脅，很可能是導(dǎo)致工業(yè)裝置內(nèi)構(gòu)件損壞的一種新機理。

對于冷態(tài)實驗室條件下一些小型簡單結(jié)構(gòu)內(nèi)構(gòu)件，實驗方法可以獲得較為準(zhǔn)確的受力特性數(shù)據(jù)，但對于工業(yè)裝置高溫高壓的苛刻條件下的大型復(fù)雜結(jié)構(gòu)內(nèi)構(gòu)件，其在流化床中的受力特性預(yù)測則難以使用實驗方法，因此，使用近年來日益成熟的計算流體力學(xué)模擬方法研究流化床中內(nèi)構(gòu)件的受力是一種有益的嘗試。Higashida等[14]采用虛擬顆粒方法（FPM 方法）模擬研究了三維鼓泡床中懸浮球體受到的動態(tài)垂直力，并將模擬計算結(jié)果與帶有拉格朗日傳感器系統(tǒng)的懸浮球體所測得的實驗值相比較。研究發(fā)現(xiàn)，球體在裝置內(nèi)床層流化過程中，其動態(tài)垂直力的波動與氣泡的運動密切相關(guān)，球體受到的合力由流體作用力與顆粒碰撞作用力兩部分組成。Yan 等[15]采用雙流體模型模擬了氣固兩相體系內(nèi)不同形狀的大顆粒受力演化規(guī)律。模擬結(jié)果顯示：隨著表觀氣速的增加，大顆粒受到的作用力也相應(yīng)增大。Nagahashi 等[16]采用FPM 模擬方法，對二維流化床中單個氣泡經(jīng)過水平圓管的受力機理進(jìn)行了模擬研究，并將模擬結(jié)果與實驗結(jié)果相比較。結(jié)果顯示：模擬與實驗中氣泡的運動軌跡吻合較好，實驗拍攝氣泡經(jīng)過圓管過程中，氣泡尾流攜帶的顆粒撞擊圓管會使圓管受到一個脈沖峰值。

本研究以啟動階段流化床層內(nèi)部構(gòu)件的非正常受力現(xiàn)象為研究對象，首先基于CFD-DEM 數(shù)值模擬方法，建立流化床擋板內(nèi)構(gòu)件表面載荷強度的統(tǒng)計方法，并通過與前期流化床啟動階段斜片擋板內(nèi)構(gòu)件受力特性的實驗數(shù)據(jù)[11]進(jìn)行對照，以驗證該載荷統(tǒng)計計算方法的正確性及合理性。進(jìn)一步，以該方法研究不同參數(shù)（碰撞模型參數(shù)、表觀氣速、顆粒粒徑）對啟動階段斜片擋板內(nèi)構(gòu)件受力特性的影響規(guī)律，旨在為下一步揭示流化床啟動階段內(nèi)構(gòu)件非正常受力現(xiàn)象的內(nèi)在機理奠定良好基礎(chǔ)，也為開發(fā)工業(yè)流化床內(nèi)構(gòu)件受力特性的預(yù)測工具奠定理論基礎(chǔ)。

1 擋板受力統(tǒng)計計算方法

氣固擋板流化床啟動階段的本質(zhì)是密相顆粒料層由未流化的固定床狀態(tài)逐漸轉(zhuǎn)變?yōu)榱骰矤顟B(tài)的一個過渡過程，相比其他多相流模擬方法，CFD-DEM 方法更適宜于準(zhǔn)確描述這一過渡過程。因此，本研究將基于CFD-DEM 方法，建立內(nèi)構(gòu)件表面載荷分布的統(tǒng)計計算方法，分析內(nèi)構(gòu)件表面受到的動態(tài)載荷信號。從統(tǒng)計力學(xué)角度分析，氣固流化床密相床層中的內(nèi)構(gòu)件受力特性是氣相分子與離散顆粒對內(nèi)構(gòu)件壁面作用特性的動態(tài)表征。CFDDEM 模擬方法可以獲得顆粒速度、位置和相互作用力(顆粒-顆粒、顆粒-壁面)等微觀顆粒相尺度信息，但是無法直接獲得工程上關(guān)注的內(nèi)構(gòu)件表面載荷強度及分布等宏觀物理量信息。

氣體壓力的物理本質(zhì)為氣固兩相流體系內(nèi)構(gòu)件受力求解提供一定參考。內(nèi)構(gòu)件氣固兩相流體系內(nèi)不同尺度下，從統(tǒng)計物理學(xué)可知固體顆粒與內(nèi)構(gòu)件表面作用的動態(tài)變化過程和氣體分子與壁面作用的動態(tài)過程相類似。此外，宏觀連續(xù)場信息通常都是借助顆粒流系統(tǒng)中每個離散顆粒和分子動力學(xué)中單個原子的作用力、位置和速度等微觀信息而獲得，故內(nèi)構(gòu)件表面受到的顆粒碰撞作用力特性可參考?xì)怏w對壁面作用形成的壓力過程及其動態(tài)作用力特性。在標(biāo)準(zhǔn)分子動力學(xué)模擬中，基于維里定理求解計算表面壓力。對于非均勻系統(tǒng)，壓力應(yīng)為張量形式，且是與位置相關(guān)的函數(shù)[17]。其中，壓力張量的公式包含兩部分：第一部分稱為動壓力，代表顆粒的運動對壓力的貢獻(xiàn)量；第二部分稱為位形壓力，代表粒子間的相互作用。

氣固兩相流體系內(nèi)構(gòu)件受力載荷的求解計算需要了解顆粒流體系內(nèi)應(yīng)力張量的組成形式。連續(xù)介質(zhì)場通常需要由離散粒子數(shù)據(jù)構(gòu)造，這些離散數(shù)據(jù)是每個原子或顆粒的位置、速度和相互作用力。在顆粒系統(tǒng)中，顆粒動理學(xué)理論是在Chapman等[18]的氣體動理學(xué)理論的基礎(chǔ)上發(fā)展起來的。一般而言，封閉固相應(yīng)力的方法主要有兩類：一類是顆粒動理學(xué)理論(KTGF)；另一類是基于CFD-DEM 方法的數(shù)值統(tǒng)計分析。KTGF 的推導(dǎo)方法主要分為以下幾種。(1) 考慮非彈性碰撞及顆粒體積影響。Gidaspow[19]基于Chapman等[18]建立的氣體動理論，進(jìn)而考慮非彈性碰撞影響和顆粒體積影響等不可忽略因素，建立起經(jīng)典的顆粒動理論。(2)考慮固相處于近平衡穩(wěn)態(tài)情況。Rao 等[20]進(jìn)一步考慮固相處于近平衡穩(wěn)態(tài)時，通過完善顆粒速度分布函數(shù)，導(dǎo)出更合理的固相本構(gòu)關(guān)系。(3)考慮體系內(nèi)介尺度結(jié)構(gòu)廣泛存在情況。Zhao等[21]考慮氣固兩相流系統(tǒng)中團聚物等介尺度結(jié)構(gòu)廣泛存在，建立了適用性更廣且更準(zhǔn)確的多尺度動理論。而對于另一類基于CFD-DEM 方法的數(shù)值計算，可通過捕捉顆粒系統(tǒng)中單顆粒的運動特性，選擇合理的統(tǒng)計平均方式，獲得合理的固相本構(gòu)關(guān)系。顆粒系統(tǒng)中應(yīng)力張量的求解方法主要分為三種。（1）體積平均法。Babic[22]最早提出了體積平均法，而后Zhu 等[23]在此基礎(chǔ)上進(jìn)一步發(fā)展，將顆粒質(zhì)量、速度和相互作用力等微觀性質(zhì)與體系密度、應(yīng)力和偶應(yīng)力等宏觀特性聯(lián)系起來。（2）平面法。Mehrabadi 等[24]采用虛功原理應(yīng)用平面法推導(dǎo)出應(yīng)力張量表達(dá)式，并證明平面法推導(dǎo)出的應(yīng)力張量形式與體積平均法等效。（3）粗?；椒╗25]。基于假設(shè)顆粒間為二元碰撞，每一對顆粒都只有單一的一個接觸點，碰撞接觸區(qū)域被一個接觸點取代，且碰撞不是瞬時的，將離散數(shù)據(jù)體積平均化并用連續(xù)性描述。此方法推導(dǎo)應(yīng)力張量表達(dá)式不需要假設(shè)顆粒是剛性的或球形的。綜上所述，顆粒體系內(nèi)的固相應(yīng)力張量的組成形式主要包含兩個部分：(1) 動應(yīng)力(kinetic stress)，由顆粒體系內(nèi)局部體積中顆粒質(zhì)心的速度相對于局部體速度的脈動產(chǎn)生；(2) 接觸應(yīng)力(contact stress)，由體系內(nèi)顆粒之間的直接接觸產(chǎn)生。

當(dāng)顆粒體系內(nèi)有擋板內(nèi)構(gòu)件存在時，內(nèi)構(gòu)件邊壁與顆粒間的碰撞也會對應(yīng)力張量有貢獻(xiàn)作用。為了建立本研究中擋板內(nèi)構(gòu)件表面載荷的統(tǒng)計計算公式，需要了解內(nèi)構(gòu)件的存在對應(yīng)力張量的貢獻(xiàn)形式。Weinhart等[26]采用Goldhirsch[25]的體積平均方法導(dǎo)出了邊界存在所貢獻(xiàn)的邊界應(yīng)力形式。研究發(fā)現(xiàn)，當(dāng)體系內(nèi)有壁面邊界存在下，應(yīng)力張量除了由顆粒速度波動產(chǎn)生的動應(yīng)力和流動顆粒間碰撞產(chǎn)生的碰撞應(yīng)力外，還有邊壁與顆粒碰撞的碰撞貢獻(xiàn)。根據(jù)牛頓第三定律可知，壁面邊界對固相顆粒的碰撞作用與顆粒施加于壁面的碰撞作用數(shù)值上相等。因此，統(tǒng)計擋板表面受到顆粒的碰撞載荷形式與邊界碰撞應(yīng)力的統(tǒng)計形式類似。

統(tǒng)計計算密相床層中的擋板內(nèi)構(gòu)件表面受到顆粒的碰撞力，旨在將微觀的顆粒壁面碰撞信息與宏觀的擋板表面受力信息聯(lián)系起來，并將信息在時間與空間上統(tǒng)計平均得到內(nèi)構(gòu)件表面受到顆粒碰撞應(yīng)力及其分布特性。本研究的擋板內(nèi)構(gòu)件形式為斜片擋板，其受到顆粒的碰撞作用主要表現(xiàn)在與氣體分布板通入氣體的流向相垂直的上、下兩表面所受的顆粒碰撞力。

根據(jù)柯西基本定律，擋板內(nèi)構(gòu)件表面上單位面積受到的顆粒碰撞力與擋板表面面積平均的碰撞應(yīng)力張量關(guān)系如式（1）所示。

式中，v為單位向量；σ為碰撞應(yīng)力張量；T為碰撞力。

基于擋板內(nèi)構(gòu)件表面受到的碰撞載荷與固相應(yīng)力中內(nèi)構(gòu)件邊壁的貢獻(xiàn)數(shù)值上相等，故需要找到一個合理的平均化方法將顆粒的微觀碰撞信息轉(zhuǎn)化為擋板表面受到的宏觀應(yīng)力信息。對體系內(nèi)邊壁存在的影響處理，Ries 等[27]提出了鏡像系統(tǒng)(the mirrored system)的方法。此方法的思想是通過類似將體系內(nèi)邊壁一側(cè)的顆粒關(guān)于邊壁鏡像對稱，來代替邊壁存在對固相應(yīng)力的影響。因此，基于已知擋板內(nèi)構(gòu)件的表面面積，為確定擋板表面面積平均的碰撞應(yīng)力張量σ，可參考此方法處理內(nèi)構(gòu)件表面。

由于本研究所采用的顆粒物料為單一粒徑的顆粒(帶粒徑分布的也可以以同樣的方法統(tǒng)計分析，但是此時兩顆粒質(zhì)心間位移連線的長度l不是常數(shù))，為了更合理有效地統(tǒng)計計算擋板表面受到的碰撞應(yīng)力，將擋板表面緊貼著的一層顆粒關(guān)于擋板表面鏡像對稱，此時與顆粒相互作用的擋板表面等效為一層粒徑與床內(nèi)顆粒相同且平鋪面積與擋板表面面積相同的平面顆粒。因此，擋板表面受到顆粒的碰撞作用等效為單層平面顆粒與虛擬顆粒之間的碰撞作用，可認(rèn)為單層平面顆粒與碰撞顆粒組成的顆粒群間存在一個法向量與擋板表面法向量相一致的“虛擬平面”，如圖1所示。

圖1 密相床層中擋板表面與顆粒接觸示意圖(a)；等效為“虛擬平面”后顆粒碰撞接觸圖示(b)Fig.1 Schematic diagram of contact between the baffle and particles in dense bed(a)and contact between an imaginary plane and particles in a dense bed(b)

水平放置的擋板主要以上、下兩個表面為碰撞受力面，故在此以擋板上表面為例（下表面處理方法一致）。單位面積的“虛擬平面”上存在足夠多數(shù)量的顆粒間碰撞，定義T為單位“虛擬平面”面積上碰撞顆粒對單層平面顆粒施加的碰撞力，則T如式（2）所示。

式中，F(xiàn)i為單位面積“虛擬平面”的上、下第i對顆粒的碰撞接觸力；Nc為單位面積上的碰撞接觸點數(shù)。

Mehrabadi 等[24,28]通過虛功原理證明了單位“虛擬平面”面積上顆粒碰撞力的統(tǒng)計求和可寫成統(tǒng)計平均的顆粒簇體積V上的求和，此采樣顆粒簇由與擋板上表面接觸的真實顆粒和鏡像虛擬擋板顆粒組成。“虛擬平面”上面積平均應(yīng)力張量σ與采樣體積內(nèi)顆粒間相互作用產(chǎn)生的碰撞應(yīng)力張量σ′相等[29-30]，即

因此，由邊壁表面受到的顆粒碰撞作用貢獻(xiàn)形式可知，統(tǒng)計平均的顆粒簇體積中顆粒間碰撞應(yīng)力張量σ的統(tǒng)計計算公式為

式中，上角標(biāo)w 代表壁面;下角標(biāo)α和β表示方向分量。

本研究中，擋板上、下表面的表面積均為S；顆粒粒徑為dp；一對顆粒的每次碰撞接觸中，兩顆粒質(zhì)心間位移連線的長度用l表示，即l的大小為顆粒質(zhì)心到碰撞接觸點間位移長度的2 倍；F為兩顆粒間的碰撞接觸力；Vw為擋板表面相對應(yīng)的均勻統(tǒng)計顆粒簇體積，如圖1(b)中虛線框所示，其值大小為以擋板表面積為底，顆粒粒徑為高的立方體體積，即Vw=Sdp；N為某一時刻擋板表面的總碰撞接觸點數(shù)。因此，擋板表面受到顆粒間碰撞應(yīng)力張量σ統(tǒng)計計算公式為

從式（5）可知，在數(shù)值模擬計算程序中，每一時間步求得擋板表面受到顆粒碰撞的應(yīng)力張量表達(dá)式為

式中，每個分量表示擋板表面受到的均布載荷；σ33為垂直擋板表面方向的碰撞均布載荷，為主要關(guān)注參數(shù)。運用此公式統(tǒng)計分析CFD-DEM 的模擬結(jié)果，即可分析出顆粒與壁面作用的應(yīng)力分布特性。

2 模擬的幾何模型及物理參數(shù)設(shè)置

本研究選取Liu 等[11]的三維方形冷模流化床實驗裝置（圖2）作為模擬研究對象。裝置的主體由方形流化床床體、預(yù)分布器、氣體分布板、單個斜片擋板和旋風(fēng)分離器等組成。方形流化床的截面面積為0.3 m×0.3 m，總高度為5 m。單個斜片由單個板條和兩個固定底座組成，固定底座通過螺栓將板條的兩端固定在流化床的床壁上，使得板條可穩(wěn)定置于密相床層中。其中，擋板板條的長為0.3 m、寬為0.05 m、厚度為0.008 m，且板條的長度可剛好橫跨床層截面，與流化床的截面邊長一致。單個斜片擋板水平安置在流化床邊壁的中心位置處，安裝高度距流化床底部分布板0.5 m。

圖2 三維方形冷模流化床實驗裝置示意圖[11]Fig.2 Schematic diagram of the cold three-dimensional fluidized bed with a square cross-section[11]

為了便于計算，本研究將上述實驗裝置進(jìn)行簡化處理，省略了流化床頂部的旋風(fēng)分離器和底部的氣體分配室。方形流化床幾何模型的長度和寬度與實驗等同，即長和寬均為0.3 m。實際實驗中流化床內(nèi)床層膨脹高度最高膨脹至2 m 左右，考慮到節(jié)省計算資源和提高計算效率，將流化床的床層稀相空間的高度縮減，故方形流化床的模型總高度定為2.2 m。模型中，擋板內(nèi)構(gòu)件的安裝位置與實驗相一致，安裝在距床層底部分布板0.5 m 的流化床壁面中心位置處，且擋板長度保持不變，仍為0.3 m。為簡化網(wǎng)格劃分并成功生成均勻六面體網(wǎng)格(六面體網(wǎng)格尺寸為0.01 m×0.01 m×0.01 m)，故將單個斜片擋板內(nèi)構(gòu)件的寬度調(diào)整為0.06 m，厚度定為0.01 m。而對于本研究中的無擋板自由床而言，除略去擋板流化床內(nèi)斜片擋板外，其余幾何結(jié)構(gòu)均與擋板流化床一致。最終得到模擬采用的方形流化床幾何結(jié)構(gòu)如圖3所示。

圖3 模擬采用的流化床幾何結(jié)構(gòu)Fig.3 The geometry of fluidized bed for simulation

本文選用劉對平[4]實驗中的粒徑為595 μm的B類非球形石英砂顆粒作為模擬體系內(nèi)的固體介質(zhì)，密度等相關(guān)物料基本性質(zhì)與實驗測量值一致。實驗前期未曾測量石英砂顆粒的球形度信息，故根據(jù)測量顆粒的起始流化氣速等參數(shù)信息，采用Hua等[31]提出的方法求解估算顆粒的球形度，最終確定模擬采用的顆粒球形度為0.86。流化床中床層物料的初始自由堆積高度為1 m。模擬體系內(nèi)流體采用常溫常壓空氣（25℃、1.2 kg/m3、1.8×10-5Pa·s）。

3 數(shù)學(xué)模型及模擬設(shè)置

3.1 數(shù)學(xué)模型

本研究所采用的粗?；疌FD-DEM 方法是基于Lu 等[32]提出的EMMS-DPM (energy-minimization multi-scale-discrete particle method，基于能量最小多尺度離散顆粒)方法。故本研究中的粗顆粒質(zhì)量mCGP=k3mp，其中，k為粗?；剩╧=dCGP/dp），mp為真實顆粒的質(zhì)量。顆粒-顆粒碰撞作用與顆粒-壁面碰撞作用的計算使用Peng 等[33]提出的方法，采用不考慮歷史碰撞影響的Hooke（胡克線彈性碰撞）模型求解計算。計算顆粒-顆粒碰撞與顆粒-壁面碰撞之間的相互作用力時，盡管非球形顆粒與顆粒和壁面之間碰撞接觸力的計算，可以通過組合顆粒球元構(gòu)建非球形顆粒幾何模型的方法實現(xiàn)[34-36]，但將大幅增加計算量。因此，本文在計算顆粒-顆粒以及顆粒-壁面間相互作用力時仍然將顆粒視為球形顆粒。顯然，這是一種提高計算精度與保證合理計算量間采取的折中方法，研究表明這種方法可以獲得較為合理的模擬結(jié)果[37-39]。氣固相間曳力的計算采用Hua 等[31]提出的非球形顆粒的曳力模型求解。在Ergun 曳力關(guān)聯(lián)式中引入球形度，并利用Ganser非球形顆粒曳力模型[40]來計算Wen 和Yu 曳力關(guān)聯(lián)式中的單顆粒曳力系數(shù)Cd。

氣固系統(tǒng)DEM 的控制方程組分別為氣相、固相的質(zhì)量守恒和動量守恒方程：

式中，ρg、ug和p分別為流體密度、速度和壓強；t為時間；g為重力加速度；εg為網(wǎng)格空隙率，即計算網(wǎng)格內(nèi)流體的體積分?jǐn)?shù)；NCGP,cell為計算網(wǎng)格中粗顆粒數(shù)目；Vcell為當(dāng)前計算流體網(wǎng)格的體積；Fdrag,i為計算網(wǎng)格內(nèi)顆粒i與氣流之間的曳力交換；τg為牛頓流體黏度應(yīng)力張量；μg為流體的剪切黏度；I為單位張量；uCGP,i為粗顆粒平動速度；dCGP,i為粗顆粒粒徑；N為某一時刻與顆粒i相互作用的顆?？倲?shù);ICGP,i為粗顆粒的轉(zhuǎn)動慣量；ωCGP,i為粗顆粒角速度；分別為法向和切向碰撞接觸力；R為顆粒質(zhì)心到碰撞接觸點的位移矢量；μt為滾動摩擦因數(shù)；RCGP,i為粗顆粒的半徑。

鑒于氣固流化床內(nèi)顆粒之間的碰撞接觸較為頻繁，本研究采用彈簧-阻尼模型對粗顆粒之間的碰撞過程進(jìn)行求解計算[33]，顆粒所受的碰撞接觸力等于法向接觸力與切向接觸力之和。

本研究中固相顆粒為非球形的石英砂顆粒，故采用Hua 等[31]提出的非球形顆粒的曳力模型來確定本研究系統(tǒng)內(nèi)氣固相間曳力。其中具體的曳力系數(shù)可表示為[41]

式中，εp為網(wǎng)格中顆粒的體積分?jǐn)?shù)；Ψ為顆粒的球形度；mi為粗顆粒質(zhì)量；ui為粗顆粒速度；k為粗?；?；Rep為顆粒Reynolds 數(shù)；K1和K2分別為顆粒斯托克斯形狀因子與牛頓形狀因子；dA為顆粒等投影面積球當(dāng)量直徑；dV為顆粒等體積球當(dāng)量直徑，且dA和dV兩參數(shù)難以通過實驗測得，故假設(shè)dA=dV=dsauter為顆粒Sauter 平均直徑；D為流化床的水力直徑；εi為網(wǎng)格內(nèi)i類顆粒的體積分?jǐn)?shù)；di為顆粒的粒徑。由于本研究體系內(nèi)顆粒為單一粒徑顆粒，故di=dave，εi=εp。

3.2 模擬設(shè)置

本研究所需方形氣固擋板流化床與無擋板自由床的三維幾何結(jié)構(gòu)由ANSYS ICEM 軟件進(jìn)行建立，模擬參數(shù)如表1 所示。采用均勻的六面體結(jié)構(gòu)化網(wǎng)格進(jìn)行劃分繪制完成對整個計算床體的網(wǎng)格表征，而后將網(wǎng)格導(dǎo)入開源計算流體力學(xué)軟件OpenFOAM 進(jìn)行求解計算。流體計算網(wǎng)格用于求解CFD 中計算流體的空隙率及其他物理量，對保證數(shù)值計算的精確性和穩(wěn)定性具有重要意義。流體計算網(wǎng)格的尺寸需足夠精細(xì)方可獲得充足的尺度信息[41-42]，且為保證計算資源與時間成本的合理性，故取邊長為10 mm的六面體網(wǎng)格作為后續(xù)模擬研究的基礎(chǔ)，網(wǎng)格尺寸約為真實顆粒粒徑的15 倍，亦是粗?；蟠诸w粒的2～5 倍，故可保證用顆粒中心算法可合理統(tǒng)計計算網(wǎng)格中的空隙率及固相速度等物理量。分別在幾何模型的X、Y和Z三個方向上繪制30×30×220 個均勻流體計算網(wǎng)格。此外，本文所采用的固相顆粒物料石英砂顆粒為Geldart B 類顆粒，劃分的流體計算網(wǎng)格足以獲得網(wǎng)格無關(guān)性的模擬結(jié)果[43]。

表1 擋板床與自由床數(shù)值模擬參數(shù)Table 1 Simulation parameters for the baffled and free fluidized bed

圖3 顯示了本文數(shù)值模擬的邊界條件設(shè)置情況。流化床床層底部作為氣相流體的入口，將其設(shè)置為速度進(jìn)口邊界條件，底部采用均勻進(jìn)氣，且速度大小依據(jù)操作氣速設(shè)定。氣相出口位于反應(yīng)器幾何模型的頂部，并設(shè)置為壓力出口邊界條件，且壓力值與大氣壓相同。為模擬前期實驗中二級旋風(fēng)分離器及濾袋的作用，在反應(yīng)器幾何模型頂部設(shè)定不允許顆粒逸出，以保證模擬體系內(nèi)的顆粒存量恒定。擋板內(nèi)構(gòu)件設(shè)置為與流化床壁面一致，設(shè)定為無滑移的邊壁。為確保計算資源與時間成本的合理性，本研究選取模擬離散顆粒的粗粒化率k=5，且流體計算網(wǎng)格的尺寸是粗?；蟠诸w粒的2～5倍，可保證模擬計算的準(zhǔn)確性[42,44]。根據(jù)流化床內(nèi)真實顆粒數(shù)目確定的粗?；w粒數(shù)為3459863個。

4 模擬結(jié)果與討論

4.1 碰撞受力統(tǒng)計計算公式驗證

擋板表面法向應(yīng)力分量σ33在顆粒與壁面的碰撞作用中占絕對主導(dǎo)地位，其統(tǒng)計結(jié)果的數(shù)量級遠(yuǎn)遠(yuǎn)大于其余8 個分量。因此，擋板表面法向應(yīng)力分量σ33可作為評價擋板表面受顆粒碰撞作用大小的參數(shù)，其余8 個分量可忽略。此外，由前面分析可知，擋板表面法向應(yīng)力分量σ33為擋板上、下表面受到顆粒碰撞的均布載荷差，即擋板受到的顆粒碰撞載荷為擋板下表面受到的碰撞均布載荷與上表面受到的碰撞載荷的差值，即

式中，q為擋板受到的顆粒碰撞均布載荷，方向豎直向上為擋板下表面受到的顆粒碰撞均布載荷；為擋板上表面受到的顆粒碰撞均布載荷。

從圖4 中可以看出，流化床啟動階段內(nèi)構(gòu)件所受碰撞載荷的CFD-DEM 模擬結(jié)果與實驗結(jié)果隨時間演變的趨勢大體上一致，數(shù)值模擬統(tǒng)計的擋板受均布載荷變化曲線顯示在流化氣體剛通入的啟動階段會產(chǎn)生一個較大的脈沖峰值，而后擋板內(nèi)構(gòu)件所受的載荷減小并進(jìn)一步波動變化。與實驗值相比，啟動階段的第一個波峰及波谷值出現(xiàn)的時間更早，這說明應(yīng)用本文所建立的載荷統(tǒng)計計算公式可以定性地成功模擬統(tǒng)計出床層由固定床向流化床轉(zhuǎn)變過程中內(nèi)構(gòu)件受到非正常受力特性這一現(xiàn)象。此外，圖中CFD-DEM 數(shù)值模擬結(jié)果的脈沖峰值與實驗值大小基本一致，說明模擬可以(半)定量復(fù)現(xiàn)實驗中觀測到的氣固流化床反應(yīng)器中密相床層流態(tài)轉(zhuǎn)變時的受力特性。最后，在后續(xù)穩(wěn)定階段載荷波峰的出現(xiàn)與氣泡尾流內(nèi)顆粒與壁面作用緊密相關(guān)。綜上所述，擋板內(nèi)構(gòu)件受顆粒碰撞載荷統(tǒng)計計算公式可以較為準(zhǔn)確地預(yù)測氣固流化床啟動階段內(nèi)構(gòu)件的受力特性。

圖4 流化床啟動階段內(nèi)構(gòu)件所受載荷實驗結(jié)果與CFD-DEM 模擬結(jié)果比較Fig.4 Comparison of the experimental and simulation results of the stress exerted on the baffle during start-up of the fluidized bed

4.2 顆粒碰撞模型參數(shù)對受力特性的影響

考慮到內(nèi)構(gòu)件表面載荷主要由固相顆粒與內(nèi)構(gòu)件表面的碰撞作用貢獻(xiàn)，在CFD-DEM 數(shù)值模擬方法中，離散顆粒與擋板壁面的碰撞作用主要由DEM 統(tǒng)計這些微觀物理量場信息。DEM 方法中，顆粒碰撞模型參數(shù)為碰撞作用的主要影響因素，顆粒楊氏模量、碰撞恢復(fù)系數(shù)、滑動摩擦因數(shù)和滾動摩擦因數(shù)這些碰撞參數(shù)對內(nèi)構(gòu)件表面與顆粒間的碰撞作用密切相關(guān)，并對擋板內(nèi)構(gòu)件表面載荷的影響尤為顯著。因此，考察固相顆粒的楊氏模量、摩擦因數(shù)和恢復(fù)系數(shù)對擋板表面受碰撞力的影響有重要意義。

顆粒楊氏模量是固相顆粒物性的重要參數(shù)，表征了顆粒表面的硬度。楊氏模量的大小取決于固相顆粒自身的性質(zhì)，其數(shù)值越大顆粒越不容易發(fā)生變形，反映了固相顆粒材料的剛性[45]。此外，在計算機數(shù)值模擬方法中，離散顆粒的楊氏模量作為模型碰撞參數(shù)是CFD-DEM 數(shù)值模擬計算的一個重要參數(shù)。顆粒楊氏模量很大程度上影響數(shù)值模擬的計算量及計算效率，其決定了DEM 中計算的時間步長。在研究氣固兩相間的時均流動特性時，顆粒的楊氏模量對統(tǒng)計計算的流動參數(shù)無顯著影響，可以將楊氏模量數(shù)值設(shè)置小些以減少計算量，提高計算效率[46]。

對于密相床層的啟動階段，體系內(nèi)顆粒-擋板壁面間的碰撞載荷統(tǒng)計則與楊氏模量的大小密切相關(guān)。基于本研究所選用的Hooke 模型，法向彈性系數(shù)kn與楊氏模量直接相關(guān)。因此，為了探究內(nèi)構(gòu)件的受力特性，需要研究楊氏模量的大小對氣固流化床啟動階段密相床層中擋板內(nèi)構(gòu)件表面受到的碰撞載荷的影響。固相顆粒物料仍為相同的B類石英砂顆粒，其楊氏模量的取值范圍介于0.1 GPa 和10 GPa 量級之間，其中前者為研究氣固流動特性統(tǒng)計時均參數(shù)常選用的經(jīng)驗值，而后者是真實石英砂顆粒的楊氏模量量級，而擋板內(nèi)構(gòu)件的楊氏模量與顆粒的楊氏模量設(shè)置相同。本小節(jié)選用的顆粒楊氏模量具體數(shù)值及對應(yīng)DEM 與CFD 中的時間步如表2 所示，模擬計算體系及其余模擬參數(shù)設(shè)置與前面小節(jié)模擬工況的設(shè)置一致。

表2 擋板床數(shù)值模擬參數(shù)Table 2 Simulation parameters for the baffled fluidized bed

圖5顯示了固相顆粒選取不同楊氏模量下密相床層中擋板內(nèi)構(gòu)件表面碰撞載荷隨時間的演化情況。發(fā)現(xiàn)流化床自啟動開始至出現(xiàn)第一個完整峰值脈沖載荷階段，因顆粒楊氏模量的不同，峰值載荷的大小與峰值出現(xiàn)的時間均受到影響發(fā)生相應(yīng)變化，而后載荷的變化不再受顆粒楊氏模量的影響。對圖中啟動階段出現(xiàn)的第一個載荷峰值位置局部放大，如圖5 內(nèi)附圖所示。（1）在顆粒初始化自由堆積成密相床層的階段(0 s 之前)，顆粒堆積過程中體系逐漸密實，新力鏈不斷生成。當(dāng)床層高度穩(wěn)定，此時體系足夠密實，力鏈充分發(fā)展幾乎不發(fā)生變化，擋板內(nèi)構(gòu)件表面載荷幾乎穩(wěn)定不變[47]。顆粒堆積的隨機性與堆積時間的差異使得未通入氣體前擋板上、下表面受到的顆粒的擠壓碰撞作用不同，致使擋板內(nèi)構(gòu)件表面載荷存在一定的差異性。但是大體上，較大的顆粒楊氏模量使得顆粒-壁面間的法向接觸力增大，致使固定床狀態(tài)時擋板受到的碰撞作用載荷更大。（2）當(dāng)體系內(nèi)通入流化空氣，啟動階段開始至第一個峰值脈沖出現(xiàn)(0～0.12 s)，峰值出現(xiàn)的時間隨楊氏模量的增加而縮減，峰值載荷的大小隨顆粒楊氏模量的增大而減小，且當(dāng)楊氏模量增大至10 GPa 量級時，此影響不再顯著，峰值大小不再發(fā)生明顯變化。說明顆粒楊氏模量的增加，加快了由固定床轉(zhuǎn)化為流化床的速度，楊氏模量大的顆粒床層內(nèi)擋板表面載荷率先達(dá)到峰值。但由于楊氏模量增大，顆粒的表面硬度更大、變形程度減小，致使顆粒的碰撞重疊量減小，碰撞接觸時間明顯縮短，使得內(nèi)構(gòu)件表面的碰撞作用載荷減小。此外，楊氏模量大的顆粒床層因顆粒硬度大、變形小，在流體通入時會更容易打破“自鎖”狀態(tài)，固相顆粒較快發(fā)生松動而流化，可能也是使得擋板表面受碰撞載荷峰值較小的原因。（3）當(dāng)擋板表面受到的第一個峰值脈沖結(jié)束后(0.12 s 之后)，內(nèi)構(gòu)件表面載荷變化基本不受顆粒楊氏模量的影響。體系進(jìn)入正常流化階段后，顆粒楊氏模量大的床層中顆粒-壁面的碰撞作用更劇烈，發(fā)生碰撞更為頻繁，以及在彈性系數(shù)增加的共同影響下，在一定程度上抵消了因碰撞重疊量減小而使得內(nèi)構(gòu)件表面載荷減小的趨勢。

圖5 不同顆粒楊氏模量下顆粒床層內(nèi)擋板受到載荷隨時間的變化Fig.5 Effect of particle Young's modulus on the stress exerted on the baffle immersed in the fluidized bed

流化床啟動階段，未流化密相床層會被氣節(jié)推起一定高度，使得擋板內(nèi)構(gòu)件受到較大的碰撞作用。而在這一短暫過程中，顆粒間的碰撞作用較為強烈，碰撞重疊量會更大，為保證模擬計算的精度，統(tǒng)計床層啟動階段顆粒最大重疊量與顆粒粒徑的比值隨時間的變化情況，如圖6 所示。從圖中可以看出，每個工況下顆粒最大重疊量與粒徑的比值都小于2%，說明啟動階段模擬計算的準(zhǔn)確性。此外，隨著顆粒楊氏模量的增大，最大重疊量與粒徑的比值減小，也進(jìn)一步說明了顆粒硬度增大，難以發(fā)生變形的特性。

圖6 不同楊氏模量下顆粒最大重疊量與粒徑比值隨時間的變化Fig.6 Variation of the ratio of the maximum particle-particle overlap to particle size with time under different Young's modulus

顆粒碰撞恢復(fù)系數(shù)是顆粒的重要物性參數(shù)之一，可以反映顆粒碰撞后恢復(fù)到變形前初始狀態(tài)的能力，即可以反映顆粒在接觸碰撞過程中能量耗散的強弱[48]。顆粒摩擦因數(shù)是固相顆粒物料的一個基本物性參數(shù)，與接觸面的粗糙程度有關(guān)。顆粒-顆粒與顆粒-壁面之間的摩擦?xí)念w粒的動能，進(jìn)而影響體系內(nèi)氣固相間的交換量以及固相顆粒在切向方向上的運動情況[49]。顆粒滾動摩擦因數(shù)作為固相顆粒物料的重要物理力學(xué)參數(shù)之一，也會影響顆粒體系內(nèi)形成的微觀力學(xué)結(jié)構(gòu)。顆粒間的滾動摩擦因數(shù)主要影響顆粒系統(tǒng)的堆積形態(tài)[50]，而在DEM 碰撞模型中，顆粒滾動摩擦因數(shù)的取值一般遠(yuǎn)小于固相顆粒的直徑[51]。因此，依據(jù)單一變量原則，分別選取0.80、0.85、0.90 這三個值作為模擬所需的顆粒碰撞恢復(fù)系數(shù)；選取顆粒的滑動摩擦因數(shù)為0.30、0.35 和0.40；取顆粒的滾動摩擦因數(shù)為0.01、0.03 和0.05。其余參數(shù)的設(shè)置與之前保持一致，模擬結(jié)果如圖7所示，可以看出碰撞恢復(fù)系數(shù)、滑動摩擦因數(shù)和滾動摩擦因數(shù)對擋板內(nèi)構(gòu)件的受力影響并不顯著。

圖7 不同顆粒碰撞恢復(fù)系數(shù)(a)、滑動摩擦因數(shù)(b)、滾動摩擦因數(shù)(c)下顆粒床層內(nèi)擋板受到載荷隨時間的變化Fig.7 Effect of particle restitution coefficient(a),particle sliding friction coefficient(b),and particle rolling friction coefficient(c)on the stress exerted on the baffle immersed in the fluidized bed

4.3 表觀氣速對受力特性的影響

流化床反應(yīng)器的流化風(fēng)速是控制流化床運行及氣固兩相流動的重要操作條件。因此，研究表觀氣速對擋板內(nèi)構(gòu)件受力特性的影響具有重要意義?；?.2 節(jié)關(guān)于碰撞模型參數(shù)對受力特性的影響研究，現(xiàn)選定一組關(guān)于顆粒屬性的最佳模型參數(shù)，即顆粒的楊氏模量為1×1010Pa，DEM 中時間步長為1×10-6s，CFD 中時間步長為1×10-5s，碰撞恢復(fù)系數(shù)為0.90，滑動摩擦因數(shù)為0.30，滾動摩擦因數(shù)為0.01。根據(jù)固相顆粒的屬性及最小流化氣速，現(xiàn)選取用于模擬研究的表觀氣速的一組值0.5、0.6、0.9 m/s。

圖8 為不同表觀氣速下，流化床啟動階段密相床層中擋板內(nèi)構(gòu)件受到的碰撞載荷瞬態(tài)變化情況。模擬結(jié)果顯示，對于初始堆積床層，表觀氣速增大，相當(dāng)于給予流化床的氣速增量增加，擋板內(nèi)構(gòu)件的受力載荷強度也相應(yīng)增大。當(dāng)表觀氣速增加的幅度越大，內(nèi)構(gòu)件表面載荷的數(shù)值增長的幅度也越大。出現(xiàn)這一現(xiàn)象的原因在于表觀氣速增大，氣速增量越大，流化床底部形成的氣節(jié)高度更高，床層中固相顆粒的動量越大，對擋板內(nèi)構(gòu)件表面的撞擊更為劇烈，因此載荷強度在數(shù)值上表現(xiàn)出更大的情況。說明在流化床床層啟動階段，表觀氣速的大小是影響擋板內(nèi)構(gòu)件受力的主要因素之一。

圖8 不同表觀氣速下顆粒床層內(nèi)擋板受到載荷隨時間的變化Fig.8 Effect of superficial gas velocity on the stress exerted on the baffle immersed in the fluidized bed

上述模擬結(jié)果說明在實際流化床啟動操作過程中，設(shè)置較小的表觀氣速可明顯減小內(nèi)構(gòu)件受到的碰撞力。因此，要避免在較大的氣速通入下完成床層啟動操作，防止擋板內(nèi)構(gòu)件因瞬時受到較大的作用力而發(fā)生強度破壞，導(dǎo)致反應(yīng)器裝置無法正常運行。

4.4 顆粒粒徑對受力特性的影響

流化床體系內(nèi)，固相顆粒粒徑的改變可使系統(tǒng)的最小流化速度發(fā)生變化，從而影響體系內(nèi)的流態(tài)轉(zhuǎn)變及顆粒碰撞動力學(xué)[49]。因此，采用不同粒徑的顆粒物料作為床料進(jìn)行數(shù)值模擬，基于單一變量原則，在表觀氣速分別按統(tǒng)一設(shè)定的工況和以最小流化氣速倍數(shù)設(shè)定的工況下，研究顆粒粒徑變化對密相床層中擋板內(nèi)構(gòu)件受力特性的影響。模擬研究所選用的固相顆粒僅粒徑不同，其余的顆粒物性參數(shù)均一致（如顆粒數(shù)目）。

基于單一變量原則，將流化床內(nèi)的流化氣速采用統(tǒng)一設(shè)定，并將其設(shè)置為0.6 m/s，以此來研究相同表觀氣速下顆粒尺寸對內(nèi)構(gòu)件受力特性的影響。本節(jié)關(guān)于顆粒屬性的模型參數(shù)設(shè)置與4.3 節(jié)一致，其余參數(shù)參考表1設(shè)置。

圖9 顯示了表觀氣速統(tǒng)一設(shè)定為0.6 m/s 工況下，不同粒徑的顆粒床層啟動階段密相床層中擋板內(nèi)構(gòu)件表面載荷強度的瞬態(tài)演變情況。從圖中可看出隨著顆粒粒徑的增大，內(nèi)構(gòu)件表面受的碰撞載荷強度整體有減小趨勢，尤其在第一個完整峰值載荷出現(xiàn)的區(qū)間，載荷強度有明顯減小的趨勢，在0.06 s 位置表現(xiàn)更為明顯。出現(xiàn)這一現(xiàn)象的原因可能在于，隨著顆粒尺寸的增加，體系內(nèi)的最小流化氣速相應(yīng)增加，而在相同表觀氣速下粒徑較大的固相顆粒床層易出現(xiàn)流化不良的現(xiàn)象，床層中氣泡的尺寸相對較小，氣固兩相的運動及相間作用程度減弱，擋板內(nèi)構(gòu)件受到的載荷作用也相應(yīng)減小，這與前期實驗中的現(xiàn)象基本一致。

5 結(jié) 論

本研究基于CFD-DEM 數(shù)值模擬方法，基于劉對平[4]前期實驗所用的三維方形冷模流化床裝置作為模擬體系，以工業(yè)上廣泛應(yīng)用的單個斜片擋板內(nèi)構(gòu)件作為研究對象，建立了擋板內(nèi)構(gòu)件表面受力載荷強度的統(tǒng)計計算方法，并借助此統(tǒng)計計算公式分析研究了氣固物性對擋板床啟動階段密相床層中單個水平斜片擋板內(nèi)構(gòu)件的受力特性的影響規(guī)律。主要結(jié)論如下。

（1）基于CFD-DEM 模擬方法建立了氣固流化床中浸沒在密相床層的擋板內(nèi)構(gòu)件表面受顆粒碰撞載荷分布的統(tǒng)計計算公式，可用于定量分析內(nèi)構(gòu)件表面受到的動態(tài)碰撞載荷信號。通過與課題組前期實驗測量數(shù)據(jù)對比，建立的內(nèi)構(gòu)件表面載荷統(tǒng)計計算公式可以定性和(半)定量復(fù)現(xiàn)實驗中流化床啟動階段內(nèi)構(gòu)件受力特性的變化，為研究內(nèi)構(gòu)件在床層啟動階段流型轉(zhuǎn)變過程中出現(xiàn)這一非正常受力現(xiàn)象的內(nèi)在機理奠定基礎(chǔ)。

（2）顆粒碰撞模型參數(shù)中，楊氏模量是影響內(nèi)構(gòu)件受力特性的關(guān)鍵性因素，而碰撞恢復(fù)系數(shù)、滑動摩擦因數(shù)和滾動摩擦因數(shù)對擋板內(nèi)構(gòu)件的受力影響并不顯著。研究發(fā)現(xiàn)流化床自啟動開始至出現(xiàn)第一個完整峰值脈沖載荷階段，不同的顆粒楊氏模量下，內(nèi)構(gòu)件受力的峰值載荷大小與峰值出現(xiàn)的時間均受到影響發(fā)生相應(yīng)變化（隨著楊氏模量的增大，峰值載荷先逐漸變小后趨于穩(wěn)定），而后階段載荷的變化幾乎不再受顆粒楊氏模量的影響。此外，在受楊氏模量影響階段，隨著楊氏模量的增大，第一個完整峰值載荷相應(yīng)減小，當(dāng)楊氏模量增大至顆粒物理真實值，峰值載荷不再發(fā)生變化。

（3）氣固流化床體系內(nèi)表觀氣速是影響內(nèi)構(gòu)件受力的主要因素之一。表觀氣速越大，密相床層中的擋板內(nèi)構(gòu)件受力載荷強度越大，與實驗結(jié)論一致。為保障工業(yè)上反應(yīng)器裝置的正常運行和內(nèi)部構(gòu)件及其支撐結(jié)構(gòu)的穩(wěn)定性，要在相對較小的氣速增幅下完成流化床的啟動操作過程。