基于兩相流理論的稀疏和致密顆粒流統一模型

何 康，施華斌，余錫平

(1.清華大學水利水電工程系，北京 100084；2.澳門大學土木與環境工程系，澳門 999078；3.港澳海洋研究中心，香港 999077；4.南方科技大學海洋科學與工程系，深圳 518055)

顆粒體和流體作為基本的兩相介質，廣泛存在于河流輸沙、海底滑坡等流動問題中[1?4]。在顆粒體和流體混合物的運動過程中，顆粒的體積分數和運動狀態均可能發生變化。如在海底滑坡過程中，滑坡體上方由于發生剪切摻混會形成低濃度的濁流區，其底部則是高濃度的致密顆粒流動區域，兩者均可能對深水工程和海底管道等設施造成重大破壞[5?8]。因此，研究多流動情形下顆粒體的運動過程具有重要意義。

數十年來，很多學者采用實驗和數值手段對稀疏顆粒體的流動過程進行了研究。早在20 世紀，研究者們就通過室內實驗對異重流等稀疏顆粒流的結構，及顆粒粒徑對流動的影響進行了分析[9?10]，也有學者建立簡化模型對顆粒體的運動過程進行預測[11]。隨著數值模擬技術的發展，近些年有學者采用直接數值模擬[12]，大渦模擬[13]和雷諾平均方程模擬[14]對稀疏顆粒流問題進行了研究。其中對顆粒體和流體的瞬時質量、動量守恒方程進行平均后可獲得兩相平均方程，進一步對各應力項、相間作用項進行封閉可建立基于連續介質假設的兩相流模型[15]。該方法在計算效率上更具優勢且保證了一定的計算精度，更適用于解決實際問題，但以往的模型大多僅考慮了稀疏顆粒體間的碰撞作用，而忽略了致密堆積狀態下顆粒體間的相互作用，因而還很難準確刻畫近底高濃度顆粒體的流動過程[16]。

近些年來，致密顆粒流問題也吸引了很多研究者的興趣。與稀疏顆粒流不同的是，致密顆粒體中顆粒體的體積分數達到自然堆積狀態，此時顆粒與顆粒間的持續接觸作用不可忽略[17]。此外，致密顆粒體和流體的相互作用，如拖曳力、孔隙水壓等均對顆粒體的運動過程有著重要影響，需要進行細致刻畫[18?19]。通過考慮顆粒體的摩擦和碰撞作用，耦合顆粒與流體的相間作用建立的兩相流模型為研究顆粒體的流動過程提供了新的視角[20]。然而，目前還未探究該兩相流模型在統一描述稀疏和致密顆粒流問題的可能性。此外，分析不同粒徑顆粒體和流體混合物的運動過程也有助于理解實際的顆粒流問題。因此，本文采用統一的兩相流模型同時對不同粒徑的稀疏和致密顆粒體流動問題開展了研究，并對流體紊動的影響，顆粒體的運動過程及其與流體的相互作用進行了分析，以加深對多流動情形下顆粒體和流體混合物流動問題的理解。

1 數學模型

1.1 基本方程

本文采用文獻[20]中的兩相流模型，平均后的顆粒體和流體控制方程為：

式中：下標 s 和 f分別為顆粒和流體相；本文僅考慮立面二維流動問題，下標i和j分別為水平和垂直方向坐標，并遵循求和約定；α為體積分數；ρ為相密度；U為速度；p為壓強；τ為剪切應力；I為相間作用力；g為重力加速度。

1.2 顆粒相本構關系

將顆粒相應力劃分為摩擦與碰撞應力，可表示為：

式中，上標 f 和 c分別為顆粒應力的摩擦與碰撞部分。

摩擦剪切應力可寫為：

摩擦區域內的顆粒相壓強采用以下公式進行計算[23]：

顆粒相碰撞應力基于顆粒動理論模化，可表示為[25]：

式中：Θs為顆粒溫度，可由顆粒溫度控制方程求解，詳見文獻 [20]；η=(1+e)/2，e為顆粒的回彈系數；R=(2?αs)/[2(1?αs)3]為顆粒的徑向分布函數[26]；和為碰撞粘性系數和碰撞體積粘性系數，有：

圖1 給出了不同區域內顆粒相壓強隨顆粒體積分數的變化。從圖1 可以看到顆粒體積分數較小時，稀疏顆粒體僅承擔碰撞壓強；當顆粒體積分數增大到可堆積狀態時，致密顆粒體將同時承擔碰撞和摩擦壓強，且顆粒體總壓強將逐漸由碰撞壓強主導轉變為由摩擦壓強主導。因此，通過綜合考慮顆粒體的摩擦與碰撞作用，本文實現了對不同濃度顆粒流的統一描述。

