姚雪蕾，袁成清，付宜風(fēng)，白秀琴

(武漢理工大學(xué)a.能源與動力工程學(xué)院可靠性工程研究所；b.船舶動力工程技術(shù)交通行業(yè)重點實驗室,武漢 430063)

管道內(nèi)壁粗糙度對沿程阻力影響的FLUENT數(shù)值模擬分析

姚雪蕾，袁成清，付宜風(fēng)，白秀琴

(武漢理工大學(xué)a.能源與動力工程學(xué)院可靠性工程研究所；b.船舶動力工程技術(shù)交通行業(yè)重點實驗室,武漢 430063)

針對在海底管道的設(shè)計階段準確獲取沿程阻力的問題，采用FLUENT軟件運用數(shù)值模擬的方法，分析粘性流體在管道中的摩阻損失與管道內(nèi)壁粗糙度的關(guān)系，將計算結(jié)果與理論結(jié)果進行對比。結(jié)果顯示，除流態(tài)從水力光滑區(qū)到混合摩擦區(qū)間有一個過渡階段存在誤差外，其余都吻合良好，證實輸油管道的內(nèi)壁粗糙度只在混合摩擦區(qū)對摩阻有較大影響，改善管壁面粗糙度可以減小阻力。若流態(tài)處于水力光滑區(qū)到混合摩擦區(qū)之間，計算摩阻系數(shù)僅靠經(jīng)驗公式是不合適的，應(yīng)考慮采用數(shù)值模擬計算的方法。

輸油管道；粗糙度；沿程阻力；FLUENT；數(shù)值模擬

隨著計算流體動力學(xué)的發(fā)展，利用FLUENT分析粘性流體在管道中的阻力特性已經(jīng)得到了很多實用性的成果，如針對粗糙的管壁[1-4]或者是產(chǎn)生局部阻力損失的典型閥件(突擴管、三通管、彎管等)[5-6]，通過數(shù)值模擬的方法可以模擬出常規(guī)實驗無法得到的不同流場的各種信息。然而，目前工程上在管道的設(shè)計階段，沿程阻力的確定一般只采用經(jīng)驗公式數(shù)學(xué)計算的方法，故考慮運用數(shù)值模擬的方法對其準確性進行驗證，為管道設(shè)計中沿程阻力的獲取提供依據(jù)。

1 沿程阻力與管道內(nèi)壁粗糙度

管道的沿程阻力損失與管長、管流的速度分布，以及管壁的粗糙度等因素有關(guān)，通常用hf表示，按達西(Darcy-Weisbach)公式計算。

(1)

式中:L——管道長度;D——管道內(nèi)徑;v——流體流速;λ——摩阻系數(shù)。

相對當量粗糙度ε=2e/D。式中e為管壁的絕對當量粗糙度，是指管內(nèi)壁凸起高度的統(tǒng)計平均值。研究表明，新的、清潔的管壁的絕對當量粗糙度僅取決于管材及制管方法，與管徑無關(guān)；使用后的管路則隨運行情況，如所輸流體性質(zhì)、腐蝕程度、運行年限、清管方法等不同而有顯著變化[7]。

根據(jù)莫迪理論，一般可以把圖分為5個區(qū)域：層流區(qū)、層流轉(zhuǎn)變?yōu)槲闪鞯倪^渡區(qū)、紊流水力光滑區(qū)、紊流混合摩擦區(qū)、紊流粗糙區(qū)[8-9]。在實際運行中，熱原油管道上最常見的流態(tài)是水力光滑區(qū)，輕油管道也多在水力光滑區(qū)，輸送低粘油品的較小直徑管道可能進入混合摩擦區(qū)。所以本研究重點在紊流水力光滑區(qū)和混合摩擦區(qū)，摩阻系數(shù)λ與相對粗糙度ε的關(guān)系。

2 經(jīng)驗公式分析

