顆粒沉降問題的高解析度CFD- DEM- IBM流固耦合數(shù)值模擬方法

底瑛棠,趙蘭浩,毛 佳

(河海大學(xué)水利水電學(xué)院,江蘇 南京 210098)

顆粒沉降現(xiàn)象在自然界中普遍存在,具有廣泛的工程應(yīng)用背景,如泥沙運(yùn)移[1]、血液流動(dòng)[2]、氣力輸送[3]等。無(wú)論從工程應(yīng)用價(jià)值或基礎(chǔ)研究角度,顆粒與流體間的相互作用機(jī)理研究都意義重大。

顆粒沉降涉及顆粒間的碰撞及流固之間的雙向強(qiáng)烈耦合,其高精度數(shù)值模擬目前仍極具挑戰(zhàn)。描述顆粒兩相流的數(shù)值方法主要有歐拉- 歐拉(Eulerian- Eulerian,E- E)方法[4- 5]和歐拉- 拉格朗日(Eulerian- Lagrangian,E- L)方法[6,8]。然而,E- E方法將流固兩相視作可以相互滲透的連續(xù)介質(zhì),無(wú)法考慮顆粒的離散特性,實(shí)際更適用于固體呈現(xiàn)顯著流態(tài)化特征的流動(dòng)。為考慮顆粒的非連續(xù)性,張濤和李紅文[6]基于E- L方法采用離散相模型(Discrete Phase Model,DPM)模擬了氣固兩相流管道內(nèi)的孔板測(cè)量工況,將顆粒視作無(wú)體積的質(zhì)點(diǎn),通過牛頓第二定律描述顆粒運(yùn)動(dòng)。DPM模型假定顆粒相非常稀疏,忽略顆粒的相互作用,不考慮其對(duì)流體運(yùn)動(dòng)的影響,實(shí)則只是單向耦合,適用于固相體積分?jǐn)?shù)小于10%～12%的情況[7]。針對(duì)稠密顆粒流,Hwang等[8]考慮了顆粒的體積效應(yīng)與粒間碰撞及相間耦合,采用密集離散相模型(Dense Discrete Phase Model,DDPM)模擬了旋轉(zhuǎn)分離器中的流動(dòng)特性。但DDPM通過顆粒動(dòng)力學(xué)理論計(jì)算顆粒間的相互作用,模化碰撞過程,并無(wú)清晰的物理碰撞面。

近年來(lái),計(jì)算流體力學(xué)與離散單元法(Computational Fluid Dynamics- Discrete Element Method,CFD- DEM)因其處理顆粒碰撞問題的優(yōu)勢(shì),被越來(lái)越多地用于顆粒兩相流的研究。Yue等[9]基于CFD- DEM方法研究了噴動(dòng)流化床中顆粒密度對(duì)噴流偏轉(zhuǎn)的影響;楊文婧等[10]采用CFD- DEM方法對(duì)固體火箭發(fā)動(dòng)機(jī)中氣固兩相流進(jìn)行了數(shù)值模擬。CFD- DEM方法能夠考慮兩相之間的耦合作用,引入顆粒接觸模型描述顆粒的物理碰撞過程,并通過空隙率考慮顆粒的體積效應(yīng),是求解顆粒兩相流宏觀特性的極佳選擇。然而,傳統(tǒng)的CFD- DEM方法通過經(jīng)驗(yàn)公式計(jì)算兩相間的耦合作用力[11],稱作非解析的CFD- DEM方法。由于缺乏對(duì)顆粒邊界的描述,顆粒周圍的流場(chǎng)信息無(wú)法得到精確解析,從而無(wú)法描述顆粒的尾流效應(yīng)[12],決定了該方法只適用于顆粒間無(wú)干擾的稀疏顆粒兩相流問題。

