基于嵌套網格的環空流體內旋轉桿柱與井筒碰撞特性研究*

岳欠杯, 王笑笑, 曹 文, 劉躍秋, 李 輝, 徐燕璐

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

0 引 言

浸沒在流體中桿管間的接觸碰撞是石油工程最常見的現象,由于碰撞引起流體域網格拓撲結構的變化,采用界面解析直接數值模擬方法研究流體中桿管間的碰撞具有很大的挑戰,因此尚未得到廣泛的研究.

傳統的流固耦合數值方法通常是基于貼體網格,最典型的例子是任意 Lagrange-Euler(ALE)界面跟蹤方法[1],該方法已得到了廣泛的應用[2].在這種方法中,由于網格服從結構形狀,可以很容易地定義流固耦合邊界,此外,還可以對固體表面附近的網格進行局部細化,在耦合界面附近的數值模擬精度較高.該方法的缺點是在處理大變形和動態邊界問題時需要不斷更新網格,并在新舊網格上交換各種數據,這不僅增加了計算量,而且降低了求解的精度和穩定性.在具有復雜幾何形狀的三維空間中,這個過程會相當復雜,特別是當流體中存在固體間接觸和碰撞問題時,會出現負體積網格,導致計算終止.為了解決這一難題,國內外學者提出了嵌套網格方法[3],采用靈活的網格耦合方法和對邊界條件的處理,將整個流體域劃分為包含整個求解區域的背景區域和一個或多個包含固體的組件區域,在每個區域生成高質量的網格,然后進行網格裝配,即建立組件網格與背景網格之間的連接,實現各個網格間的數據交換,將背景網格(或組件網格)中供體單元的值通過插值傳遞給組件網格(或背景網格)中的受體單元.目前,該方法已經廣泛成功應用于物體繞流等流固耦合領域.Miller等[4]提出了一種適用于流固耦合分析的嵌套網格方法,并通過對基準測試的數值研究證明了嵌套網格方法相對于傳統方法的優越性.基于動態結構嵌套網格方法,李映坤[5]建立了一套多物理場耦合軟件平臺,對點火過程中的燃氣沖擊特性、裝藥傳熱特性、金屬膜片機械響應和隔層變形力學特性進行了深入系統的研究,為脈沖隔離裝置的設計、雙脈沖發動機的工程研制提供了重要參考.倪同兵[6]生成了由旋翼槳葉貼體網格和背景網格組成的適用于旋翼非定常流場數值模擬的結構運動嵌套網格,從而建立了一套適用于前飛狀態直升機彈性旋翼流場與氣動特性計算的高精度CFD/CSD耦合方法.徐廣[7]對槳葉彈性變形運動的旋翼運動嵌套網格生成方法、高效貢獻單元搜索策略、旋翼流場數值模擬和槳葉動力學分析方法等進行了研究,建立了基于Navier-Stokes方程的旋翼非定常流場的數值計算方法和程序.在固體動力學分析中,有限元法是分析物體運動、幾何非線性、材料非線性和接觸非線性[8]的主要方法.通過調研發現嵌套網格方法與有限元法相結合能有效解決物體繞流、旋翼等流固耦合問題,但用于解決流體域內固體間碰撞問題的研究較少.本文將嵌套網格方法與有限元法相結合,以分域耦合方法對浸沒在流體中的旋轉桿柱與井筒的碰撞問題進行求解,建立了流體中固體與固體碰撞的數值模擬方法,并與靜止流體中球形顆粒和壁面正、斜碰撞實驗數據進行對比,證明了本文所建立的數值模擬方法的正確性,據此研究了井筒內旋轉桿柱在不同流體黏度、轉速條件下的運動與碰撞特性.

1 井筒內浸沒在流體中的旋轉桿柱力學模型

1.1 模型建立

圖1 力學模型

模型的邊界條件如下:

(a)環空流體域 (b)固體域

1.2 環空流體域控制方程

連續性方程為

(1)

式中,ρ為密度,v為速度矢量.