從莫迪圖可以定性看出，水力光滑區(qū)中相對粗糙度ε變化，摩阻系數(shù)λ不變；而混合摩擦區(qū)的摩阻系數(shù)λ會隨相對粗糙度ε的減小而減少。然而，不能明確地觀察出λ隨ε定量變化的形態(tài)和趨勢。因此，以工程上常應(yīng)用的米勒公式(式2)和Colebrook-White公式(式3)為基礎(chǔ)(見表1)，以相對粗糙度ε為變量，雷諾數(shù)Re為參量進行數(shù)學(xué)計算，研究粗糙度對輸油管道摩阻系數(shù)的影響，得出摩阻系數(shù)隨相對粗糙度變化曲線[10]。

計算結(jié)果見圖1。其中由于公式(2)為λ的隱函數(shù)，故利用埃特金法[11]進行迭代求得近似解。總體上可發(fā)現(xiàn)摩阻系數(shù)隨相對粗糙度的減小而減小，與莫迪圖結(jié)果一致。并且在相對粗糙度較大的區(qū)域，摩阻系數(shù)減小的更明顯，而在相對粗糙度較小的區(qū)域，摩阻系數(shù)減小幅度較小，并趨于穩(wěn)定。高雷諾數(shù)管內(nèi)流體要在相對粗糙度較小下才處于水力光滑區(qū)，相比低雷諾數(shù)下的流態(tài)更容易處于水力光滑區(qū)，穩(wěn)定區(qū)域更大。但水力光滑區(qū)的摩阻系數(shù)隨雷諾數(shù)減小而增大。

表1 λ的經(jīng)驗公式

圖1 摩阻系數(shù)隨相對粗糙度變化

在上述計算過程中，判斷流態(tài)為紊流水力光滑區(qū)還是混合摩擦區(qū)是以臨界雷諾數(shù)Re1為標志分界線，其值由管道層流邊界厚度與絕對當量粗糙度的比值推導(dǎo)而來。

(4)

對于某一直徑的管子，其臨界雷諾數(shù)的數(shù)值取決于管壁絕對當量粗糙度。故當管內(nèi)油流的雷諾數(shù)接近臨界值時，合理確定粗糙度數(shù)值，成為準確判別流態(tài)區(qū)和計算λ的關(guān)鍵因素。從另一方面來看，式(4)中Re1值是根據(jù)勃拉休斯(Blasius)公式等經(jīng)驗公式得來的，若采用不同的公式推導(dǎo)，將會得出不同的計算式。因此，根據(jù)臨界值來劃分流態(tài)區(qū)只是近似的，實際管道中流態(tài)從水力光滑區(qū)到混合摩擦區(qū)變化有一個過渡階段，使用經(jīng)驗公式對過渡階段的λ值計算不夠準確，需查詢莫迪圖或進行試驗來修正。

雷諾數(shù)為104～105情況下，分別查莫迪圖[9]及用公式計算得到摩阻系數(shù)λ，繪制臨近過渡階段的摩阻系數(shù)隨相對粗糙度變化曲線見圖2。

圖2 過渡階段摩阻系數(shù)隨相對粗糙度變化

由圖2可見，若根據(jù)理論公式計算，粗糙度較小時，流體處于在水力光滑區(qū)摩阻系數(shù)不隨粗糙度變化而變化，曲線是一條直線；隨著粗糙度增大，當進入混合摩擦區(qū)時，λ才開始隨相對粗糙度增大而增大。而實際上摩阻系數(shù)λ并不是在臨界雷諾數(shù)時突變，而是在Re接近Re1前就開始逐漸增大。同時還可以發(fā)現(xiàn)，在這一過渡階段之外，公式計算得出的λ值與莫迪圖查詢的值都很相近，也證明了米勒公式和Colebrook-White公式在多數(shù)工業(yè)管道流體中是適用的。

3 基于Fluent的數(shù)值模擬

研究流體流動問題的完整體系由傳統(tǒng)的理論分析方法、實驗測得方法及CFD數(shù)值方法組成，各有優(yōu)缺點，但又有很多聯(lián)系，在實際工程中將三者進行有機的結(jié)合能得到更合理的結(jié)果。用經(jīng)驗公式對管道內(nèi)流體進行摩阻計算實際上還不夠完善，特別是對于黏性流體在管中的流動，某些情況下得到的λ值并不準確，還需要通過實驗進行模擬和驗證。但實驗常受到諸多因素的限制，且經(jīng)費投入較大，相形之下，CFD更有優(yōu)勢。分析流體流動的基本方程主要是質(zhì)量守恒、動量守恒、能量守恒三種方程。如果流動涉及不同組分的混合或是相互作用，還要遵守組分守恒定律。在實際計算時，還要考慮不同的流態(tài)，對于湍流，目前工程中一般是對瞬態(tài)N-S方程做時間平均處理，同時補充反映湍流特性的其他方程，如湍動能方程和湍流耗散率方程等[12]。

