鄭曉明，張春會,,3，李雪婷，田英輝，李永靖

(1.遼寧工程技術大學土木工程學院，遼寧阜新 123000；2.河北科技大學建筑工程學院，河北石家莊 050018；3.河北省巖土與結構體系防災減災技術創(chuàng)新中心，河北石家莊 050018；4.墨爾本大學基礎設施工程系，維多利亞墨爾本 3010)

深海油氣輸送管道一般直接鋪設于海床上，在風浪流及自重等荷載作用下海底油氣管道通常貫入海床10%～50%管道直徑[1-3]。海底油氣管道通常在高溫高壓環(huán)境下運行，產(chǎn)生的熱脹應力超過海床土體所能提供的約束管道側向屈曲和軸向整體位移的摩阻力時，管道就會發(fā)生屈曲或走管等失穩(wěn)，造成管道破裂、油氣泄漏，從而污染環(huán)境。海床土體對管道的約束力與管道貫入深度密切相關，因此準確預測海底管道貫入深度，對于分析其在位穩(wěn)定性至關重要[4-7]。

國內(nèi)外對海底管道貫入深度問題開展了很多研究，RANDOLPH等[8-10]基于上限解法和下限解法研究了黏土中管道貫入深度與土體抗力的關系。AUBENY等[11]采用有限元方法研究了海底管道在黏土海床中的貫入問題，提出了基于管-土界面摩擦系數(shù)的貫入深度與貫入阻力之間的經(jīng)驗公式。王騰等[12]、DUTTA等[13]、HAWLADER等[14]、CHATTERJEE等[15]、ULLAH等[16]分別采用耦合歐拉-拉格朗日(coupled Eulerian-Lagranginan，CEL)方法和RITSS(remeshing and interpolation technique with small strain)方法研究了土體的應變率效應和應變軟化對于管道在黏土海床貫入阻力的影響。MERIFIELD等[17]采用任意拉格朗日-歐拉(arbitrary Lagranginan-Eulerian，ALE)方法研究了管道在黏土海床貫入過程中土體重度和隆起對于貫入阻力的影響。VERLEY等[18]通過開展模型試驗研究了管道在黏土中貫入阻力和位移響應，提出了預測管道黏土中貫入深度的經(jīng)驗公式。DINGLE等[19]、CHEUK等[20]開展了管道在黏土中貫入的離心機試驗，研究了管土相互作用，并通過粒子圖像測速技術觀測了管道貫入過程中土體的流動。這些研究主要聚焦于黏土海床海底管道貫入問題。在北海、墨西哥灣、印度尼西亞、中國南海和渤海等地深水區(qū)廣泛分布著砂土海床，隨著這些地區(qū)油氣資源的開發(fā)，海底管道在砂土海床貫入問題的重要性日益凸顯，然而目前對砂土海床海底管道貫入研究還較少。VERLEY等[21]通過試驗對海底管道在砂土中貫入的位移-阻力進行了研究，提出了砂土海床管道貫入深度的經(jīng)驗預測公式。ZHAO等[22]為研究砂土中的管-土相互作用，開展了相關離心機試驗和數(shù)值研究。中國現(xiàn)行的海底管道相關領域技術標準為《海底管道系統(tǒng)規(guī)范》[23]，但該規(guī)范中尚無相關砂土海床海底管道貫入深度的預測方法。目前國內(nèi)外普遍采用的相關規(guī)范為挪威的船級社發(fā)布的DNV-RP-F114[24]，該規(guī)范建議采用地基承載力公式計算海底管道在砂土海床貫入過程中所受到的阻力。這些研究主要采用經(jīng)驗公式預測海底管道貫入，僅考慮內(nèi)摩擦角對海底管道貫入深度的影響。內(nèi)摩擦角確實是描述砂土海床強度的關鍵參數(shù)，然而海底管道貫入砂土是一個復雜的管土相互作用問題，除了砂土強度外，管-土界面特性和剪脹特性等也都可能影響海底管道的貫入深度。

