高溫高壓模擬井筒應力分析與評價

侯永強, 紀 斌, 賈光政, 高 涵

(東北石油大學 機械科學與工程學院, 黑龍江 大慶 163318)

0 引 言

高溫高壓模擬井筒是用于模擬油田井下高溫高壓環境的實驗裝置．模擬井筒這類超高壓容器在加溫加壓工作過程中會受到熱應力與壓應力的耦合作用,求解此類應力耦合作用所引起的強度問題,對指導機械、化工、航空航天和核反應堆工程等領域的超高壓容器設計校核有重要意義．

針對復雜溫度場下的熱力耦合問題,國內外學者開展了大量的理論研究．通過解析法求解溫度場分布,進而采用有限元法求解熱應力和壓應力耦合作用所引起的強度問題[1-6],計算過程中溫度場與應力場獨立計算,難以準確描述耦合效應,計算準確度相對較低．Almasi等[7]和楊陽等[8]使用解析和數值方法對熱-力行為進行了研究．Manthena等[9]以Bessel函數的形式求解了瞬態二維傳導方程及其相關熱應力的解．耦合傳熱使得溫度場分布復雜且非線性,難以用解析法準確求解,導致熱應力無法準確求解．

模擬井筒加溫加壓實驗時,內部充滿液體并密封,為湍流自然對流傳熱過程．針對湍流求解過程中直接求解法(DNS)求解實際問題過程中計算量巨大,Reynolds時均方程法(RANS)求解精度不高的問題,Smagorinsky提出了大渦模擬(LES)方法．經過幾十年的發展,大渦模擬已經被實驗證明能夠較準確地求解湍流自然對流傳熱[10-14],降低對計算機資源的苛刻要求,與直接求解法相比,減少了計算的工作量．

由于湍流自然對流傳熱數學模型的N-S方程的非線性特性,它的求解十分困難,需要通過數值方法進行求解．SIMPLE系列算法求解非穩態N-S方程計算量大,效率不高．文獻[15-17]提出了投影法求解不可壓縮黏性流體的N-S方程,以提高數值求解的效率．

綜上所述,需要研究基于耦合傳熱溫度場,以求解計算熱應力的數學模型與數值求解方法．本文在此基礎上準確求解模擬了井筒的熱應力、壓應力和耦合應力的分布規律,為模擬井筒這類高溫高壓容器的設計校核提供了理論方法．

1 模擬井筒實驗裝置

圖1為模擬井筒加溫裝置,其中:①為風機,②為進風閥,③為內循環閥,④為冷卻閥,⑤為風循環管道,⑥為井式加熱爐,⑦、?為溫度傳感器,⑧為模擬井筒,⑨為加熱電阻絲,⑩為爐襯．加溫裝置采用空氣間接加熱的包覆式井式電加熱爐,為模擬井筒提供均勻穩定的高溫試驗環境．實驗裝置工作時,井式加熱爐對模擬井筒外壁進行加熱,超高壓加壓系統對模擬井筒內部加壓,模擬井筒承受高溫和高壓作用．

圖1 加溫裝置結構Fig. 1 The heating device structure

2 模擬井筒力學與熱學物理模型

模擬井筒外半徑rb=0.355 m,內半徑ra=0.175 m,高度為3.534 m,腔高度為2.50 m,rb/ra≈2.03(當rb/ra>1.2時為厚壁圓筒),模擬井筒為軸對稱厚壁圓筒．模擬井筒在理想狀態下材質均勻、形狀規整,且在三維空間內的形變也是對稱且統一的,可以不考慮其軸向彎曲．在柱坐標系下采用位移法求解其軸對稱非穩態溫度分布下的金屬熱應力．模擬井筒金屬厚壁材料為均勻的各向同性彈性體,加熱過程中變形很小,發生變形的速率很慢,因而慣性力可以忽略,可以應用線性熱彈性理論分析加熱過程的熱應力．加溫加壓過程中模擬井筒內部充滿流體,井筒上下端面采取保溫措施,可認為是絕熱面,豎直外壁面為高溫面．模擬井筒加熱溫度場徑向對稱,同時內部流體在重力場的作用下產生浮升力而形成自然對流,模擬井筒的導熱模型可用二維模型來描述．圖2為非定常高溫高壓模擬井筒熱學與力學物理模型．

