郭旭洋，金 衍，林伯韜，盧運虎，訾敬玉
(中國石油大學(xué)(北京)石油工程學(xué)院，北京 102249)
天然氣水合物儲量巨大，并具有清潔高效的特點(diǎn)，在中國國南海已經(jīng)開(kāi)展2次試采，其中第2次為水平井試采[1]。水平井降壓開(kāi)采過(guò)程中的工程擾動(dòng)會(huì )誘發(fā)水合物分解，導致井周沉積物力學(xué)性質(zhì)劣化，井周應力時(shí)空演化特征復雜，存在出砂和海底滑坡的可能性[2]。因此明確試采過(guò)程中的沉積物變形破壞機制和流-固-熱-化多場(chǎng)耦合機制對于南海水合物安全、高效鉆采具有意義。水合物沉積物膠結性較差，水合物的分解對其力學(xué)性質(zhì)影響顯著(zhù)。研究人員采用人工合成樣品和現場(chǎng)保壓取芯等手段開(kāi)展了大量三軸力學(xué)試驗，確定多種水合物沉積物試樣的強度參數和變形特性，為現場(chǎng)尺度建模分析提供關(guān)鍵力學(xué)參數。水合物和沉積物基質(zhì)共同決定水合物沉積物的力學(xué)性質(zhì)，但是基質(zhì)占據主導作用。增加水合物飽和度或減小有效圍壓均能夠使沉積物應力-應變關(guān)系由應變硬化轉向應變軟化。沉積物內聚力受水合物飽和度影響作用明顯，而內摩擦角僅在水合物飽和度較高時(shí)受水合物飽和度影響較大[3-9]。針對中國南海北部深水淺層天然氣水合物沉積物的力學(xué)試驗顯示，有效圍壓增加的情況下，沉積物性質(zhì)會(huì )呈塑性，有效圍壓和水合物飽和度綜合影響沉積物力學(xué)性質(zhì)，導致水合物沉積物鉆采過(guò)程中地層失穩風(fēng)險較大[10]。力學(xué)試驗獲取的水合物沉積物力學(xué)性質(zhì)為更大尺度的水合物沉積物鉆采過(guò)程建模分析提供了基礎數據支持，結合水合物分解動(dòng)力學(xué)、質(zhì)量守恒和能量守恒等，能夠對水合物沉積物鉆采過(guò)程中的多種物理場(chǎng)演化開(kāi)展定量建模分析。通過(guò)有限元建模分析發(fā)現，較小的井眼半徑更易導致水合物水平井發(fā)生井壁失穩，這種影響隨鉆井時(shí)間加劇[11]。水合物相變、傳熱傳質(zhì)和地層變形的耦合建模有助于確定水合物分解前緣和解釋出砂和穩產(chǎn)時(shí)間短等現象，為商業(yè)化開(kāi)發(fā)提供理論參考，而生產(chǎn)壓力、滲透率、初始飽和度和溫度等都與地層變形有一定聯(lián)系[12-13]?；赥OUGH等平臺，能夠定量表征水合物開(kāi)采過(guò)程中的傳熱傳質(zhì)問(wèn)題和水合物分解規律[14]。在此基礎上通過(guò)孔隙彈性力學(xué)理論進(jìn)行流固耦合研究，認為固體力學(xué)場(chǎng)的耦合方法對模擬結果會(huì )產(chǎn)生顯著(zhù)影響。全耦合對收斂性和準確性的提升有較大幫助，而如果采用順序耦合則需要更嚴密的耦合條件以確保數值求解的穩定性和準確性[15-17]。筆者根據南海水合物沉積物力學(xué)特征，建立水平井降壓開(kāi)采過(guò)程的流-固-熱-化多場(chǎng)耦合有限元模型，分析沉積物內的變形、破壞和海底沉降規律。
以水合物沉積物為對象構建多物理場(chǎng)耦合模型，表征沉積物內水合物分解、熱傳遞、氣水滲流和變形破壞過(guò)程。根據質(zhì)量守恒，多孔介質(zhì)內的流動(dòng)過(guò)程可以表示為
(1)
式中，φ為孔隙度；
Si為流體飽和度；
ρi為流體密度，kg/m3；
vi為流體速度，m/s；
si為匯/源質(zhì)量速度，kg/(m3·s)；
下角i代表相態(tài)類(lèi)型。
沉積物中流體滲流速度低，根據達西定律[18]，多孔介質(zhì)內流體流速表示為
(2)
式中，k為絕對滲透率，m2；
Kri為相對滲透率；
μi為黏度，Pa·s；
pi為孔隙流體壓力，Pa；
g為重力加速度，m/s2。
水合物降壓開(kāi)采過(guò)程中的傳熱問(wèn)題需要考慮沉積物固相骨架內的熱傳導、多孔介質(zhì)內流體的熱對流和水合物分解產(chǎn)生的吸熱過(guò)程[19]。該過(guò)程中的傳熱問(wèn)題可以表示為
(hC+hE)=Q.
(3)
式中，E為比內能，J/kg；
S為飽和度；
h為對流傳熱項；
Q為熱匯/源，J/(m3·s)；
下角L、G、H分別代表液相、氣相和水合物相3種相態(tài)。
降壓開(kāi)采過(guò)程中，水平井井底壓力降低，導致天然氣水合物穩定區域平衡被破壞，水合物開(kāi)始分解[20-21]。根據水合物分解動(dòng)力學(xué)模型，降壓過(guò)程中水合物分解速率表示為
RMH=-kdMMH(pe-pG)As.
(4)
式中，RMH為分解速率；
kd為反應速率，Pa·s；
MMH為水合物摩爾質(zhì)量，kg/mol；
pe為相平衡壓力，Pa；
pG為氣相壓力，Pa；
As為反應比表面積，1/m2。
由于溫度和壓力的變化會(huì )導致多孔介質(zhì)及多孔介質(zhì)內流體發(fā)生體積變化，具體可以表示為
(5)
(6)
式中，χ為壓縮系數，1/Pa；
β為膨脹系數，1/K；
T為溫度，K。
由于降壓，沉積物內水合物分解，水合物飽和度下降，導致沉積物力學(xué)性質(zhì)出現劣化[22]。
水合物分解會(huì )改變內聚力，c可以表示為水合物飽和度的函數
(7)
式中，c0為不考慮水合物部分的沉積物內聚力，Pa；
θ內摩擦角，(°)；
SMH為水合物飽和度；
α和β為系數。
為表征降壓開(kāi)采過(guò)程中的沉積物力學(xué)響應，需要根據動(dòng)量平衡對應力張量建立表達式
(8)
式中，σ為總應力張量；
ρb為體積密度，kg/m3。
根據無(wú)窮小變換，應變張量ε和位移矢量u的關(guān)系可以表示為
(9)
在多孔介質(zhì)中，本構關(guān)系表示為
δσ=C:δ(ε-εp)-αδpI.
(10)
式中，C為彈性張量；
εp為塑性應變；
α為耦合系數；
I為二階單位張量。
利用Mohr-Coulomb和Drucker-Prager模型表征彈塑性變形下的剪切破壞。屈服函數f和Drucker-Prager模型的塑性勢函數g可分別表示為
(11)
(12)
式中，I1為有效應力第一不變量；
J2為有效偏應力第二不變量；
βf、κf、βg和κg為破壞包絡(luò )線(xiàn)系數[23-25]。
Mohr-Coulomb模型可表示為
(13)
(14)
(15)
(16)
根據上述控制方程，建立水合物沉積物降壓開(kāi)采模型。由于水平井鉆采是南海試采工作的一項重要內容，也是提高產(chǎn)氣效率的關(guān)鍵手段，因此建立單個(gè)水平井模型和雙水平井模型，其中雙水平井的水平段深度不同。圖1為單個(gè)水平井筒降壓開(kāi)采模型和雙水平井筒降壓開(kāi)采模型。由于水平井筒沿x方向貫穿模型，井筒長(cháng)度較長(cháng)，可以采用平面應變假設簡(jiǎn)化數值運算，僅分析y-z截面內的二維流-固-熱-化耦合問(wèn)題。對于單井模型，井筒位于y-z截面正中(25 m，25 m)；
對于雙井模型，井筒分別位于y-z截面的(25 m，15 m)和(25 m，35 m)。對于二維y-z平面模型，y方向為SHmax，z方向為Sv，x方向施加Shmin。在厚度為50 m的沉積物層內，假設地層壓力梯度和地應力梯度的影響較小，地層壓力為14 MPa，垂向主應力為16.6 MPa，水平最大主應力為14.76 MPa，不隨深度而增大。沉積物密度為2 300 kg/m3，初始飽和度為40%，孔隙度為0.15，滲透率為1.8×10-2μm2，地層溫度為284.15 K，干濕導熱系數分別為1和3 W/(m·K-1)，降壓生產(chǎn)壓力為3 MPa。
圖1 單個(gè)水平井筒和雙水平井筒降壓開(kāi)采模型Fig.1 Sketch of models for depressurization in a single horizontal well and two horizontal wells
根據圖1所示模型，開(kāi)展水平井筒降壓開(kāi)采過(guò)程中的數值模擬，分別針對單個(gè)水平井筒和雙水平井筒開(kāi)展分析，模擬降壓開(kāi)采時(shí)長(cháng)為15 d。降壓開(kāi)采模型能夠表征該過(guò)程中井周沉積物內的水合物分解范圍、分解吸熱引發(fā)溫度變化、孔隙壓力降低和彈塑性演化，并能夠預測水平井產(chǎn)能。在二維y-z平面內，由于數值解關(guān)于z軸對稱(chēng)，僅展示y軸正半軸數值結果。
圖2為采用單個(gè)水平井和2口水平井降壓開(kāi)采5和15 d時(shí)的井周孔隙流體壓力分布。近井區域壓力降最明顯，這是由井筒降壓開(kāi)采使用的較低生產(chǎn)壓力造成的。隨著(zhù)距離向遠場(chǎng)移動(dòng)，壓力越接近初始壓力。隨著(zhù)降壓開(kāi)采時(shí)間的增加，壓降前緣的波及范圍更廣。對于單井模型，降壓開(kāi)采5 d時(shí)遠場(chǎng)壓力仍等于初始壓力14 MPa，但開(kāi)采15 d后遠場(chǎng)壓力也開(kāi)始下降。相對應地，采用雙井模型能夠增加壓降前緣的波及范圍。對比開(kāi)采5 d時(shí)的單井模型和雙井模型的孔隙流體壓力結果，發(fā)現雙井模型井周壓降范圍更大，壓降更加劇烈。開(kāi)采15 d后雙井模型的遠場(chǎng)區域壓力普遍降至低于11 MPa，單井模型原廠(chǎng)區域壓力仍可達13 MPa，說(shuō)明采用兩口水平井降壓開(kāi)采能夠更有效地建立地層內的壓降，為甲烷流動(dòng)提供促進(jìn)作用。
圖2 單個(gè)水平井筒和雙水平井筒降壓開(kāi)采5和15 d孔隙流體壓力分布Fig.2 Pore pressure distribution after 5 and 15 d of depressurization for single- and two-well scenarios
圖3為單井模型和雙井模型內降壓開(kāi)采5和15 d時(shí)的水合物飽和度分布特征，水合物飽和度為0的區域表示水合物完全分解。結果顯示水合物分解區域近似呈圓形，這是因為y-z平面內的井周壓降前緣也近似呈圓形，壓降導致水合物穩定區的相平衡遭到破壞而出現水合物分解，分解區域受到壓降前緣影響較為直接。根據圖3數據，5 d時(shí)單井水合物分解前緣距水平井筒擴展1.01 m，雙井為0.96 m；
15 d時(shí)水合物分解前緣擴展的距離在單井模型中為1.74 m，在雙井模型中則為1.58 m。說(shuō)明增加水平井數量和增加降壓開(kāi)采時(shí)間均會(huì )對水合物分解范圍產(chǎn)生影響。隨著(zhù)降壓開(kāi)采時(shí)間的增加，水合物分解前緣向遠場(chǎng)擴展，15 d內能夠擴展1.58～1.74 m。
圖3 單個(gè)水平井筒和雙水平井筒降壓開(kāi)采5 d和15 d水合物飽和度分布Fig.3 Hydrate saturation distribution after 5 and 15 d of depressurization for single- and two-well scenarios
圖4為單井模型和雙井模型降壓開(kāi)采5和15 d后的溫度分布。由于水合物分解是吸熱過(guò)程，導致分解區域及鄰近沉積物內溫度降低。降溫前緣也近似呈圓形，與壓降前緣和水合物分解前緣形狀類(lèi)似。造成這一現象的原因除了壓降與水合物分解直接相關(guān)外，也與式(1)、(3)有關(guān)，兩者均為拋物型微分方程，呈現的傳質(zhì)傳熱特征具有類(lèi)似性。圖4表明，單井模型中的降溫前緣在5 d時(shí)為1.25 m，在15 d時(shí)為2.71 m；
雙井模型中的降溫前緣在5 d時(shí)為1.21 m，在15 d時(shí)為2.67 m。這說(shuō)明開(kāi)采時(shí)間仍是影響降溫范圍的主要因素，水平井數量對降溫范圍的影響相對較小。
圖4 單個(gè)水平井筒和雙水平井筒降壓開(kāi)采5和15 d溫度分布Fig.4 Temperature distribution after 5 and 15 d of depressurization for single- and two-well scenarios
圖5 水平井累積產(chǎn)氣對比Fig.5 Comparison of cumulative gas production between studied scenarios
根據降壓開(kāi)采多物理場(chǎng)耦合模擬可以進(jìn)行水平井產(chǎn)氣預測(圖5)。對比了單井模型水平井產(chǎn)量、雙井模型井1(z=35 m)產(chǎn)量、雙井模型井2(z=15 m)產(chǎn)量和雙井模型兩口水平井平均產(chǎn)量?？傮w累產(chǎn)曲線(xiàn)顯示單井累產(chǎn)氣差異不明顯，這是因為兩口水平井井距較遠，盡管井周壓降區域互相波及，但水合物分解區域和溫度降低區域仍相隔一段距離，而水平井產(chǎn)氣主要是由水合物分解區釋放的甲烷提供，因此在井周水合物分解范圍相似的情況下，獲得的水平井累產(chǎn)氣結果也較為近似。比較單井模型和雙井模型內兩口井的累積產(chǎn)氣曲線(xiàn)，發(fā)現15 d內累產(chǎn)氣由高到低排序為雙井模型井1>單井模型>雙井模型井2。分析發(fā)現這一關(guān)系能夠和三口水平井井位深度建立起聯(lián)系，累產(chǎn)氣最高得雙井模型井1深度最小，單井模型其次，雙井模型井2深度最大。由于在y-z平面內多場(chǎng)耦合模型考慮了重力的作用，導致降壓開(kāi)采后深度更大的區域壓力相對更高、壓降相對更弱。據此，在深度更大的區域進(jìn)行水平井降壓開(kāi)采形成的壓降差異會(huì )略小于在深度較小的區域進(jìn)行降壓開(kāi)采，一定程度上減緩水合物分解和甲烷釋放，減小了對應水平井累產(chǎn)氣。此外，還計算了雙井模型兩口井累產(chǎn)的平均值，發(fā)現其數值略高于單井模型。這說(shuō)明在本研究設定的條件下，增加水平井數量能夠適當提升單井累產(chǎn)效果。4條累產(chǎn)曲線(xiàn)均呈現出開(kāi)采初期產(chǎn)氣效率很高、開(kāi)采一段時(shí)間后產(chǎn)量急劇下降的特點(diǎn)，與前期國內外部分試采工作結果具有一定的一致性。
降壓開(kāi)采誘發(fā)的沉積物變形與破壞也是需要研究的重點(diǎn)之一。圖6為采用單井模型和雙井模型降壓開(kāi)采5和15 d后的沿垂直方向的總應力分布情況。由于近井區域受到降壓開(kāi)采影響，孔隙流體壓力下降明顯，導致垂直方向總應力在井筒和近井區域處于較低值。在同一時(shí)間，增加水平井數量進(jìn)一步加劇了垂直方向總應力的降低。隨著(zhù)開(kāi)采時(shí)間增加，垂直方向總應力降低的區域逐漸增大，而且遠場(chǎng)區域也出現一定的垂向總應力降低。通過(guò)與圖2對比，發(fā)現垂直方向總應力的演化規律與孔隙流體壓力額演化規律具有一定關(guān)聯(lián)性，只是因為在孔隙彈性力學(xué)理論中，孔隙流體壓力是總應力中各正應力分量的重要組成部分。在15 d時(shí)，雙井模型內遠場(chǎng)垂向總應力比單井模型更低，也是由雙井模型遠場(chǎng)孔隙流體壓力更低造成的。
圖6 單個(gè)水平井筒和雙水平井筒降壓開(kāi)采5 d和15 d垂向應力分布Fig.6 Vertical stress distribution after 5 and 15 d of depressurization for single- and two-well scenarios
表征降壓開(kāi)采誘發(fā)的彈塑性變形能夠幫助預測氣井出砂和海底沉降。圖7為單井模型水平井、雙井模型中井1和井2同深度的塑性變形時(shí)間演化規律。由圖7(a)看出，在降壓開(kāi)采起始時(shí)，僅在井壁處存在塑性應變，儲層沉積物均無(wú)塑性應變發(fā)生。降壓開(kāi)采0.5～1 d后，由井壁至沉積物層內出現隨距離逐漸減小的塑性應變，這是由于近井區域壓降劇烈，沉積物固體骨架變形程度更大。0.5～1 d時(shí)，由于水合物分解尚不明顯，水合物飽和度對塑性變形的影響不顯著(zhù)。降壓開(kāi)采5 d后，由于水合物分解前緣的擴展，分解區向未分解區過(guò)渡的塑性變形呈現出突變特征，分解區展現的塑性應變高于未分解區。從0.5～15 d的結果顯示，水解物分解前緣和塑性應變突變點(diǎn)具有較強相關(guān)性，分解前緣外的未分解儲層內的塑性應變則由近及遠平滑降低。該結果說(shuō)明水合物分解區塑性變形特點(diǎn)明顯，更易出砂和誘發(fā)地層失穩。
圖7 不同模型井塑性變形時(shí)間演化Fig.7 Temporal evolution of plastic deformation in different model wells
由圖7(b)、(c)看出雙井模型中井1和井2各時(shí)間點(diǎn)的塑性變形分布特征與單井模型結果比較接近，均呈現出水合物分解前緣與塑性應變突變點(diǎn)吻合的情況。5、10和15 d的塑性應變在10 m以外出現小幅度升高的現象，這是受到了模型的邊界條件的影響。在y-z平面的右側邊界設置為固定位移邊界，雙井模型的壓降和變形程度和范圍更大，導致近井區域壓縮變形程度較高，由于右側邊界為固定邊界，在靠近右側邊界的部分區域出現拉伸變形。說(shuō)明在本研究中模型尺寸和邊界條件對塑性變形求解過(guò)程有影響。
(1)由于降壓開(kāi)采中的傳熱和傳質(zhì)過(guò)程控制方程均具有拋物型微分方程的特征，沉積物層內的孔隙壓力、溫度和水合物分解前緣的時(shí)空演化規律具有較強的相關(guān)性。
(2)沉積物層內的總應力時(shí)空演化規律受孔隙壓力演化的影響明顯，近井區域的總應力數值小于遠場(chǎng)受壓降影響較小區域的總應力數值。這一現象能夠通過(guò)孔隙彈性力學(xué)理論進(jìn)行解釋。
(3)在水平井降壓開(kāi)采的多場(chǎng)耦合研究時(shí)，重力作用對水平井產(chǎn)能、力學(xué)響應和孔隙壓力分布的影響較為明顯。降壓開(kāi)采后，重力作用導致深層水平井壓降程度低于淺層水平井，一定程度上抑制水合物分解前緣的擴展。本研究數模顯示，淺層水平井累積產(chǎn)氣量更高。
(4)水合物分解會(huì )明顯地劣化沉積物力學(xué)性質(zhì)，誘發(fā)較強的塑性應變。在降壓開(kāi)采初期水合物分解前緣未形成時(shí)，沉積物層塑性應變不明顯；
隨著(zhù)降壓開(kāi)采時(shí)間的增加，水合物分解區力學(xué)性質(zhì)劣化顯著(zhù)，塑性應變明顯高于未分解區,內聚力等強度參數變低,在水合物分解前緣附近的塑性應變會(huì )出現明顯的階梯式下降特征。此外，建模分析時(shí)的模型尺寸和邊界條件設置均在一定程度上可能影響塑性變形求解過(guò)程。