本文基于CEL方法，建立海底管道貫入砂土海床的大變形有限元模型，通過數(shù)值仿真模擬海底管道在風浪流及自重等荷載作用下貫入砂土海床的過程，分析管-土界面特性、砂土強度和剪脹特性等對貫入深度的影響，進而提出砂土海床海底管道貫入深度的預測公式，從而為砂土海床海底管道穩(wěn)定性分析提供基礎。

1 砂土海床海底管道貫入的CEL模型

在風浪流及自重荷載等作用下，海底管道貫入砂土海床，擠壓管道周圍土體致其向兩側流動，即海底管道貫入砂土海床是一個典型的大變形問題。CEL方法結合了拉格朗日法和歐拉法的優(yōu)勢，在小變形區(qū)域采用拉格朗日單元，在大變形區(qū)域采用歐拉單元，能夠很好地模擬巖土工程中的大變形問題，故本文采用CEL方法模擬海底管道貫入砂土海床過程。

1.1 CEL模型

在ABAQUS軟件下建立海底管道貫入砂土海床的CEL模型。海底管道通常為鋼制，模量大，變形微小，本文將管道視為剛體，使用拉格朗日單元進行模擬，管道直徑D為0.5 m；海床砂土在管道貫入過程中發(fā)生流動，采用歐拉單元模擬，選用Mohr-Coulomb本構模型，根據(jù)海床砂土常見密度、內(nèi)摩擦角、剪脹角和彈性模型范圍進行土體參數(shù)的選取[25-26]，其中土體密度ρ為1 600 kg/m3，內(nèi)摩擦角φ為39°，剪脹角ψ為15°，彈性模量E為30 MPa。為避免計算收斂問題，砂土取一小的黏聚力0.1 Pa。

本文計算模型為三維模型。砂土寬度(x軸)為10D，高度(z軸)為8D，長度方向(y軸)設定為0.1 m。在土體上方建立1D高度無材料填充的歐拉單元，允許土體自由隆起，模擬海底管道貫入過程中海床的變形過程。為確保計算結果精度，在管-土接觸區(qū)域寬度3D、高度2D范圍內(nèi)進行網(wǎng)格加密，經(jīng)過反復試算，加密區(qū)域網(wǎng)格尺寸為1/40D，所建CEL模型如圖1所示。

圖1 CEL模型Fig.1 CEL model

管-土之間的接觸面法向為硬接觸，切向為庫侖摩擦，摩擦系數(shù)為1.0，即管-土接觸完全粗糙。土體模型的左右邊界為水平方向約束，底部邊界為豎直方向約束。

在有限元計算中，先對土體進行地應力平衡計算，然后管道豎向貫入，貫入速度v=0.1 m/s，管道貫入0.5D時停止貫入，通過位移加載的方式模擬海底管道在風浪流及自重等荷載作用下貫入砂土海床的過程。

1.2 模型驗證

為驗證本文所建有限元數(shù)值模型的正確性，將CEL模擬結果與DNV-RP-F114規(guī)范計算結果進行對比。DNV-RP-F114規(guī)范中管道在砂土海床貫入的豎向抗力V計算公式為[24]

(1)

式中：γ′為土體有效重度；Nγ，Nq為土體承載力系數(shù)；B為管-土接觸寬度；dq為影響深度效應的因素；z0為修正的管道嵌入深度。

承載力系數(shù)Nq計算公式為

Nq=eπtan φtan2(45+φ/2)，

(2)

式中φ為土體內(nèi)摩擦角。

承載力系數(shù)Nγ分別采用HANSEN[27]和VESIC[28]地基極限承載力系數(shù)，其中HANSEN地基極限承載力系數(shù)Nγ為[27]

Nγ=1.5(Nq-1)tanφ。

(3)

VESIC地基極限承載力系數(shù)Nγ為[28]

Nγ=2(Nq+1)tanφ。

(4)

由式(2)—式(4)可見，DNV-RP-F114規(guī)范中豎向抗力僅與砂土的內(nèi)摩擦角有關。B，z0，dq計算方法參見DNV-RP-F114規(guī)范。