圖2 非定常高溫高壓模擬井筒熱學與力學物理模型Fig. 2 The unsteady thermal and mechanical physical model diagram for the high-temperature high-pressure wellbore hole simulator

2.1 熱應力物理方程

在柱坐標系下模擬井筒厚壁的溫度分布周向對稱,溫度的變化與軸向角度θ無關,即溫度分布是只與半徑r、高度z和時間t相關的函數T=T(r,z,t)．于是,周向位移uθ=0、徑向位移ur和軸向位移uz都只是半徑r和高度z的函數．由于溫度T的變化,模擬井筒金屬厚壁內各點的微小形變如果不受到約束,將發生線應變βT,β為模擬井筒金屬厚壁的線膨脹系數．金屬厚壁為各向同性體,線膨脹系數β不隨方向變化,因此這種線應變在所有方向均相同,因而這個過程中無切應變．模擬井筒金屬厚壁受到外在約束和其體內各部分之間相互約束,線膨脹形變并不能夠自由發生,因此產生熱應力,同時由于金屬厚壁的彈性熱應力將引起附加的形變．在變溫的情況下彈性體的應變由兩部分疊加而成: ① 熱應變分量βT; ② 在熱膨脹時由于彈性體內各部分之間的相互約束而引起的力學應變,它們和熱應力之間服從Hooke定律．根據Duhamel-Neumann法則,考慮熱應變的應力和應變在柱坐標下的關系為

(1)

(2)

(3)

γrθ=0,

(4)

γθz=0,

(5)

(6)

根據式(1)—(6)得出柱坐標下模擬井筒厚壁的熱應力為

(7)

(8)

(9)

式中,σr為模擬井筒徑向熱應力,MPa;σθ為模擬井筒周向熱應力,MPa;σz為模擬井筒軸向熱應力,MPa;r為模擬井筒任意點的半徑,m．

2.2 壓應力物理方程

模擬井筒內部充滿高壓流體,圓筒形金屬厚壁受到對稱于中心軸的均勻內壓與均勻外壓(大氣壓)作用,其受力可以等效為平面軸對稱問題．采用位移法求解,取位移分量u為基本未知函數,利用只包含應力分量的微分方程和邊界條件求出位移分量,再利用幾何方程求出應變分量,從而應用物理方程求出應力分量為

(10)

(11)

式中,pi為模擬井筒承受的內壓,MPa;p0為模擬井筒承受的外壓,MPa;Ri為模擬井筒承受的內徑,m;R0為模擬井筒承受的外徑,m;pr為模擬井筒徑向壓應力,MPa;σpθ為模擬井筒周向壓應力,MPa．

2.3 模擬井筒大渦模擬溫度場控制方程

大渦模擬通過將柱坐標系下模擬井筒耦合傳熱N-S方程進行濾波,即將變量劃分成大尺度變量與小尺度變量．濾波后流體與固體區域統一的二維非穩態模擬井筒耦合傳熱控制方程為

(12)

(13)

(14)

(15)

μeff=μ+μt,

(16)

λeff的計算公式為

λeff=λ+ρcpαt,

(17)

μ為初始設定的流場渦黏度值,Pa·s;μt為亞格子渦黏系數,Pa·s;λ為初始設定的導熱系數,W/(m·K);αt為亞格子渦擴散系數,m2/s．

大渦模擬亞格子渦黏系數采用Smagorinsky-Lilly模型,計算公式為

(18)

式中,Ls為亞格子混合長度,其定義式為

Ls=min(kd,Cs,Δ),

(19)

(20)

3 數值求解方法

3.1 溫度場控制方程的離散

