顆粒沉降問題的SPH數(shù)值模型研究

邵家儒，楊 瑜，黃志濤，鄒 麟

(1.重慶理工大學 機械工程學院， 重慶 400054； 2.慶鈴汽車(集團)有限公司， 重慶 400052)

在自然界和工業(yè)應用中存在很多顆粒兩相流動問題，如風沙泥沙輸運、污水處理、氣體除塵等。在液固兩相系統(tǒng)中，流體對固體的作用力來自壓強和剪應力兩方面，而該力通常受到流體的性質(zhì)(密度及黏性等)、相對流速、顆粒形狀等的影響[1-2]。顆粒兩相流動通常呈現(xiàn)多態(tài)性和復雜性，涉及多相流、微流體力學等多個方面，隨著研究的不斷深入，顆粒兩相流的內(nèi)涵和外延也在不斷拓展[3]。

在顆粒沉降的數(shù)值模擬研究方面，毛威等[4]采用格子Boltzmann方法模擬了熱對流條件下的顆粒沉降問題，分析了不同溫度下顆粒在流體中的沉降行為。Niu等[5]基于二維IB-LBM方法對二維圓形顆粒的沉降行為進行了數(shù)值模擬，發(fā)現(xiàn)沉降過程中領先的顆粒尾部會產(chǎn)生一個低壓渦旋，尾隨顆粒會被此尾渦捕捉，從而使其沉降速度快于領先的顆粒。Aidun等[6]使用二維格子Boltzmann方法模擬了圓柱體在充滿靜止流體的二維槽道中的沉降動力學行為，并在研究中發(fā)現(xiàn)了二維顆粒沉降的周期分岔和混沌狀態(tài)。劉漢濤等[7]應用任意拉格朗日-歐拉(ALE)方法模擬了熱對流條件下的顆粒沉降過程，指出熱對流會引起顆粒不同的動態(tài)尾跡。已有的數(shù)值模擬研究中大多沒有研究黏性、壁面與顆粒間的相互作用等因素對顆粒沉降問題的影響。

光滑粒子動力學(SPH)方法是一種拉格朗日型無網(wǎng)格粒子方法，在該方法中所模擬的介質(zhì)(流體/結構物)用一系列的粒子來代替，每個粒子具有自身的材料屬性，如質(zhì)量、密度等[8-9]。在進行物理量離散的過程中，每個粒子上的物理量僅由其支持域內(nèi)的粒子基于核函數(shù)插值決定。光滑粒子動力學方法非常適合研究涉及強非線性自由液面流動、流體與結構相互作用的問題，具有很好的自適應性[10-11]。本文基于光滑粒子動力學方法，引入雷諾平均湍流模型，構造了流體與結構的交界面算法，通過虛擬粒子技術解決了密度差異引起的數(shù)值震蕩問題。對單顆粒、雙顆粒沉降過程進行了數(shù)值模擬研究，確定了顆粒距離壁面位置、流體黏性等因素對沉降過程的影響，研究了不同顆粒間的相互作用。

1 SPH控制方程

1.1 SPH流體相控制方程

顆粒沉降問題中流體相的計算基于Navier-Stokes方程進行，描述如下[12]：

ρ▽·V

(1)

(2)

式中：ρ為密度；V為速度矢量；g為重力加速度；P為壓力；μ為動力黏性系數(shù)；R為雷諾應力張量,

(3)

δab

(4)

(5)

(6)

其中：Sab為平均應變率張量；k為湍動能；υt為渦黏系數(shù)；l為混合長度，l=Csh，cs為Smagorinsky常數(shù)，本文模擬中取值為0.12。經(jīng)過SPH離散后，控制方程中對應項可寫為如下形式：

▽iWij

(7)

(8)

(9)

(10)

其中：mj為粒子j的質(zhì)量；W為核函數(shù)；h為光滑長度；x為粒子的空間位置坐標。

任何一種不可壓縮流體理論上都可以看作是具有弱可壓縮性的，因此本文應用Monaghan提出的狀態(tài)方程[9]進行求解，描述如下：