本文CEL模擬結果與DNV-RP-F114規(guī)范計算結果對比如圖2所示。為便于結果的表示，參照VERLEY等[21]的方法將貫入阻力歸一化處理。圖2中豎坐標為w/D，其中w為管道貫入深度，橫坐標為歸一化貫入阻力V/(γD2)，其中γ為土體重度。由圖2可知，CEL模擬結果的貫入阻力隨貫入深度增加而增大，在貫入的開始階段，貫入阻力增加較為緩慢，隨著貫入深度增加，貫入阻力增長變快，這與DNV規(guī)范計算結果規(guī)律一致。從圖2還可以看出，w在0～0.1D范圍內(nèi)CEL模擬結果與Nγ取式(4)的DNV規(guī)范計算結果基本一致，w在0.1D～0.5D范圍內(nèi)CEL模擬結果與DNV規(guī)范Nγ取式(3)時計算結果基本吻合，本文CEL模擬結果整體位于DNV規(guī)范計算結果區(qū)間內(nèi)，這說明本文建立的砂土海床海底貫入的CEL模型是合理的。

圖2 CEL模擬結果驗證 Fig.2 Verification of CEL simulation results

2 砂土海床海底管道貫入深度影響因素分析

影響砂土海床海底管道貫入深度的可能因素包括管-土界面特性、土體強度、剪脹特性和彈性模量，本文以上文所述有限元模型作為基本算例，分析管-土界面摩擦系數(shù)、砂土內(nèi)摩擦角和剪脹角及砂土彈性模量對砂土海床海底管道貫入深度的影響規(guī)律。

2.1 管-土界面摩擦系數(shù)的影響

管-土界面摩擦系數(shù)α取0時，管土之間為光滑接觸；管-土界面摩擦系數(shù)α取1.0時，管土之間為完全粗糙。本文α分別取為0，0.25，0.5，0.75和1.0，利用CEL模型計算海底管道貫入的貫入阻力。圖3為不同α下的貫入深度與貫入阻力關系圖。從圖3可知，海底管道貫入深度為0.5D，α=0時歸一化貫入阻力為15.96，α分別為0.25，0.5，0.75和1.0時對應的歸一化貫入阻力為24.57，30.34，34.15和36.74，與α=0時歸一化貫入阻力相比分別增加了53.95%，90.1%，113.97%和130.2%。由此可以看出，管-土界面摩擦系數(shù)對砂土海床海底管道貫入阻力有很大影響，海底管道在貫入相同深度時，管-土界面摩擦系數(shù)越大，相應的貫入阻力越大。

圖3 不同管-土界面摩擦系數(shù)下的模擬結果 Fig.3 Simulation results under different friction coefficients of pipe-soil interface

圖4為管道貫入深度0.5D、管-土界面摩擦系數(shù)分別為0和1.0時對應的管道周圍砂土速度矢量圖。由圖4可以看出，土體的破壞區(qū)域總體表現(xiàn)為自管道底部向兩側回轉，形成類楔形破壞區(qū)，在管道兩側形成隆起，海床土體形成了連續(xù)的滑動面并自管道底部延伸到海床表面，土體破壞形式表現(xiàn)為整體剪切破壞。對比圖4 a)和圖4 b)可以看出，管-土界面摩擦系數(shù)不同時，砂土的流動狀態(tài)及影響區(qū)域存在差異。如圖4 a)所示，當管-土界面完全光滑時，管土之間無摩擦，管道向下貫入時管道正下方的土體表現(xiàn)為沿管道向管道兩側流動。由圖4 b)可以看出，管-土界面完全粗糙時，管土之間存在摩擦，管道向下貫入時，管道正下方土體在管土間摩擦力的作用下會產(chǎn)生向下的運動，管-土接觸面特性通過管-土間摩擦影響海底管道貫入時土體的流動。同時可以發(fā)現(xiàn)，當管-土接觸完全光滑(α=0)時，海床隆起的高度約為0.53D，影響深度約為0.82D，影響寬度約為2.15D；而管-土接觸完全粗糙(α=1.0)時，海床隆起高度約為0.37D，影響深度約為1.38D，影響寬度約為3.15D，相比管道光滑時土體影響區(qū)域更大，這導致管-土界面摩擦系數(shù)增大后，貫入阻力增大。

