天然氣水合物開采區內海床不均勻沉降對管道力學性能的影響

謝 鵬，常江濤，肖卓臨

(1. 中山大學 海洋工程與技術學院，廣東 珠海 519082; 2. 南方海洋科學與工程廣東省實驗室(珠海)，廣東 珠海 519082)

天然氣水合物廣泛存在于海底的沉積物中，分布范圍廣、資源量巨大、能量密度高，是一種高效的清潔能源。現行的水合物開采方法主要有加熱法、降壓法、注入外加劑法等[1]，其基本原理是通過改變水合物所處環境的溫度和壓力，打破水合物穩定存在的平衡條件，使其分解為天然氣。國際上已經有少數國家開展了水合物試采，其中產量最高的是我國于2017年5月10日實施的在我國南海神狐海域水合物試采，累計試采穩定產氣60 d，總計產氣量為30.9×104m3，其中最高日產為3.5×104m3/d，完成了驗證我國水合物開采關鍵技術的目標[2]。根據國家規劃，預計2030年實現天然氣水合物的產業化開采，建成產能達到10×108m3的大型氣田。根據最高產量3.5×104m3/d的產量進行保守估計，建造一個10×108m3的大型氣田，日均產量約為273×104m3。為實現這一目標，預計需要在氣田內布設78口同等流量的產氣井進行同時開采和聯合集輸。因此，在大規模開發的天然氣水合物開采區，必須要使用大量海底管道和管匯將各井口串聯，并進行天然氣集輸。

天然氣水合物的開采會誘發海床發生不均勻沉降，沉降會改變海底管道的受力狀態，甚至導致管道發生大變形，造成管道發生屈曲、斷裂、懸跨等失效[3-4]。大量研究表明，當天然氣水合物在儲層內分解后，水合物沉積層的結構會發生很大程度的擾動，其強度和力學特性也會發生改變，當水合物的開采量達到某一程度時，則會誘發海床發生不均勻沉降、海底滑坡、土層失穩等地質災害[5]，歷史上著名的挪威Storegga第二次海底滑坡，美國東海岸Carolina岸外的Cape Fear滑坡，西地中海的Balearic巨型濁流層和西非大陸架的海底滑坡體等事故都是由水合物分解造成的[6-7]。萬義釗等[8]基于“熱-流-固-化”四場耦合的水合物開采分析模型，研究了水合物開采早期的海床沉降規律和影響參數。研究發現，在開采第30天時，開采井位置處的最大沉降量為0.032 m，海底面沉降范圍的半徑約為166 m。開采第60天時開采井位置的最大沉降量為0.035 m，海底面沉降范圍的半徑約為232 m。金光榮等[9]通過建立南海神狐海域水合物降壓開采模型，研究了更長時間開采范圍內的海床沉降規律，獲得了海床沉降量與開采時間的關系曲線，并針對開采井500 m范圍內的海床沉降情況進行了詳細分析。結果顯示，開采第2年時，在開采井附近位置海床沉降量達到2.6 m，且隨著開采時間和開采量的增加，海床沉降量增大。

然而，現有的研究主要針對天然氣水合物開采過程中水合物的分解過程進行研究，并預測水合物開采過程引起的海床沉降規律，并未考慮布設在開采區范圍內的海底管道等工程設施，也未研究水合物開采引起的海底管道受力狀態和失效模式。為解決這一問題，基于有限元分析方法，建立天然氣水合物開采區內“海床-管道”耦合作用模型，以“生死單元”模擬了天然氣水合物開采過程中海床沉降過程，獲得了海床變形形態及其對管道應力、應變、彎矩、懸跨等力學行為的影響。研究結果可以為建立天然氣水合物開采區內海底管道的設計方法提供依據，為未來天然氣水合物的大規模產業化開采做技術儲備。

1 “管-土”相互作用的有限元分析

水合物開采過程中，位于井口附近的水合物先分解，這將會擾動開采井口附近區域的沉積層結構，改變沉積層的強度和力學特性，誘發海床發生不均勻沉降，而沉降最嚴重的區域就是開采井口，海床的沉降量以開采井口為中心向外輻射遞減。研究位于水合物開采區內場(in-field)的海底管道，該管道需要與開采井口進行連接，在大規模的水合物開發過程中預計需要在氣田內布設78口同等流量的開采井進行同時開采和聯合集輸，必須要使用大量的海底管道和管匯將各井口串聯起來，因此必然導致海底管道經過井口附近。當管道從開采井口附近穿過時，由于海床的沉降，管道在自身重力下也會發生豎向位移，位移最大處位于開采井處，由于管道發生變形，管道的應力、應變、彎矩等力學特性也會發生相應的變化。因此，有限元模型的建立以海底管道跨越單井口為研究對象，模擬在進行單井口水合物開采時引起的海床沉降規律，及其對附近管道力學性能的影響。

