黃 坤，鄧 平，李岳鵬，廖 檸，赫文博

(1.西南石油大學 石油與天然氣工程學院，四川 成都 610500； 2.舟山市港航和口岸管理局，浙江 舟山 316000； 3.重慶華潤凱源燃氣有限公司，重慶 400000； 4.浙江浙能溫州液化天然氣有限公司，浙江 溫州 325000)

0 引言

由于過濾分離設備的能力有限，天然氣在管道輸送過程中會攜帶有少量固體顆粒[1]，這些顆粒一部分來源于氣田,另一部分則是由已有固體顆粒與管壁碰撞產(chǎn)生，微粒進入到管道內(nèi)，與管輸天然氣一起形成氣固兩相流?？紤]到顆粒的運動具有不規(guī)則性，顆粒必定會不斷沖擊管壁，造成管道的沖蝕磨損[2]，尤其是彎管的磨損程度要比直管段嚴重得多[3-6]，長時間的沖蝕磨損可能會導致管道的泄漏事故[7]。因此，如何降低彎管的沖蝕磨損程度，提高彎管的耐磨性能，保證天然氣輸送管道彎管的安全運行是一個十分重要的工程課題。

近年來，國內(nèi)外學者針對如何提高管道壁面的耐磨性能進行了大量研究。Song等[8]研究了在直管內(nèi)壁上添加肋條后管壁的磨損程度；Fan等[9]在直角彎管彎曲段外側(cè)壁上添加了具有一定截面幾何形狀的肋條，并運用數(shù)值模擬技術(shù)計算了彎管壁面的沖蝕磨損速率；Pouraria等[10]針對海底管道，采用數(shù)值模擬的方法研究了T型管和標準彎頭的沖蝕磨損程度；戚勝等[11]提出了一種彎頭的替代結(jié)構(gòu)，并對這種新型裝置的流場特性和沖蝕磨損規(guī)律進行了數(shù)值模擬分析；王宇等[12]運用數(shù)值模擬技術(shù)研究了兩種不同組合彎頭內(nèi)的氣-固兩相流動特征和管壁磨損特性，并根據(jù)數(shù)值模擬結(jié)果對組合彎頭的結(jié)構(gòu)進行了優(yōu)化；Carlos 等[13]通過研究發(fā)現(xiàn)在彎頭處添加渦流腔可有效緩解彎頭的沖蝕磨損；季楚凌等[14]根據(jù)仿生學原理，對彎管內(nèi)壁面進行了改進，并運用數(shù)值模擬技術(shù)研究了改進后彎管的沖蝕磨損情況。

以上研究大多是通過改變彎管的內(nèi)部結(jié)構(gòu)或者彎管的壁面從而優(yōu)化彎管的耐磨性能，對于彎管管形優(yōu)化的研究還比較缺乏，本文針對天然氣輸送管道氣固兩相流動的90°彎管，對其彎曲段進行三段式改進，利用COMSOL軟件對三段式彎管彎曲段所對應的4個管形參數(shù)進行最優(yōu)化處理，使用Fluent軟件對優(yōu)化后的三段式彎管進行沖蝕數(shù)值模擬并與一段式彎管的模擬結(jié)果進行對比分析驗證了三段式的有效性，相關結(jié)論對今后彎管抗沖蝕結(jié)構(gòu)的研究具有一定的參考價值。

1 三段式90°彎管設計

1.1 三段式90°彎管的提出

流體介質(zhì)在流經(jīng)彎管時由于速度的突變?nèi)菀自趶澢萎a(chǎn)生較強的二次流。二次流的存在會使彎管彎曲段的流線變得更加不規(guī)則，同時彎曲段流線的不規(guī)則狀態(tài)還會傳遞到出口直管段。彎曲段流場變得更加復雜會使固體顆粒對管壁的碰撞變得更加頻繁，進而造成更嚴重的沖蝕破壞[15]。由于二次流產(chǎn)生的位置主要分布在彎曲段以及出口段，因此改變彎管彎曲段的管形可以在一定程度上減小彎管中二次流的強度進而減小管壁的沖蝕磨損程度。參考文獻[16]中的航空管路系統(tǒng)三段式彎管設計方法，本文提出了天然氣輸送管道的三段式彎管模型。

1.2 幾何模型

以90°直角彎管為原型，在幾何上保證進口直管段與出口直管段相互垂直的前提下，將原來的一段彎曲連接方式改進為三段彎曲連接方式。新的彎曲段由3個彎頭組成，與進口直管段相連接的彎頭稱為第1彎頭，其曲率半徑為R1，彎曲角度為a；與第1彎頭相連接的彎頭稱為第2彎頭，其曲率半徑為R2，彎曲角度為b；與出口直管段相連接的彎頭稱為第3彎頭，其曲率半徑為R3，彎曲角度為c。圖1為三段式彎管的二維平面圖。