動量守恒方程為

(2)

式中,f代表體積力矢量;p′為修正壓力,p′=p+(2/3)ρk;μequ為流體的有效黏度,μequ=μ+μt,μ和μt分別為流體動力黏度和湍流黏度.本文選擇使用SSTk-ω模型,其湍流黏度公式[9]為

(3)

式中,k和ω通過湍動能方程和湍動能耗散率方程求得,S是應變速率,a1為常數取值0.31,α*,F2分別為

(4)

(5)

其中,Ret=ρk/(μω),y為場點到最近壁面的距離.

1.3 環空流體域離散方程

背景網格流體域Ωf1的離散方程[10]為

(6)

式中,Af1,Bf1,Cf1,Df1,Gf1分別為背景網格流體域Ωf1的質量、對流、壓力、損耗、連續矩陣,其表達式分別為

其邊界條件為

(7)

組件網格流體域Ωf2的離散方程為

(8)

其邊界條件為

(9)

1.4 嵌套網格插值方法

由上述可知,在計算背景網格、組件網格的過程中,需要傳遞嵌套邊界條件,但是在嵌套邊界處,兩套網格尺寸往往不匹配,因此,要對組件網格與背景網格在嵌套邊界處的物理信息進行插值,本文采用3種不同的插值方法[11],具體如下.

1)嵌套邊界處網格尺寸相近

當嵌套邊界處組件網格和背景網格尺寸相近時,如圖4所示,選取邊界處任意一對單元,記為單元123,單元abc.若將單元123的信息傳遞給單元abc,則單元123記為供體單元,單元abc記為受體單元,其插值方法為

圖4 嵌套邊界處網格尺寸相近

(10)

2)供體單元尺寸大于受體單元尺寸

如圖5所示,123為供體單元,abd,bcd均為受體單元,且單元123的尺寸大于單元abd,bcd的尺寸.以供體單元123傳給受體單元abd為例,采用二階精度對其進行插值,具體過程如下:

圖5 供體單元尺寸大于受體單元尺寸

(11)

(12)

∑ψi=1, ∑xiψi=xX, ∑yiψi=yX.

(13)

求解式(13)可得各節點的線性有限元形函數.

根據Baker[13]理論,式(11)中f(X)為

f(X)=aψ1ψ2+bψ2ψ3+cψ3ψ1,

(14)

式中,a,b和c為待定系數.利用相鄰點4、5、6上的變量求解系數a,b,c,即

(15)

式中,左端矩陣下標1,2,…,6表示插值基函數的取值點;f(Xi)(i=1,2,…,6)為一階插值結果和節點已有值的誤差,

(16)

采用最小二乘法求解即可確定待定系數.確定待定系數后,就可對一階線性插值進行修正.只要確定單元的線性有限元形函數,其他單元類型的二階插值誤差修正函數計算方法類似.

3)供體單元尺寸小于受體單元尺寸

如圖6所示,134,234均為供體單元 ,abc為受體單元,且供體單元尺寸小于受體單元尺寸 ,則受體單元速度插值公式為

圖6 供體單元尺寸小于受體單元尺寸

(17)

式中,θi為供體單元中心到受體單元中心距離的倒數,n為受體單元所包含的供體單元個數.

1.5 環空流體域的嵌套網格計算框圖

基于上述模型和方法,建立環空流體域嵌套網格計算框圖,如圖7所示,其具體過程為:

圖7 環空流體域嵌套網格計算

1)建立背景網格和組件網格,并設置其邊界條件,對模型進行初始化.

2)對背景網格的固體位置處進行挖洞,根據重疊最小化原則確定背景網格和組件網格的范圍,基于式(10)、(11)、(17)對嵌套邊界處的受體單元進行插值,此過程為背景網格與組件網格的裝配.

3)將裝配后的背景網格和組件網格賦予新的邊界條件,并求解各自流體域方程,當前迭代步結束.

4)下一迭代步在當前迭代步求解的流體域中直接進行背景網格和組件網格的裝配,重復步驟2)、3).

