熊書(shū)春， 臧孟炎

(華南理工大學(xué)機(jī)械與汽車(chē)工程學(xué)院， 廣州 510640)

計(jì)算流體力學(xué)(computational fluid dynamics, CFD)與離散單元法(discrete element method, DEM)耦合模擬流-固兩相流問(wèn)題，已廣泛應(yīng)用于各種工程研究。其所需要的經(jīng)驗(yàn)參數(shù)少，可方便地考慮顆粒尺寸分布，獲得顆粒尺度的微觀信息[1-2]。目前, 按照顆粒大小與流體網(wǎng)格尺寸的相對(duì)關(guān)系,可分為解析CFD-DEM方法和非解析CFD-DEM方法[3]。在解析CFD-DEM方法中，顆粒明顯大于流體單元，一個(gè)顆粒會(huì)覆蓋多個(gè)網(wǎng)格。常用的有浸沒(méi)邊界法(immersed boundary method, IBM)[4]等，可以精確地解析每個(gè)粒子周?chē)牧鲌?chǎng)，并通過(guò)直接數(shù)值積分計(jì)算流體作用在每個(gè)粒子上的力[5]。而在非解析CFD-DEM方法中，考慮的顆粒明顯小于流體單元，因此，一個(gè)網(wǎng)格可以包含多個(gè)顆粒。這類(lèi)方法不精細(xì)求解每個(gè)顆粒周?chē)牧鲌?chǎng)，而是基于曳力模型[6]在局部平均化的流場(chǎng)網(wǎng)格中將流體-顆粒相互作用模型化。

盡管IBM等解析CFD-DEM方法可用于高精度地模擬大顆粒在流場(chǎng)中的運(yùn)動(dòng)，但由于計(jì)算效率較低，研究大多局限于小數(shù)量的球形顆粒和二維情況下的規(guī)則非球形顆粒[7-8]。對(duì)于三維非球形顆粒，目前的研究也大多使用非解析CFD-DEM方法基于多球模型(multi-sphere model)近似模擬小于流體網(wǎng)格尺寸的較小顆粒的運(yùn)動(dòng)[9-11]。而多球模型的顆粒受力點(diǎn)為整體的質(zhì)心，且有顆粒重疊時(shí)會(huì)在計(jì)算曳力時(shí)造成困難。不同于將所有組合球都作為一整個(gè)剛體運(yùn)動(dòng)的多球模型，黏結(jié)顆粒模型(bonded-particle model, BPM)可用無(wú)重疊的球形顆粒模擬各種形狀[12-13]。為實(shí)現(xiàn)高效的大顆粒耦合數(shù)值仿真，在開(kāi)源框架CFDEM中，實(shí)現(xiàn)了一種近似模擬大顆粒與流體耦合作用的新方法，將非解析CFD-DEM方法與BPM結(jié)合，模擬大顆粒在流體中的運(yùn)動(dòng)。

1 數(shù)值方法

1.1 控制方程

非解析CFD-DEM方法的控制方程是流相的連續(xù)方法和固相的離散方法的組合。不可壓縮流體的運(yùn)動(dòng)由歐拉框架下的連續(xù)方法描述，基于CFD，求解連續(xù)性方程和局部平均的N-S(Navier-Stokes)方程得到每個(gè)網(wǎng)格尺度上的流場(chǎng)參數(shù)近似解，即

(1)

(2)

式中:αf表示由“分割法”[14]計(jì)算得到的網(wǎng)格中流體所占體積分?jǐn)?shù)(孔隙率)；uf為流體速度；τ=ν?uf為黏性應(yīng)力張量；ρf和ν分別為流體密度和運(yùn)動(dòng)黏度；Rs,f為基于曳力在每個(gè)單元中計(jì)算出來(lái)的流相與顆粒相的動(dòng)量交換項(xiàng)。

固體顆粒由拉格朗日方法追蹤，基于DEM對(duì)單個(gè)顆粒的運(yùn)動(dòng)分析通過(guò)求解牛頓運(yùn)動(dòng)方程得到，顆粒間的作用采用合理的接觸模型描述。單個(gè)顆粒在基于DEM的拉格朗日框架中平移和旋轉(zhuǎn)運(yùn)動(dòng)的控制方程表示為

