基于漂移流動模型的水平井巖屑床高度瞬態計算新方法

孫曉峰 姚 笛 孫士慧 紀國棟 于福銳 孫銘浩 趙元喆 陶 亮

1.“提高油氣采收率”教育部重點實驗室·東北石油大學 2.東北石油大學三亞海洋油氣研究院3.中國石油集團工程技術研究院有限公司

0 引言

目前，水平井已是油氣資源開采的常用井型，但鉆井過程中巖屑極易在造斜段和水平段沉降形成巖屑床，從而引發鉆具高摩阻，高扭矩、鉆頭托壓、卡鉆等一系列問題[1-6]。因此，有必要深入研究大斜度與水平井段的運移規律，建立瞬態巖屑運移模型對巖屑成床與運移規律進行描述，從而準確預測井眼環空內巖屑的分布情況。

現有國內外學者所建立的巖屑瞬態運移模型主要以巖屑床高度預測分層模型為主，包括兩層巖屑運移模型與三層巖屑運移模型[7-13]。Martins等[14]以兩層巖屑穩態運移模型為基礎，首次建立了時間相關的兩層巖屑瞬態模型。Wang等[15]等提出了基于SETS方法的三層瞬態模型求解方法。Zhang等[16]基于環空固液兩相流流型建立了一種動態巖屑床分層預測模型，模型的層數可由流型確定。Guo等[17]將三層穩態模型與兩層非穩態模型相結合，建立了三層瞬態運移模型。Fallah等[18]提出了基于漂移模型的兩層巖屑瞬態模型。

就目前的研究來看，巖屑分層運移模型雖然可以計算巖屑瞬態沉積高度，但模型求解速度慢，耗費計算資源，且僅在沒有擾動時適用。在大斜度與水平井段鉆進過程中，鉆桿旋轉、震動與起下均是流動中的大擾動，鉆井液與巖屑的固液兩相流存在沙丘流、分層流、分散懸浮等多種流型，此時巖屑分層模型并不適用。因此，本文不對環空中的固液兩相流動進行區域劃分，將環空中的鉆井液與巖屑視為一個整體，僅考慮相間的滑移作用，應用漂移流動模型建立一維環空固液兩相流動瞬態計算模型，采用有限體積法對井眼環空空間進行離散，離散后應用投影法進行求解。解得各項流動參數后再通過井斜角、巖屑床孔隙度等判別條件計算大斜度井段巖屑床的高度。以期為水平井的巖屑床清除提供指導，提高水平井井眼清潔效率。

1 環空固液兩相流瞬態流動模型

在水平井鉆進過程中，巖屑以一恒定的體積濃度進入環空中。流入環空的巖屑與鉆井液混合，形成混合流體不斷向井口運移[19]，如圖1所示。筆者使用漂移流動模型對環空鉆井液—巖屑的兩相流流動進行模擬，漂移流動模型最早由Zuber等[20]提出，模型將混合流動的兩相流體看作一種流體，忽略相間的質量傳遞與能量傳遞，考慮相間的滑移速度關系[21-22]。筆者應用Walton[23]所建立的固液兩相漂移速度關系，構建井眼環空固液兩相流瞬態流動模型。

圖1 井筒環空固液兩相流動示意圖

1.1 基本假設

環空固液兩相流瞬態流動模型的基本假設如下：①固液兩相在環空中的流動是一維的，各相流動參數保存在井筒網格節點，且在同一時間步時間內流動參數值保持恒定。②模型不考慮環空中的溫度變化，且在流動過程中鉆井液與巖屑均不可壓縮，巖屑顆粒均為均勻球型顆粒，具有相同的粒徑。③巖屑顆粒在小斜度井段（井斜角小于等于30°）與直井段中始終保持均勻分散狀態，在井斜角大于30°的井段中中會沉積形成巖屑床層。④假設模型計算得到的固相巖屑全部沉積形成巖屑床，則可通過巖屑床孔隙度將巖屑床體積分數轉化為無因次巖屑床高度，本文中假設巖屑床孔隙度為52%。

1.2 控制方程

鉆井液與巖屑的連續性方程分別為：

