王健， 顧曉婷*， 張瑤瑤， 楊燕華

(1.長(zhǎng)江大學(xué)石油工程學(xué)院， 武漢 430100； 2.油氣鉆采工程湖北省重點(diǎn)實(shí)驗(yàn)室， 武漢 430100；3.新奧天然氣股份有限公司， 廊坊 065000)

中俄原油管道穿越中國(guó)多年凍土區(qū)和季節(jié)凍土區(qū)，凍土地區(qū)活動(dòng)層易受溫度影響導(dǎo)致土體凍結(jié)發(fā)生凍脹現(xiàn)象。管道敷設(shè)在凍土中受土壤活動(dòng)影響較大，土體會(huì)因水熱現(xiàn)象形成冰錐、凍脹丘等，進(jìn)而影響管道的受力變形。由于管道需要敷設(shè)在管溝中，管溝開挖和凍土回填時(shí)管道周圍土體密度不均造成管溝積水[1]，再經(jīng)凍融循環(huán)后，管道局部會(huì)出現(xiàn)應(yīng)力集中現(xiàn)象導(dǎo)致屈服彎曲。

Wu等[2]利用有限元數(shù)值模擬，研究了由凍脹丘引起的管道受力變形及應(yīng)力分布。Wen等[3]建立彈塑性有限元模型，分析管道在穿越凍土區(qū)時(shí)管道隨凍脹作用的最大變形及受力。王國(guó)麗等[4]論述了基于應(yīng)變?cè)O(shè)計(jì)方法下的管道因凍脹融沉引起的變形計(jì)算方法，并進(jìn)行了實(shí)例分析。王洪波等[5]以簡(jiǎn)化的凍脹模型分析了輸氣管道與土體的相互作用，為輸氣管道在凍土中的設(shè)計(jì)提供了相關(guān)理論。黃龍等[6]基于彈性地基梁，利用有限元模型分析了地基特性參數(shù)與溫度的關(guān)系，研究了管道在不同地基參數(shù)下的應(yīng)力變化。狄彥等[7]提出了凍土凍脹對(duì)管道的計(jì)算方法和模型以及管道基于應(yīng)變的失效判據(jù)。李超營(yíng)等[8]研究了凍土區(qū)不同斜坡角度對(duì)埋地管道的影響，分析了斜坡角度和融沉長(zhǎng)度對(duì)管道應(yīng)力變化。張旭等[9]利用非線性數(shù)值模擬方法研究漠大線管道的應(yīng)變影響規(guī)律，分析得出管道壁厚對(duì)管道受力變形影響最大。上述研究雖論述了管道及凍脹區(qū)土壤參數(shù)與管道應(yīng)力應(yīng)變的關(guān)系，但未考慮管道周圍凍土實(shí)際溫度場(chǎng)，只施加凍脹載荷進(jìn)行模擬分析，未能分析土壤溫度變化對(duì)管道應(yīng)力的影響且均未提及管溝參數(shù)。而實(shí)際上，根據(jù)《輸油管道設(shè)計(jì)規(guī)范》(GB 50253—2014)[10]，管道敷設(shè)時(shí)需進(jìn)行管溝開挖、土壤換填和回填，由于管溝換填土質(zhì)類型及含冰量等的不同，凍脹量也會(huì)有差異。尤其是在過渡區(qū)與凍脹穩(wěn)定區(qū)，會(huì)引起管道的彎曲變形，管溝參數(shù)對(duì)管道應(yīng)力應(yīng)變的影響更直接、經(jīng)濟(jì)有效。

為使有限元分析結(jié)果更符合管道敷設(shè)在一定管溝尺寸下的實(shí)際要求，針對(duì)輸油管道利用ABAQUS有限元軟件，基于耦合場(chǎng)的分析方法，考慮鋼質(zhì)管道材料的非線性特性和土壤材料隨溫度變化的非線性特性，管道內(nèi)壁承受油溫及壓力載荷，外壁承受非線性熱應(yīng)力以及經(jīng)由大氣、土壤、管道之間傳熱所獲的溫度載荷。根據(jù)中俄原油管道工程漠大線凍土工況，確定計(jì)算區(qū)域和邊界條件，對(duì)管道周圍土體溫度場(chǎng)及應(yīng)力場(chǎng)進(jìn)行耦合分析，建立了含管溝參數(shù)的熱固耦合模型。研究土壤溫度變化、管溝邊坡坡度，溝底加寬裕量等對(duì)管道應(yīng)力的影響規(guī)律，揭示了影響管道安全性的關(guān)鍵因素，旨在為今后凍土區(qū)埋地管道設(shè)計(jì)提供理論基礎(chǔ)。

