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

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

0 引言

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

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

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

1 三段式90°彎管設(shè)計(jì)

1.1 三段式90°彎管的提出

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

1.2 幾何模型

以90°直角彎管為原型，在幾何上保證進(jìn)口直管段與出口直管段相互垂直的前提下，將原來的一段彎曲連接方式改進(jìn)為三段彎曲連接方式。新的彎曲段由3個(gè)彎頭組成，與進(jìn)口直管段相連接的彎頭稱為第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。假設(shè)第1彎頭的曲率半徑R1等于原一段式彎管的曲率半徑，并且容易知道a，b，c存在如式(1)關(guān)系：

co=90°+bo-ao

(1)

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

2.2 目標(biāo)函數(shù)

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

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

(2)

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

2.3 約束條件

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

0°

(3)

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

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

1.5D≤R1,R2,R3≤6D

(4)

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

(5)

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

2.4 優(yōu)化算法

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

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

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

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

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

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

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

表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 二次流大小的對(duì)比

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

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

圖3 二次流強(qiáng)弱對(duì)比截面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)化后的三段式彎管都遠(yuǎn)遠(yuǎn)小于一段式彎管，其中二次流最大值下降幅度最大為57.58%。從圖4更能直觀的看到二次流的最大值發(fā)生在近兩頰壁面處，且三段式彎管截面處形成的二次流漩渦要小于一段式的。

4 沖蝕數(shù)值模擬

4.1 管道參數(shù)

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

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

4.2 網(wǎng)格劃分

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

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

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

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

4.4 模型有效性驗(yàn)證

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

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

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

4.5 速度分布

圖7所示為三段式彎管中心截面上的速度分布云圖。由圖7可知，入口直管段的速度分布比較均勻，當(dāng)氣流進(jìn)入彎曲段后，速度分布發(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è)的速度，同時(shí)可以看出，第2彎頭內(nèi)流場(chǎng)徑向的速度梯度最大。隨著氣流進(jìn)入出口直管段，截面上的速度梯度趨于緩和。

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

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

圖8為三段式彎管的沖蝕率云圖。觀察圖8發(fā)現(xiàn)：三段式彎管的沖蝕磨損嚴(yán)重區(qū)域?yàn)榈?彎頭末端45°～53.8°靠近內(nèi)拱兩頰處以及第2彎頭初始位置的內(nèi)拱壁處。第3彎頭的沖蝕磨損嚴(yán)重程度相比于前兩段基本可忽略不計(jì)。表4給出了各個(gè)彎管模型的最大沖蝕率，其中三段式彎管的最大沖蝕率為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的一段式彎管中的運(yùn)動(dòng)軌跡。觀察圖9(a)發(fā)現(xiàn)：當(dāng)氣體流經(jīng)三段式彎管第1彎頭時(shí)僅有少量固體顆粒與外拱壁碰撞并反彈至第1彎頭末端附近兩頰處；顆粒在第2，3彎頭內(nèi)運(yùn)動(dòng)十分順暢幾乎沒有與壁面發(fā)生碰撞，這就造成了三段式彎管嚴(yán)重沖蝕的位置是在第1彎頭末端附近的現(xiàn)象。對(duì)比圖9(a)和(b)發(fā)現(xiàn)：在兩者沖蝕最嚴(yán)重的區(qū)域，R/D=2.5的一段式彎管中的顆粒運(yùn)動(dòng)軌跡要更加混亂，可以看出三段式彎管流場(chǎng)的平穩(wěn)性得到了較大的改善，固體顆粒能更加順暢地流過彎管。

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

圖9 顆粒軌跡對(duì)比Fig.9 Comparison of particles trajectory

圖10 湍流強(qiáng)度對(duì)比Fig.10 Comparison chart of turbulence intensity

5 結(jié)論

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

2)優(yōu)化后的三段式彎管流場(chǎng)更加平穩(wěn)，彎曲段二次流的強(qiáng)度大幅降低，以一段式的截面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é)果表明，三段式彎管沖蝕嚴(yán)重的區(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%，其耐磨損性能大大提高。