砂礫形狀對(duì)儲(chǔ)氣庫井筒沖蝕的影響

彭 濤， 張菲菲， 方含之

(長江大學(xué)石油工程學(xué)院， 武漢 430100)

地下儲(chǔ)氣庫工藝作為季候調(diào)峰和能源儲(chǔ)備基礎(chǔ)設(shè)施，其完整性研究對(duì)天然氣的高效、安全、經(jīng)濟(jì)生產(chǎn)有著顯著的意義[1]。儲(chǔ)氣庫工藝有其特殊性：井筒多處于交變載荷工況；井筒中流體系持續(xù)性周期流動(dòng)、流量大、作業(yè)壓力大；儲(chǔ)氣遇水易酸化；砂巖地層(弱膠結(jié)巖層)改造時(shí)，持續(xù)拉伸作用等易導(dǎo)致砂礫彈塑性變化和儲(chǔ)層完整性下降，進(jìn)而加劇井筒內(nèi)出砂，嚴(yán)重出砂更會(huì)導(dǎo)致砂埋氣層，鉆柱磨損失效，降低生產(chǎn)壓力及產(chǎn)量[2-3]。因此，儲(chǔ)氣庫儲(chǔ)層出砂問題，亟須各研究者加以重視。

關(guān)于研究固體顆粒對(duì)筒壁和設(shè)備的沖蝕作用，中外已有大量實(shí)驗(yàn)分析和報(bào)告，主要依托Magneé[4]結(jié)合材料性質(zhì)、微觀結(jié)構(gòu)、沖擊角或硬度等建立的廣義侵蝕磨損定理或其延伸。而對(duì)一般彎管下，影響因素主要包括彎管曲率、彎管直徑、氣相流速、顆粒密度直徑等，其中最關(guān)鍵的是顆粒直徑，文獻(xiàn)[5-8]總結(jié)了一般性結(jié)論：在同一攜砂流速下，顆粒直徑越大，磨蝕作用越大。而考慮井筒沖蝕過程中，出砂基巖直徑較小且近似，但不同砂礫形狀對(duì)沖蝕影響研究較少。

因此，現(xiàn)基于實(shí)際采氣井工況，結(jié)合計(jì)算流體力學(xué)(computational fluid dynamics, CFD)思想及磨蝕工況設(shè)定方法，選取適宜湍流模型并建立沖蝕模型，提出不同儲(chǔ)層出砂形狀特性，及其對(duì)井筒磨蝕率的影響。

1 組分運(yùn)移模型

在流體動(dòng)力學(xué)計(jì)算中，應(yīng)先分析流場內(nèi)各物質(zhì)分布、流動(dòng)形態(tài)等。主要考慮儲(chǔ)氣由儲(chǔ)層進(jìn)入生產(chǎn)管并運(yùn)移至井筒造斜點(diǎn)的流動(dòng)特性，其間儲(chǔ)氣流動(dòng)機(jī)制主要是由壓力梯度所引起的黏性流動(dòng)。因此，針對(duì)造斜點(diǎn)間多組分運(yùn)移，只考慮較為簡單的氣相攜運(yùn)模型，其固相與氣相間并無能量傳遞，僅存在氣相湍流作用攜砂和砂礫沖擊壁面，并且攜砂模型下，也不需考慮各氣相間相互吸附作用。對(duì)于較小顆粒與壁面間的相互作用，分子力較慣性力占主導(dǎo)，故存在粒子由于凝聚力被壁面吸收而沒有反彈或滑移的力學(xué)現(xiàn)象。

2 數(shù)值模型建立

由于沖蝕問題作為一個(gè)長周期性行為，依托較高模擬精度的CFD方法也是分析該類問題的重要途徑[9-10]。在考慮數(shù)值方法時(shí)，不得忽略物質(zhì)連續(xù)變形特性，而針對(duì)剛體物質(zhì)，如砂礫，一般不考慮其變形。將流場內(nèi)流體做連續(xù)相介質(zhì)(流體各連續(xù)相控制方程均依托于連續(xù)性方程和N-S方程)。將砂礫作為離散項(xiàng)介質(zhì)，即以歐拉法形式定義固相微團(tuán)，通過對(duì)各個(gè)微團(tuán)研究以概括整個(gè)流場工況，并對(duì)各流固耦合力學(xué)方程進(jìn)行數(shù)值求解，模型中假設(shè)包括：

