張永星, 顧鴻博, 陳曉玲

(1.中國石油大學(xué)(北京)機(jī)械與儲運(yùn)工程學(xué)院, 北京 102249; 2.中國石油大學(xué)(北京)克拉瑪依校區(qū)工學(xué)院, 克拉瑪依 834000)

含蠟天然氣是指含有較高碳數(shù)的正烷烴,或含有由直鏈、支鏈烷烴所組成高分子碳?xì)浠衔锏葟?fù)雜混合物的天然氣。在溫度降低到析蠟溫度以下時,會有蠟晶析出,蠟晶經(jīng)過不斷運(yùn)動,最終在一定條件下,氣相中的蠟顆粒會沉積于管道內(nèi)壁。

目前國內(nèi)外天然氣集輸系統(tǒng)均面臨著不同程度的蠟沉積問題。以新疆克深2氣田處理裝置為例,原料氣在經(jīng)過J-T閥時,由于節(jié)流膨脹導(dǎo)致低溫分離器發(fā)生蠟堵,最終影響裝置的正常運(yùn)行[1]。在克拉蘇氣田大北區(qū)塊天然氣處理廠、塔河九區(qū)天然氣處理廠、克拉2氣田、塔里木油田凝析氣藏等同樣發(fā)生了不同程度的蠟堵現(xiàn)象[2-5]。伊朗南部的南帕爾斯天然氣加工廠站場設(shè)備內(nèi)也發(fā)現(xiàn)有蠟顆粒的沉積,引起站內(nèi)設(shè)備的安全問題,影響生產(chǎn)進(jìn)度[6]。因此,天然氣管道處理站管道蠟沉積不僅會減小管道的有效流通面積,降低管道的輸送能力,嚴(yán)重時還會造成蠟堵事故[7],引發(fā)安全問題。雖然含蠟天然氣在集輸系統(tǒng)各個環(huán)節(jié)都有著不同程度的蠟沉積現(xiàn)象,但由于節(jié)流閥前后會引起氣體的壓力、溫度大幅變化、從而使閥內(nèi)蠟沉積所帶來的問題顯得更為突出。

目前針對含蠟天然氣在節(jié)流設(shè)備中的蠟沉積問題的研究較少。李俊逸等[8]采用數(shù)值模擬的方法研究得到了硫沉積速率隨氣流進(jìn)口速度的增大而先增大后減小的規(guī)律。劉成川等[9]建立了元素硫在多場耦合作用下的沉積模型,認(rèn)為沉積主要取決于壓力。溫度等參數(shù)變化。劉巖等[10]采用計(jì)算流體力學(xué)與離散單元法相耦合,認(rèn)為顆粒在沉積過程中主要受二次流的影響,大部分沉積在壁面兩側(cè)附近。Sun等[11]通過Fluent進(jìn)行氣固兩相模擬,并認(rèn)為顆粒運(yùn)動與斯托克斯數(shù)緊密相關(guān)。Wang等[12]借助DPM(deformable part model)模型對天然氣內(nèi)水合物顆粒進(jìn)行模擬,得到了管道底部顆粒物濃度較大、頂部顆粒物濃度較小的趨勢規(guī)律。

然而,目前針對傳統(tǒng)孔板式節(jié)流閥的研究主要集中在焦耳湯姆遜效應(yīng)及其水合物生成等現(xiàn)象[13-16],圍繞節(jié)流閥內(nèi)蠟沉積規(guī)律的研究相對較少,對孔板式節(jié)流閥內(nèi)蠟沉積規(guī)律認(rèn)識不明確,難以有效地對實(shí)際生產(chǎn)中蠟沉積進(jìn)行主動防治。因此,開展節(jié)流閥設(shè)備蠟沉積問題的研究,揭示節(jié)流閥內(nèi)高流速下的蠟沉積規(guī)律,對提高含蠟天然氣安全高效運(yùn)行具有重要意義。為此,現(xiàn)借助計(jì)算流體力學(xué)軟件,通過UDF(user-defined function)編寫自定義函數(shù),實(shí)現(xiàn)蠟的沉積與剝離,建立孔板式節(jié)流閥的蠟沉積仿真模型,分析孔板式節(jié)流閥內(nèi)蠟沉積特性。

