高溫高壓模擬井筒流固耦合傳熱建模與數(shù)值求解

侯永強， 賈光政， 李鴻霏， 苗夏楠， 趙崇任

(東北石油大學(xué)機械科學(xué)與工程學(xué)院， 大慶 163318)

隨著油氣田勘探和開發(fā)井深的增加，井下的油層溫度和壓力也相應(yīng)增加，對射孔器材耐溫耐壓性能也提出了更高的要求[1]。研制的高壓模擬井筒加溫系統(tǒng)用于模擬井下高溫環(huán)境，對射孔器材進行耐高溫高壓性能檢測，為射孔器材性能研究和設(shè)計提供井下模擬試驗條件[2]。由于模擬井筒內(nèi)為超高壓流體，無法直接測量密封的模擬井筒中心位置的流體溫度，因此，研究并求解其升溫過程中溫度場的分布對射孔器材耐高溫高壓性能檢測具有重要意義。中外學(xué)者對腔體傳熱問題的研究大多簡化為給定內(nèi)部邊界條件的封閉腔體內(nèi)自然對流換熱過程，且針對方形腔體結(jié)構(gòu)進行建模與研究[3]，針對圓柱形厚壁井筒的厚壁與腔體內(nèi)流體的耦合換熱過程研究較少[4-5]，對模擬井筒流固耦合傳熱過程的數(shù)學(xué)模型與數(shù)值求解方法的研究更鮮有報道[6]。同時，自然對流換熱數(shù)學(xué)模型為非線性耦合問題，求解該數(shù)學(xué)模型的計算量巨大，因此高效數(shù)值求解自然對流換熱問題也是研究的重點[7-10]。針對上述問題進行模擬井筒熱流固耦合換熱研究，求解出模擬井筒的流固耦合傳熱徑向溫度分布規(guī)律，對模擬井筒加溫系統(tǒng)的設(shè)計與溫度控制有重要的工程意義。研究圓柱形厚壁井筒的厚壁與腔體內(nèi)流體的熱流固耦合換熱數(shù)值求解方法也具有重要的理論和科學(xué)意義。

1 模擬井筒耦合傳熱物理模型

超高壓模擬井筒加溫系統(tǒng)主要由井式加熱爐、模擬井筒、風(fēng)循環(huán)系統(tǒng)和支撐底座等組成，結(jié)構(gòu)示意如圖1所示。模擬井筒為內(nèi)部能夠承受射孔彈引爆產(chǎn)生的壓力沖擊的高強度空腔厚壁金屬圓柱體。

1為風(fēng)循環(huán)系統(tǒng)；2為工作介質(zhì)；3為空氣夾層；4為模擬井筒； 5為井式加熱爐；6、7、8、9、10為溫度傳感器；11為支撐底座； 12為溫度傳感器圖1 加溫裝置結(jié)構(gòu)示意圖Fig.1 Schematic diagram of the heating device

井式加熱爐與模擬井筒之間的傳熱通過風(fēng)循環(huán)系統(tǒng)的強制對流對模擬井筒外表面進行加熱；模擬井筒外表面受熱后熱量沿厚壁進行熱傳導(dǎo)到達內(nèi)壁，在內(nèi)壁處進行流固耦合傳熱；模擬井筒內(nèi)充滿水，模擬井筒選用低膨脹系數(shù)的PCrNi3MoVA Ⅳ鋼材制作，壁厚0.18 m，可承受230 MPa的壓力，加熱過程中認為模擬井筒體積增加量可以忽略不計。模擬井筒內(nèi)部水被加熱到200 ℃，加熱過程中產(chǎn)生100 MPa的壓力，為高溫高壓飽和水[11]。熱量驅(qū)動模擬井筒內(nèi)的流體(飽和水)形成自然對流，使內(nèi)部流體溫度逐漸升高。通過溫度閉環(huán)實時監(jiān)測，調(diào)節(jié)爐膛內(nèi)部的電熱元件加熱功率，控制模擬井筒外表面加熱溫度均勻。根據(jù)以上分析，加熱過程的物理模型包括圓柱型形壁容器的外壁被均勻加熱、圓柱形厚壁容器與腔體內(nèi)部流體的熱流固耦合換熱、內(nèi)部流體由于溫度不均勻而形成密度差在重力場中產(chǎn)生浮升力引起對流換熱。