(1)設(shè)置侵蝕模型下粒徑函數(shù)為常數(shù)，取常規(guī)設(shè)定值1.8×10-9。

(2)不考慮顆粒間相互作用力。

(3)不考慮沖蝕壁面時(shí)產(chǎn)生的幾何形狀變化。

取短曲率水平井造斜段方法：設(shè)計(jì)造斜率為10(°)/m(曲率半徑5.73 m)。旨在考慮攜砂氣流對(duì)管段各位置磨蝕情況的機(jī)理性分析，并且若取實(shí)際長徑比，對(duì)于計(jì)算機(jī)運(yùn)行能力要求過高，因此依據(jù)動(dòng)力學(xué)相似原理，將井筒原型定義為如圖1所示的幾何模型。

圖1 造斜幾何模型及網(wǎng)格劃分Fig.1 Deflecting geometric model and meshing

前期利用前處理軟件繪制網(wǎng)格模型，再建立沖蝕模型進(jìn)行模擬。計(jì)算時(shí)考慮穩(wěn)態(tài)、等溫及湍流條件。對(duì)各方向動(dòng)量、連續(xù)性方程等采用離散形式計(jì)算時(shí)，考慮到儲(chǔ)氣可壓及建立結(jié)構(gòu)化網(wǎng)格方法，使用二階迎風(fēng)格式。選取求解器時(shí)，考慮模型運(yùn)動(dòng)邊界條件設(shè)定方法及兩相設(shè)定形式，取壓力分離形式求解器。該求解器下壓力-速率耦合形式有以下算法：SIMPLE、SIMPLEC、PISO及COUPLED。韋魯濱等[11]研究表明，PISO算法對(duì)于各變量均進(jìn)行一次預(yù)測及兩次修正，其對(duì)計(jì)算值更為準(zhǔn)確。并且該算法對(duì)于非耦合形式的動(dòng)量方程及標(biāo)量方程(如熱邊界條件為常數(shù)時(shí))，較SIMPLEC及SIMPLE算法效率更高，故采用PISO算法。計(jì)算時(shí)對(duì)描述流場運(yùn)動(dòng)狀態(tài)各控制方程進(jìn)行迭代計(jì)算，得出解析解，進(jìn)而定性描述流固耦合問題。迭代過程中取平均磨蝕量和各方程殘差設(shè)定值為收斂依據(jù)。

2.1 湍流模型定義

定義湍流模型時(shí)，考慮到流場內(nèi)攜砂尺寸較小，并被壁面捕獲。故流場內(nèi)不易產(chǎn)生大渦街，因此不考慮大渦模擬(large eddy simulation, LES)方法，而選取一般雷諾時(shí)均法(reynolds-averaged Navier-Stokes, RNAS)方法。故選取Realizek-e模型、RNGk-e及SSTk-ω模型。

圖2為3種湍流模型不同流速下沿管壁位置最大磨蝕率驗(yàn)證結(jié)果，其中氣相速率較大時(shí)三種方法對(duì)應(yīng)沖蝕模型分析結(jié)果近似，而在流速較小時(shí)RNGk-e方法沿程振蕩幅度較大，并且進(jìn)入豎直管段時(shí)，幅度最大，因此有結(jié)論：RNGk-e及SSTk-ω方法效果均較好，其中SSTk-ω模型綜合了k-e模型在近壁區(qū)的計(jì)算特點(diǎn)，并加入橫向耗散項(xiàng)，且在定義湍流黏度時(shí)加入剪切應(yīng)力影響，故更適于黏性流動(dòng)機(jī)制模擬。因此使用SSTk-ω湍流模型。

圖2 流速6 m/s和3 m/s時(shí)不同湍流模型下的 沿程磨蝕分布Fig.2 The wear distribution of different turbulence models along the path under 6 m/s and 3 m/s

