水平90°彎管內(nèi)固液兩相流動(dòng)的數(shù)值模擬

馬曉陽(yáng)，武傳宇，陳洪立，竇華書

（浙江理工大學(xué)機(jī)械與自動(dòng)控制學(xué)院，杭州310018）

水平90°彎管內(nèi)固液兩相流動(dòng)的數(shù)值模擬

馬曉陽(yáng)，武傳宇，陳洪立，竇華書

（浙江理工大學(xué)機(jī)械與自動(dòng)控制學(xué)院，杭州310018）

為研究90°彎管內(nèi)固液兩相流動(dòng)特征，采用多相流混合模型對(duì)水平90°彎管內(nèi)水和沙粒固液兩相流動(dòng)進(jìn)行數(shù)值模擬，分析彎管典型橫截面上二次流現(xiàn)象，討論其發(fā)展變化對(duì)沙粒濃度分布的影響。模擬結(jié)果顯示：當(dāng)Re=5 ×104時(shí)，隨著入口沙粒濃度升高，彎管出口橫截面中心區(qū)域混合流體速度趨于更均勻分布，隨著入口沙粒直徑增大，沙?？焖俜e聚于管道下側(cè)，形成堆積；當(dāng)Re數(shù)增大到2×105時(shí)，在相同沙粒直徑下，彎管出口橫截面混合流體速度分布變化不大，除管道下側(cè)區(qū)域外，沙粒濃度分布變得更均勻。與實(shí)驗(yàn)結(jié)果對(duì)比表明，該模型可用于彎曲管道內(nèi)固液兩相流動(dòng)特性的有效計(jì)算。

固液兩相流；90°彎管；數(shù)值模擬；二次流；混合模型

0 引 言

固液兩相流是由固相顆粒和液相載流體組成的一種混合流體。在流動(dòng)中，固相顆粒和液相載流體有著密切的聯(lián)系，在運(yùn)動(dòng)中相互影響、相互制約。受管道幾何彎曲的影響，固液兩相流動(dòng)呈現(xiàn)出非常復(fù)雜的流動(dòng)特征［1-3］。固液兩相流管道水力輸送技術(shù)始于20世紀(jì)初，現(xiàn)已廣泛應(yīng)用于能源、化工、礦業(yè)、石油、水利、冶金及環(huán)保等各工業(yè)領(lǐng)域［4］。因此，開展管道內(nèi)固液兩相流動(dòng)特征研究對(duì)于輸送系統(tǒng)的安全運(yùn)行非常重要。

早期固液兩相流的運(yùn)動(dòng)特征研究大都基于宏觀實(shí)驗(yàn)，這些實(shí)驗(yàn)沒有從本質(zhì)上揭示兩相流運(yùn)動(dòng)規(guī)律，同時(shí)對(duì)測(cè)量技術(shù)要求苛刻。近年來(lái)，隨著固液兩相流理論及計(jì)算機(jī)模擬技術(shù)的日趨成熟，采用計(jì)算流體動(dòng)力學(xué)（computational fluid dynamics，CFD）方法研究管道內(nèi)固液兩相流運(yùn)動(dòng)規(guī)律具有成本低和周期短的優(yōu)點(diǎn)［5-6］。Taylor等［7］采用激光多譜勒測(cè)速儀對(duì)90°彎管內(nèi)發(fā)展的層流和湍流進(jìn)行測(cè)量，得到了時(shí)均速度和脈動(dòng)速度分布情況及彎曲段壁面壓力分布等重要物理量的實(shí)驗(yàn)數(shù)據(jù)。Lyne［8］首次強(qiáng)調(diào)彎管中充分發(fā)展的分層振蕩流的復(fù)雜性，提出在穩(wěn)定情況下，管道內(nèi)會(huì)出現(xiàn)一對(duì)渦。這些實(shí)驗(yàn)都是在單相流體條件下進(jìn)行。Skudarnov等［9］通過實(shí)驗(yàn)繪制了水平直管道中單種顆粒和水及雙種顆粒和水的混合流動(dòng)的壓力梯度曲線。胡曉瑋等［10］利用ASM（algebraic slip model）模型模擬了光催化劑反應(yīng)器中催化劑顆粒和液體的混合流動(dòng)，得出了混合流體速度和滑移速度分布情況。張宏兵等［11］運(yùn)用混合模型對(duì)水平直圓管內(nèi)兩相流運(yùn)動(dòng)進(jìn)行了研究，獲得挾沙油混合流體速度大于挾沙水混合流體速度的結(jié)論。目前涉及水平90°彎管內(nèi)固液兩相流動(dòng)特征的研究尚不多見，而且研究尚處于起步階段。