1 數(shù)值模型

1.1 模型假設(shè)

由于天然氣中含蠟量較低,故假設(shè)流動為稀疏的氣-固兩相流,不考慮由于溫降可能導(dǎo)致凝析油的產(chǎn)生;顆粒為球形,且符合正態(tài)分布;顆粒主要受到流體的曳力、重力、阻力、浮力等作用;由于在較低溫度的條件時,蠟的黏度變化不大,暫時忽略溫度對于黏度的影響。

1.2 數(shù)學(xué)模型

1.2.1 連續(xù)相湍流方程

k-ε(k為湍動能,ε為湍流耗散率)是目前使用最廣泛的湍流方程[17],對流場求解具有普遍的適用性。該模型通過分別求解湍動能及湍流耗散率方程來確定湍流長度和時間尺度,而湍流耗散率方程則是通過物理推理并類比相似方程得出[18],其具體控制方程見文獻(xiàn)[19]。

1.2.2 離散相運(yùn)動方程

采用DPM模型對顆粒運(yùn)動軌跡進(jìn)行追蹤,以x方向?yàn)槔?顆粒所受的平衡方程為

(1)

式(1)中:t為時間,s;u為天然氣流動速度, m/s;us為固體顆粒速度,m/s;ρ為天然氣密度,kg/m3;ρs為固體顆粒密度,kg/m3;gx為x方向的重力加速度分量,m/s2;μ為天然氣動力黏度,Pa·s;ds為蠟顆粒直徑,μm;CD為顆粒阻力系數(shù);Fx為顆粒所受x方向的外力,N。

1.3 臨界速度模型

臨界速度模型[20]通過研究顆粒的彈性性質(zhì),以顆粒物與壁面撞擊時的法向速度作為判斷依據(jù),當(dāng)顆粒的法向速度大于臨界速度,顆粒就會反彈;反之,顆粒就會發(fā)生黏附。圖1所示為單個顆粒在流體中的受力示意圖。

FL為流體作用在顆粒上的吸引力;a為與壁面接觸面積的半徑;Fpo為黏附力;FD為曳力;Dp/2-α為顆粒圓心到壁面距離圖1 單個顆粒受力示意圖Fig.1 Schematic diagram of individual particle stress

具體計(jì)算公式為

(2)

(3)

(4)

(5)

式中:Vcr為臨界速度;E為混合的楊氏模量;Es為壁面的楊氏模量;Ep為顆粒的楊氏模量;νs為壁面的泊松比;νp為顆粒的泊松比;Dp為顆粒直徑;ρp為顆粒的密度;R為運(yùn)動恢復(fù)系數(shù),這里取0.7。

Soltani等[21]通過大量實(shí)驗(yàn)對不同的剝離機(jī)制進(jìn)行研究,發(fā)現(xiàn)圓形顆粒發(fā)生剝離的主要方式為滾動。Dasetal[22]同樣發(fā)現(xiàn),滾動剝離是圓球形顆粒剝離的主要形式。故假設(shè)當(dāng)蠟晶顆粒發(fā)生剝離時,只發(fā)生滾動剝離,且為單個顆粒剝離,直接剝離與滑動剝離不會發(fā)生。滾動剝離方程[23]為

(6)

式(6)中:α為顆粒與壁面接觸后顆粒直徑所減少的部分;FL為流體作用在顆粒上的剪切力;a為顆粒與壁面的接觸半徑;針對球形顆粒,FD為所受到的曳力,可以表示為

(7)

式(7)中:uc為顆粒所在位置的流速;Acp為顆粒截面積;f為修正系數(shù);Cu為Cunningham修正系數(shù),由Sundaram[24]提出,具體表達(dá)式為

(8)

式(8)中:Kn為克努森數(shù),用來表征稀疏效應(yīng)的影響,不可忽略[25]。

Fpo具體表達(dá)式為

(9)

式(9)中:Fpo為黏附力;WA為黏著功,J/m2。黏著功取決于顆粒與壁面材料,通常通過實(shí)驗(yàn)得到,關(guān)于對蠟顆粒黏著功研究較少,在此假設(shè)與硅顆粒相同,取0.039 J/m2。

