垂直管內固-液兩相流全耦合CFD-DEM模型研究

張德勝 周 游 趙睿杰 施衛東

(1.江蘇大學國家水泵及系統工程技術研究中心， 鎮江 212013； 2.南通大學機械工程學院， 南通 226019)

0 引言

隨著我國城市化發展越來越快，基礎設施中排污設施的建設也越來越多[1-3]。污水管網作為城鎮、工業和農業區的重要排污設施，對當地生活和經濟發展起著重要的保障作用。

生活污水和工業廢水中往往含有纖維雜質、紙屑和顆粒等固體物，其中較輕的固體漂浮在水面并隨污水漂流，較重的固體分布在水流斷面上并呈懸浮狀態流動，造成固體物之間發生纏繞、摩擦和碰撞，或者與管道壁面發生摩擦碰撞，導致管道的磨損和堵塞，使得整個排污系統無法正常工作[4-6]。因此，對污水污物的兩相流動展開研究，可對排污系統的設計、安全可靠的運行以及后期的維護起到重要的指導作用。

固-液兩相輸送系統的流態極其復雜，受各工況參數影響[7-10]，其中流體的流動特性、固體物的運動特性以及流體-固體之間復雜的耦合機制受到廣泛關注。目前在固-液兩相流領域的主要研究方法包括實驗研究和數值模擬。過去，許多學者基于實驗數據將壓降或者其他過程參數與操作參數(如顆粒體積分數、輸送速度和顆粒直徑)關聯起來[11-13]，但這種方法只在特定條件下才具有適用性。此外，這種分析結果一般為流體的宏觀信息，無法通過其深入了解固-液兩相流中的微觀信息。隨著計算機技術的發展，數值模擬已經成為研究固-液兩相流的另一個重要手段。數值模擬方法可以分為Euler-Euler方法[14-17]和Euler-Lagrange[8-9,18]方法。在前者中，液相和顆粒相均被視為連續介質，求解時均Navier-Stokes方程，但該方法很大程度上依賴于本構關系，不能獲得足夠的關于顆粒-顆粒/壁面碰撞的信息。而在Euler-Lagrange方法中，流體被視為連續相，直接求解時均Navier-Stokes方程，而顆粒被視為離散相，通過計算流場中每個顆粒的運動，從而在仿真過程中能夠很好地捕捉顆粒的碰撞信息。流體和顆粒之間的相互作用通過在各個方程中添加源項實現。

在實驗方面，文獻[7]利用激光多普勒測速儀(LDA)測量了垂直管道內固-液兩相的速度分布，并且發現重顆粒趨向于集中在管道中心，而輕顆粒趨向于集中在管壁附近。在數值模擬方面，文獻[8]以NaCl顆粒和鹽水為載體進行了參數化研究，考察操作條件對水平管中NaCl顆粒兩相流動的影響。文獻[9]采用CFD-DEM (Computational fluid dynamics-discrete element method)方法，對垂直管內不同參數對顆粒水力輸送的影響進行研究，但未考慮顆粒對湍流的影響。文獻[16]利用Euler-Euler方法對顆粒-液漿流動進行模擬，結果表明，隨著顆粒尺寸增大，摩擦壓力損失減小，顆粒體積分數的徑向分布增大，且管道底部顆粒體積分數增高。文獻[19]比較具有不同葉輪型式潛水排污泵的通過性能，發現雙流道污水泵內顆粒速度更快，顆粒聚集程度更小，泵的通過性能強于雙葉片污水泵。文獻[20]對后掠式雙葉片污水泵固液兩相流場進行數值模擬，分析了固相顆粒的分布情況和泵內過流部件的磨損情況，并探討了顆粒體積分數及顆粒大小對泵水力性能的影響規律。

目前大部分模型采用傳統Saffman和Magnus升力模型，且忽略其他次要力以及顆粒-湍流的相互作用。本文提出一種全耦合CFD-DEM模型，考慮所有次要力和顆粒-湍流的相互作用，分析進口顆粒體積分數、輸送速度和顆粒直徑對固-液兩相流動的影響。