1.6 桿柱動力學方程

浸沒在流體中的桿柱碰撞的動力學方程[14]為

(18)

式中,M,K,Kc(t),C分別為固體的質量矩陣、剛度矩陣、接觸剛度矩陣和阻尼矩陣[15],其具體表達式為

Ff(t)為流固耦合界面上的流體力;Fc(t)為碰撞時的接觸力;Fs(t)為固體受到的其他力;ds是固體節點位移矢量.

1.7 流固耦合分域求解耦合界面物理量傳遞

在本文中,采用分域方法對環空流體域與固體域耦合進行求解.固體域與環空流體域在耦合界面上需傳遞力和位移信息,分別遵循力平衡和位移協調條件:

Fs=Ff,ds=df,

(19)

式中,Fs和ds為耦合界面固體側任意一點力向量和位移向量,Ff和df為流體側對應點的力向量和位移向量.

由于流體網格與固體網格在耦合邊界處不匹配,需將耦合界面處的信息進行插值,其具體插值算法見參考文獻[16].

1.8 耦合界面上歸一化的收斂準則

在流固耦合計算的任一時間步中,耦合界面傳遞的力、位移信息都需反復迭代,當結果收斂方能進入下一時間步計算,其迭代的歸一化收斂準則為

(20)

ξjk+1(t)=αξjk(t)+(1-α)ξnew(t).

(21)

1.9 分域求解流程圖

基于上述計算方法, 采用分域方法對嵌套環空流體域與桿柱固體域耦合進行求解, 其計算框圖如圖8所示.

圖8 分域求解流程

1)建立固體域模型和流體域(背景網格和組件網格)模型,將組件網格嵌入背景網格,并定義邊界條件,假設時間步t已經完成.

2)令jk=0.

3)令jk=jk+1,求解流體域方程(6)和(8),流體域內求解并保持耦合界面位移djk(t)恒定, 由更新后的流場通過式(16)進一步得到耦合界面載荷Fjk(t).

4)將Fjk(t)傳遞到固體域, 求解固體域方程(18),得到耦合界面位移dnew(t).

5)將dnew(t)傳遞到流體域,求解方程(6)和(8),得到新的耦合界面載荷Fnew(t).

6)根據收斂準則式(20),若位移dnew(t)與載荷Fnew(t)不能同時滿足收斂準則,則返回迭代步3)～5);若不同時滿足收斂準則(20),則令t=t+Δt.

7)若t≤tend,則返回迭代步2)～6);否則,迭代結束.

2 嵌套網格數值模擬方法驗證

2.1 球形顆粒自由沉降

初始時刻球形顆粒和流體都是靜止的,然后球形顆粒開始在重力和水動力作用下做自由沉降運動,流體域大小為[-5D,5D]×[-5D,5D]×[0,80D],并在[-1D,1D]×[-1D,1D]×[0,80D]區域采用相同大小的網格,Δx=Δy=Δz=0.1D,向其他區域過渡時增長率為1.2,其有限元模型如圖9所示,球形顆粒物理參數見表1.

表1 球形顆粒物理參數

圖9 流體中球形顆粒自由沉降有限元模型

基于上述方法對單個球形顆粒的自由沉降運動進行計算,并與參考文獻[15]中的數值模擬和實驗結果進行對比.其速度云圖和速度隨時間變化曲線分別見圖10、圖11.

(a)浸入邊界法(文獻[15]) (b)嵌套網格

圖11 顆粒垂向速度vp隨時間的變化

由圖10和圖11的結果對比可知,當垂向速度vp為-0.7 m·s-1時,采用本文方法得到的云圖與浸入邊界法云圖基本一致,且顆粒的垂向速度隨時間變化曲線與實驗、浸入邊界法曲線基本吻合.因此可驗證本文建立的嵌套網格數值方法對球形顆粒自由沉降過程的計算是正確的.

2.2 球形顆粒與壁面正碰撞