為解決目前非解析CFD- DEM方法的弊端,準(zhǔn)確求解流固耦合作用力,獲得解析的流場(chǎng),主要的難點(diǎn)在于流固運(yùn)動(dòng)界面的描述。傳統(tǒng)的動(dòng)網(wǎng)格方法[13]需要在每個(gè)時(shí)間步內(nèi)更新計(jì)算域網(wǎng)格,考慮到顆粒沉降運(yùn)動(dòng)過程中的大位移,采用貼體動(dòng)網(wǎng)格方法處理顆粒的運(yùn)動(dòng)邊界計(jì)算成本較高,并可能產(chǎn)生網(wǎng)格扭曲等問題。Peskin[14]提出的浸入邊界法(IBM)采用固定的結(jié)構(gòu)化網(wǎng)格,在固體邊界上覆蓋拉格朗日點(diǎn)以追蹤其運(yùn)動(dòng)邊界,回避了動(dòng)網(wǎng)格方法的缺陷。通過在Navier- Stokes動(dòng)量方程中附加體力項(xiàng)來(lái)實(shí)現(xiàn)固體邊界對(duì)流場(chǎng)的影響,浸入邊界法不僅無(wú)需實(shí)時(shí)更新網(wǎng)格,節(jié)省了網(wǎng)格生成所需的時(shí)間,也無(wú)需處理從物理平面到計(jì)算平面的坐標(biāo)和網(wǎng)格的轉(zhuǎn)換問題,正逐步被應(yīng)用到更廣泛的工程領(lǐng)域[15],對(duì)于稠密顆粒系統(tǒng)的沉降這種動(dòng)態(tài)邊界問題尤為適用。

本文基于CFD- DEM和IBM方法提出一種高解析度CFD- DEM- IBM流固耦合數(shù)值分析方法,基于自主開發(fā)的程序開展顆粒沉降問題的數(shù)值模擬研究,通過與前人工作的對(duì)比驗(yàn)證了該方法能夠描述顆粒的實(shí)際邊界,獲得完全解析的流場(chǎng),從而能夠考慮顆粒尾渦結(jié)構(gòu)的影響,為深入了解顆粒兩相流運(yùn)動(dòng)特性提供了有力工具。

1 數(shù)值方法

1.1 流體控制方程

連續(xù)的流體相由黏性不可壓縮的Navier- Stokes方程作為控制方程,為考慮顆粒邊界的影響,在動(dòng)量方程的右端附加額外的體力項(xiàng):

·un+θ1=0

(1)

(2)

式中:u為速度；t為時(shí)間；ρ為密度；p為壓力；τ=μ(u+Tu)為黏性應(yīng)力張量,μ為動(dòng)力黏度；fb為體力項(xiàng)；f為附加體力項(xiàng),表示固體邊界對(duì)流體的作用力；n為時(shí)間步；θ1和θ2的取值范圍均為[0,1],時(shí)間離散格式當(dāng)0.5≤θ1≤1且θ2=0時(shí)為顯式,當(dāng)0.5≤θ1≤1且0.5≤θ2≤1時(shí)為半隱式格式。

1.2 固體控制方程

固體的運(yùn)動(dòng)由牛頓第二運(yùn)動(dòng)方程控制:

(3)

(4)

式(3)中,接觸力包括法向接觸力Fcn和切向接觸力Fct,分別由下式計(jì)算得出:

(Fcn)i=(Fcn)i-1+ΔFcn

(5)

(Fct)i=(Fct)i-1+ΔFct

(6)

式中:下標(biāo)i表示ti時(shí)刻;Δt=ti-ti-1為時(shí)間步;力的增量滿足ΔFcn=knδn+ηnGn及ΔFct=ktδt+ηtGt;k為剛度系數(shù)；δ為接觸點(diǎn)的相對(duì)位移增量；η為阻尼系數(shù)；G為接觸顆粒的相對(duì)速度;下標(biāo)n和t分別表示法向和切向分量,其中切向力應(yīng)當(dāng)滿足庫(kù)侖摩擦定律(Fct)max=Fcn·μc,(Fct)max為最大允許靜摩擦力,μc為摩擦系數(shù)。