1 數學模型

1.1 流體控制方程

流體被視為連續相，其控制方程是質量守恒和動量守恒方程，通過在動量方程中添加源項Fp-f來考慮顆粒對流體的影響，公式為

(1)

(2)

(3)

式中ρf——流體密度，kg/m3

t——時間，s

p——壓力，Pa

μf——流體動力粘度，Pa·s

μt——流體湍流粘度，Pa·s

g——重力加速度，m/s2

Fp-f——顆粒對流體的力源項，N/m3

Ff-p——流體對顆粒的總作用力，N

Vcell——流體網格體積，m3

n——流體網格內顆粒數

αf——流體體積分數，%

顆粒被劃分為64個部分，計算流體網格中的流體體積分數，公式為

(4)

式中k——與該網格重疊的顆粒數

Vi,part——顆粒i與該網格重疊的總體積，m3

采用標準k-ε湍流模型來封閉方程，通過在方程中添加湍動能源項和湍流耗散源項來考慮顆粒對湍流的影響[21]，公式為

(5)

(6)

(7)

(8)

(9)

式中Sk——湍動能源項，kg/(m·s3)

Sε——湍流耗散源項，kg/(m·s4)

vf——流體瞬時速度，m/s

vp——顆粒速度，m/s

dp——顆粒直徑，m

νf——流體運動粘度，m2/s

fp——顆粒雷諾數有關的系數

Rep——顆粒雷諾數

C′ε3——顆粒雷諾數有關的系數

其中C′ε2=1.92,Cεp=0.085。

1.2 顆粒運動方程

顆粒運動采用離散元法(DEM)求解，粒子的運動由牛頓第二定律描述，公式為

(10)

(11)

式中mp——顆粒質量，kg

Fg——重力，N

Fdrag——Dif曳力[22]，N

Floth——Loth升力[23]，N

Fp——壓力梯度力，N

Fvm——虛擬質量力，N

Fc——顆粒-顆粒/壁面間相互作用力，N

Ip——轉動慣量，kg·m2

ωp——顆粒角速度，s-1

Tc——顆粒接觸力矩，N·m

Tf——流體作用力矩，N·m

許多因素都會導致升力的產生，而流體渦度和顆粒旋轉是產生升力的兩個主要機制。為表示流體渦度和顆粒旋轉對升力的影響，Loth升力模型將升力系數考慮為剪切和旋轉引起的升力系數的組合，這擴大了升力模型的適用范圍。根據文獻[23]，升力表示為

(12)

式中Cl——Loth升力系數

ωf——流體渦度，s-1

為比較不同升力模型的效果，本文給出Saffman和Magnus升力的表達式。流體剪切運動引起的Saffman升力[24]Fsaf表示為

(13)

式中Cls——Saffman升力系數

根據文獻[25]，升力系數Cls可表示為

(14)

(15)

顆粒旋轉引起的Magnus升力計算公式[26]為

(16)

(17)

(18)

(19)

式中Rer——顆粒旋轉雷諾數

ωf-p——相對角速度，s-1

當顆粒密度遠大于流體密度時(ρf/ρp?1)，虛擬質量力和壓力梯度力并不重要(氣-固流)。當ρf/ρp大于0.1時，虛擬質量力和壓力梯度力變得顯著，并且需要被考慮，即

(20)

(21)

式中Vp——顆粒體積，m3

Cvm——虛擬質量系數，默認值為0.5

1.3 離散隨機游走模型

(22)

(23)

式中ξi、ξj、ξk——i、j、k方向的高斯分布隨機數

u′、v′、w′ ——i、j、k方向的脈動速度分量，m/s

由于RANS模型的簡單性，假設湍流為各向同性，計算脈動速度的均方根為

(24)