應用有限差分法對模擬井筒大渦模擬溫度場的控制方程式(12)—(15)進行離散,控制方程采用交錯網格進行離散,黏性項采用Crank-Nicholson 格式離散,對流項采用Adams-Bashforth格式離散得

(21)

(22)

(23)

(24)

以上各式中差分算子的定義為

Dr(*)i, j=((*)i+1/2, j- (*)i-1/2, j)/Δr,Dz(*)i, j=((*)i+1/2, j- (*)i-1/2, j)/Δz,

(*)i+1/2, j=((*)i+1, j+ (*)i, j)/2, (*)i, j+1/2=((*)i, j+1+(*)i, j)/2,

Lh(*)i, j= ((*)i+1, j-2(*)i, j+ (*)i-1, j)/Δr2+((*)i, j+1- 2(*)i, j+ (*)i, j-1)/Δz2,

Lm(*)i, j= ((*)i+1, j- (*)i-1, j)/(2Δr),

(*)表示離散變量,i,j分別表示r,z方向的單位向量．

3.2 投影法求解溫度場控制方程

投影法求解溫度場物理方程步驟如下:

(25)

(26)

把式(23)分解為下列兩式:

(27)

(28)

第二步 壓力修正步．將離散的動量方程式(26)和(28)代入離散的連續性方程(21)中,該方法可保證連續性方程嚴格滿足．經整理得到離散的壓力Poisson方程為

(29)

式中

(30)

式中,m為迭代次數;ω為松弛因子,ω>1時為超松弛迭代．

3.3 虛擬密度法求解流固耦合傳熱

3.4 非定常熱應力數值求解方法

(31)

(32)

(33)

式(31)—(33)中的積分項采用梯形法數值積分求解．

(34)

(35)

(36)

3.5 壓應力數值求解方法

壓應力數值求解,采用與求解熱應力相同的離散化節點劃分方法,即離散化為相同的i個半徑長度ri．壓應力與熱應力在空間維度有對應的離散單元,實現壓應力與熱應力的精確耦合計算．模擬井筒加溫加壓實驗過程中壓力保持恒定,為定常壓應力,因此在時間維度不需要進行離散化．空間離散化的徑向壓應力和周向壓應力為

(37)

(38)

4 應力的耦合與強度準則

根據應力疊加原理,應力耦合為熱應力與壓應力的和,熱應力與壓應力在半徑r方向按照相同的網格離散方式進行離散,其離散化的非定常應力耦合計算公式為

(39)

(40)

(41)

根據von Mises第四強度理論,高溫高壓模擬井筒的等效應力為

(42)

(43)

(44)

5 求解結果及強度分析

5.1 力學參數及邊界條件

模擬井筒為高溫高壓厚壁容器,由PcrNi3MoVA Ⅳ材料制造,表1為PcrNi3MoVA Ⅳ的力學參數．

模擬井筒加溫加壓過程的邊界條件為:初始溫度t0=25 ℃,上、下表面為絕熱邊界條件,豎直外壁面為高溫面,加熱溫度400 ℃,中心軸為軸對稱邊界條件,邊界位置如圖2所示．壓力修正方程的邊界條件為Neumann邊界條件,即?p/?n=0(n表示外法線)．模擬井筒放置在空氣中,外部壓力為大氣壓po=0.1 MPa,內部充滿水,內部工作壓力上限為220 MPa．根據建立的模擬井筒耦合傳熱,熱應力和壓應力數學模型及數值求解方法,編寫計算程序進行數值求解．

表1 PcrNi3MoVA Ⅳ的力學參數

圖3為模擬井筒中心高度1.77 m,不同時刻徑向耦合溫度分布．在加熱過程中井筒金屬厚壁部分的溫度梯度變化趨勢為先增加后降低．越靠近r=175 mm熱耦合壁面處溫度梯度大．沿徑向方向模擬井筒腔體內壁面溫度梯度較大．

5.2 熱應力參數分析

圖4—6為徑向、周向和軸向非定常熱應力在模擬井筒內壁面(r=175 mm)、模擬井筒厚壁內部中點(r=265 mm)和模擬井筒外壁面(r=355 mm)在不同時刻t的變化規律．