圖4 不同摩擦系數(shù)下的土體速度矢量圖(w/D=0.5)Fig.4 Soil velocity vector diagram under different friction coefficients (w/D=0.5)

2.2 內(nèi)摩擦角的影響

海床砂土常見的內(nèi)摩擦角φ范圍為30°～40°，本文取φ分別為31°，33°，35°，37°和39°，開展海底管道貫入數(shù)值模擬，獲得管道貫入深度與貫入阻力之間的關系如圖5所示。由圖5可知，當貫入深度為0.5D，φ為31°時，歸一化貫入阻力為15.74，當φ為33°，35°，37°和39°時，對應的歸一化貫入阻力分別為19.79，23.62，29.40和36.74，相比于內(nèi)摩擦角φ=31°的工況分別增加了25.73%，50.06%，86.79%和133.42%，這說明內(nèi)摩擦角對管道貫入阻力有很大影響，當貫入深度相同時，砂土內(nèi)摩擦角越大，管道貫入阻力越大。

圖5 不同內(nèi)摩擦角下的模擬結果Fig.5 Simulation results under different friction angles

圖6為管道貫入深度w=0.5D、土體φ分別為31°和39°所對應的管道周圍砂土速度矢量圖。由圖6可以看出，內(nèi)摩擦角不同時，海床土體均呈現(xiàn)明顯的整體剪切破壞特征。對比圖6 a)和圖6 b)可以發(fā)現(xiàn)，在φ為31°和39°時，砂土影響深度分別約為1.13D和1.38D，內(nèi)摩擦角越大，影響深度越大，這主要是由于內(nèi)摩擦角越大，土顆粒之間表面摩擦力越大，在管道貫入過程中，砂土顆粒間存在摩擦的同時也存在相對滑動，當顆粒間摩擦越大時，相對滑動越難發(fā)生，砂土顆粒受到管道擠壓時，更易向下運動。同時，內(nèi)摩擦角越大，土顆粒之間嵌入和聯(lián)鎖作用產(chǎn)生的咬合力也更強，這導致內(nèi)摩擦角增大后，貫入阻力增大。

圖6 不同內(nèi)摩擦角下的土體速度矢量圖(w/D=0.5)Fig.6 Soil velocity vector diagram under different friction angles (w/D=0.5)

2.3 剪脹角的影響

砂土的剪脹角與松密程度有關，一般密砂剪脹角可取15°，松砂剪脹角可小于10°，本文取剪脹角ψ分別為9°，11°，13°，15°和17°，開展砂土海床海底管道貫入的CEL模擬，獲得管道深度與貫入阻力之間的關系如圖7所示。由圖7可知，貫入深度為0.5D，ψ=9°時歸一化貫入阻力為29.50，ψ分別為11°，13°，15°和17°時對應的歸一化貫入阻力為31.91，33.98，36.33和39.00，與ψ=9°時歸一化貫入阻力相比分別增加了8.17%，15.19%，23.15%和32.2%，這說明剪脹角同樣也會對海底管道在砂土海床的貫入阻力產(chǎn)生影響，當貫入深度相同、剪脹角不同時，貫入阻力隨著剪脹角的增大而增大。

圖7 不同剪脹角下的模擬結果Fig.7 Simulation results under different dilation angles

圖8為管道貫入深度為0.5D、土體剪脹角ψ=9°和15°時分別對應的管道周圍砂土速度矢量圖。由圖8可以看出，剪脹角不同時，海床土體破壞形式均表現(xiàn)為整體剪切破壞。對比圖8 a)和圖8 b)可以看出，當土體的剪脹角越大時，土體隆起高度及影響的深度和寬度范圍越大，這主要是因為隨著剪脹角增大，土體在受到管道的貫入之后產(chǎn)生的塑性體應變和塑性區(qū)范圍增大，剪脹角越大土體影響區(qū)域越大，這導致剪脹角增大后，貫入阻力增大。

