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