龐博學(xué)，楊樹人，劉麗麗，劉 超，張 瑩，馬立輝

(東北石油大學(xué)，黑龍江 大慶 163000)

基于動網(wǎng)格與UDF技術(shù)的內(nèi)桿做行星運動的環(huán)空中流場數(shù)值模擬

龐博學(xué)，楊樹人，劉麗麗，劉 超，張 瑩，馬立輝

(東北石油大學(xué)，黑龍江 大慶 163000)

應(yīng)用動網(wǎng)格技術(shù)以及UDF程序成功實現(xiàn)了內(nèi)桿在油井偏心環(huán)空中自轉(zhuǎn)同時公轉(zhuǎn)的行星運動狀態(tài)。完成了內(nèi)桿在井筒偏心環(huán)空中不同偏心率位置轉(zhuǎn)動的情況下，流體在井筒環(huán)空中流動的數(shù)值模擬。計算結(jié)果表明，內(nèi)桿偏心率增加，流體周向運動增強，湍流高黏度區(qū)域擴大，層流區(qū)域相應(yīng)增加。存在一臨界偏心率εr，在其兩側(cè)，流體作用力分別將內(nèi)桿離心外推和向心內(nèi)拉。且內(nèi)桿偏心率越大，流體內(nèi)推力就越大，流體對內(nèi)桿公轉(zhuǎn)運動方向上的阻力也越大。

動網(wǎng)格；UDF；行星運動；Fluent；偏心率

對于內(nèi)桿靜止的偏心環(huán)空中的流動，前人做了大量的研究工作。與內(nèi)桿靜止的偏心環(huán)空流動相比較，當(dāng)內(nèi)桿自轉(zhuǎn)與公轉(zhuǎn)同時進(jìn)行時，流體在此種環(huán)空中的流動要更加復(fù)雜，人們對這種流動的研究少之又少[1]。在石油鉆井工程中，尤其是對于水平井和定向井，鉆桿因受自身重力的影響，在井筒中是彎曲的，再加上電動機驅(qū)動鉆桿的自轉(zhuǎn)，在井筒與鉆桿所形成的偏心環(huán)空中的鉆井液的流動就屬于內(nèi)桿做行星運動的環(huán)空中的流動問題[2]。與此相似，在采油工程尤其是螺桿泵采油中，由于內(nèi)桿偏心，抽油桿在自轉(zhuǎn)的同時還沿外管軸線公轉(zhuǎn)，在該環(huán)空中采出液的流動也屬于內(nèi)桿做行星運動的環(huán)空中的流動問題[3]。由此可見，此種流動在工程實際中常常遇到。

本文旨在應(yīng)用計算流體動力學(xué)的方法，通過編譯UDF程序?qū)崿F(xiàn)內(nèi)管行星運動狀態(tài)，應(yīng)用Fluent動網(wǎng)格技術(shù)對內(nèi)桿行星運動的環(huán)空中的流動問題進(jìn)行數(shù)值求解。此舉可以為定向井和水平井等工況鉆井液的水力參數(shù)設(shè)計，以及螺桿泵采油生產(chǎn)工況的優(yōu)化提供理論指導(dǎo)方向，具有學(xué)術(shù)價值和工程實踐意義[4]。

1 動網(wǎng)格技術(shù)在內(nèi)桿行星運動中的應(yīng)用

1.1 模型簡化

內(nèi)桿行星運動環(huán)空流場簡化模型如圖1所示。將石油工程實際中所描述的問題做如下簡化：流體處于內(nèi)、外2個偏心的圓管構(gòu)成的環(huán)空中，外圓管靜止，內(nèi)圓管以恒定的角速度ω自轉(zhuǎn)，又以恒定的角速度Ω沿外圓管軸線公轉(zhuǎn)，并有沿外圓管軸線方向的壓力P作用在流體上。

圖1 內(nèi)桿行星運動環(huán)空流場簡化模型

1.2 動網(wǎng)格方法

為實現(xiàn)環(huán)空中內(nèi)桿的行星運動，數(shù)值模擬應(yīng)用動網(wǎng)格技術(shù)，如圖2所示。內(nèi)桿運動使得流體幾何區(qū)域變化劇烈，動網(wǎng)格更新方法采用局部網(wǎng)格重構(gòu)模型(remeshing)，動區(qū)域選擇指定剛體運動(rigid body)。模擬計算過程中，因動邊界導(dǎo)致局部網(wǎng)格超過定義的最大尺寸時網(wǎng)格一分為二，局部網(wǎng)格受擠壓緊縮低于定義的最小尺寸時合二為一，當(dāng)網(wǎng)格扭曲程度超過定義的傾斜度時亦重新劃分，以適應(yīng)運動區(qū)域的形態(tài)變化。被重新劃分的網(wǎng)格單元緊鄰動網(wǎng)格節(jié)點，運動邊界附近網(wǎng)格得到及時更新[5]。