鉆井液與巖屑的混合動量方程為：

式中αl、αs分別表示液相、固相體積分數，無量綱；ρl、ρs分別表示液相、固相密度，kg/m3；ul、us分別表示液相、固相速度，m/s；p表示壓力，Pa；q表示源項。源項考慮了質量力與鉆井液流變性等外界因素對流動的影響，包含重力項與阻力項兩部分，其表達式為：

式中Fg、Fw分別表示重力分量、阻力分量，Pa/m；θ表示井斜角，（°）；g表示重力加速度，m/s2；di、do分別表示鉆桿外徑、井眼內徑，m；f表示范寧摩擦系數，無量綱；ρm表示固液相混合密度，kg/m3；um表示固液相混合速度，m/s。

固液相混合密度與混合速度的表達式分別為：

固相體積分數與液相體積分數之間的關系滿足：

1.3 漂移速度模型

考慮巖屑與鉆井液之間的滑移關系，固相與液相間的漂移速度關系為：

式中uslide表示相間滑移速度，m/s；up表示巖屑沉降速度，m/s；c表示相間分布系數，無量綱。相間分布系數c表示速度與濃度分布的不均勻度，受井斜角、鉆桿轉速、鉆井液密度、流型等多種因素影響，通常由實驗獲得。

巖屑在冪律流體中的沉降速度up表達式為[24]：

其中阻力系數的表達式為：

顆粒雷諾數與沉降速度的關系為：

式中ds表示巖屑顆粒的當量直徑，m；Cd表示巖屑顆粒在鉆井液中的阻力系數，無量綱；Res表示巖屑的顆粒雷諾數，無量綱。

2 模型求解方法

2.1 離散格式

由于水平井井筒的軸向尺寸遠大于徑向尺寸，不考慮管柱偏心問題，可以將三維空間井筒簡化為沿井筒軸向方向的一維空間，并對井筒一維空間均勻劃分。由于本文所研究的鉆井液與巖屑的固液兩相流是標準的不可壓縮流動，為避免出現奇偶失聯型壓力場影響計算精度，選用基于交錯網格的有限體積法對井筒空間進行離散，將井筒空間分為一系列控制體，對控制體單元進行質量與動量守恒分析，離散后的井筒網格如圖2所示，交錯網格將標量儲存在控制體中心，速度參數儲存在控制體邊界。邊界處的固相體積分數與液相體積分數可通過插值獲得。

圖2 井筒環空網格離散示意圖

連續性方程的時間推進可直接采用顯格式進行時間離散，n+1時刻固相連續性方程與液相連續性方程的離散格式分別為：

對于混合動量方程，由于鉆井液與巖屑的固液兩相流為不可壓縮流動，無法通過時間推進進行求解壓力項。因此本文采用投影法求解混合動量方程的壓力泊松方程，從而得到n+1時刻各相速度與壓力的值。

2.2 投影方法

投影法最早由Chorin提出[25]，用于求解不可壓縮的納維斯托克斯方程。該方法通過計算中間時間步得到壓力泊松方程，再對壓力泊松方程求解即可得到新時刻的壓力值。相比Simple類或人工壓縮性方法求解不可壓縮流動[26-27]，投影法無需反復內迭代進行求解，極大地提高了計算效率。使用投影法對混合動量方程進行求解，首先對式（3）進行時間離散，離散化后的表達式為：

將式（16）分裂為式（17）、（18）兩部分，將時間推進分為預算步、壓力修正步和最終步三部分進行。預算步與的壓力修正步的表達式分別為：

首先計算預算步，對式（17）進行空間離散，可以得到中間時間步的計算表達式：

接下來計算壓力修正步，由于環空中的固液兩相皆為不可壓縮流體，流體的密度隨流動的變化率較低，可以認為n+1時刻守恒變量的散度為0。對式（18）求散度，可以得到混合動量方程的一維壓力泊松方程：

對上式進行空間離散化，得到離散的壓力泊松方程：

重復式（19）～（21）可以得到整個井筒空間各節點的壓力泊松方程，將其化為矩陣的形式：