本文運(yùn)用多相流混合模型，模擬了水平90°彎管內(nèi)水和沙粒兩相混合流動(dòng)特征。首先將模擬結(jié)果與實(shí)驗(yàn)結(jié)果進(jìn)行對(duì)比驗(yàn)證，然后運(yùn)用驗(yàn)證了的模型來(lái)求解水和沙粒兩相混合流動(dòng)，得出了彎管出口橫截面混合流體速度和沙粒濃度分布特征。

1 彎管模型和數(shù)值方法

1.1 水平90°彎管幾何模型

計(jì)算模型包括上游入口直管段、彎管段和下游出口直管段三部分，如圖1（a）所示。上游入口直管長(zhǎng)度L1=1 000 mm，下游出口直管長(zhǎng)度L2=1 000 mm，管道內(nèi)徑D=100 mm，彎管的曲率半徑Rc= 300 mm。坐標(biāo)系原點(diǎn)位于彎管的曲率中心，用不同的θ角度代表彎管各個(gè)典型橫截面的位置。利用GAMBIT軟件對(duì)模型計(jì)算域進(jìn)行網(wǎng)格劃分，采用六面體網(wǎng)格，橫截面上網(wǎng)格劃分采用“古錢幣”形式，如圖1（b）所示。

本文以圖1（a）所描述的水平90°三維彎管內(nèi)水和沙粒兩相混合流體流動(dòng)過程為研究對(duì)象。液相為水，密度ρ=998.2 kg/m3；固相為沙粒，密度ρ=2 650 kg/m3，沙粒直徑分0.1 mm和0.2 mm兩種，入口處沙粒初始體積分?jǐn)?shù)分4種，分別為5%、10%、15%和20%。管道z軸負(fù)向重力加速度g=9.81 m/s2。文中Re簡(jiǎn)化為單相水流的雷諾數(shù)，定義為：

式（1）中，ρ為水的密度，ν為水流速度，D為管道內(nèi)徑，μ為水的動(dòng)力黏度。

圖1 計(jì)算模型及橫截面網(wǎng)格劃分

1.2 控制方程和數(shù)值方法

本文采用無(wú)相變過程的多相流混合模型，該模型可用連續(xù)性方程、動(dòng)量方程和第二相體積分?jǐn)?shù)方程表示［12-13］。

1.2.1 連續(xù)性方程

式（2）中，ρm為混合流體密度，為第k相的密度，αk為第k相的體積分?jǐn)?shù)，vm為質(zhì)量平均速度，

1.2.2 動(dòng)量方程

式（3）中，F(xiàn)為體積力，μm為混合流體黏度，μm=為第二相k的漂移速度，vdr，k=vkvm。

1.2.3 沙粒體積分?jǐn)?shù)方程

本文利用FLUENT6.3.26軟件，采用有限體積法離散守恒方程，積分方程離散選用Second Order Upwind，流場(chǎng)數(shù)值計(jì)算選用壓力解法中的SIMPLE算法。迭代收斂標(biāo)準(zhǔn)為各個(gè)變量的相對(duì)殘差小于10-6，且監(jiān)測(cè)的進(jìn)出口流量差值小于0.1%。

為簡(jiǎn)化流動(dòng)條件，假設(shè)水平90°彎管內(nèi)為充分發(fā)展的穩(wěn)態(tài)湍流流動(dòng)，液相水不可壓縮（密度為常數(shù)），流動(dòng)為等溫流動(dòng)。入口為速度入口邊界，給定入口速度的大小和方向；出口為自由出流邊界。管壁為壁面邊界，指定壁面為無(wú)滑移邊界條件，壁面附近采用標(biāo)準(zhǔn)壁面函數(shù)來(lái)處理。湍流模型采用RNG k-ε模型，給定“湍流強(qiáng)度”和“水力直徑”。