流體域為長方體,其幾何尺寸長、寬、高分別為40 mm,40 mm,170 mm;建立流體域網格如圖12所示,在底部壁面處采用局部網格加密,且壁面處第一層網格高度為0.01 mm,邊界層增長率為1.1.玻璃板厚度為12 mm,劃分了4層單元,且玻璃板底面為固定約束.球形顆粒材料為鋼材,密度為7 800 kg·m-3,彈性模量為2.4×1011Pa,Poisson比為0.3,其計算工況見表2.

表2 球形顆粒與流體的物理參數

圖12 流體中球形顆粒與壁面碰撞有限元模型

對表中不同工況下球形顆粒與壁面碰撞和反彈過程進行計算,得到球形顆粒恢復系數e與Stokes數S的關系變化數值,并與Gondret等[17]的實驗結果進行對比,對比曲線如圖13所示.

圖13 恢復系數e與Stokes數S的關系曲線

從圖13中可以看出,當S小于10時,球形顆粒不反彈,隨著S逐漸增大,恢復系數逐漸趨向于球形顆粒干碰撞時的恢復系數,且數值模擬與實驗數值吻合較好,表明本文建立的數值方法對流體中固體與固體碰撞的模擬是可行的.

圖14為工況⑤球形顆粒與玻璃壁面碰撞和反彈過程中高度、速度隨時間的變化曲線,取相對時間t=0 s為球形顆粒與玻璃壁面第一次碰撞時間.

(a)高度 (b)速度

由圖14可知,由于流體黏性耗散和材料阻尼導致動能損失,球形顆粒后續的反彈高度越來越低,且反彈速度亦越來越小;數值模擬能夠很好地展現該工況的運動趨勢,且其數值與實驗數據吻合,表明本文建立的數值方法對流體中固體與固體正碰撞的模擬是準確的.

圖15為工況⑤球形顆粒與壁面正碰撞和反彈過程中流體渦量場隨時間的變化云圖.由圖15可知, 當球形顆粒下降時, 由于速度逐漸增大尾跡渦環也逐漸增大, 而且球形顆粒在尾跡位置開始出現與尾跡渦環方向相反的渦量.當t=-0.011 s時, 顆粒距壁面較遠, 顆粒仍處在加速過程中; 當t=-0.002 s,即球形顆粒與壁面距離約為顆粒半徑時,流體被擠出間隙,壁面處產生與顆粒尾跡渦環方向相反的渦量;當t=0 s時,顆粒第一次與壁面碰撞; 當t=0.002 s時,反彈過程中的顆粒通過主尾跡環向上運動并形成方向相反的二次渦環; 當t=0.039 s時,顆粒第一次反彈到最大高度,由于黏性耗散,顆粒周圍的初始渦環減小;當t=0.084 s時,顆粒回落并第二次與壁面碰撞.

(a)t=-0.011 s (b)t=-0.002 s (c)t=0 s

2.3 球形顆粒與壁面斜碰撞

基于上述方法對流體中球形顆粒與壁面斜碰撞和反彈過程也進行了數值模擬,流體域仍為長方體,其長、寬、高分別為60 mm,60 mm,50 mm, 建立網格如圖16所示,邊界層尺寸為0.01 mm,增長率為1.2,模型與Joseph和Hunt[18]實驗中使用的顆粒和黏性流體的物理參數相同,見表3.

表3 球形顆粒與流體的物理參數

圖16 流體中球形顆粒與壁面斜碰撞有限元模型 圖17 無量綱化入射角和反射角關系曲線

經計算,得到球形顆粒在流體中與壁面斜碰撞的無量綱化入射角與反射角變化結果,并與Joseph和Hunt[18]的實驗結果進行了對比,對比曲線如圖17所示.圖中ψin=tanθin,θin為小球碰撞點的入射角,ψout=etanθout,e為彈性恢復系數,θout為小球碰撞點的反射角.因受小球轉動速度的影響,tanθout=(vT-vω)/vN,vT為小球平動速度沿著壁面的速度分量,vω為小球轉動速度,vN為小球平動速度垂直于壁面的速度分量,即無量綱化反射角ψout與小球平動速度和轉動速度直接相關.