(3)

(4)

式中：mi、Ii、ui、wi分別為顆粒的質(zhì)量、轉(zhuǎn)動(dòng)慣量、速度和角速度；Fi,c表示顆粒間及顆粒與壁面之間的接觸力；Fi,f為周?chē)牧黧w作用在顆粒上的力，包括曳力、浮力等；其他外力(如重力、靜電力或磁力)總結(jié)為Fi,b、Ti,c為顆粒i所受轉(zhuǎn)矩。

1.2 黏結(jié)顆粒模型

顆粒可以基于BPM粘結(jié)在一起，形成大顆?；虿灰?guī)則形狀的顆粒，如圖1所示。

圖1 黏結(jié)形成的顆粒Fig.1 The bonded-particles

該模型假定兩個(gè)粒子通過(guò)彈簧-阻尼系統(tǒng)粘結(jié)在一起，表達(dá)式為

(5)

(6)

(7)

(8)

(9)

(10)

(11)

此模型采用虛擬的“黏結(jié)鍵”來(lái)黏結(jié)顆粒，這個(gè)黏結(jié)鍵可以承受切向和法向的微小位移，直到達(dá)到最大的法向和切向剪切應(yīng)力。當(dāng)法向和切向剪切應(yīng)力超過(guò)設(shè)定的黏結(jié)強(qiáng)度時(shí)，黏結(jié)破裂，表達(dá)式為

(12)

式(2)中：r為黏結(jié)半徑。

1.3 曳力模型

通過(guò)合理的耦合模型可以實(shí)現(xiàn)連續(xù)流體與離散顆粒之間的動(dòng)量交換。Gidaspow曳力模型[6]用于模擬作用在每個(gè)小球顆粒上的流固相互作用力，表達(dá)式為

(13)

(14)

式中：ds和Vs為顆粒的直徑和體積；μ為流體動(dòng)力黏度。CD計(jì)算公式為

(15)

式(15)中：Res為顆粒雷諾數(shù)，定義為

(16)

式(2)中的動(dòng)量交換項(xiàng)Rs,f計(jì)算公式為

(17)

式(17)中：n為網(wǎng)格中的顆粒數(shù)量；ΔVcell為網(wǎng)格體積；Fi,D表示網(wǎng)格中各個(gè)顆粒所受曳力。

1.4 非解析CFD-DEM方法結(jié)合BPM

提出一種近似模擬大顆粒與流體耦合運(yùn)動(dòng)的新方法，將非解析CFD-DEM方法與BPM結(jié)合，把一個(gè)大固體化整為零，可以近似地模擬占據(jù)多個(gè)網(wǎng)格的大顆粒在流體中的運(yùn)動(dòng)。

如圖2所示，大顆粒直徑為15 mm，小顆粒直徑為1.5 mm，網(wǎng)格邊長(zhǎng)為5 mm。當(dāng)顆粒大于流體網(wǎng)格尺寸時(shí)，可以用多個(gè)直徑小于網(wǎng)格尺寸的小球形顆粒黏結(jié)形成大顆粒。因此，BPM可以滿(mǎn)足非解析CFD-DEM方法的要求。在此方法中，無(wú)重疊量的小球粘結(jié)而成的大顆粒內(nèi)部存在間隙，密度相同的情況下顆粒的總體積和總質(zhì)量相較于原始大顆粒會(huì)變小。根據(jù)顆粒運(yùn)動(dòng)方程，為保證顆粒動(dòng)力學(xué)計(jì)算準(zhǔn)確，黏結(jié)顆粒各小球所受流體作用力的總和須對(duì)應(yīng)同比例減小。

圖2 流體網(wǎng)格中的大顆粒Fig.2 A large particle represented in fluid meshes

