劉唐京，王企鯤，鄒 赫

（上海理工大學 能源與動力工程學院，上海 200093）

引言

均勻的稀釋懸浮液以層流方式通過直圓管時，流體中的剛性球形顆粒最終會遷移到半徑約為0.6R（R為管道半徑）的圓環(huán)上，這個環(huán)也被稱為Segre-Silberberg 環(huán)[1].這表明顆粒在隨流體運動時，除受到主流方向的驅動力外，在徑向上還受到一個升力使得顆粒發(fā)生遷移.這個徑向上的升力是由于運動流體的慣性引起的，因此又被稱為“慣性升力”，由其“慣性升力”而引發(fā)的顆粒聚集現象稱為“慣性聚集”[2-3].

這種現象是在低Re數（Re為Reynolds 數）下管道流中發(fā)現的，之后有學者在高Re數流中也發(fā)現了顆粒的慣性聚集現象.Matas 等[4]對高Re數下懸浮粒子在Poiseuille 流中的“慣性聚集”進行實驗研究，得到了不同Re數下管道直徑與顆粒直徑比在8～42 范圍內顆粒的聚集位置.Morita 等[5]通過實驗發(fā)現，Re數在1 000范圍內的稀釋懸浮液通過直圓管時，只有一個穩(wěn)定的聚集點.由于實驗很難獲取流場中的各項力學參數，隨著數值模擬（CFD）在流體力學中的廣泛運用[6-8]，一些學者為了揭示顆粒所受慣性升力的形成特性，通過數值計算方法對顆粒的“慣性聚集”進行了研究[9-10].但由于流動是非定常的，計算模型的數學描述比較困難，Carlo 等[11]首次提出“相對運動模型”，將原本的非定常問題轉化為準定常問題.基于相對運動模型，王企鯤[12-13]對方管中顆粒的受力特性進行了數值研究，揭示了顆粒聚集的力學成因并歸納出顆粒穩(wěn)定聚集點的水動力學判據.

相對運動模型的提出，極大地簡化了計算模型.但采用相對運動模型進行數值模擬時，為了使流動達到充分發(fā)展狀態(tài)，一般需要預留L=0.058Re·D（D為管徑）[14]的管長以消除入口段的影響.這在計算高Re數工況時就需要非常長的管道，造成網格數量大，使得數值計算變得困難.考慮到管內流動屬于周期性流動，因此本文對管道進、出口設置周期性邊界條件，解決了高Re數工況下管道較長的問題，同時計算高Re數管流中顆粒慣性聚集的受力特性.

1 計算模型與計算方法

1.1 計算模型

本文考慮單個剛性球形顆粒在直圓管Poiseuille 流中運動，管道直徑為D，長度為L，小球直徑為d，如圖1所示.目前比較常用的模型是6 自由度模型，然而這種計算模型是非定常的，模型的數學描述比較復雜，并且計算過程繁瑣也很耗時，此外6 自由度模型雖然能獲得顆粒的運動軌跡，但無法獲取顆粒升力在通道內的空間分布規(guī)律.本文采用的相對運動模型是在6 自由度模型的基礎上進行簡化，它雖然無法獲得顆粒的運動軌跡，但能計算出顆粒在通道內的受力特性.通過分析受力特性來確定顆粒的聚集點，這在文獻[11-13]中給出了比較詳細的討論.

圖1 計算模型示意圖Fig.1 The sketch for the numerical model

相對運動模型是通過創(chuàng)建一個慣性坐標系將顆粒的平移速度轉移給管壁，使得計算時的流場為準定常，從而簡化計算.即將計算坐標系設置在顆粒的中心（如圖1所示）并隨顆粒一同平移，則在此慣性坐標系下，顆粒的平移速度為零，僅存在以原點為中心的旋轉速度，管壁則以顆粒的速率反向移動.

1.2 計算方法

本文基于有限體積法進行三維數值模擬，求解穩(wěn)態(tài)Navier-Stokes 方程，其控制方程如式（1）所示.流動介質為常溫常壓液態(tài)水，考慮到顆粒是懸浮顆粒，取其密度與水相同，為ρ=1 000 kg/m3.在CFD 計算中，采用雙精度來進行計算，壓力與速度的耦合采用coupled 算法，壓力方程的離散采用二階格式，動量方程的離散采用三階精度的QUICK 格式[14-15]：