模擬井筒為圓柱體，加熱溫度場徑向?qū)ΨQ，模擬井筒的導(dǎo)熱模型可以用二維模型來描述。由于井筒上下端采取保溫措施，可不考慮模擬井筒兩端散熱，柱坐標(biāo)系下模擬井筒的二維物理模型如圖2所示。圖2中，r坐標(biāo)軸為水平方向，r1、r2分別為井筒的外表面半徑與內(nèi)表面半徑，m；z坐標(biāo)軸為垂直方向，為物理模型的對稱軸；H為圓柱形封閉腔體的長度，m；水平壁面為絕熱面，豎直外壁面為高溫面，豎直內(nèi)壁面為流固耦合傳熱壁面。

2 模擬井筒耦合傳熱數(shù)學(xué)模型

2.1 模擬井筒傳熱控制方程

在柱坐標(biāo)系下，建立二維非穩(wěn)態(tài)模擬井筒傳熱控制方程。根據(jù)能量守恒定律和傅里葉定律建立模擬井筒(r2

(1)

式(1)中：TS為模擬井筒的溫度分布，K；λ1為模擬井筒的導(dǎo)熱系數(shù)，W/(m·K)；ρ1為模擬井筒的密度，kg/m3；c1為模擬井筒的比定壓熱容，J/(kg·K)；r為模擬井筒的半徑，m；z為模擬井筒的高度，m；t為時間變量，s。

模擬井筒從初始溫度(室溫25 ℃)開始加熱，其內(nèi)部流體(飽和水)的密度與溫度的函數(shù)關(guān)系在線性區(qū)域，采用Boussinesq假設(shè)，密度僅考慮動量方程中與浮升力有關(guān)的項，其余各項中的密度為常數(shù)。模擬井筒內(nèi)部流體為高溫高壓飽和水，其自然對流換熱(r

(2)

ρogβ(Tf-To)

(3)

(4)

(5)

式中：u、v分別對應(yīng)z、r方向的速度，m/s；t為時間，s；Tf為模擬井筒內(nèi)流體溫度，K；To為初始溫度，K；μ為模擬井筒內(nèi)流體動力黏度，Pa·s；p為模擬井筒內(nèi)流體壓力，Pa；ρo為初始溫度時模擬井筒內(nèi)流體密度，kg/m3；β為模擬井筒內(nèi)流體體積膨脹系數(shù)，K-1；a2為模擬井筒內(nèi)流體熱擴散率，m2/s；g為重力加速度，m2/s。

2.2 流固耦合傳熱的數(shù)學(xué)模型

模擬井筒內(nèi)壁與其內(nèi)部流體在耦合邊界上的溫度連續(xù)方程、熱流密度連續(xù)方程為[12]

Tfp=Tsp

(6)

qf=qs

(7)

式中：Tfp和Tsp分別為耦合界面上流體與固體的溫度；qf和qs分別為耦合界面上流體邊界處與固體邊界處的熱流密度。

3 控制方程的離散及數(shù)值求解方法

建立的模擬井筒耦合傳熱數(shù)學(xué)模型需要采用數(shù)值解法求解。模擬井筒厚壁與其內(nèi)部流體的流固耦合傳熱采用分區(qū)計算和邊界耦合方法求解。

3.1 控制方程的離散

3.1.1 模擬井筒厚壁區(qū)域控制方程的離散

圖3 圓柱軸對稱坐標(biāo)的網(wǎng)格系統(tǒng)Fig.3 Grid system with cylindrical axisymmetric coordinates

模擬井筒厚壁控制方程應(yīng)用有限差分法求解，采用內(nèi)節(jié)點法進行網(wǎng)格劃分。取一個弧度的中心角所包含的范圍作為研究對象，圓柱軸對稱坐標(biāo)的網(wǎng)格系統(tǒng)如圖3所示。節(jié)點P所代表的控制容積由界面N-W-S-E圍成；P、N、S、W、E表示所研究的節(jié)點及相鄰的4個節(jié)點；i和j分別對應(yīng)r和z方向的節(jié)點位置；r、z為空間坐標(biāo)，將模擬井筒厚壁計算區(qū)域半徑(r)劃分為M等份，高度z方向劃分為N等份，得到M×N個節(jié)點；半徑(r)、高度(z)、時間(t)的離散變量分別用r、z、t表示；n表示非穩(wěn)態(tài)的時間層，n+1表示時間間隔為t的下一時間層。δr、δz表示相鄰兩節(jié)點間的距離。r、z表示相鄰兩界面間的距離。采用時間項向前差分和空間項中心差分的顯示差分格式，將式(1)進行離散化得

(8)

3.1.2 內(nèi)部流體區(qū)域控制方程的離散

模擬井筒內(nèi)部流體區(qū)域采用交錯網(wǎng)格進行離散，流體區(qū)域半徑r劃分為B等份，高度z方向劃分為C等份，得到B×C個節(jié)點。空間與時間坐標(biāo)變量的表示方法與模擬井筒厚壁計算區(qū)域相同。應(yīng)用有限差分方法離散控制方程[式(2)～式(5)], 其中對流項用二階精度的Adamas-Bashforth 格式顯式處理，黏性項用Crank-Nicholson 格式離散，空間離散采用二階中心差分格式。應(yīng)用以上離散格式，將模擬井筒傳熱控制方程進行離散化得

(9)

(10)

(11)

(12)

差分算子的定義為[13]

Dr()()i-1/2,j]/Δr

(13)

Dx()()i+1/2,j-()i-1/2,j]/Δx