由圖17可知,球形顆粒反射角隨入射角增大而增大,且數值模擬結果與實驗結果基本吻合,由此可驗證采用本文方法對流體域中固體斜碰撞的數值模擬的正確性.

圖18為入射角θin=45°時球形顆粒與壁面斜碰撞和反彈過程中渦量場隨時間的變化云圖.由圖18可知,當t=-0.004 2～0 s時,流體的渦量場與球形顆粒與壁面正碰撞工況的分布相似,不同的是其渦量場傾斜了一定角度;當t=0.003 s時,球形顆粒反彈,下落時尾跡的初始渦環以相同的角度撞擊壁面,并在顆粒周圍出現一個新的渦環; 當t=0.003～0.007 9 s時,部分初始渦環隨著球形顆粒移動而逐漸脫落,直至分離.

(a)t=-0.004 2 s (b)t=-0.003 s (c)t=0 s

3 環空流體內旋轉桿柱與井筒碰撞計算結果

基于上述方法對1.1小節井筒內浸沒在流體中的旋轉桿柱與井筒的碰撞進行數值計算,計算參數見表4.

表4 旋轉桿柱和流體的物理參數

3.1 有流體與無流體計算結果對比

由表5可知固體域2 109個網格與3 234個網格第一次碰撞力變化不大,所以固體域劃分為2 019個網格,如圖19(a)所示;流體域105 396個網格和137 991個網格第一次碰撞力變化不大,所以流體域劃分為105 396個網格,邊界層尺寸為0.01 mm,增長率為1.2,如圖19(b)所示.

表5 網格無關性驗證

圖20為有環空流體和無環空流體旋轉桿柱與井筒碰撞力隨時間的變化曲線,井筒內環空流體隨桿柱運動產生的渦量如圖21所示.

圖20 旋轉桿柱與井筒碰撞力隨時間的變化

(a)t=0.03 s (b)t=0.06 s (c)t=0.071 1 s

由圖20可知,井筒與桿柱環空內無流體時,桿柱與井筒間的碰撞力最大值為6.13 N,而環空內有流體時,其碰撞力僅為0.56 N.即環空無流體時,桿柱與井筒間碰撞時的碰撞力數值明顯大于有流體的工況.因此,當井筒內有流體存在時,桿柱與井筒的碰撞劇烈程度明顯降低.

由圖21可知:當t=0.03 s時,流場受桿柱轉動影響出現渦量;到t=0.06 s時,受桿柱運動速度影響產生相反的渦量;t=0.071 1～0.074 1 s時,渦量場和小球正碰撞壁面相似;但當t=0.11 s桿柱上升到最高位置時,渦量場主要受轉速影響,桿柱平動速度的影響幾乎可以忽略;當t=0.142 s時,桿柱產生第二次碰撞.

3.2 黏度對環空流體中旋轉桿柱與井筒碰撞的影響分析

為探討流體黏度對環空流體中桿柱與井筒碰撞特性的影響,對流體黏度分別為0.01 N·s/m2,0.02 N·s/m2,0.03 N·s/m2,0.04 N·s/m2,0.05 N·s/m2進行計算,得到桿柱與井筒間碰撞力隨時間的變化曲線,如圖22所示.圖23為桿柱上第1次與井筒碰撞點的運動軌跡,圖24為桿柱中心點速度隨時間的變化.

(a)碰撞力隨時間變化 (b)碰撞力連線

圖23 不同流體黏度下桿柱上與井筒第1次碰撞點的運動軌跡

(a)0～0.5 s (b)0.4～0.5 s

由圖22(a)可知,桿柱與井筒初始碰撞時碰撞力數值較大,且最大值不受黏度影響,隨著時間推移碰撞力數值逐漸減小,且受黏度的影響越來越明顯,其數值隨著黏度增加呈下降趨勢.由圖22(b)可知,當黏度為0.01 N·s/m2,0.02 N·s/m2時,碰撞力波動下降,當黏度為0.03 N·s/m2,0.04 N·s/m2,0.05 N·s/m2時,桿柱與井筒碰撞力呈單調下降,黏度越大下降趨勢越明顯,其數值減小至0.1 N.