2.2 沖蝕模型定義

基于廣義侵蝕磨損定理不同參數(shù)設(shè)定方式，提出Oka模型[12]、Peng模型[13]和Vieria模型[14]等較好的預(yù)測模型,其中對(duì)于彎管形式?jīng)_蝕模型，Vieria等[14]和劉琦等[15]建立了與實(shí)驗(yàn)值更為接近的模型，其基本表達(dá)式為

(1)

式(1)中：vp為顆粒相對(duì)速度，m/s；b(v)為顆粒相對(duì)速度函數(shù)(或速度指數(shù)函數(shù))；α為顆粒與壁面間沖擊角,rad；C(dp)為顆粒粒徑函數(shù)，設(shè)為常數(shù)，取1.8×10-9；Aface為壁面面積，m2；mp為顆粒質(zhì)量流量，kg/s；f(α)為沖擊角函數(shù)；其中參數(shù)C、b以及f均為定義壁面的邊界條件，以常數(shù)、多項(xiàng)式、分段多項(xiàng)式或分段線性方法等定義。

可通過取適宜的碰撞模型和磨蝕參數(shù)，再定義攜砂氣流初始條件及邊界條件。

2.3 壁面函數(shù)定義

在提出對(duì)于各邊界條件設(shè)定時(shí)，主要定義沖擊角函數(shù)、反彈系數(shù)設(shè)定、粒徑函數(shù)設(shè)定(上文已討論)、速度指數(shù)函數(shù)等方法設(shè)定?？紤]沖蝕問題時(shí)，顆粒務(wù)必貫穿壁面邊界層，其邊界流動(dòng)可視為層流，因此對(duì)于整個(gè)流場設(shè)定湍流模型時(shí)，需考慮壁面區(qū)域。利用無量綱u+及以定義y+邊界層中三個(gè)子層(黏性底層、過渡層、對(duì)數(shù)律層)，劃分圖示如圖3所示。

(2)

式(2)中：u為流場時(shí)均流速,m/s；Δy為流體至壁面距離,m；ρ為流場混合密度，kg/m3;uτ為壁面摩擦速度,m/s。

黃詩嵬[16]研究表明，對(duì)近壁面區(qū)域物理量計(jì)算方法，一種是采用半經(jīng)驗(yàn)的壁面函數(shù)來過渡壁面與核心湍流區(qū)域；另一種是對(duì)湍流模型進(jìn)行修改，對(duì)邊界網(wǎng)格加密，使其適用于近壁面區(qū)域。

未加密之前，壁面y+為100～500。隨后，采用自適應(yīng)形式對(duì)壁面網(wǎng)格加密，如圖4所示。圖示各點(diǎn)為沿管各網(wǎng)格節(jié)點(diǎn)y+值最大最小范圍，將y+定義到50～400范圍值，將網(wǎng)格第一個(gè)內(nèi)節(jié)點(diǎn)位置修改至對(duì)數(shù)律層成立的范圍，使其配置到旺盛湍流區(qū)域，進(jìn)而在使用SSTk-ω模型模擬壁面區(qū)域運(yùn)動(dòng)時(shí)更為精確。

圖3 壁面三子層劃分及對(duì)應(yīng)速度分布Fig.3 Division of three sublayers of wall surface and corresponding velocity distribution

圖4 自適應(yīng)前壁面和后壁面部分y+Fig.4 y+ at adaptive front and rear wall segment

2.4 壁面反彈系數(shù)定義