圖1 三段式彎管二維平面Fig.1 Two-dimensional plan of three-section elbow

2 三段式90°彎管優(yōu)化

2.1 優(yōu)化變量

根據(jù)前面建立的三段式90°彎管幾何模型，若要確定管形那么需確定的參數(shù)包括：R1，R2，R3，a，b，c。假設第1彎頭的曲率半徑R1等于原一段式彎管的曲率半徑，并且容易知道a，b，c存在如式(1)關系：

co=90°+bo-ao

(1)

因此只需要確定其中4個管形參數(shù)a，b，R2和R3即可。

2.2 目標函數(shù)

由于二次流主要發(fā)生在彎管彎曲段及出口段，所以本文將優(yōu)化目標設置為：使彎管彎曲段的二次流平均值和彎管出口段的二次流平均值之和最小，即：

Minf(x)=Min(ξavg(w)+ξavg(c))

(2)

式中：ξavg(w)為彎曲段二次流平均值；ξavg(c)為出口段二次流平均值。

2.3 約束條件

首先，必須保證每個彎頭的彎曲角度都大于0；同時，考慮到本文是對一段式90°彎管進行優(yōu)化，為保證結(jié)構(gòu)的合理性，各彎曲段的彎曲角度還應小于90°。

0°

(3)

式中：a為第1彎頭的彎曲角度；b為第2彎頭的彎曲角度；c為第3彎頭的彎曲角度。

根據(jù)GB/T 12459-2017《鋼制對焊管件類型與參數(shù)》和某天然氣長輸管道的現(xiàn)場實際情況，本文將三段式彎管各彎曲段的曲率半徑取值限制在區(qū)間[1.5D,6D]內(nèi)。

1.5D≤R1,R2,R3≤6D

(4)

式中：D為管道外徑，mm；R1為第1彎頭的曲率半徑，mm；R2為第2彎頭的曲率半徑，mm；R3為第3彎頭的曲率半徑，mm。為了避免在進行管形參數(shù)優(yōu)化過程中出現(xiàn)某些管形計算出來的二次流最大值過大或過小的情況，設置式(5)的約束條件。

(5)

式中：ξ1為一段式彎管的二次流最大值；ξ2為優(yōu)化過程中三段式彎管的二次流最大值，其值在優(yōu)化過程中是不斷改變的。

2.4 優(yōu)化算法

COMSOL多物理場仿真軟件提供了多種優(yōu)化算法，主要包括：蒙特卡洛、Nelder-Mead，BOBYQA，COBYLA等?；趦?yōu)化對象、目標函數(shù)及約束條件，在各優(yōu)化算法中COBYLA的計算效果較好，而其他算法在優(yōu)化過程中會偏向于計算局部最優(yōu)，故本文采用COBYLA優(yōu)化算法。

COBYLA算法也被稱為線性近似約束優(yōu)化算法，該方法首先通過假設f(xi)，(i=1,2,…,m)是歐幾里得空間Rn中1個非退化單極點處的函數(shù)值，為得到下1個變量矢量，對該極值點所對應的非線性目標函數(shù)和非線性約束函數(shù)進行插值，并將其看作線性問題進行近似計算[17]。

3 優(yōu)化結(jié)果及分析

3.1 最優(yōu)管形參數(shù)的確定

以直徑400 mm彎徑比2的90°彎管的優(yōu)化為例，在COMSOL軟件中直接建立如圖2所示的三維模型(由于模型的對稱性只需建立其中一半)。

圖2 三段式彎管三維模型Fig.2 3D model of three-section elbow

定義彎管幾何參數(shù)、流體參數(shù)和邊界條件，并在軟件中輸入前面所敘述的約束條件、目標函數(shù)，采用可實現(xiàn)的k-ε湍流模型來模擬連續(xù)相的流動，最后利用COBYLA優(yōu)化算法對三段式彎管的4個管形參數(shù)進行迭代優(yōu)化，優(yōu)化結(jié)果見表1，其中初始值是通過對不同網(wǎng)格數(shù)下的模型進行試算得到的接近于最優(yōu)值的值。初始值的確定可以在一定程度上縮小迭代計算的區(qū)間，減小計算量。

表1 最優(yōu)管形參數(shù)Table 1 Parameters of optimal pipe shape

優(yōu)化后的管形尺寸為：管徑D=400 mm，R1=800 mm，a=53.8°，R2=985.5 mm，b=17.6°，R3=1 174.9 mm，c=53.8°。

3.2 二次流大小的對比

當流體介質(zhì)流經(jīng)彎曲段時會產(chǎn)生不同于主流速度方向上的二次流。二次流的形成是彎曲段管壁內(nèi)外側(cè)壓力梯度和離心力共同作用的結(jié)果，屬于主流流動引起的伴隨流動。本文將二次流近似看成主流速度在各個截面上的速度分量，這樣就可以通過比較各個截面上速度矢量的大小來判斷該截面二次流的強弱。