由圖23可知,桿柱上第1次與井筒碰撞點的運動軌跡隨黏度增大越來越趨于圓形,即桿柱運動范圍與黏度負相關,黏度越大,運動范圍越小,其軌跡越趨于圓形.

由圖24(a)可知, 桿柱中心點速度隨時間推移整體呈下降趨勢, 當桿柱與井筒發生碰撞時, 其數值產生突變, 且速度的振蕩程度和流體黏度成負相關, 黏度越低振蕩越劇烈.由圖24(b)可知,當黏度大于0.03 N·s/m2時,中心點速度在0.4 s后趨于穩定,且呈周期性變化.

3.3 轉速對環空流體中旋轉桿柱與井筒碰撞的影響分析

為探討桿柱轉速對環空流體中桿柱與井筒碰撞特性的影響,對桿柱轉速分別為12.56 rad/s,18.84 rad/s,25.12 rad/s,31.40 rad/s,37.68 rad/s進行計算,得到桿柱與井筒間碰撞力隨時間的變化曲線,如圖25所示.圖26為桿柱上第1次與井筒碰撞點的運動軌跡,圖27為桿柱中心點速度隨時間的變化.

(a)碰撞力隨時間變化 (b)碰撞力連線

圖26 不同轉速下桿柱上與井筒第1次碰撞點的運動軌跡 圖27 不同轉速下中心點速度隨時間的變化

由圖25(a)可知,桿柱與井筒間的碰撞力隨桿柱轉速增大而增大,隨時間推移整體呈下降趨勢.由圖25(b)可知, 當轉速為12.56 rad/s時, 0.1 s后碰撞力變化趨于平緩; 當轉速為18.86 rad/s和25.12 rad/s時,0.35 s后碰撞力變化趨于平緩,其數值僅為0.1 N;當轉速為31.40 rad/s和37.68 rad/s時,碰撞力始終在較大范圍內變化.

由圖26可知,桿柱運動范圍與桿柱轉速正相關,即轉速越大,其運動范圍越大.由圖27可知,桿柱中心點速度的振蕩程度與桿柱轉速成正相關,即轉速越高振蕩越劇烈.當桿柱轉速大于25.12 rad/s時,桿柱中心點速度不再隨時間推移整體呈下降趨勢,其數值在桿柱與井筒碰撞點處產生突變.

4 結 論

1)為研究環空流體內旋轉桿柱與井筒的碰撞特性,基于嵌套網格和分域求解算法,本文建立了井筒內浸沒在流體中的旋轉桿柱力學模型,將環空流體域離散為相互嵌套的背景網格和組件網格,并推導兩套網格間邊界條件信息插值公式,采用分域耦合算法對固體域與流體域耦合進行求解.該方法能夠避免流體域網格拓撲的改變,有效解決了傳統貼體網格方法在求解流體中存在固體間碰撞問題時網格易出現負體積的問題.

2)基于上述方法,對靜止流體中的球形顆粒與壁面正、斜碰撞進行數值計算,并與實驗結果進行對比,驗證了本文數值計算方法對流體域中固體與固體間碰撞模擬的正確性.

3)對不同流體黏度、不同桿柱旋轉速度下桿柱與井筒碰撞特性進行了研究,結果表明:桿柱與井筒間碰撞力隨黏度增大而降低,桿柱運動范圍與流體黏度負相關;隨著桿柱旋轉速度增大,桿柱與井筒間的碰撞力增大,桿柱運動范圍與其轉速正相關;且當桿柱與井筒發生碰撞時,桿柱上任一點速度出現突變.本文計算方法可為下一步研究石油工程領域中桿管柱偏磨、脫扣、斷裂失效提供一種行之有效的手段.