1.3 網(wǎng)格無(wú)關(guān)性驗(yàn)證

在水和沙?；旌狭黧w入口平均速度um=0.95 m/s、沙粒濃度α=10%、沙粒直徑d=0.1 mm工況下進(jìn)行網(wǎng)格無(wú)關(guān)性驗(yàn)證。分別作了總數(shù)為23萬(wàn)、31萬(wàn)、35萬(wàn)、43萬(wàn)、48萬(wàn)、53萬(wàn)6種大小網(wǎng)格，并在彎管出口橫截面水平直徑方向上（圖1（b）中y軸方向）設(shè)置9個(gè)監(jiān)測(cè)點(diǎn)，監(jiān)測(cè)其靜壓變化情況。如圖2可知，網(wǎng)格總數(shù)在31萬(wàn)以上時(shí)，各監(jiān)測(cè)點(diǎn)的靜壓變化已經(jīng)很小，誤差在1%以下。圖2結(jié)果表明即使網(wǎng)格總數(shù)繼續(xù)增加也將不會(huì)引起計(jì)算數(shù)據(jù)的顯著變化，而計(jì)算時(shí)間和計(jì)算資源將會(huì)大大的增加，所以計(jì)算采用總數(shù)為35萬(wàn)的網(wǎng)格。

圖2 網(wǎng)格無(wú)關(guān)性驗(yàn)證

1.4 模型驗(yàn)證

為確保多相流混合模型的計(jì)算可靠性，將混合模型用于求解水平直管內(nèi)水和沙粒兩相混合流動(dòng)特征。以一個(gè)水平直管路水和沙粒兩相輸送系統(tǒng)為研究對(duì)象，管道直徑D=0.103 m，管長(zhǎng)L=10 m，水和沙粒混合流體入口平均速度um=3 m/s，入口沙粒體積分?jǐn)?shù)α=24%，沙粒直徑d=90μm。模擬結(jié)果與Gillies等［14］文獻(xiàn)中實(shí)驗(yàn)結(jié)果進(jìn)行對(duì)比。本文混合模型中考慮了拖拽力、重力，忽略了升力、虛擬質(zhì)量力、Magnus力等對(duì)流動(dòng)影響很小的其他力。

圖3給出了直管出口橫截面豎直直徑方向沙粒濃度分布對(duì)比情況。由圖3可見，數(shù)值模擬結(jié)果與實(shí)驗(yàn)結(jié)果基本符合。在管道中間區(qū)域，實(shí)驗(yàn)值和模擬值非常接近，但是模型忽略了管壁粗糙度，所以在靠近彎管上側(cè)和下側(cè)區(qū)域，模擬結(jié)果與實(shí)驗(yàn)結(jié)果出現(xiàn)了偏差。

圖3 數(shù)值模擬結(jié)果與實(shí)驗(yàn)結(jié)果比較

2 數(shù)值模擬結(jié)果和分析

2.1 彎管45°橫截面混合流體速度分布分析

彎管橫截面混合流體速度分布主要受混合流體入口速度以及混合流體沙粒濃度兩個(gè)參數(shù)的影響。如圖4-6所示，與水平直管單相水流橫截面上的速度分布相比，彎管中間橫截面（θ=45°）上混合流體速度表現(xiàn)出不同的分布特征。在Re=5×104時(shí)，水平直管橫截面上水流速度極大值在圓心附近，受管道幾何彎曲的影響，混合流體速度極大值向彎管內(nèi)側(cè)偏移，且隨著沙粒濃度升高，速度極大值有向下偏移趨勢(shì)。當(dāng)Re數(shù)增大到2×105時(shí)，混合流體速度分布主要受管道幾何彎曲的影響，受沙粒濃度影響較小。如圖6所示，混合流體速度基本處于同一分布狀態(tài)。

圖4 水平直管單相水流橫截面速度分布

圖5 彎管中間橫截面（θ=45°）混合流體速度分布（Re=5×104，d=0.1 mm）

圖6 彎管中間橫截面（θ=45°）混合流體速度分布（Re=2×105，d=0.1 mm）

2.2 彎管內(nèi)二次流現(xiàn)象分析