2 模型驗(yàn)證

2.1 物理模型及網(wǎng)格劃分

以文獻(xiàn)[26]中的孔板式節(jié)流閥為研究對象,結(jié)構(gòu)見圖2。其中直管段直徑D=76 mm,孔板的孔徑d=30 mm,孔板的厚度t=10 mm,左側(cè)為節(jié)流閥入口,右側(cè)為出口。為保證湍流穩(wěn)定性,在節(jié)流閥進(jìn)口和出口分別增加10D和15D的直管段。采用ICEM軟件對計(jì)算域進(jìn)行六面體網(wǎng)格劃分,并插入膨脹層以加密壁面附近的網(wǎng)格,如圖3所示。

圖2 孔板式節(jié)流閥及路徑示意Fig.2 Orifice type throttle valve and path diagram

圖3 孔板式節(jié)流閥網(wǎng)格模型Fig.3 Orifice throttle valve mesh model

2.2 求解策略

計(jì)算過程中將二階中心差分格式用于離散擴(kuò)散項(xiàng),二階迎風(fēng)差分格式用于其他項(xiàng),時間步長為1×10-4,先對連續(xù)相進(jìn)行穩(wěn)態(tài)計(jì)算,收斂后加入顆粒相進(jìn)行非穩(wěn)態(tài)耦合計(jì)算。采用SIMPLE(semi-implicit method for pressure-linked equations)算法模擬流場特征[27],并通過UDF實(shí)現(xiàn)蠟晶顆粒的沉積、反彈、剝離等特征。

2.3 定解條件與工況設(shè)定

連續(xù)相采用壓力入口及壓力出口,壁面邊界為無滑移壁面。天然氣密度采用理想氣體模型,動力黏度為1.087×10-5kg/(m·s)。在距入口0.76 m處建立顆粒射流源平面,其中顆粒質(zhì)量流量為1.3×10-6kg/s。采用58號石蠟材料,石蠟顆粒密度為1 980 kg/m3、比熱容為986 J/ (kg·K)。同時對顆粒加載虛擬質(zhì)量力及壓力梯度力。離散相進(jìn)、出口邊界處采用逃逸邊界,壁面采用UDF自定義函數(shù)。氣流入口壓力為5～8 MPa。設(shè)定節(jié)流閥孔徑為20～60 mm,喉道長度為10～40 mm。數(shù)值模擬工況設(shè)定見表1。

表1 數(shù)值模擬工況

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

為了確保數(shù)值模擬結(jié)果獨(dú)立于網(wǎng)格,以喉道內(nèi)中軸線處速度值為參考,分別對72萬、100萬、128萬、142萬以及173萬網(wǎng)格數(shù)進(jìn)行無關(guān)性分析,結(jié)果見圖4。

圖4 網(wǎng)格無關(guān)性分析Fig.4 Mesh independence analysis

通過圖4分析可知,當(dāng)沿流動方向距離大于0.76 m后,72萬網(wǎng)格數(shù)時的速度明顯低于其他網(wǎng)格數(shù)時的速度,而對于他網(wǎng)格數(shù),速度曲線基本重合。綜合考慮計(jì)算精度及計(jì)算速度后,最終選用128萬網(wǎng)格數(shù)進(jìn)行后續(xù)數(shù)值計(jì)算。

以文獻(xiàn)[26]中五處溫度監(jiān)測點(diǎn)作為參考,將模擬結(jié)果與文獻(xiàn)中實(shí)驗(yàn)節(jié)點(diǎn)的溫度值進(jìn)行比較,具體結(jié)果見表2。

表2 文獻(xiàn)溫度與模擬溫度對比

通過表2可知,Fluent模擬結(jié)果與實(shí)驗(yàn)結(jié)果最大誤差為4.7%,誤差均不超過5%,說明所選取的模型適用于本研究。

2.5 蠟沉積模型驗(yàn)證

