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

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

0 引 言

浸沒在流體中桿管間的接觸碰撞是石油工程最常見的現(xiàn)象,由于碰撞引起流體域網(wǎng)格拓?fù)浣Y(jié)構(gòu)的變化,采用界面解析直接數(shù)值模擬方法研究流體中桿管間的碰撞具有很大的挑戰(zhàn),因此尚未得到廣泛的研究.

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

1 井筒內(nèi)浸沒在流體中的旋轉(zhuǎn)桿柱力學(xué)模型

1.1 模型建立

圖1 力學(xué)模型

模型的邊界條件如下:

(a)環(huán)空流體域 (b)固體域

1.2 環(huán)空流體域控制方程

連續(xù)性方程為

(1)

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

動量守恒方程為

(2)

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

(3)

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

(4)

(5)

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

1.3 環(huán)空流體域離散方程

背景網(wǎng)格流體域Ωf1的離散方程[10]為

(6)

式中,Af1,Bf1,Cf1,Df1,Gf1分別為背景網(wǎng)格流體域Ωf1的質(zhì)量、對流、壓力、損耗、連續(xù)矩陣,其表達(dá)式分別為

其邊界條件為

(7)

組件網(wǎng)格流體域Ωf2的離散方程為

(8)

其邊界條件為

(9)

1.4 嵌套網(wǎng)格插值方法

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

1)嵌套邊界處網(wǎng)格尺寸相近

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

圖4 嵌套邊界處網(wǎng)格尺寸相近

(10)

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

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

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

(11)

(12)

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

(13)

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

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

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

(14)

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

(15)

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

(16)

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

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

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

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

(17)

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

1.5 環(huán)空流體域的嵌套網(wǎng)格計算框圖

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

圖7 環(huán)空流體域嵌套網(wǎng)格計算

1)建立背景網(wǎng)格和組件網(wǎng)格,并設(shè)置其邊界條件,對模型進(jìn)行初始化.

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

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

4)下一迭代步在當(dāng)前迭代步求解的流體域中直接進(jìn)行背景網(wǎng)格和組件網(wǎng)格的裝配,重復(fù)步驟2)、3).

1.6 桿柱動力學(xué)方程

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

(18)

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

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

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

在本文中,采用分域方法對環(huán)空流體域與固體域耦合進(jìn)行求解.固體域與環(huán)空流體域在耦合界面上需傳遞力和位移信息,分別遵循力平衡和位移協(xié)調(diào)條件:

Fs=Ff,ds=df,

(19)

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

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

1.8 耦合界面上歸一化的收斂準(zhǔn)則

在流固耦合計算的任一時間步中,耦合界面?zhèn)鬟f的力、位移信息都需反復(fù)迭代,當(dāng)結(jié)果收斂方能進(jìn)入下一時間步計算,其迭代的歸一化收斂準(zhǔn)則為

(20)

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

(21)

1.9 分域求解流程圖

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

圖8 分域求解流程

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

2)令jk=0.

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

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

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

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

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

2 嵌套網(wǎng)格數(shù)值模擬方法驗證

2.1 球形顆粒自由沉降

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

表1 球形顆粒物理參數(shù)

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

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

(a)浸入邊界法(文獻(xiàn)[15]) (b)嵌套網(wǎng)格

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

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

2.2 球形顆粒與壁面正碰撞

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

表2 球形顆粒與流體的物理參數(shù)

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

對表中不同工況下球形顆粒與壁面碰撞和反彈過程進(jìn)行計算,得到球形顆?；謴?fù)系數(shù)e與Stokes數(shù)S的關(guān)系變化數(shù)值,并與Gondret等[17]的實驗結(jié)果進(jìn)行對比,對比曲線如圖13所示.

圖13 恢復(fù)系數(shù)e與Stokes數(shù)S的關(guān)系曲線

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

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

(a)高度 (b)速度

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

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

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

2.3 球形顆粒與壁面斜碰撞

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

表3 球形顆粒與流體的物理參數(shù)

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

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

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

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

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

3 環(huán)空流體內(nèi)旋轉(zhuǎn)桿柱與井筒碰撞計算結(jié)果

基于上述方法對1.1小節(jié)井筒內(nèi)浸沒在流體中的旋轉(zhuǎn)桿柱與井筒的碰撞進(jìn)行數(shù)值計算,計算參數(shù)見表4.

表4 旋轉(zhuǎn)桿柱和流體的物理參數(shù)

3.1 有流體與無流體計算結(jié)果對比

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

表5 網(wǎng)格無關(guān)性驗證

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

圖20 旋轉(zhuǎn)桿柱與井筒碰撞力隨時間的變化

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

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

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

3.2 黏度對環(huán)空流體中旋轉(zhuǎn)桿柱與井筒碰撞的影響分析

為探討流體黏度對環(huán)空流體中桿柱與井筒碰撞特性的影響,對流體黏度分別為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進(jìn)行計算,得到桿柱與井筒間碰撞力隨時間的變化曲線,如圖22所示.圖23為桿柱上第1次與井筒碰撞點的運動軌跡,圖24為桿柱中心點速度隨時間的變化.

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

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

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

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

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

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

3.3 轉(zhuǎn)速對環(huán)空流體中旋轉(zhuǎn)桿柱與井筒碰撞的影響分析

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

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

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

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

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

4 結(jié) 論

1)為研究環(huán)空流體內(nèi)旋轉(zhuǎn)桿柱與井筒的碰撞特性,基于嵌套網(wǎng)格和分域求解算法,本文建立了井筒內(nèi)浸沒在流體中的旋轉(zhuǎn)桿柱力學(xué)模型,將環(huán)空流體域離散為相互嵌套的背景網(wǎng)格和組件網(wǎng)格,并推導(dǎo)兩套網(wǎng)格間邊界條件信息插值公式,采用分域耦合算法對固體域與流體域耦合進(jìn)行求解.該方法能夠避免流體域網(wǎng)格拓?fù)涞母淖?有效解決了傳統(tǒng)貼體網(wǎng)格方法在求解流體中存在固體間碰撞問題時網(wǎng)格易出現(xiàn)負(fù)體積的問題.

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

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