圖2 內(nèi)桿自轉(zhuǎn)及內(nèi)桿自轉(zhuǎn)且公轉(zhuǎn)動網(wǎng)格瞬態(tài)圖

動網(wǎng)格任意控制體V，通用標(biāo)量Φ的積分守恒方程為：

時間導(dǎo)數(shù)項可由一階向后差分項用如下公式表示：

式中，n和n+1是當(dāng)前時間及下一層時間。

第n+1時間層上的體積Vn+1為：

式中，dV/dt是控制體的時間導(dǎo)數(shù)。

為滿足網(wǎng)格守恒律，控制體的體積時間導(dǎo)數(shù)為：

每個控制容積面上的點積為：

式中，δVj是整個時間步Δt上控制容積面j膨脹引起的體積改變。

2 UDF的編譯及求解模型的選擇

2.1 內(nèi)桿行星運動的用戶自定義函數(shù)(UDF)的編譯

UDF是為滿足計算個性化需求，可以動態(tài)地連接到Fluent求解器上的用戶自編的程序。通過UDF調(diào)用動網(wǎng)格宏(DEFINE_CG_MOTION)可實現(xiàn)剛體運動邊界的控制。內(nèi)桿運動UDF編譯的主體思想如模型簡化圖。自轉(zhuǎn)通過DEFINE_CG_MOTION宏中的omega[z]實現(xiàn)，公轉(zhuǎn)通過控制內(nèi)桿旋轉(zhuǎn)速度矢量vel[x]、vel[y]以及動網(wǎng)格矢徑完成。

2.2 數(shù)值求解條件及模型選擇

定義流體密度ρ為0.8×103g/m3，流變指數(shù)n為0.6，外管內(nèi)徑R為40 mm，內(nèi)桿外徑r為20 mm，內(nèi)桿自轉(zhuǎn)速度ω為30 rad/s，公轉(zhuǎn)速度Ω與ω相同，為30 rad/s，環(huán)空流體軸向壓力梯度為300 Pa/m。

此次內(nèi)桿行星運動的環(huán)空流場數(shù)值計算對控制方程的離散應(yīng)用有限體積法，采用分離式求解器求解。選擇k-epsilion 湍流模型，并開啟能量方程，用中心差分格式離散擴散項，選擇PRESTO!格式對壓力項離散，其余用一階迎風(fēng)離散，壓力速度耦合方程采用PISO算法，使用編譯型UDF文件定義剛體速度，動網(wǎng)格采用局部網(wǎng)格重構(gòu)方法，對內(nèi)桿行星運動的環(huán)空流場進(jìn)行非定常求解[6]。

3 偏心環(huán)空螺旋流與內(nèi)桿行星運動的環(huán)空中的流動

偏心環(huán)空螺旋流是僅有內(nèi)桿自轉(zhuǎn)的環(huán)空中的流動，而內(nèi)桿行星運動的環(huán)空中的流動則是在偏心環(huán)空螺旋流的基礎(chǔ)上，增加內(nèi)桿以角速度Ω繞外管軸線公轉(zhuǎn)的情況。分別模擬計算偏心環(huán)空螺旋流與內(nèi)桿行星運動的環(huán)空中的流動有助于分析內(nèi)桿公轉(zhuǎn)對流場的影響。

3.1 偏心環(huán)空螺旋流與內(nèi)桿行星運動環(huán)空中流體的速度場計算

通過流場周向速度分布圖(見圖3)可以得出結(jié)論：在僅有內(nèi)桿自轉(zhuǎn)的偏心環(huán)空螺旋流中，由于黏性力的作用，內(nèi)桿的自轉(zhuǎn)促使其壁面附近區(qū)域流體產(chǎn)生周向速度，速度標(biāo)量值以內(nèi)桿軸線為圓心呈環(huán)形區(qū)域分布，且由內(nèi)桿壁面沿徑向逐漸減小。而在內(nèi)桿行星運動的環(huán)空中，由于內(nèi)桿在環(huán)空中的偏心公轉(zhuǎn)，被轉(zhuǎn)動波及的流體區(qū)域更大，流體周向速度最大值保持在靠近內(nèi)管壁環(huán)空窄間距側(cè)。

圖3 流場周向速度分布圖

3.2 偏心環(huán)空螺旋流與內(nèi)桿行星運動環(huán)空中流體的壓力場計算