1 應(yīng)力模型

土體在凍結(jié)發(fā)生后依然會(huì)有未凍結(jié)的薄膜水，這時(shí)低溫吸力的增加和滲透率的降低會(huì)導(dǎo)致晶體的析出和冰透鏡的形成并侵入土體[11]，土壤體積增大引起地面膨脹。凍土物理學(xué)性質(zhì)會(huì)隨溫度而變化，凍土與大氣間的熱量平衡狀態(tài)會(huì)被輸油管道破壞，導(dǎo)致管道翹曲變形。輸油管道在凍土區(qū)運(yùn)行時(shí)，除受重力外，由于土壤溫度的變化還會(huì)受到管道周圍土壤的約束與壓迫而產(chǎn)生的應(yīng)力。

1.1 溫度場(chǎng)控制方程

凍土區(qū)埋地管道周圍溫度場(chǎng)分布隨時(shí)間變化而變化，其會(huì)受到管道油溫與土壤溫度的共同影響，溫度場(chǎng)的熱量平衡控制微分方程如下。

(1)在凍結(jié)區(qū)。

(1)

(2)在融化區(qū)。

(2)

(3)管壁。

(3)

式中：ρf、ρu、ρi分別為凍結(jié)區(qū)土壤、融化區(qū)土壤、管壁的密度；Cf、Cu、Ci分別為凍結(jié)區(qū)土壤、融化區(qū)土壤、管壁的熱容；λf、λu、λi分別為凍結(jié)區(qū)土壤、融化區(qū)土壤、管壁的導(dǎo)熱系數(shù)；Tf、Tu、Ti分別為凍結(jié)區(qū)土壤、融化區(qū)土壤、管壁的溫度；t為時(shí)間；r、θ為坐標(biāo)。

(4)熱固耦合方程為

(4)

(5)

式中：σx、σy為正應(yīng)力，MPa;τyx、τxy為切應(yīng)力，MPa；fx、fy為體力分量，MPa；x、y為坐標(biāo)方向。

假設(shè)初始時(shí)刻管道的溫度為Tg0,當(dāng)升溫至Tg時(shí)，管道會(huì)因?yàn)闇囟鹊淖兓l(fā)生膨脹而產(chǎn)生應(yīng)變變化，其計(jì)算公式為

εhot=-β(Tg-Tg0)

(6)

式(6)中：εhot為管線膨脹產(chǎn)生的熱應(yīng)變；Tg為管道溫度，℃；β為管道的線膨脹系數(shù)。

1.2 定解條件

初始條件:設(shè)Tw為地面溫度；Ta為地面溫度；T0為土壤初始溫度；距地面H深度的土壤溫度為TH，則大氣與土壤表面的對(duì)流邊界條件為

T|y=H=T0

(7)

(8)

式中：T為地面溫度，℃；λ為土壤導(dǎo)熱系數(shù)；α為土壤線膨脹系數(shù)；h為土壤與大氣間的對(duì)流換熱系數(shù)。

2 模型建立

2.1 管道基礎(chǔ)參數(shù)

建立X60管道模型，管道直徑813 mm,壁厚12.7 mm，密度7 800 kg/m3,埋深1.8 m(地表到管頂)屈服強(qiáng)度415 MPa，其他參數(shù)如表1所示。管線鋼材料特性由真實(shí)應(yīng)力-應(yīng)變關(guān)系來模擬，為了真實(shí)描述其塑性變形，由單向拉伸試驗(yàn)的名義應(yīng)力-應(yīng)變轉(zhuǎn)換為真實(shí)應(yīng)力-應(yīng)變，曲線如圖1[12]所示。

表1 管線鋼參數(shù)

圖1 X60管線鋼真實(shí)應(yīng)力-應(yīng)變曲線[11]