式中kf——流體湍動能，m2/s2

將上述流體的瞬時速度用于以上各種模型的計算中，所以在計算顆粒的運動時也考慮了湍流效應。

2 模型設置

2.1 數值格式

控制方程用有限體積法離散，采用標準k-ε湍流模型，采用SIMPLE算法處理速度-壓力耦合。在CFD-DEM模擬中，顆粒時間步長Δtp通常為Rayleigh時間步長ΔtR的10%～40%，并且流體時間步長需為顆粒時間步長的10～100倍。綜合考慮計算精度和時間成本，固體顆粒時間步長設為2×10-5s，流體時間步長設為2×10-4s。ΔtR計算式為

(25)

式中Ri——顆粒i半徑，m

vi——顆粒i泊松比

Yi——顆粒i楊氏模量，Pa

2.2 模型設置及邊界條件

建立與文獻[7]實驗相同的管道，并且進行網格劃分，如圖1所示。為更好捕捉壁面處流體信息，對邊界層網格進行了細化，壁面第1層網格厚度為0.26 mm，且滿足30

表1為模擬中的顆粒和流體屬性[9]。表2為本文所有的模擬情況，研究了進口顆粒體積分數Cv、輸送速度Uf和顆粒直徑dp對固-液兩相流動的影響。

表1 顆粒和流體屬性Tab.1 Particle and fluid properties

表2 本文模擬中使用的參數Tab.2 Parameters used in this simulation

3 結果與討論

3.1 模型驗證

為檢驗模型的正確性，對直徑為2.32 mm的球形顆粒進行數值模擬，并將結果與文獻[7]的實驗進行比較，圖2(圖中r表示測量點到管道中心的徑向距離，R表示管道半徑，并進行無量綱處理)所示為全耦合模型的驗證。由圖2可知，無升力的模型不能很好地預測管道中心處顆粒的分布，且管道壁面附近的顆粒體積分數略高于實驗。傳統的Saffman和Magnus升力模型預測的顆粒體積分數在管道中心和壁面處誤差很大。而考慮Loth升力的模型能夠較好地預測管道內顆粒體積分數的分布，同時各相軸向速度和流體均方根速度的模擬結果與實驗吻合較好，表明考慮Loth升力和湍流調制的全耦合模型能夠較好地預測顆粒的運動和流體的關鍵流動特性。

圖2 模型驗證Fig.2 Model validation

本文建立了網格數分別為6×104(粗網格)、9×104(中等網格)和1.1×105(細網格)的3套模型，對比3套網格預測的顆粒體積分數及各相軸向速度分布，發現粗網格和細網格與參考網格的誤差均小于3.0%。通過對比不同模擬時間(4、5、6 s)下的計算結果，發現顆粒體積分數及各相軸向速度的誤差均小于2.5%。所以，綜合考慮計算精度和成本，本文采用中等網格和總模擬時間為4 s的方案來進行數值計算。

3.2 進口顆粒體積分數對水力輸送的影響

進口顆粒體積分數Cv是固-液兩相排污系統中最重要的參數之一，其直接影響顆粒的運動以及流體的流動特性。研究Cv為1.0%、2.5%和5.0%的影響，分別對應表2中序號1、2和3。

圖3 不同Cv下顆粒分布Fig.3 Distribution of particles at different Cv values

不同Cv下顆粒的分布情況如圖3所示。從圖3中看出，當Cv為1.0%時，顆粒集中在管道中心處，而管道半徑方向的中部(r/R=0.5)和壁面附近(r/R=0.9)顆粒存在較少，甚至沒有。當Cv增加到5.0%時，顆粒同樣主要集中在管道中心，但半徑方向的中部和壁面附近均存在一定的顆粒。這是因為隨著進口顆粒體積分數的增加，管道內顆粒的數量急劇增加，顆粒更加分散，同時更多的顆粒由于升力作用向管道中心聚集。