圖1 顆粒相壓強隨顆粒體積分數的變化Fig.1 Variation of particle pressures with the particle volume fraction

1.3 流體相本構關系

考慮紊流的影響，流體相剪切應力可模化為：

基于k?ε紊流模型，有：

式中：Cμ為常數，一般取為0.09；流體相紊動能kf和紊動能耗散率 εf通過求解考慮拖曳力修正后的k?ε紊流方程得到，詳細介紹可參考文獻[20]。

1.4 相間作用力

這里主要考慮顆粒與流體的相間拖曳力，表達式為：

式中：Kfs為流體與顆粒的相間拖曳力系數；右邊第二項為顆粒由于流體紊動導致的漂移速度項，其中為流體紊動運動粘性系數；為施密特數，根據已有的研究結果取值1.0[27]。

這里采用經典的Gidaspow(1994)拖曳力公式[28]，Kfs可表示為：

式中，Cd經由以下公式確定：

式中，Res=αfρf|Uf?Us|ds/μf為顆粒雷諾數。

1.5 計算方法及邊界條件

該兩相流模型基于開源平臺MFIX(multiphase flow with interphase exchanges)離散并求解。模型的控制方程通過矩形結構化網格，以有限體積法離散，并采用改進后的SIMPLE 算法求解壓力相關方程[29]；對流項的離散采用了二階Superbee 格式，并采用了延遲校正方法處理[30]；模型的計算在滿足Courant-Friedrichs-Lewy(CFL)收斂條件下，采用自適應時間步長的方法在時空上推進。

本文以稀疏和致密顆粒流為研究對象進行數值模擬。其中稀疏顆粒流算例以文獻[10]中的實驗進行設置，如圖2(a)所示；致密顆粒流算例采用了文獻[31]中的實驗數據，其設置如圖2(b)所示。實驗中，初始時刻稀疏和致密顆粒體均放置在水槽一側，隨著閘門開啟，顆粒體在重力和密度差作用下向前發展，并最終形成穩定的沉積體和堆積體。模型應用時將所有固壁設置為無滑移邊界，對頂部采用剛蓋假定，忽略流體自由表面的波動影響。

圖2 稀疏和致密顆粒流數值模擬設置 /mFig.2 Numerical setup for dilute and dense granular flows

2 稀疏顆粒流研究

2.1 模型驗證

本文選取文獻[10]中不同粒徑的稀疏顆粒流實驗進行驗證。其中細、粗顆粒粒徑分別為ds=0.025 mm和ds=0.069 mm ，密 度 均 為ρs=3217 kg/m3，初始顆粒體積分數為=3.49×10?3，實驗設置詳見圖2(a)。在模擬設置中，顆粒摩擦角?s=30?，忽略μ0s的影響，回彈系數e=0.8；實驗中流體為水，密度取為 ρf=1000 kg/m3，粘性系數μf=0.001 kg/(m·s)。經網格敏感性分析，模擬選取5 mm×2.5 mm的均勻網格進行計算。

模擬與實測的稀疏細、粗顆粒體的前端運動距離變化結果如圖3 所示，其中以 αs/=1等值線的前端位置作為稀疏顆粒體的前端運動距離。從圖3 可以看到，該模型可以很好地刻畫稀疏細顆粒體與粗顆粒體的運動過程，證明了該模型在稀疏顆粒流問題中的有效性。此外，模擬與實驗結果均表明初始運動階段稀疏細、粗顆粒體的前端運動距離基本一致，但在t≈20 s后，稀疏粗顆粒體的前端運動距離開始明顯落后于稀疏細顆粒體的，說明不同粒徑稀疏顆粒體的運動過程存在較大的差異，需要開展進一步的分析。

圖3 模擬與實驗的稀疏顆粒體前端運動距離比較Fig.3 Comparisons of simulated and measured front distances in dilute granular flows

2.2 相間作用的影響

圖4 給出了稀疏細、粗顆粒相對濃度 αs/的分布變化。由于顆粒和流體間存在密度和速度差異，可以發現，稀疏顆粒體在與流體的剪切運動過程中前端會呈現出近似橢圓形的頭部結構，后方則是懸浮顆粒體的主體結構和狹長的尾部結構，并在其中能發現數個渦旋結構，顆粒體和流體交界面整體呈現出一系列的開爾文－亥姆霍茲不穩定性結構。此外，也能觀察到稀疏細顆粒能在較長時間內保持向前運動，并維持明顯的頭部結構；而t=15 s和t=20 s時保持懸浮狀態的粗顆粒濃度則明顯更少，其頭部結構范圍也逐漸縮小，盡管此時細、粗顆粒體的前端運動距離差異較小。