使用超松弛法對離散的壓力泊松方程進行求解，解出n+1時刻的井筒各節點處的壓力值pn+1后，即可進行最終步的計算，對式（18）進行空間離散，并將預算步與壓力修正步計算得到w*與pn+1代入離散后的式（18），即可計算得到wn+1。得到wn+1后還需進行守恒變量的還原求解n+1時刻的固相速度與液相速度。

將式（13）與式（14）進行變換，得到n+1時刻的液相體積分數與固相體積分數：

對固液兩相的漂移速度方程進行離散，得到n+1時刻的固相速度與液相速度關系式：

將式（23）、（24）、（25）代入式wn+1，進行變量還原，通量與液相速度之間的關系為：

得到n+1時刻井筒的液相速度分布與固相速度分布。至此n+1時刻井筒的壓力分布與速度分布全部求解完成，將該時刻的所有參數作為基礎，計算下一時刻井筒內各點參數值直至計算完畢。

2.3 模型計算步驟

模型的具體計算步驟如圖3所示。

圖3 算法計算步驟示意圖

①首先輸入鉆井液排量、鉆進速度、井口壓力等基本參數，結合已知n時刻井眼環空內各參數的值，計算n+1時刻邊界值。②根據n時刻的通量與源項，利用預算步計算中間時刻的守恒變量w*。③對壓力修正步求散度，得到混合動量方程的壓力泊松方程。對壓力泊松方程進行離散，將其化為矩陣形式。使用超松弛法對矩陣進行迭代求解，得到n+1時刻井筒的壓力分布pn+1。④根據n時刻的速度與體積分數，計算n+1時刻的固相體積分數與液相體積分數。將各相體積分數與計算得到的結果代入最終步，求解n+1時刻的固相速度與液相速度。

重復步驟①～④，即可計算得到全部時刻井筒內全部節點的固液相體積分數、固液速度與壓力值。

3 模型算例分析與驗證

為驗證算法的準確性與實用性，本節將應用所建立的模型模擬鉆井過程中巖屑在井眼環空中的瞬態運移規律，如圖4所示，并對數值模擬所計算結果進行分析驗證。

圖4 不同時刻環空巖屑體積分數及環空壓耗分布情況圖

3.1 模型算例分析

某小井眼水平井水平段長250 m，井眼內徑127 mm，鉆桿外徑76 mm，井斜角90°，巖屑密度2.60 g/cm3，巖屑粒徑34 mm，鉆井液密度1.18 g/cm3，鉆井液有效黏度50 mPa·s。環空排量選取8 L/s、10 L/s、12 L/s三種工況，機械鉆速30 m/h。將井筒劃分為100個網格，單位網格長度2.5 m，模擬計算總時長t=6 000 s，時間步長1 s。模擬開始后向入口處不斷注入巖屑，直至井筒中巖屑體積分數保持恒定時模擬結束。

在模擬開始0～100 s時巖屑未進入環空，此時水平段環空為純鉆井液，環空中的壓力呈線性分布，液相速度保持恒定。自100 s后開始注入巖屑，圖4-a為110 s時環空中的巖屑體積分數分布，此時巖屑剛流入環空，少量隨鉆井液向出口運移，絕大部分沉積于入口，且不同環空排量下的巖屑沉積量差異較大。由于此時三種排量下的鉆井參數一致，僅環空排量不同，且巖屑進入量較小，因此環空壓力與未進入巖屑時基本一致，環空壓耗主要為鉆井液流動所產生阻力，如圖4-b所示。

模擬時間達到280 s時，此時巖屑已大量進入環空。如圖4-c所示，高排量下巖屑隨鉆井液快速流動，自鉆頭處沿環空逐漸形成穩定巖屑床，巖屑體積分數保持恒定。低排量下巖屑仍然在不斷聚集，巖屑床仍處于發育階段，尚未完全形成，巖屑體積分數沿環空不斷降低，且環空排量越低巖屑聚集速度越慢。此時的環空壓耗與井底壓力明顯升高，環空壓耗由巖屑運移與鉆井液流動所產生阻力組成，巖屑體積分數較高的井段環空壓耗大幅增加，如圖4-d所示。