根據(jù)文獻[18]可知，彎管二次流最強的位置發(fā)生在彎曲段和出口直管段的接壤處。前文以彎曲段二次流平均值和出口段二次流平均值之和最小為優(yōu)化目標，通過優(yōu)化計算得到了最優(yōu)管形，為了驗證優(yōu)化效果，對圖3所示的A，B截面的二次流強弱進行對比。

圖3 二次流強弱對比截面Fig.3 Compared cross sections of secondary flow intensity

表2 彎曲段出口截面二次流的最大值和平均值Table 2 Maximum and average values of the secondary flow of the exit cross section of the bending section

圖4 A，B截面二次流速度矢量分布Fig.4 Distribution of secondary flow velocity vector at A and B cross sections

由表2可知，無論是二次流最大值、平均值還是最小值優(yōu)化后的三段式彎管都遠遠小于一段式彎管，其中二次流最大值下降幅度最大為57.58%。從圖4更能直觀的看到二次流的最大值發(fā)生在近兩頰壁面處，且三段式彎管截面處形成的二次流漩渦要小于一段式的。

4 沖蝕數(shù)值模擬

4.1 管道參數(shù)

為方便對比研究，本次沖蝕數(shù)值模擬共建立了4個彎管模型，如圖5所示。其中一段式90°彎管模型共3個，管徑D=400 mm，彎徑比R/D分別為1.5，2，2.5；三段式彎管的幾何尺寸根據(jù)3.1節(jié)的最優(yōu)管形參數(shù)確定，4個彎管模型的進、出口直管段長度均取為15D。常溫條件下，以天然氣作為連續(xù)相，氣體壓力8 MPa，入口速度為16 m/s，從水平直管流入，從豎直向上直管流出，離散相的顆粒密度2 600 kg/m3，粒徑10 μm，假設顆粒的初始速度與天然氣相同，質(zhì)量流量為0.001 kg/s。

圖5 彎管幾何模型Fig.5 Geometric model of elbow

4.2 網(wǎng)格劃分

以三段式彎管的網(wǎng)格劃分為例，利用ICEM軟件對三段式彎管模型進行網(wǎng)格的劃分，網(wǎng)格均采用六面體結(jié)構(gòu)化網(wǎng)格，并且對彎曲段的網(wǎng)格進行了加密處理，第一層邊界層厚度設置為0.1 mm，增長率設置為1.2，網(wǎng)格劃分情況如圖6所示。

圖6 三段式彎管網(wǎng)格劃分Fig.6 Meshing of three-section elbow

4.3 邊界條件和數(shù)值算法

整個彎管流域可被分為3個邊界：進、出口以及管壁。對于連續(xù)相，湍流計算采用RNGk-ε湍流模型，進口設置為速度進口，出口設置為壓力出口，湍流強度為5%，模擬中采用標準壁面函數(shù)法對近壁面區(qū)域進行處理。對于離散相，DPM模型中的進口和出口均設置為“Escape”，壁面設置為“Reflect”，且采用無滑移壁面。選用SIMPLE算法來計算壓力-速度耦合，利用二階中心差分格式和二階迎風差分格式分別計算擴散項和對流項，殘差均設置為10-5。本文中顆粒的碰撞模型選用碳鋼材料的碰撞模型，沖蝕預測模型選用Oka沖蝕模型[19]。由于經(jīng)過多次分離過濾后，天然氣輸送管道內(nèi)氣體所含固體量很少，所以在計算過程中忽略固體顆粒對連續(xù)相流場的影響，采用單相耦合的方法求解固體顆粒的沖蝕速率。

4.4 模型有效性驗證

為驗證本文所建立的沖蝕計算模型的有效性，對某氣田集氣站內(nèi)的彎頭進行仿真建模，結(jié)合該集氣站的實際工況條件進行參數(shù)設置，將FLUENT沖蝕模擬的計算結(jié)果與現(xiàn)場檢測報告進行了對比。所建立的90°彎頭模型的管徑D=76 mm，彎徑比R/D=1.5，流動參數(shù)設置見表3。

表3 流動參數(shù)設置Table 3 Setting of flow parameters

數(shù)值模擬計算得出的彎頭最大沖蝕率為6.541×10-8kg/(m2·s)，根據(jù)2010年該集氣站的檢測報告，彎頭的最大減薄量在1.2 mm左右，對應的最大沖蝕率為7.468×10-8kg/(m2·s)[1]，誤差在可接受范圍內(nèi)，因此，本文建立的沖蝕計算模型可以用來預測天然氣管道彎管的沖蝕磨損情況。

4.5 速度分布