圖8 不同剪脹角下的土體速度矢量圖(w/D=0.5)Fig.8 Soil velocity vector diagram under different dilation angles (w/D=0.5)

2.4 彈性模量的影響

選取5種不同的彈性模量，分別為E=10，15，20，25，30 MPa，通過有限元計算得到如圖9所示的貫入深度與貫入阻力關系圖。由圖9中可以看出，5條曲線基本重合，即彈性模量不同時貫入阻力近似相同，從而可以說明彈性模量E的變化對于貫入阻力來說沒有顯著影響。

圖9 不同彈性模量下的模擬結果Fig.9 Simulation results under different elastic modulus

3 貫入深度預測的經(jīng)驗公式

為預測砂土海床海底管道貫入深度，提出海底管道貫入砂土海床貫入深度的經(jīng)驗公式。由分析可知，海底管道在砂土海床的貫入過程受管-土界面摩擦系數(shù)、砂土內(nèi)摩擦角和剪脹角的影響較大，故本文提出考慮上述因素的海底管道貫入深度預測經(jīng)驗公式。

3.1 經(jīng)驗公式的提出

從以上數(shù)值模擬結果可以看出，當貫入阻力一定時，貫入深度受管-土界面摩擦系數(shù)、砂土內(nèi)摩擦角、剪脹角的共同影響，貫入阻力分別隨管-土界面摩擦系數(shù)、砂土內(nèi)摩擦角、剪脹角和貫入深度的增加而增大，貫入阻力與管-土界面摩擦系數(shù)、砂土內(nèi)摩擦角、剪脹角和貫入深度呈正相關關系。此外由數(shù)值模擬結果發(fā)現(xiàn)，內(nèi)摩擦角為33°，35°，37°和39°時對應的歸一化貫入阻力相比于內(nèi)摩擦角31°的歸一化貫入阻力分別增加了25.73%，50.06%，86.79%和133.42%，其增幅在3個影響因素中最大，因此可說明在一定范圍內(nèi)，內(nèi)摩擦角對于貫入阻力影響最大。

選取的4個無量綱量分別為

(5)

歸一化貫入阻力可以表示為

(6)

各無量綱量與歸一化貫入阻力關系如圖10所示，可以發(fā)現(xiàn)歸一化貫入阻力與4個無量綱量近似呈冪函數(shù)變化。冪指數(shù)函數(shù)可表示為

(7)

式中：a，b為擬合系數(shù)。

對各無量綱量及對應的歸一化貫入阻力采用式(7)進行擬合，擬合結果如圖10所示。由圖10可知，采用冪指數(shù)函數(shù)對歸一化貫入阻力與4個無量綱量進行擬合，相關系數(shù)R2平均值為0.991，擬合結果較好，說明歸一化貫入阻力隨各參數(shù)變量的變化近似呈冪函數(shù)規(guī)律。

圖10 各無量綱量與歸一化阻力關系圖Fig.10 Relationship between each dimensionless quantity and normalized resistance

在圖10基礎上，建立歸一化貫入阻力與4個無量綱量之間的函數(shù)關系即非線性乘法模型：

(8)

式中：β，x1，x2，x3，x4為擬合系數(shù)。

取自然對數(shù)，可得：

(9)

利用LINEST函數(shù)，進行多元線性回歸計算，然后得到：

β=2 966.337，x1=0.332，x2=3.676，x3=0.513，x4=0.671。

(10)

本文分別研究了管-土界面摩擦系數(shù)α范圍為0～1，砂土內(nèi)摩擦角φ范圍為31°～39°，剪脹角ψ范圍為9°～17°下的海底管道在砂土海床深度范圍為0～0.5D的貫入阻力隨貫入深度的變化關系，基于數(shù)值模擬結果，最終得出貫入深度的預測經(jīng)驗公式：

(11)

公式的適用范圍為

(12)