式中，u為平移坐標系Oxyz中流體的相對速度，p為壓強，ρ 和υ 分別為流體的密度和運動黏度.

在CFD 計算中，顆粒的平移速度UP被轉化為管壁運動的邊界條件，這里需要通過試湊的方式來獲得顆粒在該位置的最終恒定運動速度.首先假定一個速度UP0進行試算，輸出顆粒在沿流動方向（x方向）上的合力Fx0，在一定精度內，判斷合力是否為零（零指的是零量階，不是數值上為零）；若不為零則需不斷地迭代更新壁面速度，直到顆粒在流動方向上合力為零.同時，顆粒的旋轉速度也用同樣的方式進行試湊，直到顆粒在所有方向上轉矩為零.此時，顆粒在y方向上所受的力便為徑向升力.

為了提高計算效率，平移速度的初始參考值為UP=2U[1?(r+)2]，旋轉速度的初始參考值為對于UP和ω 的更新，本文利用具有超線性收斂性的“割線法”更新下一步迭代數據，分別由式（2）、（3）計算，采用這種方法通常試湊3 至4 次就能得到結果，而試湊一次也只需迭代150 次左右：

式中，UPn和ωPn分別為平移速度和旋轉速度，Fxn和MPn分別為阻力和轉矩，n=0,1,2.

顆粒的升力Fy和阻力Fx分別按式（4）、（5）計算：

式中，i和j分別為x方向和y方向的單位向量，P為應力張量，n為顆粒表面外法線單位向量，dS為面積微元，Σ 為顆粒表面.

顆粒的轉矩按式（6）計算：

式中，r為由顆粒中心指向顆粒表面的矢徑.

在CFD 計算中，當對管道進、出口采用周期性邊界條件時，計算區(qū)域內的每一處壓力分為周期性壓力和線性壓力，而在CFD 實際的計算過程中只顯現出周期性壓力少了線性壓力，因此真實的壓力應為：preal=p+β?x（p為周期性壓力，β為一個周期內的平均壓力梯度，?x=x1?x2）[14].那么，式(4)～ (6)在計算顆粒的升力、阻力以及轉矩時未包含線性壓力提供的那部分力.線性壓力對顆粒貢獻的升力、阻力及轉矩由式(7)、(8)計算，從中可以發(fā)現，線性壓力對顆粒提供的力只在流動方向上，并且對轉矩的貢獻為零，也就是說線性壓力只對顆粒的阻力有影響：

式中，I為單位二階張量.因此顆粒在流向上的實際阻力為

為了方便下文的討論與分析，本文定義如下無量綱參數：

升力系數CFL為

式中，d為顆粒直徑，D為管道直徑.顆粒的無量綱直徑為

顆粒無量綱徑向位置為

式中，r為管中心與顆粒球心的距離，R為管道半徑.Re數為

式中U為管內流體的平均速度，μ為流體的動力黏度.擾動強度為

式中，v和w分別為y方向和z方向的流速分量.

2 結果與討論

2.1 網格無關性驗證

本文采用結構化網格對計算域進行網格劃分，如圖2（a）所示.為了提高計算的穩(wěn)定性和精度，在流體與顆粒間的邊界層網格進行了加密處理，如圖2（b）所示.為了保證計算的準確性和經濟性，以及獲得顆粒的真實受力，本文先對網格無關性進行了驗證.本次驗證的計算工況為：a+=1/9，顆粒的徑向位置r+=0.1，管長L=4D.

圖2 網格示意圖：（a）管道網格；（b）顆粒周圍網格Fig.2 Schematic diagrams of the grids:a）pipeline grid;b）grids around particle

從圖3 中可看出，當網格數量達到50 萬左右時，顆粒的升力系數基本保持不變，考慮到在滿足計算精度的同時盡可能節(jié)約計算時間，本文最終確定計算域的網格數量控制在60 萬左右.對于下文采用不同管長計算的工況，網格劃分將以本次驗證的管道網格為基準，網格尺寸保持一致.

圖3 網格無關性驗證Fig.3 Grid independence verifications

2.2 周期性邊條的可行性分析