圖4 不同Cv下碰撞次數Fig.4 Number of collisions at different Cv values

碰撞次數可分為顆粒間和顆粒-壁面間碰撞次數，在z=1.7～2.2 m區域內提取碰撞次數。不同Cv下，碰撞次數隨時間變化如圖4所示?？梢钥闯觯w粒間碰撞從t=0.37 s左右才開始，這是因為顆粒需要經過一定時間運動到所研究區域內，而在t=0.5 s以后，顆粒間碰撞次數和顆粒-壁面間碰撞次數均已穩定，即顆粒兩相流已經充分發展，考慮到穩定計算和其他工況，對2～4 s之間的數據進行平均處理(虛線所示)。從圖4a、4b可以看出，顆粒間和顆粒-壁面間碰撞次數都隨著Cv增加而增加，這是因為顆粒數量的增加使得更多顆粒向管道中心聚集，管道內發生更加頻繁的碰撞。同時，管壁附近顆粒體積分數有所增加，該處的顆粒開始與壁面發生一定次數的碰撞。但顆粒-壁面間碰撞次數遠低于顆粒間碰撞次數，且顆粒-壁面間碰撞次數的波動幅度隨Cv的增加而增加。

圖5 不同Cv下各物理量分布Fig.5 Distribution of various physical quantities at different Cv values

顆粒體積分數及各相的速度分布如圖5所示。由圖5可知，顆粒趨向于聚集在管道中心，顆粒體積分數在管道中心出現最大值，并沿管徑向壁面逐漸降低，在壁面減為零。流體和顆粒的軸向速度在管道中心存在最大值，并且沿管徑向壁面降低，而流體均方根速度趨勢相反，其在管道中心處最小，并沿管徑向壁面逐漸增大，到壁面附近處增加到最大值。

由圖5a可見，隨著Cv增加，整個管道內顆粒體積分數都增大，且管道中心處顆粒體積分數增加最為明顯。這是因為Cv增大直接導致顆粒數量增多，更多的顆粒在升力作用下向管道中心聚集。而顆粒間碰撞頻率隨著顆粒數量的增加而顯著增加，更多的碰撞使得顆粒更加分散，這也就能解釋管道半徑方向的中部和壁面附近顆粒體積分數的增加。

由圖5b、5c可見，隨著Cv增加，流體軸向速度有所降低，這種變化在管道中心區域更明顯，在管壁附近忽略不計。這是因為隨著Cv增加，流動狀態從以流體為主導的流動模式轉變為流體-顆粒作用和顆粒-顆粒作用為主導的流動模式，更多數量的顆粒使得流體-顆粒間的動量交換加劇，導致流體軸向速度降低，而靠近壁面區域顆粒較少，兩相間動量交換較少使得靠近壁面區域的流體軸向速度幾乎不變。隨著Cv增加，顆粒軸向速度也出現了一定的降低，這是由于顆粒的動能被更加頻繁的顆粒間相互作用所耗散。由圖5d可以看出，顆粒對湍流的調制作用導致流體均方根速度有所降低，且調制作用隨Cv增加更加明顯。

圖6 Cv對無量綱壓降的影響Fig.6 Effect of Cv on dimensionless pressure drop

3.3 輸送速度對水力輸送的影響

輸送速度Uf與輸送系統(管道、過濾系統和泵系統等)的流量和效率密切相關，是兩相流動系統中另一重要的參數。本節研究了Uf為1.5、3.0、4.5 m/s對水力輸送特性的影響，分別對應表2中序號4、5和6。

不同Uf下顆粒分布情況如圖7所示。從圖7看出，在Uf=1.5 m/s時，顆粒在整個管道內分布最均勻，并且壁面附近處存在較多的顆粒。當Uf增加到3.0 m/s時，顆粒的分布不均勻，在管道中心處較為集中，且管道壁面附近的顆粒數減少。當Uf繼續增加到4.5 m/s時，顆粒的分布更加不均勻，大部分顆粒集中在管道中心，壁面附近幾乎沒有顆粒。因此，隨著Uf增加，管道壁面附近的顆粒逐漸減少，更多的顆粒向管道中心聚集，顆粒在管道內的分布逐漸不均勻。