1995年,美國部分學(xué)者發(fā)現(xiàn)在孔板式流量計(jì)中存在固體顆粒,會沉積到管壁,影響測量準(zhǔn)確性[28]。為驗(yàn)證臨界速度模型的準(zhǔn)確性,在2.1節(jié)孔板式節(jié)流閥模型的基礎(chǔ)上,從入口(x=0)處注入蠟晶顆粒,將數(shù)值模擬與現(xiàn)場中觀察沉積位置進(jìn)行對比。圖5為基于臨界速度模型由數(shù)值模擬得到的蠟沉積分布,而圖6為文獻(xiàn)[28]中實(shí)驗(yàn)觀察到的蠟晶沉積位置。

圖5 模擬蠟沉積位置Fig.5 Simulate wax deposition locations

圖6 現(xiàn)場固體顆粒沉積位置[28]Fig.6 On-site solid particle deposition location[28]

2020年,Peng等[29]學(xué)者對頁巖氣開采過程中孔板內(nèi)顆粒沉積進(jìn)行分析。為進(jìn)一步驗(yàn)證蠟沉積模型的準(zhǔn)確性,將數(shù)值模擬與實(shí)驗(yàn)中得到的沉積位置進(jìn)行對比。圖7為基于臨界速度模型由數(shù)值模擬得到的沉積分布,圖8為文獻(xiàn)[29]中實(shí)驗(yàn)得到的沉積位置。

圖7 模擬顆粒沉積位置Fig.7 Simulate particle deposition location

圖8 文獻(xiàn)中實(shí)驗(yàn)沉積位置[29]Fig.8 Experimental deposition location in literature [29]

從數(shù)值模擬得到的圖5可以看出,管壁四周均有沉積的蠟顆粒,但大部分蠟顆粒沉積在管道下方,與文獻(xiàn)中的實(shí)驗(yàn)結(jié)果圖6沉積趨勢相同。從模擬得到的圖7可以看出,孔板內(nèi)均出現(xiàn)顆粒沉積,隨著氣相速度增大,沉積量逐漸增加,并且在孔板入口截面處出現(xiàn)少量沉積,與文獻(xiàn)[29]中實(shí)驗(yàn)結(jié)果圖8沉積趨勢相同。通過數(shù)值模擬與實(shí)驗(yàn)結(jié)果的比較,發(fā)現(xiàn)臨界速度模型可以較為準(zhǔn)確的描述蠟顆粒沉積情況。因此,文中后續(xù)模擬分析將采用臨界速度模型對蠟沉積規(guī)律進(jìn)行研究。

3 結(jié)果分析

3.1 壓力對沉積率的影響

圖9所示為最大沉積量隨壓力變化趨勢,從中可以看出,當(dāng)管內(nèi)入口壓力從5 MPa增加到6 MPa,最大沉積量由1.06×10-6kg增加至1.67×10-6kg,增加了36.52%,當(dāng)壓力從6 MPa增加到7 MPa時,最大沉積量由1.67×10-6kg增加至1.92×10-6kg,增加了13.02%,而當(dāng)壓力從7 MPa增加到8 MPa時,最大沉積量則由1.92×10-6kg增加至2.11×10-6kg,只增加了9%。發(fā)現(xiàn)隨著進(jìn)口壓力的增加,管道內(nèi)壁的最大沉積量逐漸增加,但沉積量增加的趨勢逐漸變緩。蠟沉積量隨入口壓力呈現(xiàn)上述變化趨勢可能主要由以下兩個方面的原因。一方面由于隨著壓力的增加,管道內(nèi)的流速也隨之增加,從而使管內(nèi)湍流強(qiáng)度增加,增大了蠟顆粒撞擊壁面形成沉積的幾率,所以蠟沉積量隨著入口壓力的增大而增大;另一方面,隨著壓力增大,顆粒法相速度隨之增大,達(dá)到了部分較小粒徑的沉積條件,但由于顆粒粒徑變小,且由于粒徑呈正態(tài)分布,所以小顆粒數(shù)量相對較少,因此沉積量增加的趨勢逐漸變緩。

圖9 沉積量變化趨勢Fig.9 Variation trend of deposition rate