圖4 細顆粒和粗顆粒相對濃度分布變化Fig.4 Evolutions of the relative concentration fieldfor fine particles and coarse particles

從相間作用來看，不同粒徑顆粒體與流體的運動差異主要由其受到的相間拖曳力決定。這里給出t=5 s時細顆粒和粗顆粒與流體的相對速度Us?Uf和相對拖曳力K(Uf?Us)/ρfg的分布，如圖5所示。從圖5 可以看到，Us?Uf的方向基本垂直向下，表明顆粒體相對流體主要做垂向沉積運動，且粗顆粒與流體的相對速度要明顯大于細顆粒與流體的相對速度，即粗顆粒的沉速更大。從相間拖曳力的分布來看，其方向基本垂直向上，即對顆粒體的沉降起阻礙作用。值得注意的是，壁面處顆粒受到的拖曳力明顯更大，這是壁面附近沉積顆粒濃度較大的緣故。此外還能發現，渦旋結構內粗顆粒受到的拖曳力小于細顆粒的，而壁面附近粗顆粒受到的拖曳力作用則明顯大于細顆粒的，即粗顆粒受到的沉積作用更強。

圖5 t =5 s 時稀疏細、粗顆粒流中流體與顆粒相間速度 U s ?Uf 和相對拖曳力 K(Uf ?Us)/ρfg 比較(黃色和白色虛線分別為=1 和 =0.1的等值線)Fig.5 Comparisons of interphase velocity Uas ?nUdf relative draKg(U ff ?oUrs)c/ρefg between the fluid and solid phase in dilute fine and coarse particle flows at t(=th5e s yellow and white dash lines are contour lini =th1 and =0.1)

圖6 展示了t=20 s時細、粗顆粒體中流體和顆粒速度的比較結果。從圖中可以看到，細顆粒中流體和顆粒體的速度更大，這是因為細顆粒的重力沉降作用偏弱，更容易在運動中保持懸浮狀態，因而懸浮顆粒體及周圍流體所能維持的運動強度更大。粗顆粒的沉降速度明顯大于細顆粒，垂向沉積作用更強，更容易丟失運動能量，因而顆粒體及其周圍流體速度要更小。故而在顆粒體的向前發展過程中能發現，t≈20 s后細顆粒體的前端運動距離明顯大于粗顆粒體的，且在運動過程中處于懸浮狀態的顆粒更多。

圖6 t =20 s 時稀疏細、粗顆粒流中流體速度 Uf 與顆粒速度 Us 比較(黃色和白色虛線分別為=1 和 =0.1的等值線)Fig.6 Comparisons of fluisd velocity Uf and solid velocity Us in dilute fine and coarse particle flows at t=20 s(the yellow and white dash lines are contour lines with and =0.1)

2.3 流體紊動的影響

在稀疏顆粒體運動過程中，流體的紊動摻混可能對顆粒體和流體的不穩定結構產生影響。圖7給出了稀疏細顆粒流中流體速度 |Uf|和相對紊動粘性系數的分布變化。從t=5 s到t=20 s，可以看到流體速度 |Uf|逐漸減小，其較大值主要分布在顆粒體前端和中間區域。值得注意的是，該階段流體的紊動粘性系數仍然呈現增大的趨勢，的值甚至可達100 以上，且大都分布在顆粒體的后上方，前端和下方的值則相對較小，并在底部形成了長條形的低紊動粘性系數區域。

圖7 稀疏細顆粒流中流體速度 | Uf|和 相對紊動粘性系數的分布變化(黃色和白色虛線分別為=1 和 =0.1的等值線)Fig.7 Evolutions of the fluid velocity |aUnfd| the relative turbulent viscin dilute fine particle flows (the yellow and white dash lines are contour lines with and =0.1)

為了進一步分析紊流對流體運動的影響，圖8給出了是否考慮紊流模型時細顆粒中流體渦量|ωf|=1/2|?Uf,j/?xi??Uf,i/?xj|的比較結果。從圖8可以明顯地看到，不考慮紊流模型時的流體渦量要明顯大于考慮紊流模型時的，即不考慮紊流模型時流體的渦旋運動更為強烈，這是因為此時流體的運動沒有受到紊流的影響，因而在運動過程中耗散的能量更少。