當(dāng)水和沙粒混合流體通過彎管時(shí)，沿彎管軸線方向流動(dòng)的混合流體受到離心力作用，迫使其改變?cè)鹊牧鲃?dòng)方向，被擠壓到彎管外側(cè)。但是由于流動(dòng)的連續(xù)性，其它方位的流體迅速向管中心補(bǔ)充?；旌狭黧w剛進(jìn)入彎管時(shí)，流動(dòng)受管道幾何彎曲影響還不是很大，隨著混合流體深入到彎管以后，混合流體的總壓分布變得不均勻。如圖7所示，總壓高的區(qū)域逐漸向彎管外側(cè)移動(dòng)。

圖7 彎管典型橫截面上混合流體總壓分布（右側(cè)為彎管外側(cè)）

圖8給出了在Re=5×104時(shí)彎管典型橫截面上水和沙?；旌狭黧w流動(dòng)情況。從圖8可以清晰地看到二次流的發(fā)生，這和單相液流通過彎管時(shí)得到的特征是相似的。彎管中總壓高的外側(cè)流體沿著彎管上側(cè)壁面和下側(cè)壁面向總壓低的內(nèi)側(cè)流動(dòng)，同時(shí)管中心流體受離心力作用不斷向外側(cè)流動(dòng)，二次流由此形成。隨著θ角度的增加，彎管典型橫截面上下兩個(gè)渦心位置不斷變化。在θ=0°到θ=60°之間，上渦心不斷向彎管內(nèi)側(cè)移動(dòng)，下渦心逐漸向管中心移動(dòng)。在傳熱傳質(zhì)領(lǐng)域，二次流可以提高換熱系數(shù)；但在管道輸送中，二次流的產(chǎn)生會(huì)造成混合流體的總壓和能量損失。

圖9 彎管典型橫截面上沙粒濃度分布（Re=5×104，α=10%，d=0.1 mm）

圖9給出了在Re=5×104時(shí)彎管典型橫截面上沙粒濃度的分布情況。在彎管入口橫截面（θ= 0°）上沙粒濃度呈不均勻分布，分層現(xiàn)象非常明顯。彎管內(nèi)二次流加強(qiáng)了水和沙粒的混合作用，但是由于此時(shí)流速較低，二次流混合作用效果并不明顯。

2.3 彎管出口橫截面混合流體速度和沙粒濃度分布

圖10給出了不同入口速度和沙粒濃度條件下，彎管出口橫截面（θ=90°）混合流體速度和沙粒濃度分布情況。圖10中y/D表示彎管出口橫截面水平直徑方向無(wú)量綱距離，z/D表示豎直直徑方向無(wú)量綱距離，見圖1（b）所示。圖10（a）和（b）顯示出在彎管出口橫截面水平直徑方向上混合流體速度分布不再關(guān)于管軸線對(duì)稱分布。在較小入口速度um=0.5 m/s（Re=5×104）時(shí)，靠近彎管內(nèi)側(cè)的速度隨著入口沙粒濃度的升高而增大，此時(shí)出現(xiàn)了類似“塞狀流”特征，即管中心某個(gè)半徑圓柱內(nèi)流速保持不變。當(dāng)入口速度增大到um=2 m/s（Re=2×105）時(shí)，混合流體速度分布隨沙粒濃度升高變化不明顯，基本處于同一分布狀態(tài)。只在靠近彎管內(nèi)側(cè)區(qū)域，混合流體速度隨著入口沙粒濃度升高有略微增大。

圖10（c）和（d）給出了彎管出口橫截面豎直直徑方向上沙粒濃度分布情況。在較小入口速度um=0.5 m/s（Re=5×104）時(shí)，彎管出口橫截面沙粒濃度分布大致可以分為3個(gè)區(qū)域：管道上側(cè)附近的低濃度區(qū)域，管道中心的均勻分布區(qū)域和管道下側(cè)附近的高濃度區(qū)域。由于固液兩相存在密度差，受重力的影響，因而在管道上側(cè)區(qū)域濃度低，下側(cè)區(qū)域濃度高。隨著入口速度增大到um=2 m/s（Re=2× 105），除靠近管道下側(cè)區(qū)域外，其它區(qū)域沙粒濃度趨于均勻分布。此時(shí)管道內(nèi)水流速度足以托舉沙粒，使其懸浮于管道中，使得沙粒能夠很好地跟隨載流體水流動(dòng)。