2.2 土壤基礎(chǔ)參數(shù)

凍土區(qū)土壤與普通土壤不同之處在于凍土區(qū)的土壤對(duì)溫度具有較高的敏感性，土壤的溫度會(huì)影響其彈性模量及泊松比。水分凍結(jié)成冰引起的體積膨脹是凍脹的主要原因[13]，凍土的凍脹量可以由土壤的線膨脹系數(shù)來確定，粉質(zhì)黏土與砂土的含水率為29.82%、12%，根據(jù)式(9)可得凍結(jié)土壤的線膨脹系數(shù)α1=-0.000 375、α2=-0.000 138[14]。

(9)

式(9)中：α為土壤線膨脹系數(shù)；ΔT為溫度變化量，℃；ω為土壤含水率。

土體凍融的導(dǎo)熱系數(shù)和比熱會(huì)由較大變化，熱參數(shù)的變化可以由不同溫度下的焓值來定義。焓與溫度的對(duì)應(yīng)關(guān)系如式(10)所示,焓值單位：J/m3。管溝內(nèi)土壤用砂土換填，土壤參數(shù)如表2[15]所示。

表2 回填砂土土壤參數(shù)[15]

(10)

對(duì)時(shí)間進(jìn)行微分得

(11)

再將式(1)代入式(11)得

(12)

式中：H(T)為土壤焓值；ρ為土壤密度；c(T)為比熱容；λ(T)為熱傳導(dǎo)系數(shù)。

2.3 有限元模型

利用ABAQUS軟件建立含管溝的穿越凍土區(qū)輸油管道三維管土模型，利用Drucker-Prager彈塑性本構(gòu)模型模擬土壤，管溝剖面如圖2所示，管溝邊坡坡度根據(jù)土壤及其物理學(xué)性質(zhì)確定，包括土壤黏聚力，密度等。管溝溝底寬度為加寬裕量和管道直徑之和，加寬裕量一般不小于0.5 m。3D模型采用8節(jié)點(diǎn)縮減積分C3D8R熱固耦合六面體單元對(duì)管道及土壤模型進(jìn)行結(jié)構(gòu)離散，每個(gè)節(jié)點(diǎn)具有一個(gè)溫度自由度，網(wǎng)格互不重疊。考慮到管道周圍溫度梯度大離散過程中對(duì)管道及管溝土壤區(qū)域進(jìn)行網(wǎng)格加密，在管道壁厚方向劃分單元層數(shù)，遠(yuǎn)離管道及管溝的土壤網(wǎng)格劃分較粗略。管土模型最終的網(wǎng)格劃分方案網(wǎng)格單元數(shù)為28 522，節(jié)點(diǎn)數(shù)43 578。管-土間的相互作用由接觸面法定義，允許管道表面與土體表面發(fā)生接觸與分離，屬于有限滑動(dòng)接觸，摩擦系數(shù)為0.3。模型上邊界地表與大氣會(huì)以對(duì)流傳熱的方式進(jìn)行熱交換，對(duì)土壤表面施加對(duì)流換熱系數(shù)模擬真實(shí)情況，對(duì)流換熱系數(shù)取17.8 W/m2·K[16]。管道與周圍土壤為面-面接觸，熱交換形式主要為熱傳導(dǎo)，模型左右邊界距離管道足夠遠(yuǎn)，且管道水平方向的溫度傳播有限，因此將模型兩側(cè)設(shè)置為絕熱邊界。天然地表下地溫年變化深度以下的土壤溫度梯度且超出了輸油管道熱傳導(dǎo)范圍，認(rèn)為計(jì)算區(qū)域的下邊界為恒溫。管土作用模型如圖3所示。

K為加寬裕量，mm；D為管道外徑，mm；A為管溝坡度；α為管溝坡度角，(°)

圖3 管土相互作用模型軸向方向

2.4 模型驗(yàn)證