圖7 不同Uf下顆粒分布Fig.7 Distribution of particles at different Uf values

不同Uf下，碰撞次數隨時間的變化如圖8所示?？梢钥闯?，在不同Uf下，開始出現顆粒間碰撞和顆粒流達到穩定的時間不同。在較低的Uf下，開始出現顆粒間碰撞的時間和顆粒流達到穩定的時間都較長，分別為1.1 s和1.5 s。這是因為顆粒從管道入口運動到出口的時間由顆粒速度決定，較低的Uf導致顆粒速度較低，顆粒到達管道出口的時間就較長。從圖8還可以看出，顆粒間碰撞次數和顆粒-壁面間碰撞次數都隨Uf的增加而增加，這是因為隨著Uf的增加更多的顆粒向管道中心聚集，中心區域內發生更為頻繁和劇烈的碰撞，導致少部分顆粒在碰撞的作用下向管壁運動，并與管壁發生碰撞。同樣，顆粒間碰撞次數遠大于顆粒-壁面間碰撞次數,且顆粒-壁面間碰撞次數的波動幅度隨Uf增加而增加。

圖8 不同Uf下碰撞次數Fig.8 Number of collisions at different Uf values

圖9 不同Uf下各物理量分布Fig.9 Distribution of various physical quantities at different Uf values

圖9為不同Uf下各物理量分布。由圖9a可見，隨著Uf增加，管壁附近顆粒體積分數逐漸降低，顆粒體積分數峰值逐漸向管道中心靠近且逐漸增加，當Uf為4.5 m/s時顆粒體積分數在中心存在最大值。這是因為隨著Uf增加，壁面附近滑移速度增加，顆粒受到的升力變大，推動更多的顆粒向管道中心聚集。

由圖9b、9c可見，隨著Uf增加，流體和顆粒的軸向速度也隨之增加，其徑向梯度也隨輸送速度的增加變得更加陡峭，這在管壁附近更加明顯。輸送速度對流體均方根速度的影響如圖9d所示。流體均方根速度隨Uf增加而增加，且陡峭程度逐漸增加。這是因為隨著Uf增大，管道內雷諾數增大，從而湍流強度總體變大。在壁面附近，流體均方根速度的增加更加明顯。

圖10 Uf對無量綱壓降的影響Fig.10 Effect of Uf on dimensionless pressure drop

由圖10可見，隨著Uf增加，管內壓降也快速增加，這與文獻[29]中的實驗一致。這是因為隨著Uf增加，流體-壁面的摩擦損耗加劇，以及更加頻繁和劇烈的粒子間相互作用，導致更多的能量被耗散，從而造成更大的壓降。

3.4 顆粒直徑對水力輸送的影響

在排污設施中，輸送的固體污物往往大小各異。因此，研究不同粒徑的顆粒對兩相流動的影響具有重要的意義。對相對較粗的顆粒(dp為1.0、1.5、2.0 mm)進行數值計算，研究顆粒直徑對固-液兩相流動的影響。

不同dp下顆粒的分布情況如圖11所示。從圖11 看出，當dp較小(dp為1.0 mm)時，顆粒在管道內分布非常均勻，并且壁面附近也存在一定量的顆粒。當增大顆粒直徑到1.5 mm時，顆粒較為集中，管道壁面附近的顆粒減少。當dp繼續增加到2.0 mm時，顆粒分布極不均勻，大部分顆粒集中在管道中心，壁面附近幾乎沒有顆粒。因此，從顆粒的分布可以看出，小顆粒在管道內分布均勻，而大顆粒在管道中心聚集，且壁面附近幾乎不存在顆粒，表現出極不均勻性。