在低Re數下，文獻[2-13]基于相對運動模型，邊界條件為進口給定均勻的相對速度，出口為壓力出口.而在高Re數下，這種邊界條件通常需要很長的管道才能計算出可靠的結果，因此考慮采用周期性邊界條件來減小管長.為了驗證周期性邊界條件的可行性，本文對Re=350，a+=1/9 的工況進行了數值模擬.本次模擬分別選取管長L=4D和L=50D的管道，L=4D的管道對進、出口采用周期性邊界條件，L=50D的管道則與文獻[13]采用相同的邊界條件，模擬結果如圖4所示.

由圖4 可知，兩種邊界條件得到的升力系數曲線完全重合，升力系數為零且該點一階導數小于零的點即為顆粒的穩(wěn)定聚集點，因此Re=350 時a+=1/9 的顆粒主要聚集在r+≈0.76 處.而文獻[4]的實驗結果為r+≈0.77，兩者比較吻合，這說明周期性邊界條件是可行的.當Re=350 時，對進、出口采用周期性邊界條件只需4 倍管徑長度的管道便可計算出結果，而用文獻[13]中的邊界條件卻需要50 倍管徑長度.因此在求解高Re數流中顆粒的慣性聚集時，采用周期性邊界條件可以有效地減小管長，降低計算量.

圖4 不同邊界條件的模擬結果Fig.4 Simulation results under different boundary conditions

2.3 周期長度的確定

當計算高Re數工況采用周期性邊界條件時，需知道管長L（周期長度）為多少時，獲得的計算結果才真實可靠.針對這個問題，本文采用不同周期長度的管道對無量綱直徑a+=1/9 的顆粒進行了模擬分析，Re=350.

從圖5 中可看出，當L≥3D時，隨著周期長度的增加，計算結果將保持不變，并且與實驗結果是吻合的（2.2 小節(jié)中已做對比），這從力的角度來看L=3D的管道便可以得到穩(wěn)定的計算結果.在相對運動模型中，顆粒是相對靜止的，這會對管中的流體產生擾動，而較短的周期長度有可能會使這種擾動延伸到邊界上，影響計算結果.而且這有可能會出現不同周期長度得到相同的升力分布，但流場卻不一定是相同的情況.為了確保計算結果的可靠性，本文對不同周期長度顆粒附近以及靠近進、出口處的擾動強度進行了對比，如圖6所示.考慮到越靠近壁面，顆粒對流體造成的擾動越強，因此本文選取顆粒靠近壁面（r+=0.8）的工況進行對比.

圖5 不同周期長度下的升力分布Fig.5 The lift distribution under different period lengths

從圖6 中可以發(fā)現，隨著周期長度的增加，靠近管道進、出口處的擾動強度是不斷減小的，當L≥4D時，靠近進、出口處的擾動基本可以忽略不計.這說明當周期長度大于4D時，流場將不會在發(fā)生變化，結合上文升力分布結果，對于Re=350 的工況，L=4D是可行的計算周期.而對更高的Re數，本文也進行了驗證，結果如圖7所示.當Re數達到800 時，計算結果也符合：1）靠近管道進、出口處擾動強度為零；2）顆粒的聚集點與文獻[4]的實驗結果相符.因此，對于Re<1 000、a+≤1/9 的工況，本文確定4D為計算周期.

圖6 不同周期長度下的擾動強度對比：（a）L=2D；（b）L=3D；（c）L=4D；（d）L=5DFig.6 Comparisons of disturbance intensities under different period lengths:a）L=2D;b）L=3D;c）L=4D;d）L=5D

圖7 Re=800 的計算結果：（a）升力分布；（b）擾動強度Fig.7 Calculation results for Re=800:a）lift distribution;b）disturbance intensity

2.4 周期性邊條的應用

2.4.1Re<1 000

前文對周期性邊條的可行性進行了驗證，并對Re<1 000 工況給出了一個可行的計算周期為L=4D.在此基礎上，本文在這里對不同粒徑的顆粒進行了數值模擬，研究不同Re數下顆粒的受力特性以及對顆粒聚集點的影響.