設(shè)定壁面反彈系數(shù)時(shí)，可取葛令行等[17]以依托粒子跟蹤測速技術(shù)(particle tracking velocimetry，PTV)對(duì)碰撞前后速度分解，能夠獲得切向和法向反彈系數(shù)。或采用經(jīng)典反彈模型(已知入射角和入射速率、反彈速率，利用數(shù)學(xué)模型計(jì)算)定義反彈系數(shù)，由于需配套設(shè)備要求高，且操作復(fù)雜。鑒于實(shí)驗(yàn)平臺(tái)有限，考慮采取指數(shù)多項(xiàng)式定義各反彈系數(shù)。由于小粒徑砂礫法向反彈系數(shù)與入射角正相關(guān)，且變化明顯，故顆粒與延展性壁面碰撞時(shí)，各反彈系數(shù)較大。取李潔瓊[18]對(duì)280～400 μm粒子與壁面發(fā)生彈性碰撞實(shí)驗(yàn)所得各彈性系數(shù)多項(xiàng)式，并考慮砂礫入射角度和通過防砂篩管逃逸砂礫直徑，進(jìn)行插值計(jì)算，得出反彈系數(shù)多項(xiàng)式如下。

切向反彈系數(shù)：

et=0.997 6+0.942 1α-1.642α2+0.706 7α3

(3)

法向反彈系數(shù)：

en=0.874 1+0.964 1α-1.752 4α2+0.961 8α3

(4)

2.5 速度指數(shù)定義

對(duì)于一般沖蝕模型，其速度指數(shù)函數(shù)的設(shè)定，主要與沖蝕壁面材質(zhì)、壁面溫度梯度有關(guān)，由于一般套管鋼質(zhì)壁面導(dǎo)熱系數(shù)較大，各處溫差幅度小，故壁面溫度梯度影響可忽略。對(duì)于延展性材料，參考文獻(xiàn)[19-20]，沖蝕模型中速率指數(shù)為2；Stack等[21]認(rèn)為速度指數(shù)應(yīng)介于2～2.5，甚至后期已開發(fā)模型的速度指數(shù)超過3。依托沖蝕模型，并考慮與Stack方法下的模型存在幾何相似，取速度指數(shù)為2.8。

2.6 沖擊角函數(shù)定義

對(duì)于不同沖蝕模型下沖擊角函數(shù)定義方法，可通過調(diào)研各有關(guān)文獻(xiàn)進(jìn)行賦值。例如Tusla Angle Dependent模型[22]提出經(jīng)典砂礫沖蝕鋼制壁面各沖蝕率改寫公式和相應(yīng)沖擊角函數(shù)。對(duì)于侵蝕性磨損如海底管道磨蝕情況，也提出利用Tabak-off和Grant模型[23]。對(duì)地層砂礫和井筒生產(chǎn)套管沖蝕模型各參數(shù)，取經(jīng)典模型，以多個(gè)線性段方法定義砂礫沖擊角函數(shù)：

f(θ)=Aθ+Bθ2+Cθ3+Dθ4+Eθ5+Fθ6+Gθ7+Hθ8

(5)

表1 砂礫沖擊角取值

3 模型分析

考慮用圓球度表征對(duì)不同儲(chǔ)氣庫出砂形狀變化，如基巖地層，其圓球度較大，而針對(duì)頁巖屑，其圓球度較小。圓球度(SF)為描述研究對(duì)象邊界復(fù)雜程度：

(6)

式(6)中：Ap為研究對(duì)象表面積，m2；Vp為研究對(duì)象體積，m3。

SF值為0～1，越趨近于1，邊界過度越光順，SF=1時(shí)，邊界區(qū)域?qū)賵A形；若邊界起伏程度越大，其形狀系數(shù)越小。計(jì)算圓球度表面積時(shí)，可取樊騏鋮[24]提出的像素投影法。對(duì)預(yù)測模型影響主要依據(jù)顆粒邊界與壁面碰撞前后運(yùn)動(dòng)方程，將各撞擊點(diǎn)前后軌跡整合到對(duì)應(yīng)磨損模型，再進(jìn)行磨蝕預(yù)測。

取砂礫密度2 500 kg/m3，防砂篩管網(wǎng)孔尺寸0.23 mm(忽略持續(xù)過流沖蝕對(duì)網(wǎng)孔尺寸變化率影響)，直徑小于0.2 mm的砂礫從進(jìn)口(inlet)被攜入，其質(zhì)量流量取1 kg/s，進(jìn)口流速6 m/s；在設(shè)置砂礫圓球度時(shí)，考慮基巖砂礫、頁巖屑設(shè)定等，其他設(shè)定參數(shù)不變，參考實(shí)驗(yàn)數(shù)據(jù)如表2所示。