基于每個(gè)小球顆粒計(jì)算流體的作用力并將其分別施加在每個(gè)組成球上，不需要計(jì)算整體質(zhì)心和轉(zhuǎn)動(dòng)慣量，這與將力施加到組合顆粒的質(zhì)心的多球模型方法相反。由于每個(gè)黏結(jié)組成球所處流場(chǎng)的位置不同，流體作用在每個(gè)小粒子上的力可能不同，但是由于黏結(jié)鍵的作用，所有粒子宏觀上都可以保持幾乎一致的運(yùn)動(dòng)狀態(tài)?？傮w而言，黏結(jié)力和力矩是內(nèi)力，所有粒子運(yùn)動(dòng)參數(shù)的平均值可用于描述整個(gè)黏結(jié)顆粒的運(yùn)動(dòng)狀態(tài)。在非解析CFD-DEM方法中，流體對(duì)球形顆粒無(wú)力矩作用，只對(duì)黏結(jié)的小球顆粒施加力在其質(zhì)心上。但是，各個(gè)小球所受流體作用力相對(duì)于整個(gè)黏結(jié)顆粒的質(zhì)心會(huì)產(chǎn)生力矩作用，且能通過(guò)黏結(jié)鍵傳遞，因此該方法可近似地表示流體對(duì)非球形大顆粒的力矩作用，模擬在流體中旋轉(zhuǎn)的非球形顆粒。

由于非解析CFD-DEM方法和解析CFD-DEM方法計(jì)算流體對(duì)顆粒作用力的原理不同，且黏結(jié)顆粒內(nèi)部存在間隙，因此必須修正曳力模型，使黏結(jié)顆粒所受流體曳力總和與大顆粒所受曳力之比等于質(zhì)量之比，以保證黏結(jié)顆粒整體運(yùn)動(dòng)的準(zhǔn)確性。引入校正因子K，即

(18)

(19)

式(19)中：A為大顆粒的迎風(fēng)面積。

盡管使用更多的組成球可以更準(zhǔn)確地表示大顆粒的形狀，但是計(jì)算效率會(huì)大幅降低。文中將小球的直徑都設(shè)置為大顆粒(等效)直徑的1/10左右，這可以同時(shí)滿(mǎn)足計(jì)算方法對(duì)網(wǎng)格尺寸、計(jì)算準(zhǔn)確性和效率的要求。根據(jù)顆粒的楊氏模量和泊松比來(lái)設(shè)置黏結(jié)鍵的法向和切向剛度，以保證鍵的精確變形。法向和切向黏結(jié)強(qiáng)度設(shè)置為足夠大，可以防止黏結(jié)斷裂。

該方法在開(kāi)源框架CFDEM中實(shí)現(xiàn)，通過(guò)開(kāi)源軟件OpenFOAM的CFD求解器計(jì)算流場(chǎng)，LIGGGHTS軟件基于DEM計(jì)算、更新顆粒的位置、速度、受力等信息，二者通過(guò)耦合接口進(jìn)行質(zhì)量、動(dòng)量的傳遞和數(shù)據(jù)交換，實(shí)現(xiàn)并行雙向CFD-DEM耦合。為保證耦合的穩(wěn)定性和準(zhǔn)確性，DEM時(shí)間步長(zhǎng)應(yīng)始終小于CFD時(shí)間步長(zhǎng)。本文中CFD時(shí)間步長(zhǎng)都設(shè)置為DEM時(shí)間步長(zhǎng)的10倍，這意味著CFD和DEM的數(shù)據(jù)交換在耦合接口中每10個(gè)DEM時(shí)間步長(zhǎng)進(jìn)行一次。

2 研究方法驗(yàn)證

2.1 單個(gè)球形大顆粒沉降

分別使用結(jié)合BPM的非解析CFD-DEM方法和IBM方法模擬重力作用下單個(gè)球形大顆粒在黏性流體中的沉降運(yùn)動(dòng)，并與Cate等[15]的實(shí)驗(yàn)結(jié)果進(jìn)行比較，驗(yàn)證數(shù)值方法的準(zhǔn)確性。

如表1所示，設(shè)計(jì)了四組與Cate試驗(yàn)對(duì)應(yīng)的算例，雷諾數(shù)由試驗(yàn)得到的顆粒穩(wěn)態(tài)速度算出。大顆粒的直徑為15 mm，密度為1 120 kg/m3，球體起始位置距底部高120 mm，重力加速度為9.81 m/s2。計(jì)算域?yàn)?00 mm×100 mm×160 mm，兩種方法的網(wǎng)格數(shù)都為25×25×40(IBM方法使用局部網(wǎng)格細(xì)化[3])，CFD和DEM的計(jì)算時(shí)間步長(zhǎng)分別為1×10-5s和1×10-6s，BPM的黏結(jié)參數(shù)如表2所示。