3.2 經(jīng)驗公式驗證

為驗證經(jīng)驗公式的準確性，將ZHAO等[22]關于管道在砂土貫入的離心機試驗土體參數(shù)代入至本文經(jīng)驗公式中，與其離心機試驗結果進行對比驗證，對比結果如圖11所示。由圖11可以看出，離心機試驗數(shù)據(jù)具有較大的波動性，但整體呈現(xiàn)出貫入阻力隨貫入深度的增加而增大的趨勢，與經(jīng)驗公式結果規(guī)律一致。經(jīng)驗公式計算結果位于離心機試驗結果波動范圍之間，與離心機試驗結果具有良好的相關性，這說明本文所提出的預測海底管道在砂土海床貫入深度的經(jīng)驗公式是合理的。

圖11 經(jīng)驗公式驗證Fig.11 Verification of empirical formula

3.3 經(jīng)驗公式與DNV規(guī)范對比

為驗證擬合公式的準確性，在既有模型參數(shù)的基礎上重新選取了6種不同工況進行計算，計算工況如表1所示。將CEL模擬結果分別與DNV-RP-F114規(guī)范計算結果和所提出的經(jīng)驗公式計算結果進行對比，結果如圖12所示。

表1 計算工況

圖12 經(jīng)驗公式與DNV規(guī)范對比結果Fig.12 Comparison results between empirical formula and DNV specification

圖12 a)為CEL模擬結果與DNV-RP-F114規(guī)范計算結果對比，由圖12 a)可以看出，工況Ⅰ-4～Ⅰ-6的CEL模擬結果位于DNV-RP-F114規(guī)范計算結果區(qū)間，而工況Ⅰ-1～Ⅰ-3的CEL模擬結果與DNV-RP-F114規(guī)范計算結果相關性較差，說明DNV-RP-F114規(guī)范只考慮砂土內(nèi)摩擦角的影響，雖在一定程度上可預測海底管道在砂土海床的貫入深度，但具有一定的局限性。圖12 b)為CEL模擬結果與本文所得經(jīng)驗公式計算結果對比，如圖12 b)所示，CEL模擬結果與經(jīng)驗公式計算結果相差較小，擬合效果較好，由此說明，本文基于管-土界面摩擦系數(shù)、砂土內(nèi)摩擦角、剪脹角所提出的砂土海床海底管道貫入深度預測經(jīng)驗公式適用性更強，所提出的經(jīng)驗公式在一定條件范圍內(nèi)能夠較好地預測海底管道在砂土海床的貫入深度。

4 結 論

本文基于ABAQUS中的CEL方法建立了海底管道在砂土海床豎向貫入的有限元模型，模擬了海底管道在豎向貫入下砂土海床的大變形過程，研究了管-土界面摩擦系數(shù)、砂土內(nèi)摩擦角、剪脹角及彈性模量對貫入深度的影響，主要得到結論如下。

1)在海底管道貫入過程中，管道兩側的海床發(fā)生隆起，在海床土體中形成了連續(xù)的滑動面且自管道底部延伸至海床表面，土體破壞表現(xiàn)為整體剪切破壞形式。

2)管-土界面摩擦系數(shù)、砂土內(nèi)摩擦角及剪脹角均會影響砂土海床海底管道貫入深度，砂土彈性模量則對海底管道貫入深度沒有顯著影響。海底管道貫入阻力隨管-土界面摩擦系數(shù)、砂土內(nèi)摩擦角和剪脹角增加而呈冪函數(shù)增大。

3)所提出的考慮管-土界面摩擦系數(shù)、砂土內(nèi)摩擦角及剪脹角影響的貫入深度預測經(jīng)驗公式能夠較好地預測海底管道在砂土海床的貫入深度。

本文主要針對砂土海床海底管道貫入進行了數(shù)值研究，后續(xù)研究中還應補充管道在砂土中豎向貫入相關試驗，對所提出的貫入深度預測經(jīng)驗公式進行修正和完善，為下一步開展的海底管道在位穩(wěn)定性研究提供基礎。