1.1 海底管道模型

基于ABAQUS有限元軟件進行數值模擬，管道模型全長600 m，管道半徑0.31 m, 厚度0.02 m，如圖1所示。管道模型采用X65鋼，主要材料參數如表1所示。ABAQUS中適用于接觸問題的單元須為減縮積分單元[10]，因此管道采用減縮積分的4節點曲面殼單元(S4R)，該單元是一種通用的殼單元類型，具有很好的適應性，既可以用于厚殼問題的模擬，也可以用于薄殼問題的模擬。減縮積分是相對于完全積分來說的，減縮積分單元比完全積分單元在每個方向上少使用一個積分點，在模型受彎曲時可以消除計算中的剪切鎖閉問題且對計算結果的精度影響不大。

圖1 管道模型Fig. 1 Pipeline model

表1 X65 鋼材主要參數Tab. 1 Main parameters of X65

采用Romberg-Osgood方程建立管道材料的本構關系，其表達式為：

(1)

式中：ε是材料應變；σ是材料應力；E為彈性模量；σs是屈服強度；A和B是描述材料塑性應變的硬化參數，A=1.29，B=25.58 。

圖2 海床模型Fig. 2 Seabed model

1.2 海床模型

將海床的不均勻沉降作為海底管道力學分析的邊界條件，基于有限單元法的實體單元模擬水合物儲層被開挖后海床的沉降變形過程，構造出“等效的”不均勻海床表面變形形態，不需要對水合物開采過程中的熱量傳遞、固氣轉變等過程進行分析。分析時，海床選擇減縮積分的8節點六面體單元(C3D8R)，該單元對模擬非常大的網格扭曲問題有很好的適應性，由于模型中海床會發生較大的變形，因此選用該單元進行網格劃分。海床模型全長600 m，寬300 m，高200 m，其中天然氣水合物層深度100 m，厚度50 m，海床上表面設置寬度為20倍管道直徑的條帶作為和管道接觸的區域，如圖2所示。海床模型材料簡化為彈塑性材料，模型主要參數如表2所示。該模型的分析參數基于我國南海北部陸坡神狐海域天然氣水合物的成藏規律確定，研究結果可以代表未來在該海域內海底管道沉降及受力分析的一般情況。

表2 海床主要參數Tab. 2 Main parameters of seabed

1.3 接觸設置

管道外壁與海床上表面之間設置接觸面，采用單純的“主控-從屬”接觸算法。在計算過程中，為了獲得較好的模擬結果，選擇管道的外表面作為主接觸面，海床的上表面取20倍管徑寬度的條帶區域作為從接觸面。接觸作用屬性設置為法向“硬”接觸，切向采用罰函數摩擦公式定義切向摩擦力。

2 基于“生死單元”模擬海床不均勻沉降過程

“生死單元”功能是為了克服有限元本身缺陷提出的一種分析方法，被廣泛應用于礦產開采、隧道鉆進、基坑開挖等巖土工程領域，用來模擬被開挖部分“消失”的情況[11-12]。該單元彌補了有限單元法常規單元的不足。由于有限元本身是基于連續介質力學的，所研究的物體需要是連續的，在這樣的理論框架下，單元本身是不會消失的。然而在水合物開采過程中，由于水合物分解并經開采井運輸后，勢必在原區域產生短暫的環空，并造成海床沉降。因此，在進行數值分析時，需要模擬水合物儲層單元消失或者完全失效的情況。為了模擬這種工況，采用ABAQUS“生死單元”功能進行研究。在移除單元之前，ABAQUS自動計算和儲存待移除單元和余下單元交界面上的節點力，這些節點力在移除步中逐漸減少為零。在有限元分析過程中，通過分步移除水合物層一定范圍內的單元，模擬水合物逐步開采的過程，在此基礎上得到海床的不均勻沉降模型，如圖3所示。具體分析過程可分為13個分析步。

圖3 海床的不均勻沉降模擬Fig. 3 Simulation of uneven subsidence of seabed

第1分析步：管道左端施加x、y、z三個方向的約束；右端施加x方向的約束，約束管道的橫向位移。海床底面施加固定約束，海床左右兩端施加x、y、z三個方向的約束。

第2分析步：為管道施加重力載荷，重力加速度取一較小值0.1 m/s2,使得管道與海床之間的接觸作用平穩的建立起來。同時在管道的右端施加一水平張力T0來模擬管道內部的軸向拉力，其值為200 kN，此值主要根據工程經驗得到。

第3分析步：為海床施加體力載荷，海床的重度為20 kN/m3,將管道受到的重力加速度改為10 m/s2。