圖10 彎管出口橫截面水平直徑方向混合流體速度分布和豎直直徑方向沙粒濃度分布

由圖10（c）和（g）可知，在um=0.5 m/s（Re=5 ×104）時(shí)，隨著沙粒直徑增大到d=0.2 mm，彎管出口橫截面豎直直徑方向上沙粒濃度分布變得越來(lái)越不均勻。從管中心到管上側(cè)某些區(qū)域出現(xiàn)無(wú)沙粒區(qū)，而靠近管道下側(cè)區(qū)域沙粒濃度則迅速增大。這種輸送條件很可能引起沙粒堆積，形成不能移動(dòng)的沙粒層床，進(jìn)而減小有效流動(dòng)面積，引起管道堵塞，對(duì)沙粒輸送非常不利。如圖10（h）可知，當(dāng)入口流速增大到um=2 m/s（Re=2×105）時(shí)，上述情況消失，管中心區(qū)域沙粒濃度又趨于相對(duì)均勻分布。

同時(shí)由圖10（h）可以發(fā)現(xiàn)在離管底附近區(qū)域（z/D=0.15）的沙粒濃度突然發(fā)生了變化，而在圖10（d）中則沒有出現(xiàn)。圖10（h）和（d）的對(duì)比說明在沙粒較小直徑和較低濃度時(shí)，該區(qū)域沙粒濃度突變不會(huì)發(fā)生，這可能是由于該區(qū)域急劇變化的速度梯度和較高的沙粒濃度引起了沙粒間的碰撞而造成，這個(gè)解釋的正確性需要后續(xù)工作進(jìn)一步驗(yàn)證。

3 結(jié) 論

本文采用無(wú)相變過程的多相流混合模型，研究了水平90°彎管內(nèi)水和沙粒兩相混合流動(dòng)特征。數(shù)值模擬結(jié)果給出了水和沙?；旌狭鲃?dòng)過程中一些重要的流動(dòng)參數(shù)，結(jié)果表明：混合流體速度在90°彎管典型橫截面上的分布，與單相水流明顯不同。在彎管θ=45°橫截面上，速度極大值向彎管內(nèi)側(cè)偏移，且在Re較小時(shí)，如Re=5×104，混合流體速度的極大值隨著沙粒濃度的升高，呈下移趨勢(shì)。受管道幾何彎曲的影響，所受離心力大的流體向彎管外側(cè)運(yùn)動(dòng)，導(dǎo)致外側(cè)總壓升高，使得混合流體沿著管道上下兩側(cè)從彎管外側(cè)向總壓低的內(nèi)側(cè)流動(dòng)，呈現(xiàn)出類似于單相流體經(jīng)過彎管時(shí)出現(xiàn)的二次流現(xiàn)象。在Re=5 ×104和d=0.1 mm時(shí)，彎管出口中心區(qū)域混合流體速度在水平直徑方向上分布隨著沙粒濃度的升高而趨于更均勻分布；豎直直徑方向上沙粒濃度分布呈明顯的分層現(xiàn)象。隨著沙粒直徑增大到d=0.2 mm，沙粒濃度分布變的極不均勻。當(dāng)Re增大到2 ×105時(shí)，混合流體速度分布受沙粒濃度影響變小，沙粒濃度在彎管出口橫截面中心區(qū)域分布更均勻。

研究結(jié)果可為管道水力輸送系統(tǒng)局部設(shè)計(jì)提供參考，也可以對(duì)水力輸送過程中固相顆粒在管道中的堆積程度進(jìn)行預(yù)測(cè)。實(shí)際工程中很多管道都是運(yùn)動(dòng)著的，后續(xù)工作將對(duì)旋轉(zhuǎn)彎管內(nèi)固液兩相流動(dòng)特征進(jìn)行研究。

［1］倪晉仁，王光謙.固液兩相流研究的兩種方法之比較［J］.泥沙研究，1992（3）：95-102.

［2］Ling J，Skudarnov P V，Lin C X，et al.Numerical investigations of liquid-solid slurry flows in a fully developed turbulent flow region［J］.International Journal of Heat and Fluid Flow，2003，24（3）：389-398.

［3］Lin C X，Ebadian M A.A numerical study of developing slurry flow in the entrance region of a horizontal pipe［J］.Computers&Fluids，2008，37（8）：965-974.

