底瑛棠,趙蘭浩,毛 佳

(河海大學(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),通過(guò)牛頓第二定律描述顆粒運(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通過(guò)顆粒動(dòng)力學(xué)理論計(jì)算顆粒間的相互作用,模化碰撞過(guò)程,并無(wú)清晰的物理碰撞面。

近年來(lái),計(jì)算流體力學(xué)與離散單元法(Computational Fluid Dynamics- Discrete Element Method,CFD- DEM)因其處理顆粒碰撞問(wèn)題的優(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方法能夠考慮兩相之間的耦合作用,引入顆粒接觸模型描述顆粒的物理碰撞過(guò)程,并通過(guò)空隙率考慮顆粒的體積效應(yīng),是求解顆粒兩相流宏觀特性的極佳選擇。然而,傳統(tǒng)的CFD- DEM方法通過(guò)經(jīng)驗(yàn)公式計(jì)算兩相間的耦合作用力[11],稱作非解析的CFD- DEM方法。由于缺乏對(duì)顆粒邊界的描述,顆粒周?chē)牧鲌?chǎng)信息無(wú)法得到精確解析,從而無(wú)法描述顆粒的尾流效應(yīng)[12],決定了該方法只適用于顆粒間無(wú)干擾的稀疏顆粒兩相流問(wèn)題。

為解決目前非解析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)過(guò)程中的大位移,采用貼體動(dòng)網(wǎng)格方法處理顆粒的運(yùn)動(dòng)邊界計(jì)算成本較高,并可能產(chǎn)生網(wǎng)格扭曲等問(wèn)題。Peskin[14]提出的浸入邊界法(IBM)采用固定的結(jié)構(gòu)化網(wǎng)格,在固體邊界上覆蓋拉格朗日點(diǎn)以追蹤其運(yùn)動(dòng)邊界,回避了動(dòng)網(wǎng)格方法的缺陷。通過(guò)在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)換問(wèn)題,正逐步被應(yīng)用到更廣泛的工程領(lǐng)域[15],對(duì)于稠密顆粒系統(tǒng)的沉降這種動(dòng)態(tài)邊界問(wèn)題尤為適用。

本文基于CFD- DEM和IBM方法提出一種高解析度CFD- DEM- IBM流固耦合數(shù)值分析方法,基于自主開(kāi)發(fā)的程序開(kāi)展顆粒沉降問(wèn)題的數(shù)值模擬研究,通過(guò)與前人工作的對(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 顆粒沉降問(wèn)題的CFD- DEM- IBM方法示意

1.3 浸入邊界法

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

將動(dòng)量方程式(2)沿特征線展開(kāi)可得到時(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)可見(jiàn),流體的速度、壓力與附加體力項(xiàng)相互耦合,需要采用迭代方案求解流場(chǎng)。根據(jù)分步投影法,定義不包含壓力項(xiàng)與附加體力項(xiàng)的中間速度場(chǎng)u*,k=un+Kn,流場(chǎng)的求解過(guò)程如圖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ù)以下流程展開(kāi)第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é)果展開(kāi)對(duì)比。為凸顯高解析度CFD- DEM- IBM方法的優(yōu)勢(shì),非解析CFD- DEM方法的計(jì)算結(jié)果也在文中給出。其中,非解析方法中流固相互作用力依據(jù)式(13)計(jì)算:

(13)

式中:Cd為阻力系數(shù)；ρf為流體密度；di為顆粒pi的直徑；uf,i為顆粒pi周?chē)W(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耦合模型用于分析顆粒沉降問(wèn)題的正確性與可靠性,首先選擇經(jīng)典的單顆粒沉降問(wèn)題[19- 20]展開(kāi)計(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è)置為開(kāi)邊界外,其余三邊均設(shè)置為無(wú)滑固壁邊界。

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

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

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

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

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

為初步驗(yàn)證本文提出的方法對(duì)于稠密顆粒兩相流問(wèn)題的適用性、高解析度描述流場(chǎng)與顆粒行為的能力,進(jìn)行2個(gè)相同圓形顆粒自由下落問(wèn)題研究。計(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)處,由靜止受到重力作用開(kāi)始下落。

圖4 雙顆粒沉降流速矢量

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

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

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

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

2.3 顆粒群沉降

為進(jìn)一步驗(yàn)證本文所提出的CFD- DEM- IBM方法處理稠密顆粒兩相流問(wèn)題的能力,本節(jié)展開(kāi)了504個(gè)圓形顆粒在密閉容器內(nèi)的沉降過(guò)程[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è)顆粒沉降過(guò)程中的速度矢量圖?？梢?jiàn),沉降開(kāi)始時(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)。

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

3 結(jié) 論

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

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

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

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