中俄原油管線漠河-大慶段輸油工程穿越中國(guó)東北大興安嶺，該區(qū)域凍土類型較多，凍土層厚度平均在2～5 m，土壤巖性如表3[17]所示，管道的建設(shè)和運(yùn)行會(huì)影響周圍土壤的水熱性能，使其發(fā)生凍土凍脹現(xiàn)象[18]。為驗(yàn)證本文模型的可靠性，以該地區(qū)為例根據(jù)現(xiàn)有資料[19],選取凍土模型，尺寸為50 m(長(zhǎng))×10 m(寬)×10 m(高)，凍土層厚度5 m，凍脹區(qū)長(zhǎng)度10 m。區(qū)域管道埋深1.8 m，油壓為4 MPa，輸油溫度隨季節(jié)變化，管道垂直位移驗(yàn)證如圖4(a)所示。通過漠大地區(qū)地溫監(jiān)測(cè)數(shù)據(jù)來模擬管道周圍土壤溫度場(chǎng)變化[20]，地溫沿縱深變化驗(yàn)證曲線如圖4(b)所示。由圖4可以看出，本文管道垂直位移變化與文獻(xiàn)[3]中位移變化趨勢(shì)一致，模擬的土壤溫度變化與實(shí)測(cè)地下土壤溫度梯度相似，因此可以驗(yàn)證本文有限元模擬結(jié)果的可靠性。

表3 土壤熱物理性質(zhì)[17]

圖4 模型驗(yàn)證對(duì)比

3 結(jié)果與分析

為了獲得輸油管道在寒季運(yùn)行對(duì)其周圍土體的熱力影響規(guī)律，利用有限元模型模擬含管溝土壤的凍脹，對(duì)埋地輸油管道進(jìn)行了熱固耦合，分析了管溝參數(shù)、土壤溫度場(chǎng)多種因素對(duì)管道的應(yīng)力影響。有限元應(yīng)力云圖如圖5所示。

圖5 有限元應(yīng)力云圖

3.1 溫度對(duì)管道應(yīng)力的影響

保持管溝邊坡坡度(高∶寬)A為1，溝底加寬裕量K為1 m，根據(jù)地溫觀測(cè)數(shù)據(jù)研究地表溫度對(duì)管道應(yīng)力的影響。取地表溫度較低的1、2、3、11、12月，分別得到圖6所示的管頂和管底沿軸向Von Mises應(yīng)力隨溫度變化的曲線。

圖6 溫度影響下沿軸向Von Mises應(yīng)力

由圖6可知，土體發(fā)生凍脹會(huì)在管道底部和頂部出現(xiàn)應(yīng)力集中現(xiàn)象，在管道頂部有3個(gè)高風(fēng)險(xiǎn)失效區(qū)域，分別處于凍脹區(qū)域及非凍脹區(qū)與過渡區(qū)的交界處，這是因?yàn)橥寥腊l(fā)生凍脹后產(chǎn)生垂直位移，管土相互作用力會(huì)使管道受力變形，管底承受壓力，管頂承受拉力。管道頂部最大Von Mises應(yīng)力位于非凍脹區(qū)與過渡區(qū)的交界處，隨地表溫度的降低而增大，即在1月份管道Von Mises應(yīng)力達(dá)到峰值264.09 MPa。管道底部由于土壤的凍脹作用會(huì)在凍脹區(qū)域和非凍脹區(qū)與過渡區(qū)的交界處產(chǎn)生3個(gè)峰值點(diǎn)，最大Von Mises應(yīng)力位于凍脹區(qū)域中部，隨地表溫度的降低而增大，即在1月份管道Von Mises應(yīng)力達(dá)到峰值368.98 MPa。

圖7(a)中顯示了當(dāng)溝底加寬裕量K=1 m，管溝坡度A為1、1.5和2時(shí)，管道頂部與底部在1、2、3、11、12月份中受到的最大Von Mises應(yīng)力。在管道頂部，當(dāng)A=2時(shí)，管道Von Mises應(yīng)力受溫度影響最大，管道受溫度影響變化量為85.38 MPa，當(dāng)A=1時(shí)，管道Von Mises應(yīng)力受溫度影響最小，管道受溫度影響變化量為56.84 MPa。在管道底部，當(dāng)A=1時(shí)，管道Von Mises應(yīng)力受溫度影響最大，管道受溫度影響變化量為112.17 MPa，當(dāng)A=2時(shí)，管道Von Mises應(yīng)力受溫度影響最小，管道受溫度影響變化量為76.52 MPa。