由圖4可知,模擬井筒內壁面(r=175 mm)的徑向非定常熱應力隨著加熱時間的增加緩慢增加,模擬井筒厚壁內部中點(r=265 mm)與外壁面(r=355 mm)位置,加熱初始時間段徑向熱應力較大,之后隨著加熱時間的增加熱應力減小,并趨于定值．

圖3 不同時刻徑向耦合溫度分布

圖5 周向非定常熱應力的時間分布

由圖5、6可知,周向與軸向的非定常熱應力變化規律相同,在模擬井筒內壁面(r=175 mm)與外壁面(r=355 mm)處加熱初始時間段周向與軸向熱應力較大,之后隨著加熱時間的增加周向與軸向熱應力減小,并趨于定值,模擬井筒厚壁內部中點(r=265 mm)處周向熱應力趨于零．

圖7—9為徑向、周向和軸向熱應力沿徑向在不同時刻的變化規律．由圖7—9可知,只考慮溫度載荷時,徑向熱應力相較于周向和軸向的熱應力小1個數量級;周向與軸向熱應力沿徑向由內到外按先正值后負值分布．

5.3 壓應力參數分析

圖10為只考慮模擬井筒內部流體壓力作用時徑向和周向壓應力沿徑向的變化規律．由圖10可知,只考慮內壓時,周向壓應力為正值(拉應力),徑向壓應力為負值(壓應力)．

圖7 徑向非定常熱應力沿半徑的分布

圖9 軸向非定常熱應力沿半徑的分布

5.4 等效應力求解結果參數分析

圖11、12為等效熱應力和等效應力耦合在厚壁井筒內壁面(r=175 mm)、厚壁井筒厚壁內部中點(r=265 mm)和厚壁井筒外壁面(r=355 mm)處在不同時刻t的變化規律．非定常等效熱應力與等效應力耦合隨著加熱時間的增長,其值先增加到最大值,之后逐漸減小．

圖11 非定常等效熱應力沿徑向的時間分布

5.5 模擬井筒的設計校核

圖13為加溫加壓條件下,模擬井筒不同壁厚與等效應力最大值的關系．隨著模擬井筒壁厚的增加,等效熱應力最大值增大,等效壓應力的最大值減小,等效應力耦合最大值基本保持不變．壁厚越大,模擬井筒的結構強度越大,因此最大等效壓應力值越小．但隨著壁厚的增加,模擬井筒在加熱過程中,其金屬厚壁的溫度梯度會增加,因此最大等效熱應力值增大．等效應力耦合為等效壓應力與等效熱應力的耦合作用結果, 在設計壁厚較大時(0.12～0.2 m)時, 壁厚發生變化時等效應力耦合最大值的變化率較小．在設計壁厚較小時(0.04～0.12 m),壁厚發生變化時等效應力耦合最大值的變化率較大．模擬井筒壁厚與最大等效壓應力、等效熱應力和等效應力耦合為非線性關系．因此在模擬井筒的設計過程中要綜合考慮等效熱應力,等效壓應力和等效應力耦合的大小,三者均需滿足強度要求．

圖13 模擬井筒壁厚與等效應力關系Fig. 13 The relationship between the wall thickness and the equivalent stress in the wellbore hole simulator

模擬井筒為超高壓容器,根據《超高壓容器》(GB/T34019—2017)的規定[18],理想彈塑性模型的屈服強度取材料在設計溫度下屈服強度的0.9倍,屈服強度與等效壓應力的比不小于2.20．設計模擬井筒的技術參數為:井筒壁厚為0.18 m,最高加熱溫度為400 ℃,最高工作壓力為220 MPa．應用建立的模擬井筒熱應力與壓應力數學模型與數值求解方法進行理論分析計算,計算結果為等效應力耦合計算結果為1 022 MPa,最大等效壓應力為477 MPa,等效熱應力為766 MPa．模擬井筒PcrNi3MoVA Ⅳ材料在設計溫度400 ℃下屈服強度為1 150 MPa,設計技術參數下等效應力耦合計算結果為1 022 MPa,小于0.9倍的屈服強度1 035 MPa,屈服強度與等效壓應力的比為2.30,均滿足標準規定．屈服強度與等效熱應力的比為1.50,滿足材料的強度要求．設計的最小壁厚為0.18 m,模擬井筒能夠在溫度400 ℃、工作壓力220 MPa的參數下安全工作．