如圖8所示，在低Re數下，顆粒在徑向上升力分布是類拋物線，不同大小的顆粒主要聚集在r+=0.6 ～0.7 之間，與管壁有一定的距離.而隨著Re數的增大，顆粒的升力分布以及聚集點都出現了明顯的變化.主要表現如下：高Re數下顆粒的升力不再呈類拋物線分布，顆粒受到的升力具有一定的波動，在r+=0.5 ～ 0.7 之間出現一段升力相對較小區(qū)域，而在這個區(qū)域內有出現第二個聚集點（內環(huán)）的趨勢.此外，隨著Re數的增大，顆粒主要聚集在r+=0.75 ～ 0.85 之間（外環(huán)），向著壁面靠近.即Re數越大，顆粒的聚集位置越靠近壁面；而與顆粒粒徑的關系則與之相反，粒徑越大，聚集位置越向著管中心遷移.

圖8 不同Re 下顆粒的升力分布：（a）Re=50；（b）Re=350；（c）Re=500；（d）Re=800Fig.8 The lift distribution of particles under different values of Re:a）Re=50;b）Re=350;c）Re=500;d）Re=800

在本次的計算結果中，只發(fā)現了一個穩(wěn)定聚集點，這與Morita 等[5]的實驗結果是一致的，當Re<1 000 時，如果管道足夠長，內環(huán)將消失，所有的粒子都將聚集在外環(huán)上.而在Matas 等[4]的實驗中，他們在管道的上游區(qū)域發(fā)現了內環(huán)的存在，但在下游區(qū)域觀察到內環(huán)上的顆粒有向外環(huán)遷移的趨勢.本文認為這可能是Matas 等[4]實驗的管道不夠長，顆粒的遷移未完全發(fā)展，顆粒要脫離r+=0.5 ～ 0.7 這個升力相對較小的區(qū)域需要較長的時間.

2.4.2Re≥1 000

從圖8 的升力分布來看，在更高Re數下有可能出現第二個穩(wěn)定聚集點.為了探究Re≥1 000 時是否有第二個穩(wěn)定聚集點的出現，本文選用了a+=1/17 的顆粒進行模擬.

本次模擬仍是用L=4D的管道進行計算，首先對周期長度的可靠性進行驗證，結果如圖9（a）所示.從其擾動強度來看，4 倍管徑周期長度符合計算精度，且顆粒的升力分布與上文中高Re數的分布特征一樣.因此本文認為對于Re≤1 600，a+=1/17 的工況，L=4D的管道依然是可行的.

從圖9（b）顆粒的升力分布可以發(fā)現，當Re數達到1 200 時，a+=1/17 的顆粒在徑向上有三個聚集點，其中有兩個是穩(wěn)定聚集點，分別在r+≈0.63，0.87 處.這說明當Re>1 000 時，對于小粒徑的顆粒是有可能存在兩個穩(wěn)定聚集點的.由于大粒徑的顆粒在Re數達到1 000 時計算不穩(wěn)定，所以對于更大粒徑的顆粒本文沒有繼續(xù)進行深入研究.

圖9 a + =1/17，Re≥1 000 的計算結果：（a）Re=1 600 時的擾動強度；（b）升力分布Fig.9 The calculation results of a+ =1/17,Re≥1 000:a）the disturbance intensity at Re=1 600;b）the lift distribution

2.5 流場分析

為了探究低Re數和高Re數下管內顆粒慣性升力分布不同的原因，本文對顆粒所在橫截面（即x=0 截面）的流場進行了分析，以a+=1/9 的顆粒為例，如圖10 ～ 12所示.圖10 ～ 12 分別為r+=0.4，0.6，0.8 時不同Re數下z方向的速度云圖和該截面上的速度矢量圖，對z方向的速度無量綱化為.

圖10 x=0 截面的速度云圖和矢量圖（r+ =0.4）：（a）Re=50；（b）Re=350；（c）Re=500；（d）Re=800Fig.10 Velocity contours and velocity vectors of section x = 0r+ =0.4）:a）Re=50;b）Re=350;c）Re=500;d）Re=800

圖12 x=0 截面的速度云圖和矢量圖（r+ =0.8）：（a）Re=50；（b）Re=350；（c）Re=500；（d）Re=800Fig.12 Velocity contours and velocity vectors of section x = 0r+ =0.8）:a）Re=50;b）Re=350;c）Re=500;d）Re=800