圖7所示為三段式彎管中心截面上的速度分布云圖。由圖7可知，入口直管段的速度分布比較均勻，當氣流進入彎曲段后，速度分布發(fā)生了較大的改變：在第1彎頭的內(nèi)拱側(cè)附近和第2彎頭的外拱側(cè)附近出現(xiàn)了速度的最大值，在第2彎頭的內(nèi)拱側(cè)附近則出現(xiàn)了低速區(qū)，速度減小了很多，尤其是在第1彎頭和第2彎頭的交界處；第1彎頭內(nèi)拱側(cè)的速度大于外拱側(cè)的速度，而第2，第3彎頭外拱側(cè)的速度大于內(nèi)拱側(cè)的速度，同時可以看出，第2彎頭內(nèi)流場徑向的速度梯度最大。隨著氣流進入出口直管段，截面上的速度梯度趨于緩和。

圖7 三段式彎管速度云圖Fig.7 Velocity nephogram of three-section elbow

4.6 沖蝕模擬結(jié)果分析

圖8為三段式彎管的沖蝕率云圖。觀察圖8發(fā)現(xiàn)：三段式彎管的沖蝕磨損嚴重區(qū)域為第1彎頭末端45°～53.8°靠近內(nèi)拱兩頰處以及第2彎頭初始位置的內(nèi)拱壁處。第3彎頭的沖蝕磨損嚴重程度相比于前兩段基本可忽略不計。表4給出了各個彎管模型的最大沖蝕率，其中三段式彎管的最大沖蝕率為2.41×10-10kg/(m2·s)，該數(shù)值甚至比R/D=2.5的一段式彎管的最大沖蝕率3.08×10-10kg/(m2·s)還小0.67×10-10kg/(m2·s)，因此，三段式彎管的局部耐磨性能得到改善。

圖8 三段式彎管沖蝕率云圖Fig.8 Erosion rate nephogram of three-section elbow

管形三段式彎管R/D=1.5R/D=2R/D=2.5最大沖蝕速率(×10-10kg·m-2·s-1)2.4111.14.613.08

圖9展示了顆粒在三段式彎管和R/D=2.5的一段式彎管中的運動軌跡。觀察圖9(a)發(fā)現(xiàn)：當氣體流經(jīng)三段式彎管第1彎頭時僅有少量固體顆粒與外拱壁碰撞并反彈至第1彎頭末端附近兩頰處；顆粒在第2，3彎頭內(nèi)運動十分順暢幾乎沒有與壁面發(fā)生碰撞，這就造成了三段式彎管嚴重沖蝕的位置是在第1彎頭末端附近的現(xiàn)象。對比圖9(a)和(b)發(fā)現(xiàn)：在兩者沖蝕最嚴重的區(qū)域，R/D=2.5的一段式彎管中的顆粒運動軌跡要更加混亂，可以看出三段式彎管流場的平穩(wěn)性得到了較大的改善，固體顆粒能更加順暢地流過彎管。

圖10是三段式彎管和R/D=2.5的一段式彎管在Z=0 mm切面處的湍流強度云圖，對比發(fā)現(xiàn)：三段式彎管內(nèi)拱側(cè)附近大部分區(qū)域的湍流強度在90%左右，而R/D=2.5的一段式彎管內(nèi)拱側(cè)附近的湍流強度達到了115%左右。局部湍流強度高會增強顆粒運動的無序性，加劇顆粒對壁面的碰撞。因此，三段式彎管的最大沖蝕率要小于R/D=2.5的一段式彎管的最大沖蝕率。

圖9 顆粒軌跡對比Fig.9 Comparison of particles trajectory

圖10 湍流強度對比Fig.10 Comparison chart of turbulence intensity

5 結(jié)論

1)對D=400 mm，R/D=2的一段式90°彎管進行三段式改進，利用COMSOL對4個管形參數(shù)進行優(yōu)化，最終得到4個管形參數(shù)的數(shù)值分別為：a=53.8°；b=17.6°；R2=985.5 mm；R3=1 174.9 mm。

2)優(yōu)化后的三段式彎管流場更加平穩(wěn)，彎曲段二次流的強度大幅降低，以一段式的截面A和三段式的截面B為例，二次流的最大值、最小值以及平均值分別從5.128 65，0.673 39，2.199 2 m/s降低到了2.175 47，0.357 72，1.133 5 m/s，下降幅度分別為57.58%，46.88%，48.46%。

3)Fluent沖蝕模擬結(jié)果表明，三段式彎管沖蝕嚴重的區(qū)域主要位于第1彎頭末端45°～53.8°靠近內(nèi)拱兩頰處以及第2彎頭初始位置的內(nèi)拱壁處。三段式彎管的最大沖蝕率為2.41×10-10kg/(m2·s)，比起D=400 mm，R/D=1.5，2，2.5的一段式90°彎管分別下降了約80%，50%，20%，其耐磨損性能大大提高。