侯永強, 紀(jì) 斌, 賈光政, 高 涵

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

0 引 言

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

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

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

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

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

1 模擬井筒實驗裝置

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

圖1 加溫裝置結(jié)構(gòu)Fig. 1 The heating device structure

2 模擬井筒力學(xué)與熱學(xué)物理模型

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

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

2.1 熱應(yīng)力物理方程

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

(1)

(2)

(3)

γrθ=0,

(4)

γθz=0,

(5)

(6)

根據(jù)式(1)—(6)得出柱坐標(biāo)下模擬井筒厚壁的熱應(yīng)力為

(7)

(8)

(9)

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

2.2 壓應(yīng)力物理方程

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

(10)

(11)

式中,pi為模擬井筒承受的內(nèi)壓,MPa;p0為模擬井筒承受的外壓,MPa;Ri為模擬井筒承受的內(nèi)徑,m;R0為模擬井筒承受的外徑,m;pr為模擬井筒徑向壓應(yīng)力,MPa;σpθ為模擬井筒周向壓應(yīng)力,MPa．

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

大渦模擬通過將柱坐標(biāo)系下模擬井筒耦合傳熱N-S方程進行濾波,即將變量劃分成大尺度變量與小尺度變量．濾波后流體與固體區(qū)域統(tǒng)一的二維非穩(wěn)態(tài)模擬井筒耦合傳熱控制方程為

(12)

(13)

(14)

(15)

μeff=μ+μt,

(16)

λeff的計算公式為

λeff=λ+ρcpαt,

(17)

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

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

(18)

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

Ls=min(kd,Cs,Δ),

(19)

(20)

3 數(shù)值求解方法

3.1 溫度場控制方程的離散

應(yīng)用有限差分法對模擬井筒大渦模擬溫度場的控制方程式(12)—(15)進行離散,控制方程采用交錯網(wǎng)格進行離散,黏性項采用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)代入離散的連續(xù)性方程(21)中,該方法可保證連續(xù)性方程嚴(yán)格滿足．經(jīng)整理得到離散的壓力Poisson方程為

(29)

式中

(30)

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

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

3.4 非定常熱應(yīng)力數(shù)值求解方法

(31)

(32)

(33)

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

(34)

(35)

(36)

3.5 壓應(yīng)力數(shù)值求解方法

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

(37)

(38)

4 應(yīng)力的耦合與強度準(zhǔn)則

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

(39)

(40)

(41)

根據(jù)von Mises第四強度理論,高溫高壓模擬井筒的等效應(yīng)力為

(42)

(43)

(44)

5 求解結(jié)果及強度分析

5.1 力學(xué)參數(shù)及邊界條件

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

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

表1 PcrNi3MoVA Ⅳ的力學(xué)參數(shù)

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

5.2 熱應(yīng)力參數(shù)分析

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

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

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

圖5 周向非定常熱應(yīng)力的時間分布

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

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

5.3 壓應(yīng)力參數(shù)分析

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

圖7 徑向非定常熱應(yīng)力沿半徑的分布

圖9 軸向非定常熱應(yīng)力沿半徑的分布

5.4 等效應(yīng)力求解結(jié)果參數(shù)分析

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

圖11 非定常等效熱應(yīng)力沿徑向的時間分布

5.5 模擬井筒的設(shè)計校核

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

圖13 模擬井筒壁厚與等效應(yīng)力關(guān)系Fig. 13 The relationship between the wall thickness and the equivalent stress in the wellbore hole simulator

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

6 模擬井筒加溫加壓實驗

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

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

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

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

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

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

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

圖17 實驗與計算升溫曲線

7 結(jié) 論

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