流體對(duì)顆粒的作用力Ff為施加于浸入邊界點(diǎn)上的力,與作用于流體網(wǎng)格點(diǎn)上的附加體力項(xiàng)f含義不同。Ff可由參考文獻(xiàn)[16]計(jì)算得到:

(7)

式中:Ω為固體邊界SP外、計(jì)算域邊界Sout內(nèi)的流體域；Ω′為整個(gè)計(jì)算域；Ω″為SP所圍固體域,如圖1所示。

圖1 顆粒沉降問題的CFD- DEM- IBM方法示意

1.3 浸入邊界法

為處理顆粒的移動(dòng)邊界,引入直接力浸入邊界法[17],如圖1所示。浸入邊界法采用基于歐拉框架的正交笛卡爾網(wǎng)格離散整個(gè)計(jì)算域,并在顆粒的物理邊界上布置浸入邊界點(diǎn),流體網(wǎng)格點(diǎn)與浸入邊界點(diǎn)之間的信息傳遞通過插值函數(shù)I與分布函數(shù)D完成。

將動(dòng)量方程式(2)沿特征線展開可得到時(shí)間離散后的動(dòng)量方程:

un+1=un+Kn+Rn+θ2+fn+1Δt

(8)

fn+1Δt=D[Vn+1-I(un+Kn+Rn+θ2)]

(9)

1.4 數(shù)值求解方案

1.4.1 流體相的迭代求解

由式(9)可見,流體的速度、壓力與附加體力項(xiàng)相互耦合,需要采用迭代方案求解流場(chǎng)。根據(jù)分步投影法,定義不包含壓力項(xiàng)與附加體力項(xiàng)的中間速度場(chǎng)u*,k=un+Kn,流場(chǎng)的求解過程如圖2所示。

圖2 流體域迭代求解步驟

其中,內(nèi)循環(huán)中附加體力項(xiàng)根據(jù)下式求解:

fn+1,k,iΔt=fn+1,k,i-1Δt+D[Vn+1-I(un+Kn+Rn+θ2+fn+1,k,i-1Δt)]

(10)

內(nèi)循環(huán)的收斂判斷則依據(jù)為‖Vn+1-I(un+Kn+Rn+θ2+fn+1,k,i-1Δt)‖<ε,其中‖·‖為任一范數(shù),容差ε取為1×10-5。外循環(huán)的收斂判斷依據(jù)為‖I(un+Kn+Rn+θ2,k)-I(un+Kn+Rn+θ2,k-1)‖<ε。如圖2所示,采用迭代的方式隱式求解附加體力項(xiàng),在保持流體的速度、壓力與附加體力項(xiàng)相互耦合的同時(shí),滿足無(wú)滑移邊界條件,并在附加體力與壓力項(xiàng)的迭代計(jì)算內(nèi)進(jìn)行速度校正,確保結(jié)果滿足散度自由。

1.4.2 流固耦合系統(tǒng)的迭代求解

顆粒與流體間存在雙向相互作用,形成強(qiáng)流固耦合系統(tǒng),采用分區(qū)方案求解分屬不同控制方程的流固兩相。根據(jù)標(biāo)準(zhǔn)Galerkin加權(quán)余量法進(jìn)行空間離散,流固兩相的控制方程迭代格式可記作:

(Δun,k,Δpn,k)=RHSf(un,pn,δn+1,k-1)

(11)

Δδn,k=RHSs(δn,un+1,k,pn+1,k)

(12)

式中:RHSf與RHSs分別為流體方程與固體方程的右端項(xiàng)。

現(xiàn)假設(shè)第n步計(jì)算已完成,依據(jù)以下流程展開第n+1步的第k次迭代:

(1) 求解流體方程式(11),更新速度場(chǎng)un+1,k及壓力場(chǎng)pn+1,k；

(2) 求解固體方程式(12),更新固體位置δn+1,k；

