秦浩,汪道兵*,楊凱,張紹良,孫東亮,宇波
1 北京石油化工學院機械工程學院,北京 102617
2 中國石油大學(北京)非常規(guī)油氣科學技術研究院,北京 102249
3 清華大學體育部運動表現(xiàn)與數(shù)據(jù)科學實驗室,北京 100084
干熱巖作為一種新興的清潔可再生能源,可在不受環(huán)境氣候制約的條件下,廣泛用于發(fā)電、供暖等領域[1-2]?,F(xiàn)階段干熱巖地熱能開發(fā)技術的關鍵在于通過水力壓裂技術人工造縫在地下形成人工熱交換系統(tǒng),然而由于干熱巖埋藏很深,基巖溫度較高且地應力各向異性較強,巖石塑性特征增強,常規(guī)水力壓裂技術形成的人工裂隙結構較為單一,地下?lián)Q熱面積不充分因而導致采熱效率不理想。為了解決以上難題,通過向已有人工裂縫內加入暫堵劑對裂縫實施臨時封堵,逼迫裂縫轉向形成多條新裂縫,提升縫網(wǎng)復雜程度的暫堵轉向壓裂技術成為提高干熱巖采熱能力的最具潛力技術之一[3-6]。因此,研究和掌握暫堵劑在干熱巖人工裂隙內的運移規(guī)律是干熱巖暫堵轉向壓裂成敗的關鍵所在。
暫堵劑在干熱巖人工裂隙內的運移過程為典型的固液兩相流動過程,涉及到顆粒運移、聚集及相互碰撞等復雜物理過程,目前的實驗研究因受實驗條件限制難以對暫堵劑在裂隙的這一復雜流動過程進行精準的刻畫。Zhang等借助基于3D打印技術重現(xiàn)的裂縫模型,研究了纖維長度、濃度、注入速度等對暫堵過程的影響,然而實驗過程中無法觀察暫堵劑縫內運移過程[7]。為了實現(xiàn)觀察裂縫中的纖維和顆粒暫堵劑的動態(tài)封堵過程,相關學者開發(fā)了配有高速攝像機的暫堵實驗系統(tǒng),但該系統(tǒng)無法實現(xiàn)高壓阻力,難以獲得真實工況的暫堵劑流動狀態(tài)[8]。因此,數(shù)值模擬是研究暫堵劑在人工裂隙內的運移規(guī)律的有效方法,相比實驗成本也較為低廉。目前針對兩相流的數(shù)值模擬方法主要包括歐拉-歐拉法和歐拉-拉格朗日法兩類[9]。歐拉-歐拉法在計算顆粒-流體兩相流過程中將顆粒作為擬流體與真實流體共同占據(jù)流體單元空間,計算量相對歐拉-拉格朗日法較小,計算過程中可不受顆粒數(shù)量的限制,但基于歐拉-歐拉方法的模型用來模擬裂縫內流體中顆粒的運移機制,因顆粒擬流體化處理造成關鍵顆粒信息丟失,無法真實捕捉顆粒在裂隙內的真實運動及其間相互作用[10-12]。在歐拉-拉格朗日法中,將固體顆粒視為離散相,顆粒離散相計算采用離散單元法,顆粒占據(jù)流體體積真實體現(xiàn)顆粒特性,并且具有豐富的接觸力模型,可精確捕捉暫堵劑在干熱巖人工裂隙內的運動及顆粒間相互作用過程;同時,將流體視為連續(xù)相,計算采用Navier-Stokes方程求解,過程中考慮顆粒與流體間相互作用。目前應用歐拉-拉格朗日法研究暫堵劑在人工裂縫中的運移過程尚處于探索階段,僅有少量學者如Dahi Taleghani等,但其采用的接觸模型較為簡單,難以刻畫顆粒間的相互作用過程,且針對顆粒間相互作用力對運移過程的影響研究相對較少[13-16]。
基以上述分析,本文通過歐拉-拉格朗日描述體系,即流體流動采用歐拉方法計算,顆粒運動通過拉格朗格朗日方法追蹤,建立模擬暫堵劑顆粒在干熱巖人工裂縫內運移過程的CFD-DEM雙向耦合計算模型,該耦合算法可精準捕捉暫堵劑在裂隙內運移過程中的位置、運動速度、接觸力及其他相互作用力等信息,分析暫堵劑攜帶液黏度、顆粒間摩擦系數(shù)、暫堵劑質量濃度以及攜帶液流動狀態(tài)對于暫堵劑縫內運移過程的影響。本論文的研究成果對指導干熱巖暫堵轉向壓裂的暫堵劑用量優(yōu)化具有重要理論及指導意義。
由于CFD-DEM計算量巨大,為不失一般性,本文采用了小尺度的人工裂縫表征單元來模擬顆粒型暫堵劑在人工裂縫內的運移過程,如圖1所示,裂縫長度和高度分別為200 mm和50 mm,入口和出口處的裂縫寬度分別為1 mm和0.5 mm。為了保證CFD-DEM耦合計算結果的準確性和收斂性,將裂縫幾何模型劃分為長、寬(入口/出口)、高單元長度為1 mm×0.25 mm/0.125 mm×0.5 mm的六面體網(wǎng)格單元,共有80 000個網(wǎng)格節(jié)點。
圖1 裂縫網(wǎng)格劃分Fig. 1 Fracture meshing
本模型假設裂縫內的流體為不可壓縮牛頓流體,在裂縫入口與出口處的邊界條件均為恒定壓力條件,裂縫入口壓力和出口壓力分別為30 MPa和29 MPa。裂縫壁面為無滑移邊界,即假設暫堵劑流體在裂縫壁面處的速度為零。
控制單元體的流動滿足質量守恒和動量守恒,不可壓縮流體流動分別滿足流體力學的連續(xù)性方程和Navier-Stokes方程,具體表達式如下[17]:
作為一種典型的拉格朗日方法,離散元法主要根據(jù)牛頓第二定律來描述顆粒的運動速度、加速度及顆粒間相互作用力等,跟蹤粒子的轉動和平動的控制方程分別為[18]:
計算顆粒間法向碰撞力與切向接觸力采用Cundall和Strack提出的線彈性阻尼器(the linear spring-dashpot,LSD)模型[13],該模型適用于多顆粒接觸碰撞是最常用的粘彈性力-位移模型[19-20]。在LSD模型中法向碰撞力是由法向彈性力與法向黏滯力組成。其中彈性力由胡克定律計算,在計算過程中節(jié)約了碰撞動能。而粘滯力與顆粒的相對運動速度成正比,計算過程中產生了碰撞動能耗散。因此,法向碰撞力定義為:
式中:kn為法向彈性剛度因子;ηn為法向阻尼系數(shù);νrn為法向相對速度。
同樣,切向碰撞力是由切向彈性力與切向黏滯力組成,計算過程中假設在具有恒定切向彈性剛度的碰撞顆粒接觸區(qū)域沒有微滑移,將切向碰撞力定義為:
式中:kt為切向彈性剛度因子;ηt為切向阻尼系數(shù);νrt為切向相對速度。
由于顆粒在裂隙內的運動方式與流體流動過程二者相互影響,困此本文采用流體流動與顆粒運動的雙向耦合算法對固液兩相模型進行精準求解計算。如圖3所示,為CFD-DEM耦合算法的程序流程圖,主要包括耦合計算,DEM循環(huán)計算以及CFD計算3個模塊。計算開始前首先對耦合計算模塊、DEM模塊及CFD模塊進行初始化;模型參數(shù)初始化后進入耦合計算模塊,根據(jù)顆粒位置和流體網(wǎng)格信息計算每個流體單元孔隙度,并通過顆粒速度、流體壓力及速度等得到每個流體單元中的流體-顆粒相互作用力;然后,將耦合計算模塊的數(shù)據(jù)傳遞給DEM循環(huán)計算模塊,通過DEM模塊循環(huán)計算后得到下一個流體時間步長中所有粒子的新位置、平動速度和轉動速度;最后,將DEM計算結果傳遞給CFD計算模塊,求解流體相質量方程、動量方程,從而得到流體的速度場與壓力場。判斷是否達到計算時間,若未完成則將計算所得流體及顆粒數(shù)據(jù)傳遞至耦合模塊進行下一循環(huán)計算;若完成計算,則輸出最終計算結果。
圖3 CFD-DEM耦合總體框架Fig. 3 The overall framework for CFD-DEM coupling
為了確保上述CFD與DEM雙向耦合同步進行和數(shù)值解的交換,在DEM循環(huán)計算中迭代次數(shù)m為流體相計算時間步長Δtf與顆粒相計算時間步長Δtp的應滿足如下比值關系式:
圖2 顆粒i, j相互碰撞作用力示意圖Fig. 2 Schematic diagram of collision force between particles i and j
為了驗證本文中模型解的可靠性,將本文獲得的數(shù)值模擬所得到的出口流量變化的數(shù)據(jù)與文獻[14]中的結果進行對比,參考文獻中模型數(shù)據(jù)如下表1所示,其中文獻中數(shù)據(jù)包括暫堵前和暫堵后兩部分。通過圖4中數(shù)據(jù)對比可以看出,Vinj/Vfrac小于0.1的范圍內,參考數(shù)據(jù)與數(shù)值模擬中出口流量基本相同,最大誤差小于4.6%,而當Vinj/Vfrac大于0.1后文獻中模擬裂縫內堵劑顆粒逐漸出現(xiàn)粘結局部封堵,使得出口流量明顯小于本文模擬結果,但兩者總體變化趨勢基本相同。為了進一步驗證模型的可靠性,本文通過該模型根據(jù)文獻[15]中支撐劑顆粒在裂縫內的運移及沉積實驗過程構建了相應的數(shù)值模型,如圖5所示為顆粒運移過程后期,顆粒沉積及分布數(shù)值模擬結果與實驗結果對比,對比結果表明,數(shù)值模擬結果與實驗中顆粒分布結果基本吻合。綜上所述,數(shù)值模型模擬結果可靠,可用作模擬暫堵劑顆粒在干熱巖人工裂縫內的運移過程。
表1 參考文獻相關屬性Table 1 References related attributes
圖4 數(shù)值模擬結果與參考模型出口流量對比圖Fig. 4 The comparison of the numerical simulation results with the reference model
圖5 顆粒沉積及分布數(shù)值模擬結果與實驗結果對比圖Fig. 5 The Comparison of numerical simulation results and experimental results of particle deposition and distribution
本文針對暫堵劑攜帶液黏度、顆粒間摩擦系數(shù)、暫堵劑質量濃度等因素對于暫堵劑運移過程的影響,開展了相關的數(shù)值模擬研究。首先對暫堵劑攜帶液黏度分別在0.03 mPa·s、30 mPa·s、60 mPa·s、90 mPa·s及120 mPa·s條件下的暫堵劑運移過程進行分析;隨后通過對靜摩擦系數(shù)及滾動摩擦系數(shù)分別為1/0.5、0.8/0.4、0.6/0.3、0.4/0.2的情況下裂縫內暫堵劑顆粒間的相互作用進行研究,分析摩擦系數(shù)對運移過程的影響;通過對比暫堵劑顆粒質量濃度在4%、5%、6%、7%及8%條件下顆粒間相互作用力變化及顆粒運動狀態(tài),分析顆粒質量濃度對縫內暫堵劑輸運過程的影響。最后對湍流作用下裂縫內暫堵劑顆粒運動及相互作用進行分析。其中暫堵劑顆粒、攜帶液具體參數(shù)及干熱巖壁面相關參數(shù)如下表2所示。
表2 相關屬性Table 2 The relevant properties(a)攜帶液相關參數(shù)
由于雷諾數(shù)大于4000,本模型采用RNGk-ε湍流模型。考慮到流體流動及顆粒運移時間會因流體性質、顆粒尺寸以及裂縫幾何形狀尺寸等因素的變化而產生急劇變化,為此我們采用無量綱時間對相關結果進行分析來避免其影響。定義特征時間τ=d/g對計算時間進行無量綱化處理,其中d為顆粒直徑,g為重力加速度。為了確保CFD與DEM雙向耦合同步進行和數(shù)值解的交換,本文在模擬中選用CFD時間步長為5×10-5s,DEM時間步長為2.5×10-6s,每個流體時間步內,DEM循環(huán)計算中迭代20次。
干熱巖人工裂隙內存在明顯的干熱巖壁面-攜帶液流體-暫堵劑顆粒間換熱過程。如圖6所示為干熱巖壁面溫度為473 K時,裂隙內攜帶液流體與壁面換熱后縫內流體溫度分布及暫堵劑顆粒在與攜帶液流體換熱后顆粒溫度分布??p內攜帶液流體在與干熱巖壁面換熱后溫度逐漸提升,但由于換熱距離較小,導致出口處流體溫度僅升高4.5 ℃;而縫內暫堵劑在受到流體溫度及壁面溫度共同作用后,溫度也出現(xiàn)相應提升,同樣因換熱面積較短且暫堵劑導熱系數(shù)相對較小,使得顆粒溫度未出現(xiàn)較大幅度提升,僅為流體溫度增大量的11.1%。在該模型的基礎上分別對帶液黏度、暫堵劑質量濃度、攜帶液流動狀態(tài)和顆粒間摩擦系數(shù)等因素對干熱巖人工裂隙內暫堵劑運移規(guī)律進行分析。
圖6 攜帶液流體及暫堵劑顆粒溫度分布Fig. 6 The Temperature distribution of carrying fluid and temporary plugging agent particles
4.2.1暫堵劑攜帶液黏度
圖7為顆粒平均速度隨暫堵劑攜帶液黏度變化圖。從圖中可以看出,隨攜帶液黏度的增大,縫內暫堵劑顆粒平均運動速度逐漸減??;顆粒進入裂縫初期,受縫內流體作用,其平均速度迅速提升,當流體攜帶暫堵劑顆粒充滿整個流域后,暫堵劑顆粒運動速度趨于穩(wěn)定;當攜帶液黏度由0.03 mPa·s提升至120 mPa·s的過程中,暫堵劑在縫內的平均運動速度降低超過77%。
圖7 顆粒平均速度隨攜帶液黏度變化圖Fig. 7 The diagram of average particle velocity changing with viscosity of carrier fluid
為了進一步研究攜帶液黏度對于暫堵劑顆粒速度的具體影響,對暫堵劑穩(wěn)定運動狀態(tài)下攜帶液流速與暫堵劑顆粒運動進行分析。從圖8中可以看出,當縫內流動趨于穩(wěn)定,因裂縫形狀及縫內攜帶液壓力變化影響,由入口至出口流動過程中,攜帶液流速不斷上升。受縫內攜帶液流動的作用,暫堵劑顆粒在入口至出口運動過程中速度不斷增大。隨攜帶液黏度的變化,其顆粒速度變化幅度發(fā)生相應的減小。在黏度為0.03 mPa·s時,顆粒速度由入口時的1.2 m/s增大至28.84 m/s左右,速度增大超過27 m/s;而當流體黏度增大至30 mPa·s時,出口處顆粒速度僅為19.6 m/s左右,增大18.4 m/s左右。同樣,當黏度增大至60 mPa·s、90 mPa·s時,暫堵劑顆粒在出口處速度分別增大12 m/s、8.5 m/s左右;而當流體黏度為120 mPa·s時,出口處顆粒速度僅為6.523 m/s??梢缘贸?,隨流體黏度的增大,暫堵劑顆粒在裂縫內流動的過程中,速度增幅將逐漸降低,出口處顆粒速度減小。
圖8 穩(wěn)定運動狀態(tài)下流體流速與暫堵劑顆粒運動速度云圖Fig. 8 The cloud diagram of fluid velocity and velocity of temporary plugging agent particles in stable motion state
攜帶液黏度在影響暫堵劑顆粒運動的同時,對攜帶液的流動過程也存在一定的影響。如圖9所示,出口流量由黏度為0.03 mPa·s時的0.68 m3/s下降至120 mPa·s時0.14 m3/s左右。出口流量隨攜帶液黏度增大不斷減小的主要原因是黏度增大造成裂縫內攜帶液流動阻力增大,流速明顯下降,相同時間內通過裂縫出口流出的攜帶液總量減少。同時可以看出攜帶液黏度在0.03 mPa·s、30 mPa·s變化時裂縫出口流量在注入流體體積與裂縫體積之比Vinj/Vfrac=2附近達到峰值,而后保持相對穩(wěn)定出口流量;但隨著黏度逐漸增大,峰值區(qū)間不斷向前推移,黏度在60 mPa·s、90 mPa·s 、120 mPa·s時,出口流量峰值分別出現(xiàn)在Vinj/Vfrac=1.5、Vinj/Vfrac=1.3、Vinj/Vfrac=1。
圖9 出口流量隨注入流體體積與裂縫體積之比變化圖Fig. 9 The diagram of outlet flow rate as the ratio of injected fluid volume to fracture volume
4.2.2顆粒間摩擦系數(shù)
暫堵劑顆粒在裂縫內隨攜帶液運移的過程中,顆粒相互之間因接觸而存在相互作用,顆粒與顆粒之間的相互作用主要受摩擦系數(shù)的影響,為了研究顆粒間摩擦系數(shù)對于顆粒型暫堵劑在裂縫內運移過程的影響,對不同暫堵劑質量濃度下的運移過程進行了分析。圖10(a)中顯示的是在暫堵劑質量濃度為6%時顆?!w粒間靜/滾動摩擦系數(shù)分別為1/0.5、0.8/0.4、0.6/0.3、0.4/0.2、0.2/0.05情況下暫堵劑顆粒間相互作用力的變化,從圖中可以看出,裂縫內顆粒處于穩(wěn)定運動狀態(tài)時靜/滾動摩擦系數(shù)變化的過程中,顆粒間相互作用力基本保持在1.4~1.6×10-4N,說明顆?!w粒間摩擦系數(shù)對顆粒間相互作用力的影響較小。
為了進一步研究顆粒間摩擦系數(shù)對暫堵劑顆粒運動的影響,對各靜/滾動摩擦系數(shù)下顆粒平均運動速度進行分析,如圖10(b)顆粒平均運動速度結果表明在靜/滾動摩擦系數(shù)由1/0.5減小至0.2/0.05的過程中,裂縫內顆粒平均運動速度均保持在5.8 m/s左右,變化較小。從圖11中可以看出隨著靜/滾動摩擦系數(shù)的減小,顆粒間相互作用力變化相對較小的同時,裂縫內暫堵劑顆粒的角速度逐漸增大,由325 rad/s左右提升至600 rad/s左右,因為顆粒的角速度變化與顆粒間的滾動摩擦相關,如式9所示,角速度與滾動摩擦系數(shù)成反比。在相同的滾動摩擦力的作用下,滾動摩擦系數(shù)越小,其角速度則越大。
圖10 不同靜/滾動摩擦系數(shù)下顆粒間相互作用力及顆粒平均運動速度圖Fig. 10 Interaction force and average velocity of particles under different static/rolling friction coefficients
圖11 顆粒角速度與顆粒間相互作用力隨摩擦系數(shù)變化圖Fig. 11 The diagram of angular velocity and interaction force of particles varying with friction coefficient
其中:μr為滾動摩擦系數(shù);R為接觸點到質心的距離;Fn為法向力;ω為顆粒在接觸點的單位角速度矢量。
通過上述分析可以看出,顆粒間摩擦系數(shù)對暫堵劑顆粒在裂縫內運動速度的影響較小,僅滾動摩擦系數(shù)較小時,顆粒角速度發(fā)生改變,但對裂縫全段內顆粒的運移過程并未產生影響,因此顆粒間摩擦系數(shù)對于裂縫內顆粒運動的影響較小。
4.2.3暫堵劑質量濃度
圖12顯示的是4%、5%、6%、7%以及8%顆粒質量濃度下的運移過程中縫內顆粒作用力的變化,從圖中可以看出隨著質量濃度的增大,其作用力也將隨之更大。暫堵劑最初進入裂縫入口時,由于入口段裂縫寬度大,暫堵劑顆粒初始速度相對較低,顆粒無法快速向前運動,出現(xiàn)前端堆積,導致顆粒間的接觸碰撞增加,顆粒間相互作用力較大,如圖13所示,隨后暫堵劑顆粒不斷隨攜帶液向裂縫出口運動,當t/τ=3.5后裂縫內暫堵劑顆粒隨攜帶液運動充滿整個裂縫,達到較為穩(wěn)定的運移狀態(tài);與暫堵劑顆粒最初進入裂縫階段相比,暫堵劑達到穩(wěn)定運動狀態(tài)后,顆粒間相互作用力減小超過12.64%。
圖12 不同質量濃度下顆粒間相互作用力變化Fig. 12 The interaction force between particles changes at different mass concentrations
圖13 暫堵劑顆粒裂縫前端堆積顆粒運動速度與相互作用力Fig. 13 The velocity and interaction force of the stacked particles at the front of the fracture of temporary plugging agent particles
圖14 暫堵劑穩(wěn)定運動狀態(tài)下相互作用力降幅Fig. 14 The decrease of interaction force under the stable movement of temporary plugging agent
暫堵劑顆粒質量濃度由4%增大至8%的過程中,縫內暫堵劑顆粒間相互作用力不斷增大,當暫堵劑運動達到相對穩(wěn)定的狀態(tài)后,質量濃度增大一倍,顆粒間的相互作用力增大超過43.8%,可以說明質量濃度對于暫堵劑縫內運移過程具有重要的影響。如圖15所示,穩(wěn)定運動狀態(tài)下,暫堵劑質量濃度的增大,顆粒運動速度降低,顆粒間相互作用力增大,主要是由于質量濃度增大后,相同裂縫尺度內顆粒間相互接觸碰撞增多,一定程度上減緩了顆粒運動速度,有助于后期暫堵劑在縫內的形成有效封堵。但過高的濃度會致使實際施工困難并帶來材料的浪費,因此在實際施工過程中應根據(jù)實際的工況,確定相應的暫堵劑質量濃度。
圖15 暫堵劑穩(wěn)定運動狀態(tài)下顆粒間相互作用力與運動速度的關系Fig. 15 The relationship between particle interaction force and velocity under the stable movement of temporary plugging agent
4.2.4攜帶液流動狀態(tài)
攜帶液在裂縫內流動時易出現(xiàn)明顯的湍流現(xiàn)象,與正常流動情況下相比,湍流區(qū)域內攜帶液流速、暫堵劑顆粒運動及顆粒間相互作用等方面產生較大變化。如圖16所示,在紅色圓圈內分別標注暫堵劑質量濃度為6%情況下,裂縫內局部區(qū)域出現(xiàn)湍流作用時,顆粒速度、湍流動能及顆粒間相互作用力的變化。與裂縫內正常流動區(qū)域相比,湍流區(qū)內湍流動能超過9.5 m2/s2,出現(xiàn)明顯增大,局部區(qū)域湍流導致攜帶液流動速度的增大使得該區(qū)域內顆粒速度增大超過50倍。如圖16(c)所示,湍流造成顆粒局部旋渦流動,造成大量顆粒間相互接觸碰撞,出現(xiàn)顆粒間相互作用力突增,達到1.29 N,較正常流動情況下增大超過150倍。由上述分析可以得出,在裂縫內攜帶液的流動方式對暫堵劑的運移過程具有較大的影響。
圖16 湍流情況下顆粒速度、攜帶液湍流動能及顆粒間相互作用力云圖Fig. 16 The cloud diagram of particle velocity, turbulent kinetic energy of carrying fluid and interaction force between particles in turbulent condition
(1)本文建立了顆粒型暫堵劑在人工裂隙內運移過程的CFD-DEM耦合模型,數(shù)值模擬結果表明影響干熱巖人工裂隙中暫堵劑顆粒運移的主要因素為暫堵劑攜帶液黏度、暫堵劑質量濃度以及攜帶液的流動狀態(tài)。
(2)顆粒間摩擦系數(shù)對運移過程中暫堵劑顆粒間相互作用力及顆粒運動速度影響較小,分析造成這種現(xiàn)象的原因可能是顆粒間因摩擦產生的作用力相對較小,難以克服流體流動對顆粒速度的影響。
(3)攜 帶 液 黏 度 由0.03 mPa·s提 升 至120 mPa·s的過程中,暫堵劑在縫內的平均運動速度降低超過77%。攜帶液黏度增大,暫堵劑顆粒在裂縫內運動速度增幅逐漸降低,出口處顆粒速度及攜帶液流量明顯減小。
(4)隨暫堵劑質量濃度的增大,顆粒間相互作用力不斷增大,暫堵劑運動達到相對穩(wěn)定的狀態(tài)后,顆粒質量濃度增大一倍,顆粒間的相互作用力增大超過43.8%。當裂縫內攜帶液出現(xiàn)局部湍流現(xiàn)象時,暫堵劑顆粒運動速度與顆粒間相互作用力明顯發(fā)生變化,顆粒間相互作用較力正常流動情況下增大超過150倍。