圖10所示為管道內(nèi)蠟沉積的位置。從圖10可以看出,蠟顆粒雖然在節(jié)流閥后的整個管道內(nèi)均有沉積,但主要沉積位置發(fā)生在喉道出口(x=0.77～0.775 m)處以及距喉道出口0.17 m(x=0.94 m)左右的管壁處。第二處沉積位置隨著壓力增大逐漸向管道出口移動,壓力增加1 MPa,沉積位置向后移動約0.005 m,并且沉積寬度為0.01 m左右。當(dāng)壓力為5 MPa時,沉積位置主要在0.935～0.945 m;當(dāng)壓力為6 MPa時,沉積位置主要在0.938～0.948 m;當(dāng)壓力為7 MPa時,沉積位置主要在0.945～0.955 m;當(dāng)壓力為8 MPa時,沉積位置主要在0.948～0.958 m。管道內(nèi)蠟顆粒主要沉積在上述兩個位置可能是湍動能強(qiáng)度和渦流共同作用的結(jié)果。

圖10 模擬蠟顆粒沉積位置示意圖Fig.10 Schematic diagram of simulated wax particle deposition location

圖11所示為沿流動方向中軸線處湍動能分布情況,從中可以看出湍動能在喉道出口(x=0.77 m)和距喉道出口0.17 m處附近,即兩處主要的沉積位置,出現(xiàn)了兩個明顯的峰值。說明在該處速度的波動幅度更大,并且這種劇烈程度隨著壓力的增加而增大,從而導(dǎo)致上述兩處位置蠟沉積量隨著壓力增加而增加。同時發(fā)現(xiàn)第二處湍動能峰值位置隨著壓力的增大逐漸后移,該趨勢與主要沉積位置隨著壓力的增大相應(yīng)后移相一致。另外還發(fā)現(xiàn)在喉道出口(x=0.77 m)位置處,隨著壓力的增加,湍動能的增大幅度逐漸減小,并在7～8 MPa時十分接近,而在距喉道出口0.17 m處附近,隨著壓力的增加,湍動能的增大幅度差別不大,從而導(dǎo)致隨著壓力的增大后者的湍動能大于前者的湍動能。這也是導(dǎo)致在距喉道出口0.17 m處蠟沉積量大于喉道出口處沉積量的原因之一。

圖11 沿流動方向湍動能變化Fig.11 Variation of turbulent kinetic energy along the flow direction

圖12為兩個主要沉積截面處垂直于流動方向的湍動能分布。從圖12(b)中可以看到,湍動能出現(xiàn)中間下凹兩端凸起的駝峰形主要是因?yàn)檩S線處速度較為穩(wěn)定,隨著氣流向徑向方向流動,形成漩渦,不穩(wěn)定程度加劇,湍動能逐漸增大,達(dá)到峰值。圖13所示為所研究最低壓力(5 MPa)和最高壓力(8 MPa)下管道內(nèi)的流場圖。氣體經(jīng)過節(jié)流閥后,速度明顯增加,并在沉積處產(chǎn)生兩處渦流。隨著壓力增大,渦流影響范圍沒有發(fā)生明顯變化,但變化程度更為劇烈。當(dāng)顆粒隨流體運(yùn)動過程中,在渦流位置由于流體方向發(fā)生突然改變,顆粒在慣性作用下撞擊到壁面,發(fā)生沉積。因此,隨著壓力增大,顆粒與壁面發(fā)生撞擊的幾率越大,最終導(dǎo)致管內(nèi)沉積量增加。

圖12 沉積截面處垂直于流動方向的湍動能變化Fig.12 Variation of turbulent kinetic energy at the deposition section perpendicular to the flow direction

3.2 蠟顆粒直徑的影響

圖14所示為顆粒粒徑與最大沉積量之間的關(guān)系?？梢钥闯?在節(jié)流直徑、喉道長度一致時,最大沉積量隨粒徑增大而增大,并在粒徑為50 μm時發(fā)生陡增。以入口壓力為5 MPa時為例,粒徑由30 μm增加到50 μm時,最大沉積量由0.067 7×10-6kg增加到0.142 6×10-6kg,增加了52.52%;粒徑由80 μm增加到100 μm時,最大沉積量由0.739 4×10-6kg增加到2.145 4×10-6kg,增加了65.53%。這是因?yàn)殡S著粒徑的增大,蠟顆粒具有更大的慣性。因此,當(dāng)渦流方向突然改變時,在慣性的作用下,較大粒徑的蠟顆粒更易保持原來的方向運(yùn)動,脫離原流場,從而撞擊至節(jié)流閥管壁的幾率增大,導(dǎo)致更大的沉積量。