圖9 給出了不考慮紊流模型時細顆粒和粗顆粒的相對濃度分布。與圖3 中考慮紊流模型的計算結果相比，不考慮紊流模型后懸浮顆粒體的前端會呈現明顯的“抬起”現象，即顆粒體的頭部會被水流抬起，且前端運動距離要更小，與實驗結果不符。此外，不考慮紊流模型后懸浮顆粒體的后端結構發展得更為無序，可見紊流對顆粒體和流體交界面處的結構有著重要影響。

圖9 不考慮紊流模型時細顆粒和粗顆粒濃度分布演變Fig.9 Evolutions of the concentration field for fine particles and coarse particles without turbulence model

為了進一步分析流體紊動對流體和顆粒體的影響，圖10 給出了不考慮紊流t=20 s時細、粗顆粒體前端流體速度Uf和顆粒速度Us。與圖5 相比，不考慮紊流模型時渦旋結構處流體和顆粒的速度均要比考慮紊流時的更大，這是因為在運動過程中沒有受到流體紊動粘性的抑制作用。從運動方向上來看，不考慮紊流模型時流體和顆粒的運動有明顯的上下起伏變化，而考慮紊流后流體和顆粒的流動方向則基本一致向前。數值模擬結果表明，顆粒體與流體在剪切運動過程中會在顆粒體后上方產生較大的紊動粘性系數，其分布結構抑制了流體渦旋結構的發展。在流體紊動耗散作用下，流體的運動方向傾向于向前發展，進而也促進了顆粒體的運動。因而紊流盡管使得流體的渦量和速度減少了，但在運動方向上有利于流體和顆粒體的發展。

圖10 不考慮紊流模型 t =20 s 時刻稀疏細、粗顆粒流中流體速度 Uf 和顆粒速度 Us 比較(黃色和白色虛線分別為 =1 和 =0.1的等值線)Fig.10 Comparisons of the fluid velocity Uf and the solid velocity Us in dilute fine and coarse particle flows at t=20 swithout the turbulence model (the yellow and white dash lines are contour lines with and =0.1)

圖11 模擬與實驗中的顆粒堆積表面比較Fig.11 Comparisons of simulated and measured granular deposit surfaces

3 致密顆粒流研究

3.1 模型驗證

本文同時對不同粒徑顆粒堆積體的坍塌過程進行了模擬，以對該模型在致密顆粒流問題中的有效性進行驗證。其中細、粗顆粒粒徑分別為ds=0.84 mm和ds=6.02 mm，顆粒摩擦角分別為?s= 2 6°和 ?s= 3 1°，密度為 ρs=3.6×103kg/m3，初始體積分數為=0.59，實驗設置詳見圖2(b)。在模擬中，細、粗顆粒體中的取值分別為2.5 kg/(m·s)和0.5 kg/(m·s) ，回彈系數均取為e=0.8；實驗中流體為水，密度取為 ρf=1000 kg/m3，粘性系數μf=0.001 kg/(m·s)。經網格敏感性分析，模擬選取5 mm×5 mm的均勻網格進行計算。

3.2 拖曳力的影響

為了對細、粗顆粒體的坍塌過程開展進一步分析，圖12 給出了不同時刻細、粗顆粒與流體的相對速度Us?Uf分布，其中虛線表示顆粒體的堆積表面。從中可以看到，細顆粒坍塌體中顆粒與流體的相對速度較小，且最大值僅集中于顆粒體表面；粗顆粒坍塌體中顆粒與流體的相對速度則較大，且在顆粒體內部也會存在較大的相對速度。

圖12 顆粒體與流體相對速度 U s ?Uf 的變化Fig.12 Evolutions of the relative velocity between solid and fluid phases Us ?Uf

進一步地，對流體作用在顆粒體上的相對拖曳力K(Uf?Us)/ρfg進行分析，所得結果如圖13所示。在t?=1時刻，流體拖曳力的較大值均分布在顆粒體堆積表面附近，其中細顆粒堆積表面附近的流體拖曳力朝向外側，有利于顆粒體的向外發展；而粗顆粒體受到的流體拖曳力則朝向左上方，與顆粒體坍塌方向相反；細、粗顆粒體在內部受到的流體拖曳力則均朝向左上方，對顆粒體的運動起阻礙作用。在t?=3 和t?=5時刻，細顆粒堆積表面附近的仍存在較大的流體拖曳力，其方向仍然朝向外側；而粗顆粒堆積表面附近的流體拖曳力則明顯更小，其中=3時流體拖曳力對顆粒體前端仍然起阻礙作用；細、粗顆粒體內部的拖曳力方向則均呈現較為復雜的局部變化，既可能促進顆粒體流動也可能阻礙顆粒體的運動。