第4分析步：從第4分析步(即第1步開采)開始模擬水合物的開采過程，開采井口設置在海床模型上表面的中心點處，開采位置設置在開采井口下100 m深的水合物層處，開采厚度50 m，開采范圍20 m。開采過程通過ABAQUS中的單元生死功能實現，通過移除水合物層一定范圍內的單元，來模擬逐步開采的過程。

第5～13分析步：重復第4分析步，開采范圍以20 m為梯度逐步增大，共開采10步，最終開采范圍為200 m。通過這種逐步開采的方式獲得海床的沉降以及管道的力學響應。

3 天然氣水合物開采區內海床沉降規律及驗證

3.1 天然氣水合物開采過程中海床不均勻沉降量

通過ABAQUS進行數值模擬，對水合物進行分步開采，得到開采過程中的每一步海床沉降量及形態分布結果，如圖4、圖5所示。提取第5、第7、第9和第10步開采的位移云圖，開采影響范圍分別為100 m、140 m、180 m和200 m，分別如圖5所示。從圖中可以明顯的發現，隨著開采范圍的增大，海床的位移逐漸增大，第10步開采時，海床沉降量達到最大值3.543 m，位移云圖顯示的較大位移區域逐漸向開采井處收縮。

圖4 海床位移云圖Fig. 4 Displacement nephogram of seabed

圖5 不同開采步下模型的位移云圖Fig. 5 Displacement nephogram of models under different mining steps

3.2 海床沉降量規律的驗證

以“生死單元”法模擬實體單元消失引起海床沉降變形，等效模擬真實開采過程中的海床沉降量和海床表面形態分布規律，結果如圖6所示。x=0處為開采井，圖6(a)為海床沉降量在模型長度方向上隨開采井距離的變化曲線，圖6(b)為海床沉降量在模型寬度方向上隨開采井距離的變化曲線；曲線之間的變化趨勢相同，海床沉降量隨著與開采井距離的增大而減小。開采過程中開采井處的海床沉降量如表3所示。

圖6 海床沉降量在模型長度和寬度方向上的變化Fig. 6 Variation of seabed subsidence with the length and width of the model

表3 開采過程中開采井處的海床沉降量Tab. 3 Subsidence of seabed at the mining well in mining process

為驗證該模型模擬結果的準確性和適用性，與基于“熱-水動力-機械耦合”模型模擬的水合物開采引起的海床沉降情況進行對比和驗證。金光榮等[9]通過使用TOUGH2Biot模擬器，擴展了Biot合并模型，從而得到了用于水合物開采的“熱-水動力-機械耦合”模型，通過改變開采井底部壓力，獲得了神狐海域天然氣水合物開采過程中海床不均勻沉降的規律。

將“生死單元”法獲得的海床沉降量與文獻進行對比，結果如圖7所示，圖7(a)為“生死單元”法海床沉降計算結果，圖7(b)為金光榮等[9]通過降壓法開采水合物得到的海床沉降曲線。圖7(a)分別為第5、第7和第9步開采時海床沉降的結果，海床最大沉降量分別為1.621 m、2.068 m和2.767 m，沉降發生的主要區域位于開采井400 m范圍內；圖7(b)分別為開采時間為第0.5、第1和第2年時海床沉降的結果，海床最大沉降量分別為1.52 m、1.80 m和2.50 m，沉降發生的主要區域位于開采井500 m范圍內。兩種結果較為接近，海床的最大沉降都位于開采井處，當海床與開采井的距離增大時，沉降量都會逐漸減小。將圖7(a)中的曲線擬合到圖7(b)中，可以看出文中“生死單元”法模擬的水合物開采引起的海床表面沉降形態與“熱-水動力-機械耦合”法模擬的海床表面變形形態能夠很好的吻合，海床沉降的數值也很接近前人的研究成果。驗證了通過使用ABAQUS中單元生死功能模擬水合物開采進而引發海床沉降這一方法是可行的，以該模型獲得的海床沉降規律作為管道受力分析的邊界條件是合理的。

圖7 海床位移曲線Fig. 7 Displacement curve of seabed

4 海床不均勻沉降引起的管道力學性能變化

基于獲得的海床沉降規律，研究海床沉降對管道造成的力學行為變化，分析管道的豎向位移、曲率和彎矩、應力應變分布、管道懸跨等問題。

4.1 管道豎向位移

在水合物開采過程中，管道豎向位移隨著海床的沉降發生變化。圖8(a)為管道位置變化示意，圖8(b)為每一開采步的管道豎向位移量，表4為開采過程中開采井處管道的沉降量。結果顯示，管道在開采井附近的豎向位移最大，為3.407 m，當管道位置遠離開采井時，由于海床的沉降量逐漸減小，管道的豎向位移也逐漸減小。管道的位移曲線關于開采井位置對稱且沉降主要發生在開采井周圍400 m范圍內，沉降的最大值總是出現在開采井處。這一形態與海床的表面沉降形態是一致的。

圖8 管道的豎向位移Fig. 8 Vertical displacement of pipeline

