底瑛棠,趙蘭浩,毛 佳

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

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

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

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

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

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

1 數(shù)值方法

1.1 流體控制方程

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

·un+θ1=0

(1)

(2)

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

1.2 固體控制方程

固體的運動由牛頓第二運動方程控制:

(3)

(4)

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

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

(5)

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

(6)

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

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

(7)

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

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

1.3 浸入邊界法

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

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

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)可見,流體的速度、壓力與附加體力項相互耦合,需要采用迭代方案求解流場。根據(jù)分步投影法,定義不包含壓力項與附加體力項的中間速度場u*,k=un+Kn,流場的求解過程如圖2所示。

圖2 流體域迭代求解步驟

其中,內(nèi)循環(huán)中附加體力項根據(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所示,采用迭代的方式隱式求解附加體力項,在保持流體的速度、壓力與附加體力項相互耦合的同時,滿足無滑移邊界條件,并在附加體力與壓力項的迭代計算內(nèi)進行速度校正,確保結(jié)果滿足散度自由。

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

顆粒與流體間存在雙向相互作用,形成強流固耦合系統(tǒng),采用分區(qū)方案求解分屬不同控制方程的流固兩相。根據(jù)標(biāo)準(zhǔn)Galerkin加權(quá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àn)假設(shè)第n步計算已完成,依據(jù)以下流程展開第n+1步的第k次迭代:

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

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

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

2 數(shù)值算例

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

(13)

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

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

2.1 單顆粒沉降

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

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

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

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

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

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

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

圖4 雙顆粒沉降流速矢量

在兩同軸顆粒的沉降過程中,初始位于上方的拖尾顆粒將會超越下方引導(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ū)⒘鲌龈叨冉馕?正確再現(xiàn)了DKT現(xiàn)象,與文獻[21]中的結(jié)果具有很好的一致性。

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

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

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

2.3 顆粒群沉降

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

圖6 504個顆粒沉降速度矢量

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

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

3 結(jié) 論

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

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

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

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