(3) 計(jì)算‖Δζn,k-Δζn,k-1‖判斷收斂,若小于容差εζ,則滿足收斂條件,結(jié)束當(dāng)前時(shí)間步的計(jì)算,否則令k=k+1,重復(fù)(1)—(3)直至滿足收斂條件。式中ζ可取為變量u、p或位移δ,εζ取為1×10-5。

2 數(shù)值算例

為驗(yàn)證本文方法的準(zhǔn)確性與有效性,進(jìn)行單、雙、群顆粒的沉降行為模擬,并與前人數(shù)值計(jì)算結(jié)果展開對(duì)比。為凸顯高解析度CFD- DEM- IBM方法的優(yōu)勢(shì),非解析CFD- DEM方法的計(jì)算結(jié)果也在文中給出。其中,非解析方法中流固相互作用力依據(jù)式(13)計(jì)算:

(13)

式中:Cd為阻力系數(shù)；ρf為流體密度；di為顆粒pi的直徑；uf,i為顆粒pi周圍網(wǎng)格點(diǎn)插值得到的流體速度；vi為顆粒pi的速度；εi為pi處的流體體積分?jǐn)?shù)；χ=3.7-0.65exp[-(1.5-lgRei)2/2],為與顆粒雷諾數(shù)Rei有關(guān)的參數(shù)。

非解析方法要求流體網(wǎng)格的尺寸遠(yuǎn)大于顆粒尺寸[18],考慮到本文算例的具體參數(shù),采用非解析方法計(jì)算時(shí),網(wǎng)格尺寸取為約2.5倍顆粒直徑。

2.1 單顆粒沉降

為驗(yàn)證本文提出的高解析度CFD- DEM- IBM耦合模型用于分析顆粒沉降問題的正確性與可靠性,首先選擇經(jīng)典的單顆粒沉降問題[19- 20]展開計(jì)算。計(jì)算域?yàn)? cm×6 cm的矩形域,域內(nèi)充滿密度為ρf=1 g/cm3、黏度為μf=0.1 g/(cm·s)的牛頓流體。初始時(shí)刻顆粒位于計(jì)算域中線位置,從距離底面4 cm處由靜止受重力下落,顆粒直徑dp=0.25 cm,密度ρp=1.25 g/cm3,其法向及切向剛度均取為2×107N/m。整個(gè)計(jì)算域采用規(guī)則的矩形網(wǎng)格離散,網(wǎng)格尺寸為1/256 cm,顆粒邊界上布置402個(gè)浸入邊界點(diǎn),重力加速度g=9.81 m/s2,除上部邊界設(shè)置為開邊界外,其余三邊均設(shè)置為無(wú)滑固壁邊界。

表1 最大雷諾數(shù)對(duì)比

圖3對(duì)比了單顆粒沉降的垂向位置與垂向速度。由圖3(b)可見,顆粒沉降過程的加速沉降、勻速沉降、觸底減速沉降3個(gè)階段都得到了很好的再現(xiàn)。值得注意的是,采用傳統(tǒng)非解析的CFD- DEM方法的計(jì)算結(jié)果也呈現(xiàn)在圖中,可以看出非解析方法的結(jié)果雖與其他方法的結(jié)果具有一定差異,但顆粒沉降過程中3個(gè)階段的典型流態(tài)、顆粒的穩(wěn)定沉降速度等均呈現(xiàn)出較好的一致性。

圖3 單顆粒沉降垂向位置與垂向速度時(shí)程

實(shí)際上,非解析CFD- DEM方法要求流體網(wǎng)格的尺寸需遠(yuǎn)大于顆粒尺寸[18],更適合于規(guī)模更大的顆粒沉降問題。而在此算例中,由于顆粒尺寸一定而計(jì)算域有限,流體網(wǎng)格尺寸難以達(dá)到理想要求,勢(shì)必造成計(jì)算結(jié)果與前人研究的差異?？偟膩?lái)說(shuō),本文方法與以往研究結(jié)果吻合度較高,證明了本文提出的CFD- DEM- IBM方法處理顆粒沉降問題的準(zhǔn)確性與可靠性；另外,對(duì)于稀疏的顆粒兩相流系統(tǒng),無(wú)論是非解析還是解析的CFD- DEM方法均可適用,能夠描述顆粒的沉降行為。

