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

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

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

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

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

1 SPH控制方程

1.1 SPH流體相控制方程

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

ρ▽·V

(1)

(2)

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

(3)

δab

(4)

(5)

(6)

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

▽iWij

(7)

(8)

(9)

(10)

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

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

(11)

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

1.2 SPH固體相模型

剛體上任意一點(diǎn)a的運(yùn)動(dòng)可以用式(12)進(jìn)行描述。

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

(12)

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

(13)

(14)

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

(15)

(16)

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

2 液固交界面算法

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

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

在標(biāo)準(zhǔn)的SPH邊界算法中，邊界粒子被當(dāng)成流體粒子直接參與相關(guān)物理量的計(jì)算。由于該方法在被截?cái)嗟闹С钟蛑兄苯佑?jì)算，導(dǎo)致相關(guān)物理量震蕩非常劇烈，不能得到高精度的計(jì)算結(jié)果。為了恢復(fù)核函數(shù)的歸一化性質(zhì)，提高邊界粒子相關(guān)物理量的計(jì)算精度，本文應(yīng)用如下公式計(jì)算邊界粒子上的相關(guān)物理量：

(17)

(18)

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

(19)

(20)

(21)

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

(22)

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

3 數(shù)值算例

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

圖2 顆粒沉降模型示意圖

3.1 單顆粒中心沉降

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

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

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

3.2 單顆粒偏心沉降

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

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

為了進(jìn)一步研究流體黏性阻尼力對顆粒的影響，這里稍微增大動(dòng)力黏性至0.001 Pa·s，對單顆粒偏心沉降過程進(jìn)行了新的數(shù)值模擬?？梢钥吹剑郝晕⒃龃箴ば院?，顆粒的整個(gè)運(yùn)動(dòng)形態(tài)發(fā)生了很大的變化，當(dāng)顆粒越過平衡位置后，由于黏性力及顆粒的兩體相互作用力，迫使顆粒趨近平衡位置，使得顆粒沿對稱軸搖擺著沉降。顆粒沉降路徑的差異必然導(dǎo)致流場中渦旋結(jié)構(gòu)的變化，黏性顆粒后的尾渦更加明顯，作用區(qū)域也更大，如圖5所示。

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

3.3 雙顆粒沉降

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

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

圖6 雙顆粒沉降過程

4 結(jié)束語

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