圖13 流體與顆粒體相對拖曳力 K(Uf ?Us)/ρfg 的變化Fig.13 Evolutions of the relative drag force between solid and fluid phases K(Uf ?Us)/ρfg

3.3 流體動壓的影響

如圖14 所示，受到顆粒坍塌體的驅動，周邊流體也會產生運動。在=1時刻，顆粒體和流體的最大速度均集中于右前方，并且可以看到流體在運動過程中發展出了明顯的渦旋結構，即在顆粒體右上方形成了流動中心，該處流體速度較小而周圍流體速度較大。在t?=3時刻，顆粒和流體的流動速度明顯變大，其最大值仍集中于右前方，且在顆粒堆積表面上方形成了2 個渦旋區域。在t?=5時刻，顆粒體的速度明顯減小，但顆粒體表面上方的渦旋區域仍維持較大的運動強度，可繼續影響顆粒體的堆積形態。

圖14 細顆粒體坍塌過程中顆粒速度 Us 和流體速度 Uf的變化Fig.14 Evolutions of the particle velocity Us and the fluid velocity Uf during the collapse process of fine particles

考慮到流體壓強的變化可能對顆粒體的運動產生影響，這里去除流體靜壓的影響，定義流體相 的 動 壓 強=pf?ρfg(Hw?y)。圖15 顯 示 了細、粗顆粒體坍塌過程中流體動壓及相對動壓梯度力的變化。從圖15 可以看到，在t?=1時流體正壓大多集中于顆粒體前端，負壓大多集中于顆粒體的上方，且細顆粒中的動壓絕對值和動壓梯度力明顯更大，這是流體受到的顆粒阻力越大，產生的動壓越不容易消散的緣故。t?=3 和t?=5時顆粒體內部的正壓逐漸消散，但在顆粒堆積表面上方則發展出了2 個負壓區域，相應的相對動壓梯度力大小可達1 以上，即與靜水壓梯度力相比不可忽略。值得注意的是，t?=1時坍塌體前端正壓區域內動壓梯度力促進了顆粒體的向前發展；t?=3 和t?=5時發展出的負壓區域則主要對顆粒表面形態的塑造起作用，對顆粒堆積體前端的影響較小。

圖15 坍塌過程中流體動壓 (云圖)和動壓梯度力 ? ?/ρfg(箭頭)的變化Fig.15 Evolutions of the dynamic fluid pressure (cloud charts) and the gradient of dynamic fluid pressure??/ρfg(arrows) during the collapse process

4 結論

本文應用綜合考慮不同流態下顆粒體本構關系、流體與顆粒體相間作用以及流體紊動影響的統一兩相流模型，對稀疏和致密顆粒流問題開展了數值模擬研究，揭示了不同粒徑顆粒體與流體的相互作用機制。得到以下結論：

(1) 該兩相流模型可以很好地模擬稀疏和致密顆粒體的運動過程，對稀疏顆粒體的前端運動距離和致密顆粒體的堆積表面均能很好地刻畫，證明了該模型在稀疏和致密顆粒流問題中的有效性。

(2) 對稀疏顆粒流的數值研究結果表明，不同粒徑顆粒體與流體的相間作用有明顯的差異。其中，細顆粒體與流體的相間速度明顯更小，受到的沉積作用更弱，更容易在流體中保持懸浮狀態，因而稀疏細顆粒體及周圍流體所能維持的運動強度要更大，顆粒體的前端運動距離也要更長。

(3) 紊流在稀疏顆粒體的運動過程中起著重要作用。模擬結果表明：紊動粘性系數與流體自身粘性系數的比值可達100 以上，考慮紊流后流體的運動明顯受到紊動粘性系數的耗散作用，流體渦旋結構的發展受到了抑制，流體則傾向于向前運動，進而促進了稀疏顆粒體的發展。

(4) 流體與顆粒的相間拖曳力在不同粒徑致密顆粒體中有著不同的分布。其中細顆粒體堆積表面在較長時間內均有著較大的拖曳力，其方向朝向堆積表面外側，促進了顆粒體的向外發展；粗顆粒體堆積表面的拖曳力則在初始時刻較大，阻礙了顆粒體的運動，但在后續坍塌過程中逐漸減小。

(5) 流體動壓對致密顆粒體的運動起著重要作用，且在細顆粒體中更為明顯。顆粒坍塌過程中產生的流體動壓梯度力與靜水壓梯度力相當，其中，正壓在初始時刻促進了顆粒體前端的發展，負壓則主要影響了顆粒體堆積表面的塑造。