(14)

(15)

(16)

Lh()i,j=[()()i,j+()i-1,j]/Δr2+[()

(17)

Lm()i,j=[()()i-1,j]/2Δr

(18)

3.2 控制方程的數(shù)值求解

3.2.1 模擬井筒控制方程數(shù)值求解

令r=ir，z=jz，F(xiàn)0=1t/ρ1c1rz，模擬井筒厚壁區(qū)域的控制方程如式(8)所示，整理得

(19)

3.2.2 模擬井筒內(nèi)部流體區(qū)域的控制方程數(shù)值求解

流體區(qū)域采用投影法求解，模擬井筒內(nèi)部流體區(qū)域的控制方程為不可壓縮流動的Navier-Stokes方程組，是一個非線性很強的方程組，直接耦合數(shù)值求解該方程組計算量巨大[14]，常用的計算方法都是將速度和壓力解耦求解[15]。投影法最初由Chorin提出，其原理主要是Helmholtz-Hedge矢量分解定理[16]，任意矢量場可以分解為一個無散度矢量場和一個無旋矢量場之和。投影法把一個非線性系統(tǒng)化成一系列橢圓問題，解耦的同時也解決了非線性項的問題，解耦后的Navier-Stokes方程分成三步進行求解，求得速度和壓力值，免去了在各步之間反復(fù)迭代過程。因而能非常高效的求解模擬井筒內(nèi)部流體區(qū)域的控制方程。

采用Helmholtz-Hedge矢量分解方法，通過引入中間速度u*、v*，把式(10)分解為式(20)、式(21)：

(20)

(21)

把式(11)分解為式(22)、式(23)：

(22)

(23)

求解分解后的動量方程[式(20)、式(22)]，求解中間速度場u*、v*。

離散的壓力Poisson方程為

ap′i,j+bp′i+1,j+cp′i-1,j+dp′i,j+1+ep′i,j-1+f=0

(24)

(25)

式(25)中：m為迭代次數(shù)；ω為松弛因子，ω>1時為超松弛迭代。

由求得的u*、v*、pn+1，應(yīng)用式(21)、式(23)，求得n+1時刻的un+1、vn+1，應(yīng)用式(12)，求得下一時刻的溫度T。

3.3 流固耦合傳熱的數(shù)值求解

在計算過程中為了節(jié)省計算資源，提高計算速度，井筒厚壁控制區(qū)域與井筒內(nèi)流體區(qū)域分別滿足各自計算區(qū)域穩(wěn)定性的網(wǎng)格數(shù)。耦合界面處固體與液體的網(wǎng)格數(shù)不同的情況下，熱流密度采用線性插值方法進行數(shù)據(jù)的傳遞，保證空間上的數(shù)據(jù)傳遞正確。

熱流固耦合熱流密度連續(xù)方程采用熱平衡法求解。熱平衡法直接將能量守恒原理以及傅里葉導(dǎo)熱定律應(yīng)用于節(jié)點所代表的控制容積。熱流固耦合邊界坐標(biāo)網(wǎng)格系統(tǒng)如圖4所示。

圖4 熱流固耦合邊界坐標(biāo)網(wǎng)格系統(tǒng)Fig.4 Fluid-solid coupled heat transfer boundary coordinate grid system

圖4中對于節(jié)點P，在t時刻所代表的控制容積(由虛線N-W-S面圍成)建立熱平衡關(guān)系，各項熱流量都以進入單位體積元體P的方向為正。N、S、W節(jié)點分別對應(yīng)節(jié)點P相鄰的上部、下部和左邊節(jié)點。根據(jù)能量守恒定律，非穩(wěn)態(tài)節(jié)點的能量守恒表達形式為

(26)

從節(jié)點N，通過界面n傳導(dǎo)到節(jié)點P的熱流量可表示為

(27)

從節(jié)點S，通過界面s傳導(dǎo)到節(jié)點P的熱流量可表示為

(28)

從節(jié)點W，通過界面s傳導(dǎo)到節(jié)點P的熱流量可表示為

(29)

根據(jù)式(7)耦合邊界上熱流密度連續(xù)，節(jié)點P與內(nèi)部流體換熱的熱流量可表示為

(30)

(31)

內(nèi)邊界單位體積元體P′的能量守恒方程為

φN+φS+φW+φP=ΔEP′

(32)

將式(27)～式(31)代入式(32)，整理得

a′PT′P=aPTP+aNTN+aSTS+aWTW+afTf

(33)

3.4 邊界條件

模擬井筒耦合傳熱數(shù)學(xué)模型的邊界條件為：上、下表面z=0、z=H處絕熱，r=r1壁面為高溫面，r=0為對稱軸，r=r2壁面為流固耦合傳熱面，邊界位置如圖2所示。壓力修正方程的邊界條件為Neumann邊界條件，即p/n=0(n為外法線)[17]。其他邊界條件如表1所示。

表1 模擬井筒耦合傳熱邊界條件Table 1 Simulated wellbore coupling heat transfer boundary conditions

3.5 流固耦合傳熱求解

式(1)～式(7)組成了模擬井筒流固耦合非穩(wěn)態(tài)傳熱數(shù)學(xué)模型，采用流體與固體瞬態(tài)耦合算法進行求解。即在每個[t,t+t]時間步長內(nèi)，按照以下步驟進行求解。

(1)假定耦合邊界上的溫度分布，作為流體區(qū)域的邊界條件。

(2)對流體區(qū)域進行求解，得出耦合邊界上的局部熱流密度和溫度梯度，作為固體區(qū)域的邊界條件。

(3)按照式(7)求解固體區(qū)域，得出耦合邊界上新的溫度分布，再以此分別作為流體區(qū)域的邊界條件。重復(fù)上述計算步驟直到收斂。

4 分析及計算

模擬井筒及其內(nèi)部流體的物理參數(shù)為：模擬井筒外壁半徑r1=0.355 m；內(nèi)壁半徑r2=0.175 m；高度H=2.5 m。初始溫度T0=25 ℃；外壁加熱溫度400 ℃恒定；模擬井筒鉻鎳鋼的密度ρ1=7 830 kg/m3；比熱容C1=463 J/(kg·K)。其他物性參數(shù)，如模擬井筒鉻鎳鋼的導(dǎo)熱系數(shù)(1)、飽和水的導(dǎo)熱系數(shù)(2)、水的運動黏度ν、水的熱擴散率(a)均為對應(yīng)溫度下的物性參數(shù)。固體區(qū)域半徑r方向節(jié)點數(shù)為M=12；高度z方向節(jié)點數(shù)為N=157，時間步長dt1=1 s。傅里葉數(shù)F=0.051 6即(F<0.25)，滿足穩(wěn)定性條件。流體區(qū)域經(jīng)網(wǎng)格無關(guān)性驗證后取半徑r方向節(jié)點數(shù)為B=32；高度z方向節(jié)點數(shù)為C=450，時間步長dt2=0.001 s。半徑與高度方向CFL數(shù)均小于0.1，滿足流場計算的穩(wěn)定條件。

根據(jù)建立的模擬井筒耦合傳熱數(shù)學(xué)模型和求解方法，按照流固耦合傳熱求解步驟編寫計算程序進行數(shù)值求解。加熱的初始時間t<0.14 h，此過程中模擬井筒金屬厚壁的升溫速度遠大于其內(nèi)部水的升溫速度。模擬井筒金屬厚壁的比熱容C1和密度ρ1較大，對模擬井筒外壁恒溫加熱過程中，金屬厚壁需要吸收一定量的熱量，模擬井筒內(nèi)壁與外壁產(chǎn)生較大的溫度梯度。模擬井筒內(nèi)壁面溫度升高很慢，基本維持在初始溫度25 ℃。模擬井筒內(nèi)壁面與水之間無溫差，沒有進行流固熱耦合熱量交換。t=0.14 h時，模擬井筒內(nèi)壁面與水之間產(chǎn)生溫差，開始進行流固熱耦合熱量交換。t=0.14 h模擬井筒內(nèi)水的速度場如圖5(a)所示。

由圖5、圖6可知，加熱時間在0.14 h

加熱時間在1 h

流體的自然對流由圖5(a)變化到圖5(b)的狀態(tài)，水由沿井筒內(nèi)壁面的快速向上運動，發(fā)展到整個流體區(qū)域的環(huán)狀運動。因自然對流比較微弱，模擬井筒高徑比較大，在上部和下部形成兩個微弱環(huán)狀自然對流運動。

模擬井筒中心高度位置(高度H=1.1 m)處的徑向溫度分布計算結(jié)果如圖8所示，r=0處為模擬井筒徑向中心，r=175 mm為流固熱耦合處(模擬井筒內(nèi)壁面)，r=355 mm處為模擬井筒外壁面。圖8中，0≤r≤175 mm為模擬井筒內(nèi)部飽和水的溫度分布，175 mm

圖5 模擬井筒內(nèi)水的速度場圖Fig.5 Simulate the velocity field of water in the wellbore

圖6 流固耦合溫度場分布Fig.6 Flow-solid coupling temperature field distribution

圖7 流體流線圖Fig.7 Fluid flow diagram

圖8 流固耦合徑向溫度分布Fig.8 Fluid-solid coupling radial temperature distribution

5 實驗研究

為了驗證計算結(jié)果的正確性，采用模擬井筒加溫實驗裝置進行模擬井筒流固耦合傳熱升溫實驗。加溫實驗裝置為井式電加熱爐，采用晶閘管調(diào)節(jié)器對電加熱功率進行控制調(diào)節(jié)，實現(xiàn)溫度精確控制。圖9為模擬井筒加溫實驗裝置結(jié)構(gòu)。通過風(fēng)循環(huán)系統(tǒng)的空氣循環(huán)對模擬井筒外表面進行均勻加熱，主要技術(shù)如表2所示。

實驗裝置爐溫調(diào)整范圍為0～500 ℃，選用WFN型鎧裝熱電偶溫度傳感器，測量范圍為0～1 100 ℃，測試位置為厚壁井筒外壁高度方向均布4個。模擬井筒內(nèi)溫度范圍為0～300 ℃，選用WRN型鎧裝熱電偶溫度傳感器測量范圍是0～1 000 ℃，頂部熱電偶(圖1中的溫度傳感器6)通過模擬井筒的上端蓋插入其內(nèi)部，測試位置為模擬井筒內(nèi)液體的高度為H=2.28 m，徑向為r=0.1 m處；底部熱電偶(圖1中的溫度傳感器12)通過模擬井筒的下端插入其內(nèi)部，測試位置為高度H=0.22 m，徑向位置r=0.1 m處。

仿真計算與實驗裝置采用相同的物理參數(shù)，實測與仿真計算升溫曲線如圖10所示。仿真計算結(jié)果與實驗實測數(shù)據(jù)基本一致，都是初始加熱時間段升溫速度緩慢，之后溫度快速上升，頂部溫度略高于底部溫度，自然對流效應(yīng)很弱。圖11給出了實測與仿真計算升溫曲線的相對誤差，初始加熱時間段0～2 h相對誤差最大，頂部相對誤差最大值為18%，低部相對誤差最大值為15%。隨著加熱時間的不斷推進相對誤差不斷減小，頂部和低部相對誤差均在5%以內(nèi)。頂部溫度的計算值與實驗測試值的均方根誤差為4.66 ℃，低部溫度的計算值與實驗測試值的均方根誤差為5.19 ℃。實驗測試值與計算值的均方根誤差較小。通過實驗驗證了所建立的數(shù)學(xué)模型和數(shù)值求解方法求解超高壓模擬井筒流固耦合傳熱問題的正確性。

圖9 模擬井筒加溫實驗裝置Fig.9 Simulated wellbore heating experiment device

表2 加溫實驗裝置的主要技術(shù)參數(shù)Table 2 Main technical parameters of the heating experimental device

圖10 實測與仿真計算升溫曲線Fig.10 Measured and simulated temperature rise curve

圖11 仿真計算升溫曲線的相對誤差Fig.11 Simulation of the relative error of the heating curve

6 結(jié)論

研究了模擬井筒的厚壁腔體與腔體內(nèi)流體(飽和水)之間的動態(tài)耦合傳熱過程，得出以下結(jié)論。

(1)基于模擬井筒耦合傳熱過程的分析，建立了模擬井筒耦合傳熱數(shù)學(xué)模型；應(yīng)用有限差分法離散模擬井筒耦合傳熱數(shù)學(xué)模型，得出井筒厚壁和腔體內(nèi)流體(飽和水)數(shù)學(xué)模型的離散格式。

(2)采用投影法對腔體內(nèi)流體的非穩(wěn)態(tài)自然對流傳熱數(shù)學(xué)模型進行數(shù)值求解，提高了非線性耦合數(shù)學(xué)模型的計算效率和求解速度；采用迭代法分別計算井筒厚壁與腔體內(nèi)流體區(qū)域的傳熱過程，在厚壁與流體的交界處應(yīng)用熱平衡法進行耦合傳熱計算，實現(xiàn)了對模擬井筒加熱過程的求解。

(3)計算結(jié)果表明，模擬井筒厚壁升溫速度遠大于其內(nèi)部流體升溫速度，內(nèi)部流體具有邊界層型流動和傳熱的特點；模擬井筒整體的升溫速度主要受其內(nèi)部高溫高壓飽和水的熱傳導(dǎo)影響，從理論上解釋了厚壁模擬井筒加熱升溫速度緩慢的原因。

(4)實驗數(shù)據(jù)與仿真計算結(jié)果的對比分析，驗證了所建立的數(shù)學(xué)模型求解高溫高壓模擬井筒流固耦合傳熱過程的正確性。