通過流體動態(tài)壓力分布圖(見圖4)可以得出結(jié)論：在偏心環(huán)空螺旋流中，流體作用力峰值出現(xiàn)在內(nèi)桿壁表面，且內(nèi)桿周向受力基本均勻；而在內(nèi)桿行星運動的環(huán)空中，流體作用力峰值保持在內(nèi)桿公轉(zhuǎn)迎流面窄間距處，即流體有阻止內(nèi)桿公轉(zhuǎn)并將內(nèi)桿向環(huán)空內(nèi)側(cè)推回的趨勢。

圖4 流體動態(tài)壓力分布圖

4 內(nèi)桿偏心程度對流場的影響

在內(nèi)桿自轉(zhuǎn)且公轉(zhuǎn)的環(huán)空流動中，流體與內(nèi)桿存在相互作用力，內(nèi)桿在流體中轉(zhuǎn)動的同時伴有橫向擾動， 這種擾動使得內(nèi)桿公轉(zhuǎn)半徑即相對于外管的偏心距在時刻變化，流體對內(nèi)桿的徑向作用力與內(nèi)桿偏心程度相互作用。由于內(nèi)桿偏磨、失穩(wěn)等情況與偏心率密切相關(guān)，因此有必要針對不同偏心率對內(nèi)桿行星運動流場的影響計算分析。

定義內(nèi)桿偏心率ε=e/emax，其中e是內(nèi)桿相對于外管軸線的偏心距，emax是最大偏心距(即內(nèi)桿外壁接觸外管內(nèi)壁的完全偏心)，此模型中，emax=20 mm。分別模擬計算內(nèi)桿偏心率ε=0.25、ε=0.5以及ε=0.75時內(nèi)桿行星運動的環(huán)空流體場。

4.1 不同偏心率環(huán)空流體周向速度場分布

通過不同偏心率環(huán)空流場周向速度分布圖(見圖5)可以得出結(jié)論：由于內(nèi)桿在環(huán)空中的偏心公轉(zhuǎn)，隨著偏心距的增加，被轉(zhuǎn)動波及的流體區(qū)域擴大，流體周向運動增強，速度標(biāo)量值增大。

圖5 流場周向速度分布圖

4.2 不同偏心率環(huán)空流體湍流強度分布

通過不同偏心率環(huán)空流場湍流強度分布圖(見圖6)可以得出結(jié)論：內(nèi)桿偏心距對流場湍流強度影響明顯，湍流強度峰值集中出現(xiàn)在環(huán)空寬間距側(cè)，且隨著偏心距的增加，湍流強度減弱，層流區(qū)域增加。

圖6 流場湍流強度分布圖

4.3 不同偏心率環(huán)空流體高黏度區(qū)域分布

通過不同偏心率環(huán)空流場湍流黏度分布圖(見圖7)可以得出結(jié)論：由于內(nèi)桿在環(huán)空中旋轉(zhuǎn)導(dǎo)致漩渦，渦流擴散使得流體產(chǎn)生遠(yuǎn)大于分子水平的表觀黏度，即湍流黏度，使組分黏度增加，湍流黏度峰值區(qū)域集中出現(xiàn)在寬間距側(cè)，但繞內(nèi)桿轉(zhuǎn)動方向出現(xiàn)旋轉(zhuǎn)偏移。且偏心距越大，湍流高黏度區(qū)域越大。

圖7 流場湍流黏度區(qū)域分布圖

4.4 偏心率對內(nèi)桿在環(huán)空中所受流體力的影響

內(nèi)管在環(huán)空流體中轉(zhuǎn)動的過程中，流體作用力將影響內(nèi)管受力平衡，此作用力由流體壓強和周向剪應(yīng)力綜合作用產(chǎn)生[7]。將流體作用力F分解為沿公轉(zhuǎn)路線徑向上的流體力Fj和沿公轉(zhuǎn)路線切向上的流體力Fq。定義Fj正方向指向公轉(zhuǎn)軌道軸心，F(xiàn)q正方向與內(nèi)桿公轉(zhuǎn)運動速度矢量方向一致。分別模擬內(nèi)桿偏心率ε為0.05、0.20、0.35、0.50、0.65、0.80及0.95情況下行星運動時的流體場，計算得出內(nèi)管所受流體沿公轉(zhuǎn)路線徑向的作用力Fj和沿公轉(zhuǎn)路線切向的作用力Fq，繪制如圖8所示曲線。

圖8 內(nèi)桿在不同偏心率ε下所受流體作用力曲線