(11)

其中：B=c2ρ0/γ，c為人工聲速，其取值通常為流場最大特征速度的10倍，這樣既可以使流體的可壓縮性保持在1%以下，又可以節(jié)約計算時間；ρ0為參考密度；γ=7。

1.2 SPH固體相模型

剛體上任意一點a的運動可以用式(12)進行描述。

ua=uc+uca=uc+ωc×rca

(12)

式中：uc為剛體質(zhì)心的平動速度；ωc為剛體自轉(zhuǎn)的角速度；rca為從剛體質(zhì)心到點a的位移向量。剛體質(zhì)心的平動和剛體自轉(zhuǎn)的角速度可由如下方程求解：

(13)

(14)

其中：F為整個剛體受到的總的流體作用力；M為剛體的質(zhì)量；Lc為剛體受到的外力對剛體質(zhì)心的力矩；Ic為剛體相對于其質(zhì)心的轉(zhuǎn)動慣量。對固體相方程進行SPH離散得：

(15)

(16)

式中：aj為粒子j的加速度；rj為粒子j的空間位置；R0為旋轉(zhuǎn)中心的位置。

2 液固交界面算法

如圖1所示，固壁邊界由3類固定的虛粒子構成，均勻分布在邊界外，其層數(shù)由所應用核函數(shù)的支持域半徑來確定，從而保證在真實粒子的支持域內(nèi)不會出現(xiàn)粒子缺失[12]。

圖1 流固耦合交界面模型示意圖

在標準的SPH邊界算法中，邊界粒子被當成流體粒子直接參與相關物理量的計算。由于該方法在被截斷的支持域中直接計算，導致相關物理量震蕩非常劇烈，不能得到高精度的計算結果。為了恢復核函數(shù)的歸一化性質(zhì)，提高邊界粒子相關物理量的計算精度，本文應用如下公式計算邊界粒子上的相關物理量：

(17)

(18)

β0+βx(xi-xj)+βy(yi-yj)]Wij

(19)

(20)

(21)

對于顆粒沉降問題而言，固體顆粒與流體間往往存在一定的密度差異。為了消除密度差異對計算帶來的影響，本文采用了虛擬粒子技術，即當支持域內(nèi)存在密度不同的粒子時，采用如下公式進行計算：

(22)

該計算格式可以保證核函數(shù)的歸一化性質(zhì)得到滿足。

3 數(shù)值算例

本文基于矩形槽道建立了顆粒沉降問題的數(shù)值模型，如圖2所示。矩形槽道寬為d，高為h，內(nèi)部液體初始時刻保持靜止。槽道內(nèi)固體顆粒半徑為R，初始保持靜止，因密度差異在自身重力的作用下在槽道內(nèi)開始沉降。內(nèi)部液體設定為水，其密度為ρf=1 000 kg/m3，顆粒的密度為ρs=1 250 kg/m3。

圖2 顆粒沉降模型示意圖

3.1 單顆粒中心沉降

本文首先對單顆粒的中心沉降問題進行了數(shù)值模擬。模型中動力黏性為0.01 Pa·s，顆粒半徑R=0.001 25 m，槽道寬度d=0.02 m，高度h=0.06 m。

從圖3中可以看到:在單顆粒中心位置沉降過程中，整個流場基本保持對稱。顆粒密度與流體密度的差異使得顆粒沿重力方向運動。運動過程中產(chǎn)生的拖曳力帶動顆粒附近流體向下運動，同時顆粒下降過程中對兩側流體的擠壓使其具有向上的速度矢量，從而形成渦旋結構。隨著時間的增加，顆粒的速度逐漸增加，后期在黏性阻尼力的作用下將達到一個穩(wěn)定的速度。這一過程中，顆粒附近的漩渦結構逐漸變化，由于邊界的影響，最后形成2個細長的漩渦結構。

圖3 槽道中心位置處顆粒沉降過程

3.2 單顆粒偏心沉降