本研究中基于FLUENT軟件進行三維數(shù)值模擬。模擬對象為一段直管內(nèi)油品的湍流流動。通過改變?nèi)肟诙肆魉賤來控制雷諾數(shù)Re，并分別在不同管壁粗糙度下進行數(shù)值模擬，得出處于管中試驗段某兩端面的壓力的差值，根據(jù)伯努利方程和達西公式聯(lián)立計算即可得到沿程損失阻力系數(shù)λ，再將所得的值與莫迪圖曲線或經(jīng)驗公式計算得出的數(shù)值對比，判斷其是否正確。

3.1 物理模型

在GAMBIT中建立直徑300 mm，長4 m的圓截面直管，其中前3 m是前置段，用來讓湍流充分發(fā)展，后1 m為試驗段。管內(nèi)流體選取我國中原地區(qū)原油，基本參數(shù)：50 ℃密度846.6 kg/m3，運動粘度10.32 mm2/s。

3.2 劃分網(wǎng)格

管內(nèi)流體的阻力主要是流體與壁面間的作用產(chǎn)生的，在壁面附近會有較大的速度梯度。故在端面上創(chuàng)建邊界層網(wǎng)格加密。具體設(shè)置為：first per centage(第一層邊界層網(wǎng)格的高度關(guān)于寬度的百分比)為20,Growth factor取1.2，rows取5層。端面網(wǎng)格劃分采用默認方式，單元格為四邊形，網(wǎng)格個數(shù)count設(shè)為50。得到mesh faces為473，體網(wǎng)格mesh volumes為100 276。網(wǎng)格劃分情況見圖3，圖中坐標X方向為流向，Y、Z方向為徑向。

圖3 模型網(wǎng)格劃分

后續(xù)模擬計算中，經(jīng)過測試，在此網(wǎng)格密度上即使再提高其密度，對計算的精確性也無法進一步明顯的提高，而計算量反而大大增加。因此選用此密度進行網(wǎng)格劃分可達到計算要求。

3.3 邊界條件與初始化

將在GAMBIT中劃分好的網(wǎng)格文件導(dǎo)入FLUENT中。采用標準的k-ε湍流模型與標準壁面函數(shù)法(Standard wall function)，變化參數(shù)主要有兩個:①入口速度及相關(guān)的湍流參數(shù);②壁面粗糙度。具體設(shè)置如下。

1)模型求解器(Solver)選擇壓力基、顯式、定常流。

2)粘性模型(viscous model)選擇雙方程模型(k-epsilon[2 eqn]Standard)。

3)流體材料設(shè)置為原油。

4)速度進口(velocity-inlet)邊界條件設(shè)置為入口端速度，變化范圍為2～6 m/s。湍流方法指定為湍流強度I與水力直徑DH,充分發(fā)展的管流核心湍流強度I可按經(jīng)驗公式來估算。

(5)

式中:ReDH——以水力直徑為特征長度求出，水力直徑定義為流通面積與濕周比值的4倍。

5)出流(outflow)邊界條件保持默認值。

6)固壁(wall)邊界條件設(shè)置壁面粗糙度(wall roughness)，在此默認粗糙度厚度值KS(roughness height)等同于絕對當量粗糙度e。變化范圍為0.05～0.30 mm；粗糙度常數(shù)CS(roughness constant)保持默認值0.5，認為管壁粗糙度均勻分布。

7)其他條件如求解方法(Solutions)等默認FLUENT的初始設(shè)定，使用進口的量對全場完成初始化。

3.4 迭代計算