從圖11、12 可以明顯地看出，在顆粒的周圍有二次流產生，而二次流可能會對顆粒的升力造成影響.從速度云圖及矢量圖來看，當Re=50 時，顆粒周圍并沒有明顯的二次流動，尤其是顆粒更靠近通道中心時，但隨著顆粒向壁面靠近，可以發(fā)現有微弱的二次流產生.然而當Re≥350 時，即使顆粒更靠近通道中心也會有微弱的二次流產生，且隨著Re數的增大以及顆粒向壁面靠近，二次流變得越來越強烈.此外從速度矢量圖可以發(fā)現：在低Re數時，二次流主要向著顆粒的下方流動；而在高Re數時，二次流在顆粒靠近通道中心時先是向顆粒下方流動，而后隨著顆粒向壁面靠近以及Re數的增大，其逐漸向顆粒左右兩側流動.

圖11 x=0 截面的速度云圖和矢量圖（r+ =0.6）：（a）Re=50；（b）Re=350；（c）Re=500；（d）Re=800Fig.11 Velocity contours and velocity vectors of section x = 0（r+ =0.6）:a）Re=50;b）Re=350;c）Re=500;d）Re=800

對此本文認為，由于顆粒周圍二次流的影響，顆粒在徑向上的升力分布才出現圖8所示的變化.在低Re數時，二次流主要向顆粒下方流動且其強度隨著顆?？拷诿娑龃?，這會給顆粒一個向下的力，而在圖8中也能明顯地看到在靠近壁面時升力下降得更快，這說明二次流對顆粒的升力是有影響的.在高Re數時，顆粒所受升力在r+=0.4 ～ 0.7 之間相對平緩且升力系數較小，這與前文所說的二次流強度隨顆粒徑向位置變化相對應.而在Re≥500，r+=0.75 時升力出現回升，本文認為這是由于二次流流向變化所引起的，在r+=0.8 時，二次流主要向顆粒兩側成對稱流動，這使得其對升力的影響減弱.

3 結論

Carlo[2]提出的相對運動模型在求解低Re數下顆粒的慣性聚集是比較成熟的.但在高Re數下，如果仍對進口給定均勻流，為了消除入口段的影響需要很長的管道，造成網格數量大，計算成本高，因此本文嘗試對管道進、出口采用周期性邊界條件以減小計算域管長.本文主要對周期性邊界條件的可行性進行了驗證并求解了高Re數下顆粒受力特性，得到了以下結論：

1）在求解高Re數流的慣性聚集，周期性邊界條件的使用可以有效地減小管長，這很大程度上提高了數值計算的效率以及經濟性.當Re<1 000 時，a+≤1/9 的顆粒用4D周期就可以計算出可靠的結果.對于粒徑細小的顆粒，如文中a+=1/17 的顆粒，4D周期可計算的Re數高達1 600.

2）在低Re數下，顆粒在徑向上的升力呈拋物線分布，且顆粒主要聚集在離壁面較遠的區(qū)域.隨著Re數的不斷增大，顆粒的聚集位置向著壁面靠近，且其升力分布出現了較大的波動，它將不再呈類拋物線分布，在r+=0.5 ～ 0.7 之間出現了一段升力相對較小的區(qū)域，而在這個區(qū)域內有出現新聚集點的趨勢.

3）當Re≤1 000 時，本文只發(fā)現了一個聚集點，新的聚集點并沒有出現.但當Re>1 000 時，本文用a+=1/17 的顆粒進行計算得到了兩個穩(wěn)定的聚集點，這說明在高Re數流中小粒徑的顆粒有可能出現兩個穩(wěn)定的聚集點.

4）顆粒周圍有二次流的產生，其強度隨著Re數的增大而增大，且隨著顆粒越靠近壁面，二次流的強度也會增加.在低Re數時，二次流主要向著顆粒的下方流動；而在高Re數時，二次流在顆?？拷ǖ乐行臅r先是向顆粒下方流動，而后隨著顆粒向壁面靠近以及Re數的增大，其逐漸向顆粒左右兩側流動.而受二次流的影響，顆粒所受升力在高Re數和低Re數呈現不同的空間分布規(guī)律.