表1 實(shí)驗(yàn)和仿真參數(shù)Table 1 Experimental and simulation parameters

表2 黏結(jié)顆粒參數(shù)Table 2 Bonded-particle parameters

圖3顯示了在不同Re的條件下顆粒沉降速度、顆粒距底部高度與顆粒直徑比值的時(shí)間歷程。其中虛線表示非解析CFD-DEM方法結(jié)合BPM獲得的數(shù)值計(jì)算結(jié)果(顆粒沉降速度與顆粒高度值都是615個(gè)顆粒的平均值)，實(shí)線表示IBM的模擬結(jié)果，而散點(diǎn)符號(hào)表示實(shí)驗(yàn)結(jié)果。對(duì)比可知本文提出的仿真方法所得分析結(jié)果與IBM和實(shí)驗(yàn)結(jié)果一致性都很好，驗(yàn)證了該方法的有效性。

圖3 仿真結(jié)果與實(shí)驗(yàn)對(duì)比Fig.3 Comparison of simulation results and experiments

2.2 計(jì)算效率驗(yàn)證

為了驗(yàn)證新方法的效率及其與大顆粒數(shù)量之間的關(guān)系，分別使用該方法和IBM模擬了重力作用下一個(gè)、四個(gè)和九個(gè)大球形顆粒在黏性流體中沉降的運(yùn)動(dòng)。流體密度為1 000 kg/m3，動(dòng)力黏度為0.1 N·s/m2；大顆粒的直徑為15 mm，密度為2 000 kg/m3，現(xiàn)象時(shí)間為0.5 s。

兩種方法模擬所用CPU時(shí)間對(duì)比如圖4所示，顯然非解析CFD-DEM方法結(jié)合BPM的計(jì)算效率顯著高于IBM。因?yàn)镮BM需要顆粒周?chē)浅＞?xì)的網(wǎng)格，且需要通過(guò)高精度的直接數(shù)值積分來(lái)計(jì)算曳力，盡管使用局部網(wǎng)格細(xì)化[3]減少了整體網(wǎng)格數(shù)量，但耦合計(jì)算過(guò)程仍耗時(shí)非常多。而本文提出的新方法，盡管BPM用多個(gè)小顆粒近似表達(dá)大顆粒，增加了離散元部分的計(jì)算量，但非解析CFD-DEM方法使用較粗的網(wǎng)格和標(biāo)準(zhǔn)化曳力模型確保了耦合計(jì)算效率。從圖4中還可以看到，大顆粒數(shù)量越多，該方法的高效性越明顯。

圖4 兩種方法計(jì)算耗時(shí)對(duì)比Fig.4 Comparison of calculation time between the two methods

3 非球形大顆粒的運(yùn)動(dòng)模擬

3.1 立方體大顆粒沉降

使用BPM表示立方體大顆粒，基于非解析CFD-DEM方法模擬其在黏性流體中的沉降，并且通過(guò)最大沉降速度仿真結(jié)果與表3所示的三組實(shí)驗(yàn)結(jié)果[16]進(jìn)行比較，驗(yàn)證該方法對(duì)于模擬非球形大顆粒流場(chǎng)運(yùn)動(dòng)的有效性。

表3 實(shí)驗(yàn)和仿真參數(shù)Table 3 Experimental and simulation parameters

計(jì)算域設(shè)置為100 mm×100 mm×400 mm，其高度足夠保證顆粒能達(dá)到穩(wěn)態(tài)速度且不會(huì)碰到底部，網(wǎng)格數(shù)為40×40×160。顆粒的起始位置在計(jì)算域的中心線上，以避免旋轉(zhuǎn)運(yùn)動(dòng)。顆粒尺寸為8 mm×8 mm×8 mm，其體積相當(dāng)于直徑為9.93 mm的球體。使用615個(gè)直徑為1 mm的小球黏結(jié)成該正方體大顆粒，黏結(jié)參數(shù)同上。