計算何時收斂對結(jié)果的準確性有重要影響，本研究中，模型的收斂采用兩個判斷標準，一是監(jiān)視每個守恒值的殘差，X、Y及Z方向速度的殘差分別要求低于1×10-7，連續(xù)性的殘差低于1×10-2，湍動能k和耗散率ε的殘差低于1×10-6。另一個是出口流量值，只有當出口流量值穩(wěn)定，且入口流量與出口流量近似相等時，才可以判定為收斂。完成以上所有設(shè)定后，開始迭代計算。

4 FLUENT計算結(jié)果

4.1 試算情況

以入口速度設(shè)置為2 m/s，粗糙度為0.05 mm的邊界條件為例，使用標準的k-ε湍流模式對該模型進行試算。迭代停止后各方向速度及連續(xù)性殘差監(jiān)測見圖4。

圖4 殘差監(jiān)控

由圖4可見，其收斂條件已經(jīng)符合上述要求。另外，質(zhì)量流量通量報告表明，出口與入口的質(zhì)量流量誤差為4.577 636 7×10-5，也表明計算收斂，應(yīng)用FLUENT軟件對本例的計算是可行的。

管道出口速度在徑向上的變化見圖5，軸線壓降見圖6。

圖5 出口速度徑向分布

圖6 軸線壓降

由圖5、6可見，管道后端速度呈充分發(fā)展狀態(tài)，壓降呈線性?？芍四Ｐ陀嬎阌蜻x取足夠長，可應(yīng)用FLUENT進行模擬計算。

4.2 計算結(jié)果分析

通過FLUENT的Report菜單中的surface Integrals命令可以獲得圓管任一端面的面積加權(quán)平均壓強值。模擬選擇x=3.4 m和x=3.9 m兩端面的計算值，分別記為pa1和pa2，通過其壓強差得到沿程阻力損失進而算得摩阻系數(shù)λ。

(6)

式中:D——管徑0.3 m;L——距離0.5 m;ρ——密度,846.6 kg/m3;v——流速,m/s。

為了驗證FLUENT仿真的準確度，將模擬計算值λ1與通過按表1中經(jīng)驗公式計算得出的數(shù)值λ2進行對比。結(jié)果對比見到表2。

由表2可見，若流態(tài)為水力光滑區(qū)，基于FLUENT仿真的模擬計算結(jié)果與經(jīng)驗公式計算值基本接近，最大誤差為2.33%。在此狀態(tài)下對三維管道的數(shù)值模擬取得了合理的結(jié)果，并具有較高的計算精度。而流態(tài)為混合摩擦區(qū)時，基于FLUENT的模擬計算結(jié)果與經(jīng)驗公式計算值誤差較大，最大誤差絕對值達15.65%，數(shù)值模擬沒有達到較好的效果，需要分析誤差原因。

1)流態(tài)為水力光滑區(qū)時，相對粗糙度ε變化，摩阻系數(shù)λ幾乎不變。模擬計算結(jié)果符合理論，證明了摩阻系數(shù)λ只和雷諾數(shù)有關(guān)，與相對粗糙度ε無關(guān)。

2)當雷諾數(shù)Re很接近臨近雷諾數(shù)Re1的值時，模擬計算的摩阻系數(shù)λ開始逐漸增大，但不顯著。例如流速為3 m/s，粗糙度為0.25 mm時λ值較之前不隨粗糙度變化的λ的值增大率為1.8%，與公式計算值的誤差為2.23%，此時Re比Re1只近似小于1 800。分析認為FLUENT的數(shù)值模擬結(jié)果可以證明，出現(xiàn)此變化的原因應(yīng)該是管道中的流態(tài)變化有一個過渡階段，實際上摩阻系數(shù)λ并不是在臨界雷諾數(shù)時突變，而是在Re接近Re1前的過渡階段就開始逐漸增大。

圖7 數(shù)值模擬計算的λ隨ε變化

3)流態(tài)為混合摩擦區(qū)時，數(shù)值模擬的值與公式計算值有較大的誤差，使用標準k-ε湍流模型進行數(shù)值模擬的值小于經(jīng)驗公式計算值。但變化趨勢相同，如圖7，和前面的分析一致，符合莫迪圖。證明λ與Re和ε均有關(guān)，在雷諾數(shù)相同時，λ隨ε增大而增大。