圖7 溫度影響下最大Von Mises應(yīng)力

圖7(b)中顯示了當(dāng)管溝坡度A=1一定，溝底加寬裕量K為1、1.5、2 m時(shí)，管道頂部與底部在1、2、3、11、12月份中受到的最大Von Mises應(yīng)力。在管道頂部，當(dāng)K=1 m時(shí)，管道Von Mises應(yīng)力受溫度影響最大，管道受溫度影響變化量為63.48 MPa，當(dāng)K=2 m時(shí)，管道Von Mises應(yīng)力受溫度影響最小，管道受溫度影響變化量為51.06 MPa。在管道底部，當(dāng)K=1 m時(shí)，管道Von Mises應(yīng)力受溫度影響最大，管道受溫度影響變化量為73.43 MPa，當(dāng)K=2 m時(shí)，管道Von Mises應(yīng)力受溫度影響最小，管道受溫度影響變化量為50.11 MPa。

3.2 管溝坡度對(duì)管道應(yīng)力的影響

在最冷月并保持溝底加寬裕量K=1 m不變的情況下，管溝邊坡坡度(高∶寬)A取1、1.5、和2，分別得到了管頂與管底Von Mises應(yīng)力沿軸向分布的曲線，如圖8所示。

圖8 管溝坡度影響下沿軸向Von Mises應(yīng)力

由圖8可知，土體發(fā)生凍脹會(huì)在管道底部和頂部出現(xiàn)應(yīng)力集中現(xiàn)象，在管道頂部有3個(gè)高風(fēng)險(xiǎn)失效區(qū)域，分別處于凍脹區(qū)域及非凍脹區(qū)與過渡區(qū)的交界。在管道底部，由于土體對(duì)管道的壓力作用，管道在凍脹區(qū)域中部出現(xiàn)應(yīng)力集中現(xiàn)象，形成高風(fēng)險(xiǎn)失效區(qū)域。在非凍脹區(qū)與過渡區(qū)的交界處同樣出現(xiàn)應(yīng)力集中，但隨著管溝坡度的增大，應(yīng)力集中的現(xiàn)象得到了緩解。在溝底加寬裕量一定時(shí)，管頂凍脹區(qū)域中部最大Von Mises應(yīng)力隨管溝坡度的增大而增大，管溝坡度為2時(shí)，最大Von Mises應(yīng)力為285.46 MPa。管頂非凍脹區(qū)與過渡區(qū)交界處的最大Von Mises應(yīng)力隨管溝坡度的增大而減小，管溝坡度為1時(shí)，最大Von Mises應(yīng)力為264.09 MPa。管底凍脹區(qū)域中部最大Von Mises應(yīng)力隨坡度的增大而增大，管溝坡度為2時(shí)最大Von Mises應(yīng)力為372.07 MPa。當(dāng)管溝坡度為1時(shí)，在非凍脹區(qū)與過渡區(qū)的交界處出現(xiàn)應(yīng)力集中現(xiàn)象，此時(shí)Von Mises應(yīng)力為185.99 MPa。

由圖9可知，管道底部同一溝底加寬裕量下，管溝坡度為1.5時(shí)的最大Von Mises應(yīng)力相比管溝坡度為1時(shí)的最大增量為0.34%，最小增量為0.04%。管溝坡度為2時(shí)的最大Von Mises應(yīng)力相比管溝坡度為1.5時(shí)的最大增量為0.79%，最小增量為0.26%。管道頂部同一溝底加寬裕量下，管溝坡度為1.5時(shí)的最大Von Mises應(yīng)力相比管溝坡度為1時(shí)的最大增量為6.20%，最小增量為2.81%。管溝坡度為2時(shí)的最大Von Mises應(yīng)力相比管溝坡度為1.5時(shí)的最大增量為15.97%，最小增量為8.47%。

圖9 管溝坡度影響下最大Von Mises應(yīng)力

3.3 管溝加寬裕量對(duì)管道應(yīng)力的影響

在最冷月保持管溝坡度比1∶1不變，加寬裕量K取1、1.5、2 m，分別得到了管頂和管底的Von Mises應(yīng)力沿軸向分布的曲線，如圖10所示。