圖14 粒徑對沉積量的影響Fig.14 Effect of particle size on deposition amount

3.3 節(jié)流孔徑的影響

圖15所示為入口壓力為5 MPa時,節(jié)流孔直徑與最大沉積量的關(guān)系。分析可知,節(jié)流孔徑增大,管壁上蠟顆粒的最大沉積量逐漸下降。當(dāng)孔徑由20 mm增加至30 mm時,管壁最大沉積量由1.293×10-6kg降至0.952 7×10-6kg,下降26.37%;當(dāng)孔徑繼續(xù)增大至40 mm時,管壁沉最大沉積量由0.952 7×10-6kg降至0.493 9×10-6kg,下降48.15%;當(dāng)管徑為50 mm和60 mm時,沉積量接近為0。從不同孔徑下流線示意圖16可以看出,由于管內(nèi)為超音速流動,隨著孔徑的增大,速度也隨之增大。并且渦流明顯變小,顆粒大部分隨氣流向前流動,顆粒撞擊壁面幾率減小,節(jié)流孔出口壁面處對蠟顆粒的阻礙作用減小,因此在大孔徑下,蠟顆粒更容易被氣流帶走,而不發(fā)生沉積。

圖15 孔徑對沉積量的影響Fig.15 Effect of pore size on deposition amount

3.4 喉道長度的影響

圖17所示為管道內(nèi)最大沉積量隨喉道長度的變化規(guī)律。入口壓力為5 MPa,隨著喉道長度增加,管內(nèi)最大沉積量逐漸上升,但增加幅度逐漸減小。由10 mm增加至20 mm時,最大沉積量增加了20.83%,湍動能增加15.02%;當(dāng)增加至30 mm時,沉積率增加了13.99%,湍動能增加3.38%;當(dāng)增加至40 mm時,最大沉積量增加了6.43%,湍動能增加0.038%。

圖17 喉道長度對沉積量的影響Fig.17 Effect of throat length on deposition amount

圖18為沿流動方向中軸線處湍動能變化,隨著喉道長度的增加,喉道出口截面處湍動能略有增大,增幅逐漸減小,而距喉道出口0.17 m處左右,湍動能逐漸減小,并且峰值出現(xiàn)后移。這是因?yàn)殡S著喉道長度的增加,氣流經(jīng)過節(jié)流后更容易達(dá)到穩(wěn)定的狀態(tài),導(dǎo)致湍動能變化不大。這表明此時流動已經(jīng)逐漸穩(wěn)定,更接近于圓管內(nèi)湍動能分布。由于湍動能變化幅度減小,最終導(dǎo)致最大沉積量增長率也隨之減小。

圖18 沿流動方向湍動能變化Fig.18 Variation of turbulent kinetic energy along the flow direction

4 結(jié)論

(1)傳統(tǒng)孔板式節(jié)流閥的蠟沉積主要出現(xiàn)在2個位置:節(jié)流閥節(jié)流出口(x=0.77～0.775 m)處壁面處以及距節(jié)流閥喉道出口0.17 m左右的管壁處,并且管道內(nèi)顆粒最大沉積量隨壓力的增加而增加,但增加趨勢逐漸變緩。沉積寬度為0.01 m,當(dāng)壓力增加1 MPa時,第二處沉積位置后移約0.005 m。

(2)在節(jié)流直徑、喉道長度一定時,管道內(nèi)最大沉積量隨著粒徑的增大而增大,并且當(dāng)粒徑大于50 μm后,最大沉積量明顯增加。

(3)管道內(nèi)壁最大沉積量與節(jié)流孔徑呈負(fù)相關(guān)。其中,當(dāng)節(jié)流孔徑增加到50 mm時,最大沉積量幾乎為0;最大沉積量與喉道長度呈正相關(guān),但隨著喉道長度的增加,管內(nèi)最大沉積量的增長趨勢逐漸變緩。