4.3 誤差原因

1)經(jīng)驗公式使用的是科爾布魯克-懷特(Colebrook-White)公式，通過試算迭代得到的近似解是不精確的。

2)由于湍流模式和數(shù)值計算中都有不可避免的誤差，同時模擬的管道中流體可能處于水力光滑區(qū)至混合摩擦區(qū)間的過渡階段，不完全適用于科爾布魯克-懷特公式或數(shù)值模擬。從數(shù)據(jù)上來看，隨著雷諾數(shù)增大，λ誤差整體上從15%下降為14%。由此可以推測，流體雷諾數(shù)比臨界雷諾數(shù)Re1大得越多，較完全地處于混合摩擦區(qū)時，模擬的計算值越有符合經(jīng)驗公式值的趨勢。

表2 模擬計算與公式計算結(jié)果對比

注：流態(tài)Ⅱ表示為水力光滑區(qū)，Ⅲ表示為混合摩擦區(qū)。

3)標準k-ε湍流模型在壁面區(qū)使用了不夠精確的近壁函數(shù)的半經(jīng)驗公式，需要修正模型或精細網(wǎng)格。

為了驗證猜想2)，選取幾組參數(shù)進行數(shù)值模擬，讓Re值適中地處于Re1和Re2值之間。其中Re2為判斷混合摩擦區(qū)和粗糙區(qū)的臨界雷諾數(shù)，計算式見表1。數(shù)據(jù)處理同前。結(jié)果見表3。

結(jié)果表明，數(shù)值模擬得到的λ與經(jīng)驗公式值誤差很小，符合猜想2)。通過FLUENT的數(shù)值模擬發(fā)現(xiàn)，水力光滑區(qū)至混合摩擦區(qū)間的過渡階段(Re值接近于Re1的值)的流態(tài)，不完全適用于科爾布魯克-懷特公式，或是使用標準k-ε湍流模型進行數(shù)值模擬。當流態(tài)較完全地處于混合摩擦區(qū)時的結(jié)果更吻合，能適用于工程。

為了驗證猜想3)，對模型劃分了更精細的網(wǎng)格。具體設(shè)置為：first percentage取18,Growth factor取1.1，rows取6層。端面網(wǎng)格劃分采用默認方式，單元格為四邊形，網(wǎng)格個數(shù)count設(shè)為5。得到mesh faces為616，體網(wǎng)格mesh volumes為146 608。使用標準k-ε湍流模式中非平衡壁面函數(shù)法(Non-Equilibrium Wall Function)。以流速為6 m/s，粗糙度為0.3 mm為例，結(jié)果見表4。

表3 計算結(jié)果

表4 計算結(jié)果比較

結(jié)果表明，計算精確性沒有明顯提高，而計算量卻大大增加。由此判斷猜想3)不是產(chǎn)生誤差的最主要原因。使用更精細的模型要付出更大的模擬計算量和時間，不適用于工程實際應(yīng)用。

5 結(jié)論

1)數(shù)值模擬證實了輸油管道的內(nèi)壁粗糙度只在流態(tài)為混合摩擦區(qū)對摩阻有較大影響，即摩阻系數(shù)λ隨相對粗糙度ε的增大而增大。且摩阻系數(shù)的變化幅度與粗糙度成正比，在相對粗糙度較大的區(qū)域，摩阻系數(shù)增大的更明顯，在相對粗糙度較小的區(qū)域，摩阻系數(shù)增大的幅度較?。?/p>

2)經(jīng)驗公式推知流態(tài)區(qū)的劃分只是近似的，當流態(tài)處于變化臨界時摩阻系數(shù)用經(jīng)驗公式計算的準確性有待驗證。仿真分析表明，管道中的流態(tài)從水力光滑區(qū)到混合摩擦區(qū)間有一個過渡階段，在這一階段摩阻系數(shù)λ隨與相對粗糙度ε增大而小幅度增大。因此，若想獲取準確的沿程阻力，在此階段計算摩阻系數(shù)僅靠經(jīng)驗公式是不合適的，數(shù)值模擬計算是應(yīng)考慮的方法；

