1 數值計算方法

濕顆粒碰撞過程的直接數值模擬通過流體體積模型（volume of fluid，VOF）模型追蹤氣液界麵，連續表麵張力模型（continuum surface force，CSF）模型計算界麵的表麵張力，兩者共同作用，控製液膜的形變過程。顆粒運動則采用六自由度模型（six degrees of freedom，6 DOF）結合重疊網格技術進行描述，為了模擬顆粒的反彈運動並提取分析數據，本文開發並應用了用戶自定義函數。

1.1 VOF模型

采用VOF模型來模擬氣液兩相流體的流動問題，通過求解體積分數的對流方程實現相界麵的瞬態追蹤。在VOF模型中，通過跟蹤不同相的體積分數來描述不同相的位置和界麵形狀[式(1)]。

∂α∂t+∇⋅(αu)=0(1)

式中，α為液體的體積分數，%，α=1處的網格充滿液體，α=0處的網格充滿氣體，0<α<1的網格區域含有自由麵；t為時間，s；u為流體的速度矢量，m/s；∇為梯度算子。

假設流體為不可壓縮流體，且氣液兩相的密度保持不變，涉及質量與動量守恒。

連續性方程為式(2)。

∇⋅u=0(2)

動量方程為式(3)。

∂(ρu)∂t+u⋅∇(ρu)=-∇p+∇⋅[μ⋅(∇u+∇uT)]+g+FSV(3)

式中，ρ為流體的局部平均密度，kg/m3；p為流體的壓力，Pa；μ為動力黏度，Pa·s；g為流體自身重力，N；FSV為表麵張力的等價體積力，N。

局部平均密度ρ和動力黏度μ表達式為式(4)和式(5)。

ρ=αρl+(1-α)ρg(4)

μ=αμl+(1-α)μg(5)

式中，下角標l為液體相；下角標g為氣體相。

1.2 連續表麵張力模型

連續表麵張力模型能夠將液體界麵上的力轉換為周圍的體積內連續分布的力，被廣泛應用於液滴流動模擬中。界麵處表麵張力可以表示為式(6)。

Fsf=σκ∇φ(6)

式中，Fsf為表麵張力，N；σ為表麵張力係數；∇φ為界麵切向的梯度算子；κ為表麵曲率。

界麵處的單位法向量計算如式(7)、式(8)。

n=-∇φ|∇φ|(7)

n=nwcosθ+twcosθ(8)

式中，θ為氣-液界麵在壁麵處的接觸角，(°)；nw和tw分別為與壁麵法向和切向的單位向量。得到界麵的曲率為式(9)。

κ=∇⋅∇φ|∇φ|(9)

1.3 顆粒運動模型

控製顆粒運動位移的方程為式(10)。

dupdt=1mpFp+Fv+Fsf+Fg(10)

式中，Fp為壓差阻力，N；Fv為黏性阻力，N；Fg為重力，N；mp為顆粒質量，kg；up為顆粒速度，m/s。

為進一步分析濕顆粒碰撞過程中的能量損失，可以將其分為以下四部分：壓差阻力引起的能量損失比Ep、黏性阻力引起的能量損失比Ev、由濕顆粒間非完全彈性碰撞引起的能量損失比Ec、由表麵張力引起的動能耗損比Esf。並將顆粒1剛接觸液膜時刻的初始動能定義為E0。計算見式(11)~式(13)。

Ep=(∫h0htFpdh+∫hthmFpdh)/E0Fp=-∫Spiz⋅ndS(11)

Ev=(∫h0htFvdh+∫hthmFvdh)/E0Fv=μ∫S∇uz⋅ndS(12)

Ec=(12mut2-12mus2)/E0(13)

式中，h0為顆粒1接觸液膜時質心高度，μm；ht為顆粒1與顆粒2碰撞時質心高度，μm；hm為顆粒1碰撞後反彈到最高點的質心高度，μm；ut為顆粒1碰撞前的速度，m/s；us為顆粒1碰撞後的速度，m/s。

1.4 重疊網格

本文使用的重疊網格方法廣泛應用於模擬移動網格的邊界運動問題。重疊網格技術的工作原理包括挖孔、重疊最小化和供體搜索。首先，生成互不約束的背景網格和重疊網格，滿足N-S求解器的三階離散化。在求解的迭代過程中，流場數據通過內嵌邊界進行傳輸。其次，確定執行各自職責的單元，包括挖孔單元、重疊邊界單元（接收單元）和計算單元（供體單元）。挖孔單元是指通過挖掘操作丟棄的一些煩瑣單元。重疊邊界單元主要圍繞挖孔單元，並被視為每個重疊網格中的邊界。重疊邊界單元與計算單元之間的關係通過插值確定。具體實現流程如圖1所示。

圖1   重疊網格計算求解流程