2.2 雙顆粒沉降的DKT現(xiàn)象

為初步驗(yàn)證本文提出的方法對(duì)于稠密顆粒兩相流問題的適用性、高解析度描述流場(chǎng)與顆粒行為的能力,進(jìn)行2個(gè)相同圓形顆粒自由下落問題研究。計(jì)算域與顆粒大小及網(wǎng)格尺寸等參數(shù)與2.1節(jié)中取相同,不同之處為本節(jié)中流體黏度μf取為0.01 g/(cm·s),顆粒密度ρp取為1.5 g/cm3。如圖4(a),初始時(shí)刻,兩顆粒分別位于((1+0.001) cm,5 cm)及((1-0.001) cm,4.5 cm)處,由靜止受到重力作用開始下落。

圖4 雙顆粒沉降流速矢量

在兩同軸顆粒的沉降過程中,初始位于上方的拖尾顆粒將會(huì)超越下方引導(dǎo)的顆粒,發(fā)生經(jīng)典的追趕—接觸—翻滾(Drafting- Kissing- Tumbling,DKT)現(xiàn)象[21]。如圖4所示,下方顆粒沉降過程中產(chǎn)生的低壓尾渦結(jié)構(gòu)作用于上方的顆粒,造成上方顆粒受到的阻力降低從而加速下落(圖4(b)),直至與下方顆粒接觸(圖4(c))，最終兩顆粒由于接觸結(jié)構(gòu)的不穩(wěn)定而分開(圖4(d))。可見,本文方法能夠?qū)⒘鲌?chǎng)高度解析,正確再現(xiàn)了DKT現(xiàn)象,與文獻(xiàn)[21]中的結(jié)果具有很好的一致性。

為定量對(duì)比計(jì)算結(jié)果,圖5將雙顆粒沉降垂向位置與水平位置時(shí)程圖與前人研究結(jié)果展開對(duì)比。圖中,p1為初始時(shí)刻位于上方的顆粒,p2為初始位于下方的顆粒。由圖可見,在0.2 s前,采用本文所提出的CFD- DEM- IBM方法進(jìn)行計(jì)算得到的水平與垂向位置與前人[21]的結(jié)果都保持高度一致。考慮到翻滾現(xiàn)象主要取決于擾動(dòng)的發(fā)展,受到不同方法實(shí)施過程中數(shù)值誤差的影響,極具不穩(wěn)定性,在翻滾階段產(chǎn)生的現(xiàn)象可能會(huì)有所不同,這一流場(chǎng)特性將影響顆粒的位置,因此,圖5中0.2 s后的結(jié)果差異的產(chǎn)生是合理的。

圖5 雙顆粒沉降垂向位置與水平位置時(shí)程

另外,由圖5可見,將非解析方法的計(jì)算結(jié)果與本文方法的計(jì)算結(jié)果對(duì)比發(fā)現(xiàn),非解析CFD- DEM方法計(jì)算所得結(jié)果完全屬于兩顆粒單獨(dú)沉降現(xiàn)象的累加,未能反映出兩顆粒間的相互作用,不能考慮到顆粒尾流的影響。對(duì)于此類顆粒間存在相互干擾的問題,通過經(jīng)驗(yàn)公式計(jì)算流固耦合作用力做法是不準(zhǔn)確的,因此,非解析的CFD- DEM方法不適用于稠密顆粒兩相流問題的研究。

2.3 顆粒群沉降