通過不同偏心率環(huán)空流場內(nèi)桿所受流體作用力曲線可以得出結(jié)論：內(nèi)桿在公轉(zhuǎn)軌道徑向上(Fj方向)所受流體作用力一般情況下不平衡。存在臨界偏心率εr(約為0.65)，使得當(dāng)ε<εr時，F(xiàn)j始終小于零，流體作用力有將內(nèi)桿向離心方向外推的趨勢；ε>εr時，F(xiàn)j始終大于零，流體作用力有將內(nèi)桿向公轉(zhuǎn)軌道軸心方向內(nèi)推的趨勢，且偏心率ε越大，內(nèi)推作用力越大。而在公轉(zhuǎn)軌道切向上(Fq方向)，內(nèi)桿所受到的流體作用力始終小于零，說明流體力一直在阻止內(nèi)桿沿公轉(zhuǎn)軌道的運動，內(nèi)桿偏心率越大，流體阻力越大。

5 結(jié)語

應(yīng)用Fluent動網(wǎng)格技術(shù)和UDF可以成功實現(xiàn)偏心環(huán)空中內(nèi)桿的行星運動。對內(nèi)桿在不同偏心率下行星運動的環(huán)空中的流動研究得出如下結(jié)論。

1)隨著內(nèi)桿偏心率ε的增加，被轉(zhuǎn)動波及的流體區(qū)域擴大，流體周向運動增強。

2)內(nèi)桿旋轉(zhuǎn)導(dǎo)致渦流擴散使得流體產(chǎn)生遠(yuǎn)大于分子水平的湍流黏度，其峰值區(qū)域在環(huán)空寬間距側(cè)，但繞內(nèi)桿轉(zhuǎn)動方向出現(xiàn)旋轉(zhuǎn)偏移。隨著偏心距的增加，湍流高黏度區(qū)域擴大，湍流強度減弱，層流區(qū)域增加。

3)存在臨界偏心率εr：當(dāng)ε<εr時，流體作用力將內(nèi)桿向離心方向外推；當(dāng)ε>εr時，流體作用力將內(nèi)桿向軸心內(nèi)推，偏心率越大，內(nèi)推作用力越大，流體對內(nèi)桿公轉(zhuǎn)運動方向上的阻力越大。

[1] 崔海清,張淑云,修德艷,等. 流體在內(nèi)管做行星運動的環(huán)空中流動的壓力梯度[J]. 中國石油大學(xué)學(xué)報：自然科學(xué)版,2008(6):76-78，87.

[2] 裴曉含. 冪律流體在內(nèi)管做行星運動的環(huán)空中流動時內(nèi)管壁的受力分析[D].大慶：大慶石油學(xué)院,2007.

[3] 季海軍. 冪律流體在內(nèi)管做行星運動的環(huán)空中的流動[D].大慶：大慶石油學(xué)院,2005.

[4] 楊樹人. 粘彈性流體偏心環(huán)空非定常流的數(shù)值計算[D].大慶：大慶石油大學(xué),2006.

[5] 張來平,鄧小剛,張涵信. 動網(wǎng)格生成技術(shù)及非定常計算方法進(jìn)展綜述[J]. 力學(xué)進(jìn)展,2010(4):424-447.

[6] 楊樹人,龐博學(xué),劉麗麗. 天然氣井應(yīng)用渦流工具排水采氣的流場分析[J]. 石油礦場機械,2014(10):13-16.

[7] 馮順新,李啟兵,符松. 內(nèi)管自轉(zhuǎn)且公轉(zhuǎn)時環(huán)管冪律流的受力分析[J]. 計算力學(xué)學(xué)報,2008(2):212-217.

責(zé)任編輯彭光宇

NumericalSimulationofFlowFieldinEccentricAnnulusWhoseInnerTubehasPlanetaryMotionbasedonDynamicMeshandUDF

PANG Boxue，YANG Shuren，LIU Lili，LIU Chao，ZHANG Ying，MA Lihui

(Northeast Petroleum University, Daqing 163000，China)

Movement of the inner tube rotation and revolution in eccentric annulus was achieved successfully by means of dynamic mesh and UDF. The flow in the annulus whose inner tube at different eccentricity rotates and revolves simultaneously was simulated. The results showed that: with the eccentricity increase, circumferential movement of fluid enhance, aera of high turbulent viscosity and laminar region expanded. There was a critical eccentricity on each side of which fluid force pushed inner tube outward and pulled it inward. The larger the eccentricity was, the greater fluid force that pulled inner tube inward was. Furthermore, resistance to the motion of the inner tube in the direction of revolution by fluid becomed larger.

dynamic mesh, UDF, planetary motion, Fluent; eccentricity

TH 311

：A

龐博學(xué)(1991-)，男，碩士研究生，主要從事復(fù)雜流體流動與數(shù)值模擬等方面的研究。

2014-12-28