3)在過渡階段之外，流態(tài)可以認為是完全處于水力光滑區(qū)或混合摩擦區(qū)。在此種情況下，經(jīng)驗公式計算、莫地圖查詢與數(shù)值模擬得出的λ值都很相近，證明了米勒公式和Colebrook-White公式在多數(shù)工業(yè)管道流體沿程阻力計算中具有一定的適用性和準確性。

[1] VIJIAPURAPU，S， JIE C. Numerical Simulation of Turbulent Flow in a Ribbed Pipe Using Large Eddy Simulation[J]. Numerical Heat Transfer, Part A: Application(ISSN:1040-7782), 2007,51(12):1137-1165.

[2] VIJIAPURAPU，S，.& JIE C. Performance of turbulence models for flows through rough pipes[J], Applied Mathematical Modelling(ISSN: 0307-904X), 2010,34(6):1458-1466.

[3] ZHANG,Y B & Che,D F. Numerical simulation of turbulent flow and wall mass transfer in a rectangular channel roughened by v-shaped grooves[J],Numerical Heat Transfer,Part A:Application(ISSN:1040-7782),2014,66(5):551-581.

[4] 張建卓, 李康康, 高 猛,等. 粗糙度對流體泄漏和摩擦阻力的仿真研究[J].微計算機信息,2011,27(6):238-240.

[5] 孫 琳.圓管突擴(縮)局部阻力試驗及流場數(shù)值模擬[D].哈爾濱:東北農(nóng)業(yè)大學(xué),2010.

[6] 趙 月.基于CFD的管道局部阻力的數(shù)值模擬[D].大慶:東北石油大學(xué),2011.

[7] 楊筱蘅.輸油管道設(shè)計與管理[M].北京:中國石油大學(xué)出版社,2001.

[8] 楊樹人,汪志明,何光渝,等.工程流體力學(xué)[M].北京:石油工業(yè)出版,2006.

[9] 詹德新,王家楣.工程流體力學(xué)[M].武漢:湖北科學(xué)技術(shù)出版社,2001.

[10] 楊曉鴻,朱薇玲.粗糙度對輸氣管道摩阻系數(shù)的影響[J].石油化工設(shè)備,2005,34(1):25-28.

[11] 羅祖帖,Excel在管路阻力計算中的應(yīng)用[J].化肥設(shè)計,2004,42(2):15-17.

[12] 王福軍.計算流體動力學(xué)分析:CFD軟件原理與應(yīng)用[M].北京：清華大學(xué)出版社,2004.

Numerical Simulation of Effect of Wall Roughness on the On-way Resistance Based on FLUENT

YAO Xue-lei, YUAN Cheng-qing, FU Yi-feng, BAI Xiu-qin

(a. Reliability Engineering Institute, School of Energy and Power Engineering; b. Key Laboratory of Marine Power Engineering and Technology (Ministry of Transport), Wuhan University of Technology, Wuhan 430063, China)

Aiming at the problem of obtaining accurately the on-way resistance of submarine pipeline at the design stage, numerical simulation method is adopted to analyze the relationship between loss of viscous fluid friction and wall roughness of the pipeline based on FLUENT. Comparing the calculated results with theoretical results and the actual situations, the results show there is a good agreement, except for the flow pattern transiting from the hydraulic smooth region to the mixed friction region. It is confirmed that the friction coefficient is strongly influenced by pipeline wall roughness only in the turbulent mixing friction region. In this case, reducing the surface roughness of the wall can cut down the resistance. If the flow pattern transits from the hydraulic smooth region to the mixed friction region, it is inappropriate to calculate friction coefficient only by the empirical formulas, and numerical simulation is the method that should be considered.

oil pipeline; roughness; on-way resistance; FLUENT; numerical simulation

10.3963/j.issn.1671-7953.2015.06.024

2015-07-17

國家自然科學(xué)基金優(yōu)秀青年基金項目 (51422507)

姚雪蕾(1992-)，女，碩士生

TE832；P756.2

A

1671-7953(2015)06-0101-07

修回日期：2015-09-10

研究方向:摩擦學(xué)系統(tǒng)及表面工程

E-mail: yaoxuelei2013@outlook.com