為進(jìn)一步驗(yàn)證本文所提出的CFD- DEM- IBM方法處理稠密顆粒兩相流問題的能力,本節(jié)展開了504個(gè)圓形顆粒在密閉容器內(nèi)的沉降過程[22]的數(shù)值模擬。封閉的矩形域大小為2 cm×2 cm,其內(nèi)部流體的密度為ρf=1 g/cm3,黏度為μf=0.01 g/(cm·s)。504個(gè)顆粒具有相同的直徑dp=0.062 5 cm,密度為ρp=1.01 g/cm3,計(jì)算域內(nèi)顆粒的體積分?jǐn)?shù)為38.66%,初始時(shí)刻504個(gè)顆粒的位置分布如圖6(a)所示。計(jì)算域采用1/256 cm的網(wǎng)格進(jìn)行劃分,每個(gè)顆粒表面布置150個(gè)浸入邊界點(diǎn),計(jì)算域的4個(gè)邊界均設(shè)置為無(wú)滑固壁。

圖6 504個(gè)顆粒沉降速度矢量

圖6給出了504個(gè)顆粒沉降過程中的速度矢量圖?？梢?沉降開始時(shí),在重力作用的影響下所有顆粒都均勻下落；由于壁面的阻滯效應(yīng),靠近兩側(cè)壁面逐漸形成一對(duì)渦結(jié)構(gòu),且在t=1 s到t=3 s間可以觀測(cè)到Rayleigh- Taylor不穩(wěn)定性的發(fā)展；隨后,在t=3 s到t=5 s間,顆粒被渦結(jié)構(gòu)帶向計(jì)算域底部,計(jì)算域下方持續(xù)能夠觀察到渦結(jié)構(gòu)的變化；繼而,一部分粒子被帶往計(jì)算域上方,但最終,計(jì)算域內(nèi)的渦結(jié)構(gòu)趨于消散,整個(gè)稠密顆粒兩相流系統(tǒng)近乎穩(wěn)定下來(lái)。

通過定性分析可知,504個(gè)顆粒的沉降行為與文獻(xiàn)[22]中的結(jié)果具有較高的一致性。通過壁面阻滯效應(yīng)及渦結(jié)構(gòu)的形成與演化的成功模擬、稠密顆粒系統(tǒng)主要運(yùn)動(dòng)特征的正確捕捉,反映出本文所提出的CFD- DEM- IBM方法處理顆粒間碰撞、計(jì)算流固兩相間耦合作用力、處理顆粒與壁面碰撞行為的準(zhǔn)確性與可靠性。總的來(lái)說(shuō),對(duì)該問題的精確再現(xiàn)也進(jìn)一步證明了本文所提出的方法處理稠密顆粒兩相流問題的能力以及對(duì)流場(chǎng)的高精度解析。

3 結(jié) 論

(1) 高解析度CFD- DEM- IBM流固耦合數(shù)值模擬方法能夠方便地處理顆粒的移動(dòng)邊界并精確計(jì)算流固之間的相互作用力,獲得詳細(xì)的流場(chǎng)信息,深入研究顆粒行為。

(2) 在稀疏顆粒兩相流問題中,本文方法及非解析的CFD- DEM方法均可用于模擬顆粒行為；但對(duì)于稠密顆粒兩相流系統(tǒng),非解析的CFD- DEM方法不能準(zhǔn)確計(jì)算流固之間的相互作用力,無(wú)法反映顆粒尾流作用的影響,不適用于此類問題。

(3) 本文提出的CFD- DEM- IBM方法能夠?qū)崿F(xiàn)對(duì)流場(chǎng)的高精度解析,準(zhǔn)確、可靠地模擬顆粒沉降行為,對(duì)解決稠密顆粒兩相流問題尤其具有優(yōu)越性。

(4) 為獲得高解析度流場(chǎng),本文方法對(duì)網(wǎng)格要求較高,且流固之間的強(qiáng)耦合需要通過多次迭代實(shí)現(xiàn),導(dǎo)致本文方法計(jì)算成本較高;CFD- DEM- IBM方法的并行技術(shù)將是下一步的研究重點(diǎn)。