6 模擬井筒加溫加壓實驗

通過模擬井筒加溫加壓實驗,對建立的熱力學數學模型與數值求解方法的數值求解結果的正確性進行實驗驗證．圖14為模擬井筒加溫裝置,通過安裝在井式電加熱爐爐膛內及模擬井筒內部的熱電偶傳感器,實時測試爐膛內與模擬井筒內部流體的溫度．圖15為加壓液控系統及加壓泵,通過安裝在模擬井筒加壓液控系統的壓力傳感器,實時測試加壓液體的壓力．

圖16為加溫加壓實驗曲線,在實驗過程中,模擬井筒通過井式電加熱爐進行加熱升溫,爐溫達到400 ℃后,保持爐溫400 ℃加熱．模擬井筒通過加壓液控系統進行升壓,模擬井筒內部壓力達到220 MPa后開始保壓,保壓約15 h后開始卸荷,并停止加熱降溫．保壓期間最高壓力為222.47 MPa,最低壓力為219.16 MPa．

模擬井筒安裝有溫度傳感器,圖1中的溫度傳感器?通過模擬井筒的下端接口插入其內部,測試模擬井筒腔高度為0.22 m,徑向位置為r=0.12 m處內部流體的溫度．實驗與計算升溫曲線如圖17所示,升溫過程的計算值與實驗測試值的最大誤差絕對值為4.79 ℃,計算結果與實驗測試數據基本一致．圖18給出了計算升溫曲線相對實驗升溫曲線的相對誤差,初始加熱時間段0～2 h內相對誤差最大,相對誤差絕對值最大值為4.18%．隨著加熱時間的增加,相對誤差不斷減小,相對誤差絕對值在2.5%以內．驗證了建立的模擬井筒大渦模擬溫度場數學模型及投影法求解流固耦合傳熱過程的準確性．

實驗結果表明,在最高400 ℃的加熱環境與內部加壓最高220 MPa的實驗參數條件下,設計最小壁厚為0.18 m的模擬井筒,能夠安全且無泄漏的工作．通過加溫加壓實驗,驗證了建立的模擬井筒熱應力與壓應力數學模型與數值求解方法的正確性．

圖14 模擬井筒加溫裝置Fig. 14 The wellbore hole simulator heating device

圖15 加壓液控系統及加壓泵Fig. 15 The pressurized hydraulic control system and the pressure pump

圖16 加溫加壓實驗曲線Fig. 16 Heating and pressurization experimental curves

圖17 實驗與計算升溫曲線

7 結 論

本文基于熱力學及大渦模擬理論,建立了模擬井筒大渦模擬溫度場物理方程．基于熱彈性力學理論建立了模擬井筒熱應力物理方程．給出了投影法數值求解溫度場控制方程的算法,梯形法數值積分求解熱應力控制方程的算法,給出了溫度場與應力控制方程的離散格式．通過虛擬密度法對流固耦合傳熱過程耦合求解,根據應力疊加原理對熱應力與壓應力耦合求解．對模擬井筒加溫加壓過程中的非穩態溫度分布、熱應力、壓應力及其耦合作用進行了分析,對模擬井筒進行了校核計算．在最高400 ℃的加熱環境下與內部加壓最高220 MPa的實驗參數條件下,設計最小壁厚為0.18 m的模擬井筒,能夠安全且無泄漏的工作．通過模擬井筒加溫加壓實驗,驗證了所建立的模擬井筒這類高溫高壓容器熱力學數學模型及數值求解結果的正確性,對指導高壓容器設計有重要意義．