［4］白曉寧，胡壽根.漿體管道的阻力特性及其影響因素分析［J］.流體機(jī)械，2000，28（11）：26-29.

［5］宋文吉，肖 睿，馮自平，等.潛熱輸送介質(zhì)顆粒沉降速度的固-液兩相流模擬［J］.工程熱物理學(xué)報(bào)，2010，31（10）：1693-1696.

［6］Chen L，Duan Y，Pu W，et al.CFDsimulation of coal-water slurry flowing in horizontal pipelines［J］.Korean Journal of Chemical Engineering，2009，26（4）：1144-1154.

［7］Taylor A，Whitelaw J H，Yianneskis M.Curved ducts with strong secondary motion-velocity measurements of developing laminar and turbulent flow［J］.ASME Transactions Journal of Fluids Engineering，1982，104：350-359.

［8］Lyne W H.Unsteady viscous flow in a curved pipe［J］. Journal Fluid Mech，1971，45：13-31.

［9］Skudarnov P V，Kang H J，Lin C X，et al.Experimental investigation of single-and double-species slurry transport in a horizontal pipeline［C］//The Proceedings of ANS 9th International Topical Meetingon Robotics and Remote Systems.Seattle，WA，2001：4-9.

［10］HU X，GUO L.Numerical investigation of catalystliquid slurry flow in the photocatalytic reactor for hydrogen production based on algebraic slip model［J］. International Journal of Hydrogen Energy，2010，35（13）：7065-7072.

［11］張宏兵，陳露露，謝榮華，等.水平圓管固液兩相穩(wěn)態(tài)流動(dòng)特性數(shù)值模擬［J］.化工學(xué)報(bào)，2009，60（5）：1162-1168.

［12］Wang J，Wang S，Zhang T，et al.Numerical investigation of ice slurry isothermal flow in various pipes［J］. International Journal of Refrigeration，2012，36：70-80.

［13］Raisee M，Alemi H，Iacovides H.Prediction of developing turbulent flow in 90°-curved ducts using linear and non-linear low-Re k-εmodels［J］.International Journal for Numerical Methods in Fluids，2006，51（12）：1379-1405.

［14］Gillies R G，Shook C A，Xu J.Modelling heterogeneous slurry flows at high velocities［J］.The Canadian Journal of Chemical Engineering，2004，82（5）：1060-1065.

Numerical Simulation of Solid-liquid Two-phase Flow in a Horizontal 90°Elbow Pipe

MA Xiao-yang，WU Chuan-yu，CHEN Hong-li，DOU Hua-shu
（School of Mechanical Engineering&Automation，Zhejiang Sci-Tech University，Hangzhou 310018，China）

To study solid-liquid two-phase flow feature in 90°elbow pipe，this paper uses multiphase flow mixture model for numerical modeling of water and sand solid-liquid two-phase flow in horizontal 90° elbow pipe，analyzes secondary flow phenomenon on typical cross section of elbow pipe and discusses the influence of its development and change on concentration distribution of sand.The simulation result shows that，when Re=5×104，with the increase of sand concentration at the entrance，the speed of mixed fluid in the central area of cross section of elbow pipe at the exit tends to more uniform distribution；with the increase of sand diameter at the entrance，sand accumulates in the underside of pipe rapidly and forms stacking；when Re increases to 2×105，under the same sand diameter，the distribution of speed of mixed fluid on the cross section of elbow pipe at the entrance does not change a lot；except the underside area of pipe，sand concentration distribution becomes more uniform.The comparison with the experimental result shows that this model can be used for effective calculation of solid-liquid two-phase flow characteristics in bent pipe.

solid-liquid two-phase flow；90°elbow pipe；numerical simulation；secondary flow；mixture model

TH31

A

（責(zé)任編輯：康 鋒）

1673-3851（2014）03-0228-07

2013-12-03

國(guó)家自然科學(xué)基金（51205363）；浙江理工大學(xué)科研啟動(dòng)基金（11130032241201）

馬曉陽(yáng)（1990-），男，浙江東陽(yáng)人，碩士研究生，研究方向?yàn)閮上嗔骱土黧w機(jī)械。

竇華書，電子郵箱：huashudou@yahoo.com