圖5顯示了顆粒沉降速度曲線，可見(jiàn)仿真結(jié)果的最大值與實(shí)驗(yàn)結(jié)果吻合良好。與球體的沉降運(yùn)動(dòng)類(lèi)似，立方體顆粒的速度在沉降過(guò)程中逐漸增加，直到達(dá)到穩(wěn)態(tài)為止。圖6為算例C2中密度為4 450 kg/m3的顆粒達(dá)到穩(wěn)態(tài)時(shí)流場(chǎng)速度等高線圖。結(jié)果表明本文提出的新方法可成功模擬IBM等解析CFD-DEM方法難以實(shí)現(xiàn)的三維情況下非球形大顆粒在流體中的耦合運(yùn)動(dòng)。

圖5 顆粒沉降速度曲線與實(shí)驗(yàn)測(cè)量值的比較Fig.5 Comparison of the simulated particle settling velocities curves with the experimental values

圖6 算例C2的流場(chǎng)速度等高線圖Fig.6 Contour map of flow field velocity of case C2

3.2 橢球形大顆粒沉降

為了驗(yàn)證該方法適用于可旋轉(zhuǎn)的非球形大顆粒，進(jìn)一步使用非解析CFD-DEM方法結(jié)合BPM模擬了三維情況下單個(gè)橢球形大顆粒在充滿(mǎn)黏性流體的細(xì)長(zhǎng)豎直管中的沉降運(yùn)動(dòng)。橢球的長(zhǎng)軸長(zhǎng)23.8 mm，兩短軸長(zhǎng)11.9 mm，其體積約等于直徑為15 mm的球體。橢球顆粒通過(guò)黏結(jié)610個(gè)直徑為1.5 mm的小球而形成，黏結(jié)參數(shù)同上。顆粒和流體密度分別為1 100 kg/m3和1 000 kg/m3；流體動(dòng)力黏度為0.1 N·s/m2；計(jì)算域?yàn)?5.2 mm×95.2 mm×1 800 mm，網(wǎng)格數(shù)量為18×18×360。橢球顆粒最初重心位于計(jì)算域的中軸線，長(zhǎng)軸沿水平方向傾斜45°。

圖7顯示了四個(gè)時(shí)刻的顆粒運(yùn)動(dòng)位置和流場(chǎng)速度分布，其運(yùn)動(dòng)過(guò)程與Xia等[17]模擬的結(jié)果類(lèi)似：開(kāi)始時(shí)，橢球顆粒緩慢下落并逆時(shí)針擺動(dòng)，逐漸靠近右壁面；然后，粒子擺動(dòng)到長(zhǎng)軸水平狀態(tài)，但是它有繼續(xù)轉(zhuǎn)動(dòng)的趨勢(shì)。接著顆粒繼續(xù)下降并逆時(shí)針擺動(dòng)，并逐漸接近左壁面，且隨著沉降過(guò)程擺動(dòng)幅度逐漸減小至穩(wěn)定的水平狀態(tài)。

圖7 不同時(shí)刻橢球粒子的位置和流場(chǎng)速度云圖Fig.7 The position of the ellipsoidal particle and the flow field velocity contours at different time

4 結(jié)論

為高效率分析大顆粒在流場(chǎng)中的運(yùn)動(dòng)響應(yīng)，提出將非解析CFD-DEM與BPM結(jié)合的方法，并通過(guò)大顆粒在流體中的沉降運(yùn)動(dòng)仿真分析進(jìn)行計(jì)算精度和效率的評(píng)價(jià)，研究結(jié)論如下。

(1)仿真與實(shí)驗(yàn)結(jié)果的對(duì)比分析表明，使用非解析CFD-DEM方法結(jié)合BPM可較準(zhǔn)確地模擬球形大顆粒在流體中的沉降運(yùn)動(dòng)，且計(jì)算效率顯著高于IBM方法。

(2)非解析CFD-DEM與BPM結(jié)合的方法可有效模擬三維情況下非球形大顆粒在流體中包括平動(dòng)和轉(zhuǎn)動(dòng)的沉降運(yùn)動(dòng)，使非球形大顆粒在流場(chǎng)中的高效率運(yùn)動(dòng)響應(yīng)分析成為可能。