基于非解析計算流體力學和離散單元法的大顆粒在流場中的高效率運動模擬

熊書春， 臧孟炎

(華南理工大學機械與汽車工程學院， 廣州 510640)

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

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

1 數值方法

1.1 控制方程

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

(1)

(2)

式中:αf表示由“分割法”[14]計算得到的網格中流體所占體積分數(孔隙率)；uf為流體速度；τ=ν?uf為黏性應力張量；ρf和ν分別為流體密度和運動黏度；Rs,f為基于曳力在每個單元中計算出來的流相與顆粒相的動量交換項。

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

(3)

(4)

式中：mi、Ii、ui、wi分別為顆粒的質量、轉動慣量、速度和角速度；Fi,c表示顆粒間及顆粒與壁面之間的接觸力；Fi,f為周圍的流體作用在顆粒上的力，包括曳力、浮力等；其他外力(如重力、靜電力或磁力)總結為Fi,b、Ti,c為顆粒i所受轉矩。

1.2 黏結顆粒模型

顆?？梢曰贐PM粘結在一起，形成大顆粒或不規則形狀的顆粒，如圖1所示。

圖1 黏結形成的顆粒Fig.1 The bonded-particles

該模型假定兩個粒子通過彈簧-阻尼系統粘結在一起，表達式為

(5)

(6)

(7)

(8)

(9)

(10)

(11)

此模型采用虛擬的“黏結鍵”來黏結顆粒，這個黏結鍵可以承受切向和法向的微小位移，直到達到最大的法向和切向剪切應力。當法向和切向剪切應力超過設定的黏結強度時，黏結破裂，表達式為

(12)

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

1.3 曳力模型

通過合理的耦合模型可以實現連續流體與離散顆粒之間的動量交換。Gidaspow曳力模型[6]用于模擬作用在每個小球顆粒上的流固相互作用力，表達式為

(13)

(14)

式中：ds和Vs為顆粒的直徑和體積；μ為流體動力黏度。CD計算公式為

(15)

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

(16)

式(2)中的動量交換項Rs,f計算公式為

(17)

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

1.4 非解析CFD-DEM方法結合BPM

提出一種近似模擬大顆粒與流體耦合運動的新方法，將非解析CFD-DEM方法與BPM結合，把一個大固體化整為零，可以近似地模擬占據多個網格的大顆粒在流體中的運動。

如圖2所示，大顆粒直徑為15 mm，小顆粒直徑為1.5 mm，網格邊長為5 mm。當顆粒大于流體網格尺寸時，可以用多個直徑小于網格尺寸的小球形顆粒黏結形成大顆粒。因此，BPM可以滿足非解析CFD-DEM方法的要求。在此方法中，無重疊量的小球粘結而成的大顆粒內部存在間隙，密度相同的情況下顆粒的總體積和總質量相較于原始大顆粒會變小。根據顆粒運動方程，為保證顆粒動力學計算準確，黏結顆粒各小球所受流體作用力的總和須對應同比例減小。

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

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

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

(18)

(19)

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

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

該方法在開源框架CFDEM中實現，通過開源軟件OpenFOAM的CFD求解器計算流場，LIGGGHTS軟件基于DEM計算、更新顆粒的位置、速度、受力等信息，二者通過耦合接口進行質量、動量的傳遞和數據交換，實現并行雙向CFD-DEM耦合。為保證耦合的穩定性和準確性，DEM時間步長應始終小于CFD時間步長。本文中CFD時間步長都設置為DEM時間步長的10倍，這意味著CFD和DEM的數據交換在耦合接口中每10個DEM時間步長進行一次。

2 研究方法驗證

2.1 單個球形大顆粒沉降

分別使用結合BPM的非解析CFD-DEM方法和IBM方法模擬重力作用下單個球形大顆粒在黏性流體中的沉降運動，并與Cate等[15]的實驗結果進行比較，驗證數值方法的準確性。

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

表1 實驗和仿真參數Table 1 Experimental and simulation parameters

表2 黏結顆粒參數Table 2 Bonded-particle parameters

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

圖3 仿真結果與實驗對比Fig.3 Comparison of simulation results and experiments

2.2 計算效率驗證

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

兩種方法模擬所用CPU時間對比如圖4所示，顯然非解析CFD-DEM方法結合BPM的計算效率顯著高于IBM。因為IBM需要顆粒周圍非常精細的網格，且需要通過高精度的直接數值積分來計算曳力，盡管使用局部網格細化[3]減少了整體網格數量，但耦合計算過程仍耗時非常多。而本文提出的新方法，盡管BPM用多個小顆粒近似表達大顆粒，增加了離散元部分的計算量，但非解析CFD-DEM方法使用較粗的網格和標準化曳力模型確保了耦合計算效率。從圖4中還可以看到，大顆粒數量越多，該方法的高效性越明顯。

圖4 兩種方法計算耗時對比Fig.4 Comparison of calculation time between the two methods

3 非球形大顆粒的運動模擬

3.1 立方體大顆粒沉降

使用BPM表示立方體大顆粒，基于非解析CFD-DEM方法模擬其在黏性流體中的沉降，并且通過最大沉降速度仿真結果與表3所示的三組實驗結果[16]進行比較，驗證該方法對于模擬非球形大顆粒流場運動的有效性。

表3 實驗和仿真參數Table 3 Experimental and simulation parameters

計算域設置為100 mm×100 mm×400 mm，其高度足夠保證顆粒能達到穩態速度且不會碰到底部，網格數為40×40×160。顆粒的起始位置在計算域的中心線上，以避免旋轉運動。顆粒尺寸為8 mm×8 mm×8 mm，其體積相當于直徑為9.93 mm的球體。使用615個直徑為1 mm的小球黏結成該正方體大顆粒，黏結參數同上。

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

圖5 顆粒沉降速度曲線與實驗測量值的比較Fig.5 Comparison of the simulated particle settling velocities curves with the experimental values

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

3.2 橢球形大顆粒沉降

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

圖7顯示了四個時刻的顆粒運動位置和流場速度分布，其運動過程與Xia等[17]模擬的結果類似：開始時，橢球顆粒緩慢下落并逆時針擺動，逐漸靠近右壁面；然后，粒子擺動到長軸水平狀態，但是它有繼續轉動的趨勢。接著顆粒繼續下降并逆時針擺動，并逐漸接近左壁面，且隨著沉降過程擺動幅度逐漸減小至穩定的水平狀態。

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

4 結論

為高效率分析大顆粒在流場中的運動響應，提出將非解析CFD-DEM與BPM結合的方法，并通過大顆粒在流體中的沉降運動仿真分析進行計算精度和效率的評價，研究結論如下。

(1)仿真與實驗結果的對比分析表明，使用非解析CFD-DEM方法結合BPM可較準確地模擬球形大顆粒在流體中的沉降運動，且計算效率顯著高于IBM方法。

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