表2 顆粒形狀系數(shù)表

在相同進(jìn)口條件和邊界條件下，對(duì)于井筒內(nèi)出砂，流場對(duì)砂礫攜砂作用最大位置為圖5位置1、2處。位置1主要由于砂礫進(jìn)入造斜管段，受力方向改變，砂礫沿管內(nèi)部偏移；對(duì)于位置2處，曳力完全沿豎直方向且由于渦街作用，砂礫進(jìn)一步沿內(nèi)部偏移，導(dǎo)致在進(jìn)入豎直管段后對(duì)壁面磨蝕作用減小。圖6和圖7為對(duì)于不同攜砂形狀系數(shù)設(shè)定(即砂礫邊界光順度)沖蝕行為。

圖5 砂礫攜運(yùn)模型圖示Fig.5 Gravel transport model diagram

圖6 不同圓球度沖蝕及砂礫軌跡模型Fig.6 Model of erosion and gravel track with different sphericity

圖7 管壁各處最大磨蝕分布Fig.7 Maximum wear distribution on pipe wall

由圖6、圖7可知，由于進(jìn)口直管段沙礫平行于壁面流動(dòng)，進(jìn)而對(duì)于進(jìn)口段，其沖蝕率較??；并且儲(chǔ)氣紊流波動(dòng)對(duì)于顆粒路徑也有影響，且在彎管處，多存在紊流增強(qiáng)及二次渦街等，因此加大了流場運(yùn)動(dòng)復(fù)雜性，進(jìn)而彎管外側(cè)中心處磨蝕最嚴(yán)重；考慮不同形狀系數(shù)對(duì)磨蝕影響，由于邊界越粗糙，其流線軌跡分布越復(fù)雜，即對(duì)造斜段外側(cè)沖蝕分布更廣，受曳力較小，在貫穿黏性底層時(shí)能量耗損較大或被底層吸收，即顆粒與壁面間相互作用較弱，壁面磨蝕程度較輕。

顆粒圓球度越小，不規(guī)則顆粒在流場中受曳力越大，并且最大磨蝕位置也會(huì)向上偏移。形狀越不規(guī)則的顆粒(如形狀系數(shù)為0.32、0.56時(shí))，其儲(chǔ)氣攜砂能力越強(qiáng)，顆粒在進(jìn)入豎直管段時(shí)，仍存在一定磨蝕率。

對(duì)圓球度與磨蝕率進(jìn)行線性回歸分析時(shí)，如圖8所示，可決系數(shù)R2=0.903 5，因此在其他參數(shù)變化時(shí)，如氣相流速，可考慮只取兩個(gè)參考數(shù)據(jù)做出相關(guān)擬合曲線方程，對(duì)其他不同形狀系數(shù)粒子求解時(shí)只需代入單一值進(jìn)行求解。

圖8 圓球度與最大磨蝕率關(guān)系Fig.8 The relationship between SF and maximum erosion rate

4 結(jié)論

通過以上研究，得出結(jié)論及建議如下。

(1)對(duì)于一般氣固兩相流，NRAS湍流模型均適宜。若不過分考慮計(jì)算能力，其中SSTk-ω模型更多考慮近壁面計(jì)算特性，且結(jié)合自適應(yīng)網(wǎng)格形式時(shí)，該模型更適宜作為沖蝕問題湍流計(jì)算方法。

(2)形狀系數(shù)越小的顆粒，由于邊界越粗糙，受曳力較小，其流線軌跡分布越復(fù)雜，在貫穿黏性底層時(shí)能量耗損較大或被底層吸收，即顆粒與壁面間相互作用較弱，壁面磨蝕程度越輕，且該情況下其儲(chǔ)氣能力越強(qiáng)。

(3)商業(yè)軟件對(duì)于各邊界調(diào)節(jié)設(shè)定只能通過設(shè)定形式，故對(duì)模型修改、二次開發(fā)等具有局限性，后期研究可考慮利用如OpenFoam等開源軟件進(jìn)行設(shè)定。