表4 開采過程中開采井處的管道沉降量Tab. 4 Subsidence of pipeline at the mining well in mining process

4.2 管道的曲率及彎矩分布

圖9 管道曲率Fig. 9 Curvature of pipeline

圖10 管道彎矩Fig. 10 Bending moment of pipeline

4.3 管道的應力應變分布

由于天然氣水合物的開采導致管道發生較大的沉降，盡管沉降的跨度較大，但是依然會對管道的應力應變造成影響。由于海床沉降主要發生在開采井附近區域，距離開采井附近150 m區域內的沉降量較大，對管道應力應變的影響更為顯著，因此只需選擇該區域進行分析即可。如圖11(a)所示，管道應力和應變的最大值均出現在管道中點的底部。每一步開采后，管道的最大應力和應變值如圖11(b)所示，隨著水合物開采范圍的擴大，管道的最大應力和應變都隨之增加，應力值的范圍為9.97～24.37 MPa，應變值的范圍為4.38×10-5～1.17×10-4。

圖11 管道的最大應力、應變Fig. 11 Maximum stress and strain of pipeline

為了探究該150 m范圍內管道上應力應變隨水平位置的變化趨勢，須對每一步開采所對應的管道的應力應變進行數據提取，得到的結果如圖12所示。

圖12 開采過程中管道的應力、應變分布Fig. 12 Stress and strain distribution of pipeline in mining process

如圖12所示，每一步開采后，管道的應力應變隨水平位置的變化趨勢是一致的，當管道位置靠近開采井時，管道的位移和變形逐漸增大，導致管道的應力應變逐漸增大；當管道位置遠離開采井時，管道的位移和變形逐漸減小，導致管道的應力應變逐漸減小。管道的應力應變曲線關于開采井位置對稱且最大值總是出現在開采井處。從圖11和圖12可以看出，天然氣水合物開采導致海床的不均勻沉降會改變管道的應力應變狀態，開采范圍的擴大會導致管道內部的應力和應變逐漸增大，但是其數值都還在安全范圍以內，因此水合物的開采對管道的應力應變影響相對較小，在應力應變的角度上，管道是安全可靠的。

4.4 管道的懸跨

在水合物開采過程中，管道雖然隨著海床的沉降發生豎向位移，但由于管道本身具有彎曲剛度且受到軸向拉力的約束，在豎向位移發生到一定程度后，則不再發生變化。因此，當海床的沉降量達到一定值時，管道不再繼續沉降，而是與海床發生分離，在海床表面形成懸跨，且逐步增大，如圖13所示。在開采過程中海床和管道的相對位移變化量如圖14所示。如圖14(a)所示，天然氣水合物開采在前4步即開采范圍在80 m以內時，海床與管道的最大沉降量相等，這表示海床與管道是始終接觸的，當天然氣水合物開采到第5步即開采范圍達到100 m時，海床與管道的最大沉降量分別為1.621 m和1.576 m，出現0.045 m的懸空，此時海床和管道已發生分離。隨著開采范圍的繼續增大，海床與管道之間的懸空距離整體呈上升趨勢，如圖14(b)所示。當海水流經懸跨管道時，在一定流速條件下，懸跨管道兩側會出現旋渦，并以一定的頻率交替脫落，從而在管道表面形成周期載荷，將會使懸跨管道在順流方向及橫流方向上發生振動，即渦激振動，很可能對海底管道造成疲勞損傷乃至破壞[13-15]。

圖13 管道與海床分離產生懸跨(開采步10)Fig. 13 Suspension span caused by separation of pipeline and seabed (Mining step 10)

圖14 開采過程中管道和海床的沉降量及懸跨高度變化Fig. 14 Change of seabed subsidence and height of suspension span of pipeline during mining process

5 結 語

經過研究和計算分析，可以得到以下結論：

1) 有限單元中“生死單元”法可以等效模擬天然氣水合物開采過程中海床的變形規律，為海底管道的沉降分析提供邊界條件，對于模擬海底管道在海床沉降邊界下的安全校核有重要的指導意義。

2) 天然氣水合物開采引起的海床沉降區域主要集中在開采井附近400 m的范圍內，其對管道的應力、應變、位移、曲率和彎矩均產生影響。由海床沉降引起的管道應力應變基本處于彈性范圍內，不會造成管道本身的強度破壞；但是當海床沉降達到一定程度時，鋪設在開采區內的海底管道將與海床發生分離，形成懸跨，并引發渦激振動風險，對管道造成疲勞破壞。

3) 基于天然氣水合物“單井開采”模型進行海床沉降研究，與實際工程中的“多井聯合開采”模式仍有一定差距。后續研究中可進一步考慮多井聯合開采工況下海床的變形規律及其影響，為天然氣水合物的大規模開發提供技術支撐。