基于嵌套網(wǎng)格的環(huán)空流體內(nèi)旋轉(zhuǎn)桿柱與井筒碰撞特性研究*

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

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

0 引 言

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

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

1.1 模型建立

圖1 力學(xué)模型

模型的邊界條件如下:

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

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

連續(xù)性方程為

(1)

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

動(dòng)量守恒方程為

(2)

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

(3)

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

(4)

(5)

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

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

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

(6)

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

其邊界條件為

(7)

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

(8)

其邊界條件為

(9)

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

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

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

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

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

(10)

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

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

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

(11)

(12)

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

(13)

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

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

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

(14)

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

(15)

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

(16)

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

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

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

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

(17)

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

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

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

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

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

2)對(duì)背景網(wǎng)格的固體位置處進(jìn)行挖洞,根據(jù)重疊最小化原則確定背景網(wǎng)格和組件網(wǎng)格的范圍,基于式(10)、(11)、(17)對(duì)嵌套邊界處的受體單元進(jìn)行插值,此過(guò)程為背景網(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 桿柱動(dòng)力學(xué)方程

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

(18)

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

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

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

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

Fs=Ff,ds=df,

(19)

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

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

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

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

(20)

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

(21)

1.9 分域求解流程圖

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

圖8 分域求解流程

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

2)令jk=0.

3)令jk=jk+1,求解流體域方程(6)和(8),流體域內(nèi)求解并保持耦合界面位移djk(t)恒定, 由更新后的流場(chǎng)通過(guò)式(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)不能同時(shí)滿(mǎn)足收斂準(zhǔn)則,則返回迭代步3)～5);若不同時(shí)滿(mǎn)足收斂準(zhǔn)則(20),則令t=t+Δt.

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

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

2.1 球形顆粒自由沉降

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

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

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

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

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

圖11 顆粒垂向速度vp隨時(shí)間的變化

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

2.2 球形顆粒與壁面正碰撞

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

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

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

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

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

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

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

(a)高度 (b)速度

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

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

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

2.3 球形顆粒與壁面斜碰撞

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

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

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

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

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

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

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

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

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

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

3.1 有流體與無(wú)流體計(jì)算結(jié)果對(duì)比

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

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

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

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

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

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

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

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

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

(a)碰撞力隨時(shí)間變化 (b)碰撞力連線(xiàn)

圖23 不同流體黏度下桿柱上與井筒第1次碰撞點(diǎn)的運(yùn)動(dòng)軌跡

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

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

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

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

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

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

(a)碰撞力隨時(shí)間變化 (b)碰撞力連線(xiàn)

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

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

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

4 結(jié) 論

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

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

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