圖10 管溝加寬裕量影響下沿軸向Von Mises應(yīng)力

由圖10可知，管頂與管底會(huì)因凍脹產(chǎn)生應(yīng)力集中，分別處于凍脹區(qū)域及非凍脹區(qū)與過渡區(qū)的交界處，最大Von Mises 應(yīng)力出現(xiàn)在管底凍脹區(qū)域中部。當(dāng)管溝坡度比一定時(shí)，最大應(yīng)力隨溝底加寬裕量的增加而增大，凍脹區(qū)內(nèi)管道最大應(yīng)力為372.8 MPa，位于管底凍脹區(qū)中部；非凍脹區(qū)與過渡區(qū)的交界處也有應(yīng)力集中出現(xiàn)，最大Von Mises應(yīng)力為192.4 MPa。管道頂部最大應(yīng)力為277.7 MPa，位于非凍脹區(qū)與過渡區(qū)的交界，凍脹區(qū)內(nèi)管道最大應(yīng)力為272.4 MPa，位于凍脹區(qū)域中部。

由圖11可知，同一管溝坡度下，管道底部加寬裕量為1.5 m時(shí)的最大Von Mises應(yīng)力相比加寬裕量為1 m時(shí)的最大增量為1.15%，最小增量為0.73%。加寬裕量為2 m時(shí)的最大Von Mises應(yīng)力相比加寬裕量為1.5時(shí)的最大增量為0.18%，最小增量為0.05%。管道頂部加寬裕量為1.5 m時(shí)的最大Von Mises應(yīng)力相比加寬裕量為1 m時(shí)的最大增量為8.72%，最小增量為2.84。加寬裕量為2 m時(shí)的最大Von Mises應(yīng)力相比加寬裕量為1.5 m時(shí)的最大增量為5.46%，最小增量為0.94%。

圖11 管溝加寬裕量影響下最大Von Mises應(yīng)力

4 結(jié)論

(1)管溝坡度一定時(shí)，管道頂部由于凍脹會(huì)出現(xiàn)3個(gè)峰值，分別位于凍脹區(qū)及非凍脹區(qū)與過渡區(qū)的交界處，且隨溝底加寬裕量的增加而增大，最大Von Mises應(yīng)力出現(xiàn)在非凍脹區(qū)與過渡區(qū)的交界處。管道底部會(huì)出現(xiàn)一個(gè)高風(fēng)險(xiǎn)失效區(qū)，位于凍脹區(qū)域中部，在K=2 m時(shí)達(dá)到最值。

(2)管溝溝底加寬裕量一定時(shí)，管道頂部會(huì)因凍脹的發(fā)生出現(xiàn)3個(gè)峰值，分別位于凍脹區(qū)及非凍脹區(qū)與過渡區(qū)的交界處。在凍脹區(qū)域中部，最大Von Mises應(yīng)力隨管溝坡度的增大而增大。在非凍脹區(qū)與過渡區(qū)的交界處，最大Von Mises應(yīng)力隨管溝坡度的增大而降低。管道底部由于凍脹的發(fā)生出現(xiàn)高風(fēng)險(xiǎn)失效區(qū)，位于凍脹區(qū)域中部，且隨溝底加寬裕量的增加而增大。在非凍脹區(qū)與過渡區(qū)的交界處會(huì)出現(xiàn)峰值，但隨溝底加寬裕量的增加，應(yīng)力集中現(xiàn)象會(huì)緩解。

(3)管道穿越凍土區(qū)時(shí)，由于土體在寒季的凍脹現(xiàn)象受地表溫度影響較大，管道受土體凍脹的影響在管頂和管頂出現(xiàn)應(yīng)力集中的現(xiàn)象，分別位于凍脹區(qū)域中部及非凍脹區(qū)與過渡區(qū)的交界處。管道所受應(yīng)力隨地表溫度的降低而增大，在最冷月達(dá)到最值。當(dāng)溝底加寬裕量一定時(shí)，管溝坡度的增加會(huì)增大溫度對(duì)管道的受力影響。當(dāng)管溝坡度一定時(shí)，溝底加寬裕量的減小會(huì)增大溫度對(duì)含管道的受力影響。