圖4-e、f為模擬時間達到1 300 s時環空巖屑分布情況與壓耗分布情況，此時高排量下的環空巖屑體積分數已趨于穩定，固液兩相的流動狀態不再隨時間發生變化。低排量下的巖屑也開始在入口處形成穩定巖屑床，并沿環空向出口不斷運移。此時低排量下的井底壓力已高于高排量下的井底壓力，這是因為巖屑運移所產生阻力遠高于鉆井液流動阻力，環空巖屑體積分數成為影響環空壓耗的主要因素。

模擬時間達到4 000 s，此時排量為10 L/s、12 L/s的環空中已形成穩定巖屑床，環空巖屑體積分數保持不變，排量為8 L/s的環空中部分井段的巖屑體積分數仍未穩定，直至模擬達到5 400 s時才形成穩定巖屑床，此時三種工況下的井底壓力與環空壓耗也趨于穩定。

圖5-a、b為排量為8 L/s時環空固、液相速度的分布情況。排量一定的情況下巖屑運移速度主要受巖屑體積分數影響，巖屑體積分數越大則流動產生的阻力越大，運移速度越慢。且巖屑進入環空后不斷聚集并沉積形成巖屑床，導致鉆井液的流道減小，由于巖屑運移速度較低且密度高于鉆井液密度，因此巖屑體積分數越高的井段液相速度越高。

圖5 不同時刻環空固、液相速度分布情況圖

圖6為不同環空排量下形成巖屑床時間，環空排量較低時巖屑運移速度也較小，從鉆頭處運移至出口的時間更長。且由于低排量時鉆井液流動速度緩慢，導致大量巖屑聚集，環空中的巖屑體積分數增高，巖屑產生堆積效應，極大地降低了巖屑運移速度。因此低排量下形成穩定巖屑床所需時間要遠遠大于高排量下形成穩定巖屑床所需時間，巖屑的堆積效應在長水平段鉆進過程中尤為明顯，隨著水平段長度的增加，巖屑堆積導致井下摩阻呈指數性增長，嚴重影響了水平段的安全鉆進。利用本文所建立模型，可以計算得到形成穩定巖屑床所需時間，從而合理設計停鉆循環時間或改變泵排量與鉆桿轉速，在巖屑床影響正常鉆進前及時進行清除，預防高摩阻所產生的卡鉆、鉆具磨損等問題，提升大位移井鉆井過程中的井眼清潔效率[28-30]。

圖6 不同環空排量下巖屑床形成時間圖

3.2 模型驗證

為驗證本文所建立模型的準確性，將模型計算結果與穩態巖屑床高度計算模型[31]進行了對比。取本文算例中的鉆井參數進行計算，結果如圖7所示。

圖7 模型對比驗證圖

由圖7可知，在環空排量較大時，本文模型與穩態模型計算結果較為接近，在環空排量較小時本文模型所計算巖屑床高度低于穩態模型，由于穩態模型所預測巖屑床高度略高于室內實驗的實測值1%～8%，因此可以認為本文所建立模型可以基本滿足實際工況中對巖屑床高度進行預測的需求。

4 結論

1）提出了一種一維瞬態模擬算法，根據鉆井攜屑實際工程特點將環空固液兩相流動簡化為一維，采用基于離散網格有限體積法與投影方法，對模型進行離散化求解，得到了巖屑床高度在井筒中的變化規律。通過與穩態巖屑床高度模型進行對照，驗證了本文模型與計算公式的準確性。

2）應用投影方法對固液混相動量守恒的納維斯托克斯方程進行數值求解，首先利用中間時間步求解壓力泊松方程得到環空中的壓力分布，再將連續性方程與漂移流動關系代入混合動量方程，求解得到固相速度、液相速度、固相體積分數與液相體積分數。算法將壓力與速度解耦分別進行求解，且無需反復進行內迭代計算，極大地提升了計算效率。

3）所建立模型可用于研究巖屑床在環空中的形成規律，從而實現水平井環空巖屑床高度的實時監測，且模型求解效率較高，可以滿足長水平井段巖屑床計算和監測的需求。模型計算結果可以為水平井的巖屑床高效清除與卡鉆事故風險預測提供一定的理論指導。