圖11 不同dp下顆粒分布Fig.11 Distribution of particles at different dp values

在不同dp下，碰撞次數隨時間的變化如圖12所示。從圖12a看出，顆粒間碰撞次數隨dp增加而減少。隨著dp增加，相同濃度下顆粒數量的急劇減少會降低顆粒間碰撞次數，而大顆粒在管道中心更加聚集會增加中心處顆粒間碰撞次數，由此可見，顆粒數量的急劇減少對顆粒間碰撞次數的影響大于顆粒向中心聚集帶來的影響。從圖12b中看出，在顆粒直徑較小時，顆粒分布較為分散，顆粒-壁面間碰撞次數比顆粒間碰撞次數少，但未到可忽略的程度。在顆粒直徑較大時，顆粒主要集中在管道中心處，顆粒-壁面間碰撞次數遠小于顆粒間碰撞次數。

圖12 不同dp下碰撞次數Fig.12 Number of collisions at different dp values

圖13為不同dp下各物理量的分布。由圖13a可知,同濃度下，顆粒直徑越小，在管道內的分布越分散。這是因為小顆粒在管道內的顆粒數量較大，管道中心處顆粒間的接觸次數增加，使得顆粒更分散，導致近壁面處存在一定顆粒。當輸送速度不變時，流體軸向速度徑向梯度基本一致，但是大顆粒與流體滑移速度更大，顆粒受到的升力就更大，管道壁面處顆粒向中心聚集，這也能解釋大顆粒更加向管道聚集的現象。

由圖13b、13c可知，流體軸向速度隨著dp增加呈現略微降低的趨勢，這種趨勢在管道中心較為明顯。同時發現，隨著dp增加，顆粒速度也存在一定量的降低。這是因為隨著dp增加，流體對大顆粒的夾帶作用變小，即大顆粒的隨動性較低，因此顆粒速度較低。由圖13d可知,dp對流體均方根速度的影響很小。

由圖14可知，隨著dp增加，壓降的變化不大，這在文獻[30]的實驗中也有相同的報道。因此，輸送速度對壓降的影響大于進口顆粒體積分數，而顆粒直徑的影響最小。

4 結論

(1)升力驅使顆粒向管道中心聚集，并且Loth升力比傳統Saffman和Magnus升力更加適合固-液兩相流的研究?？紤]湍流調制的全耦合模型能夠較好地預測顆粒的運動以及流體的流動特性。

圖13 不同dp下各物理量分布Fig.13 Distribution of various physical quantities at different dp values

圖14 dp對無量綱壓降的影響Fig.14 Effect of dp on dimensionless pressure drop

(2)隨著進口顆粒體積分數增加，更多顆粒聚集在管道中心，顆粒分布更加分散，同時顆粒間碰撞次數急劇增加。在管道中心處，流體軸向速度和顆粒軸向速度隨進口顆粒體積分數增加而逐漸減小，而管壁附近變化很小，并且流體均方根速度在整個管道范圍內有一定的下降趨勢。此外，壓降隨著進口顆粒體積分數的增加而增加。

(3)隨著輸送速度增加，管道中顆粒更加聚集，顆粒間碰撞次數逐漸增加，顆粒體積分數峰值逐漸向管道中心移動且峰值逐漸變大，同時流體軸向速度、顆粒軸向速度以及流體均方根速度均增加，且速度梯度在壁面附近更加陡峭。此外，壓降隨著輸送速度的增加而增加，且影響比進口顆粒體積分數大。

(4)在所研究的顆粒直徑范圍內，隨著顆粒直徑增加，顆粒分布逐漸不均勻，顆粒體積分數在管道中心處存在的峰值逐漸升高且沿管徑向壁面的下降幅度更大。顆粒直徑對液體軸向速度和均方根速度分布的影響較小。顆粒軸向速度隨顆粒直徑的增加而降低。此外，顆粒直徑對壓降的影響最小。