如果顆粒不從中心位置沉降，而是中心偏離對稱軸一段距離，那么沉降過程中將會產(chǎn)生很多新的現(xiàn)象。這里對模型進行修改，模擬了顆粒的偏心沉降問題。模型中d=0.01 m，顆粒中心距離左側壁面0.007 5 m。首先忽略流體的黏性作用，如圖4所示。由于計算域的不對稱性，剛體與剛體壁面存在兩體相互作用，驅(qū)動顆粒向平衡位置運動，在顆粒運動過程中驅(qū)動周圍流體從而形成尾渦，尾渦隨顆粒的運動與時間的推移逐漸發(fā)展并最終脫落。由于不考慮黏性的作用，顆粒在流場中沒有黏性阻尼作用，流場中復雜的渦結構使得顆粒越過平衡位置，最終與剛體壁面接觸。

圖4 無黏性單顆粒偏心沉降過程

為了進一步研究流體黏性阻尼力對顆粒的影響，這里稍微增大動力黏性至0.001 Pa·s，對單顆粒偏心沉降過程進行了新的數(shù)值模擬。可以看到：略微增大黏性后，顆粒的整個運動形態(tài)發(fā)生了很大的變化，當顆粒越過平衡位置后，由于黏性力及顆粒的兩體相互作用力，迫使顆粒趨近平衡位置，使得顆粒沿對稱軸搖擺著沉降。顆粒沉降路徑的差異必然導致流場中渦旋結構的變化，黏性顆粒后的尾渦更加明顯，作用區(qū)域也更大，如圖5所示。

圖5 有黏性單顆粒偏心沉降過程

3.3 雙顆粒沉降

增加沉降顆粒的數(shù)目，流場中又會出現(xiàn)一些新的現(xiàn)象。本算例對雙顆粒沉降問題進行了數(shù)值模擬，分析了固體顆粒間的相互作用及引起的流場變化情況。模型中顆粒的半徑、密度均與本文之前算例保持一致，槽道寬度保持為0.01 m，水的動力黏性為0.001 Pa·s。該模型中兩顆粒初始時刻的中心分別距離左側壁面0.007 5、0.002 5 m，距離槽道底部分別為0.04 m和0.035 m。

圖6為雙顆粒沉降過程。由于初始時刻兩顆粒的中心位于偏離槽道中心線±0.025 m處，顆粒在初始沉降階段仍然會趨向于槽道中心線運動，此階段兩顆粒的相互作用并不明顯，然而隨著時間的增加，上方顆粒A開始進入下方顆粒B的尾渦區(qū)，出現(xiàn)了拖拽現(xiàn)象，此時兩顆粒的垂直速度都明顯增大，顆粒A的速度將大于顆粒B，兩顆粒間的相對距離逐漸減小，最終相互接觸。接觸后在扭轉(zhuǎn)力矩的作用下，顆粒的旋轉(zhuǎn)作用十分明顯，親和后的顆粒開始分開，同時顆粒A超過顆粒B繼續(xù)向下沉降，而顆粒B在扭轉(zhuǎn)力矩和邊界力的作用下產(chǎn)生了一個向上的速度。當顆粒B的動能由黏性力和漩渦作用消耗掉后又會向下繼續(xù)沉降。SPH方法可以很好地模擬出雙顆粒沉降問題中的DKT過程，即拖拽(drafing)、接觸(kissing)與翻轉(zhuǎn)(tumbling)。

圖6 雙顆粒沉降過程

4 結束語

本文基于光滑粒子動力學方法建立了矩形槽道內(nèi)的顆粒沉降數(shù)值模型。應用雷諾平均湍流模型以及涉及密度差異的交界面算法，模擬了流場中的渦旋結構以及流體與固體顆粒間的相互作用。通過SPH模型分析可以發(fā)現(xiàn)：在單顆粒沉降中，顆粒的初始位置、流體的黏性對沉降軌跡有較大的影響。在雙顆粒沉降中，SPH模型很好地模擬出了顆粒在黏性力、邊界力作用下出現(xiàn)的拖拽、接觸與翻轉(zhuǎn)過程。