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

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

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

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

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

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

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

動(dòng)網(wǎng)格；UDF；行星運(yùn)動(dòng)；Fluent；偏心率

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

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

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

1.1 模型簡化

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

圖1 內(nèi)桿行星運(yùn)動(dòng)環(huán)空流場(chǎng)簡化模型

1.2 動(dòng)網(wǎng)格方法

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

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

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

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

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

第n+1時(shí)間層上的體積Vn+1為：

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

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

每個(gè)控制容積面上的點(diǎn)積為：

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

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

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

UDF是為滿足計(jì)算個(gè)性化需求，可以動(dòng)態(tài)地連接到Fluent求解器上的用戶自編的程序。通過UDF調(diào)用動(dòng)網(wǎng)格宏(DEFINE_CG_MOTION)可實(shí)現(xiàn)剛體運(yùn)動(dòng)邊界的控制。內(nèi)桿運(yùn)動(dòng)UDF編譯的主體思想如模型簡化圖。自轉(zhuǎn)通過DEFINE_CG_MOTION宏中的omega[z]實(shí)現(xiàn)，公轉(zhuǎn)通過控制內(nèi)桿旋轉(zhuǎn)速度矢量vel[x]、vel[y]以及動(dòng)網(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)桿行星運(yùn)動(dòng)的環(huán)空流場(chǎng)數(shù)值計(jì)算對(duì)控制方程的離散應(yīng)用有限體積法，采用分離式求解器求解。選擇k-epsilion 湍流模型，并開啟能量方程，用中心差分格式離散擴(kuò)散項(xiàng)，選擇PRESTO!格式對(duì)壓力項(xiàng)離散，其余用一階迎風(fēng)離散，壓力速度耦合方程采用PISO算法，使用編譯型UDF文件定義剛體速度，動(dòng)網(wǎng)格采用局部網(wǎng)格重構(gòu)方法，對(duì)內(nèi)桿行星運(yùn)動(dòng)的環(huán)空流場(chǎng)進(jìn)行非定常求解[6]。

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

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

3.1 偏心環(huán)空螺旋流與內(nèi)桿行星運(yùn)動(dòng)環(huán)空中流體的速度場(chǎng)計(jì)算

通過流場(chǎng)周向速度分布圖(見圖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)桿行星運(yùn)動(dòng)的環(huán)空中，由于內(nèi)桿在環(huán)空中的偏心公轉(zhuǎn)，被轉(zhuǎn)動(dòng)波及的流體區(qū)域更大，流體周向速度最大值保持在靠近內(nèi)管壁環(huán)空窄間距側(cè)。

圖3 流場(chǎng)周向速度分布圖

3.2 偏心環(huán)空螺旋流與內(nèi)桿行星運(yùn)動(dòng)環(huán)空中流體的壓力場(chǎng)計(jì)算

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

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

4 內(nèi)桿偏心程度對(duì)流場(chǎng)的影響

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

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

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

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

圖5 流場(chǎng)周向速度分布圖

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

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

圖6 流場(chǎng)湍流強(qiáng)度分布圖

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

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

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

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

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

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

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

5 結(jié)語

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

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

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

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

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

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

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

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

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

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

[7] 馮順新,李啟兵,符松. 內(nèi)管自轉(zhuǎn)且公轉(zhuǎn)時(shí)環(huán)管冪律流的受力分析[J]. 計(jì)算力學(xué)學(xué)報(bào),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ù)雜流體流動(dòng)與數(shù)值模擬等方面的研究。

2014-12-28