劉小麗 李喬花 廖周全

(①山東省海洋環(huán)境地質(zhì)工程重點實驗室， 青島 266100， 中國) (②中國海洋大學， 環(huán)境科學與工程學院， 青島 266100， 中國)

0 引 言

海底管線常被稱為海洋油氣資源開發(fā)中的“生命線”，其安全問題歷來受到眾多研究者和海洋工程設計人員的關注。極端氣候下風暴潮掀起的巨浪會導致管線周圍海床的液化進而造成海底管線的失穩(wěn)破壞(張衍濤， 2009; 付長靜等， 2014)，因此較多文獻對波浪作用下海底埋置管線周圍海床的響應及液化特征進行了計算研究。當海床土滲透系數(shù)相對較小而波浪荷載相對較大時，波浪循環(huán)荷載下海床土的塑性變形會引起超孔隙水壓力的不斷累積，當某一深度的累積孔壓超過其上覆有效應力時，即引起該深度位置海床土的液化，進而威脅管線的安全。

目前波浪作用下海床累積響應的計算方法主要有兩種，第1種方法是直接利用彈塑性本構模型來模擬循環(huán)荷載下孔壓的累積過程(Dunn et al.，2006;Ye et al.，2015; 王小雯等， 2018)，該種方法相對較為復雜，本構模型的描述需要大量不易測試的參數(shù)，因而其應用受到一定限制； 另一種方法是將循環(huán)剪切試驗得到的半經(jīng)驗孔壓累積源項表達式與滲流連續(xù)方程相結(jié)合以計算累積孔壓的發(fā)展過程(Seed et al., 1978)，該種方法因計算簡便參數(shù)較少而得到了較多應用(McDougal et al.，1989; 李小均等， 2009;Sumer et al.，2012)?；诳讐豪鄯e源項的方法雖然簡單，但由于其未考慮累積孔壓與海床應力的耦合效應，可能會導致計算結(jié)果出現(xiàn)較大偏差，因此，Liu et al. (2019)基于該方法的計算思想，提出了一種能夠考慮累積孔壓與海床應力耦合發(fā)展的計算模型，且未增加額外的計算參數(shù)。

基于孔壓累積源項計算海床累積響應的方法已經(jīng)被應用于波浪作用下管線-海床系統(tǒng)的響應分析中(Zhao et al.，2014; Foo et al.，2019)，但均未考慮累積孔壓與海床應力的耦合效應。波浪作用下孔壓累積與海床應力的耦合效應對埋置管線周圍海床累積響應及其液化特征的影響如何，目前還未見相關的分析研究。為此，基于Liu et al. (2019)提出的計算模型，考慮累積孔壓與海床應力的耦合效應，對非線性行進波作用下埋置管線周圍海床土的累積響應特征進行模擬計算，并通過與不考慮孔壓累積與海床應力耦合計算結(jié)果的對比分析，揭示累積孔壓與應力的耦合效應對海底埋置管線周圍海床累積響應特征的影響及其作用機制。

為了闡述方便，文中將未考慮孔壓累積與海床應力耦合的累積響應計算模型稱為“非耦合模型”，而將能夠考慮孔壓累積與海床應力耦合發(fā)展的海床響應計算模型稱為“耦合模型”。

1 數(shù)值模型

波浪作用下含埋置管線的海床響應計算示意圖如圖 1 所示，其中：h為海床厚度，d為水深，H為波高，管線埋置在海床水平方向的中間位置，D為管徑，e為埋深(指管線中心到海床表面的距離)。波浪在海床表面沿x軸正方向傳播，海床土視為均勻各向同性的多孔彈性介質(zhì)。

圖 1 波浪作用下含埋置管線海床計算示意圖Fig. 1 Sketch of the wave-induced seabed response around a buried pipeline

1.1 海床響應控制方程

海床響應耦合計算模型的控制方程推導詳見Liu et al. (2019)的文獻，此處給出其控制方程為：

(1)

(2)

(3)

(4)

式中：Kf為孔隙水的彈性體積模量，Kf=2×109N·m-2；Sr為海床土體的飽和度；P0為絕對靜水壓力。

mv為海床土的體積壓縮系數(shù)，可表示為：

(5)

ψ為累積孔壓增長函數(shù)(孔壓累積源項)，可表示為：

(6)

式中： |τmax|為海床周期最大剪應力；T為波浪周期；α=0.34Dr+0.084，β=0.37Dr-0.46，均為與海床相對密度有關的孔壓累積計算參數(shù)；Dr為土體相對密度；σ′0為初始平均有效應力。

對于含埋置管線的海床，管線自重會影響周圍海床的初始平均有效應力，因此，σ′0表達式為(Zhao et al.，2014)：

(7)

式中：σ′x0和σ′z0分別為重力作用下海床土的水平向和豎向有效正應力。

根據(jù)相關文獻(Jeng et al.，2015)，非耦合模型海床的控制方程為：

(8)

式中：pres為波致累積孔隙水壓力，其余符號同上述。

1.2 管線計算控制方程

管線為彈性體，其平衡方程為：

(9)

(10)

式中：up、wp分別為管線的水平位移和豎向位移；Gp和μp分別為管線的剪切模量和泊松比。

1.3 邊界條件

考慮波浪的非線性特征，采用二階Stokes行進波。在海床上表面，波致孔隙水壓力為海床表面的波壓力Pb(x，t)，其表達式為(Jeng et al.， 1997)：

Pb(x，t)=P0(C11cos(λx-ωt)+C02+

C20cos 2(λx-ωt)+C22cos 2(λx-ωt)

(11)

C11=1

(12)

(13)

(14)

(15)

海床底面為不透水剛性邊界，土體位移和孔隙水壓力豎向梯度均為0，即?p/?z=u=w=0； 海床兩側(cè)為周期性邊界； 管線表面為不透水邊界，即?p/?n=0； 設管線與周圍海床土共同位移，即管土邊界的位移保持連續(xù)。

已有研究表明，海床液化后對波浪的傳播具有一定的衰減效應(Tong et al.，2018)，因而當海床液化后，液化區(qū)海床土的性質(zhì)以及海床表面的邊界條件都會發(fā)生一定的變化，如何考慮這種變化，尚待更深入的理論與試驗研究，本文同目前絕大部分文獻的處理方式相同(Ye et al.，2015; Foo et al.，2019)，暫對此不作考慮。

1.4 模型的求解

利用Comsol Multiphysics多場耦合有限元分析軟件建立數(shù)值計算模型并求解，耦合模型在求解過程中孔壓累積源項中的最大剪應力項隨著波浪作用時間而進行周期性更新。由于管線的存在會影響其周圍海床土體響應，因此管線及其周圍土體區(qū)域為重點研究區(qū)域，此部分的網(wǎng)格進行了局部加密，如圖 2 所示。

2 模型驗證

為了驗證數(shù)值計算模型的可靠性，將本文的計算結(jié)果同Turcotte et al. (1984)的波浪水槽模型試驗結(jié)果進行了對比，試驗中海床和管道參數(shù)見表 1，3組波浪參數(shù)分別標注在圖 3 中； Cheng et al.(1986)也利用與此試驗參數(shù)進行了數(shù)值計算結(jié)果的對比驗證，對比結(jié)果如圖 3 所示。

表 1 試驗參數(shù)Table 1 Parameters of the flume test

從圖 3 中可以看出，本文數(shù)值模型計算得到的結(jié)果與試驗結(jié)果(Turcotte et al.，1984)以及文獻中數(shù)值計算結(jié)果(Cheng et al.， 1986)吻合較好，表明了本文數(shù)值計算模型的可靠性。

圖 3 波浪作用下埋置管線周圍孔壓對比圖Fig. 3 Comparison of wave-induced pore pressure around a buried pipeline

3 結(jié)果與討論

考慮管線自重對海床初始應力場的影響，利用耦合模型和非耦合模型分別對波浪作用下含埋置管線的海床響應進行模擬計算，并通過對比分析，揭示管線周圍海床土響應特征與機制。根據(jù)近海區(qū)域可能出現(xiàn)的極端波浪條件和海床土特征(常方強， 2009； 李佳峻等， 2017)，數(shù)值模型計算參數(shù)見表 2。

表 2 數(shù)值模型計算參數(shù)Table 2 Parameters of the numerical model

圖 4 海床初始豎向有效應力分布圖Fig. 4 Isosurface of the initial vertical effective stress

圖 5 初始平均有效應力分布圖Fig. 5 Isosurface of the initial mean effective stress

3.1 管線影響下的海床初始應力場

為考慮管線自重對海床初始應力場的影響，首先計算重力作用下含埋置管線海床的初始應力場，結(jié)果如圖 4 和圖 5 所示，分別為管線周圍海床的初始豎向有效應力與初始平均有效應力等值面分布圖(應力以拉為正，壓為負)。

由圖 4 和圖 5 可以看出，受管線自重的影響，管線周圍海床的初始應力不再是簡單的線性分布，與無管線時相比，管線頂部和底部周圍海床的初始豎向有效應力和平均有效應力均有所增加，而在管線兩側(cè)靠近管線處則有所減小，例如，管線底部初始豎向有效應力由無管線時的25 kPa增加為現(xiàn)在的30 kPa，增加量為20%。海床的初始應力場對于累積孔壓的發(fā)展和海床液化都有重要影響，經(jīng)數(shù)據(jù)分析可知，在距離管線邊界約1.8 m處，管線自重對海床應力的影響小于5%，可知管線自重的影響范圍約為距管線中心2倍管徑(2D)的區(qū)域。

海床初始應力場會影響孔壓累積源項，進而會對累積孔壓的發(fā)展產(chǎn)生影響(Zhao et al.，2014)，因此，本文基于上述考慮管線自重影響的海床初始應力場進行累積孔壓的發(fā)展計算。

3.2 管線周圍海床累積孔壓響應特征與機制

3.2.1 海床累積孔壓發(fā)展特征

如圖 6 所示，為管線頂面、底面和側(cè)面各點以及相對應的遠端各點(遠點，即不受管線影響的區(qū)域)累積孔壓隨時間的發(fā)展曲線。耦合模型的計算結(jié)果表明，管線周圍各點海床的累積孔壓發(fā)展速度均遠快于同一深度遠端各點的累積孔壓發(fā)展速度； 而非耦合的計算結(jié)果表明，管頂接觸點(點1)位置累積孔壓的發(fā)展速度小于其對應遠點(點2)累積孔壓的發(fā)展速度； 管中接觸點(點3)與其對應遠點(點4)的累積孔壓發(fā)展速度差異不大； 管底接觸點(點5)累積孔壓的發(fā)展則快于其遠點(點6)處，至100個波浪作用周期(100T)時，兩者的累積孔壓相差約4 kPa。同時，可以發(fā)現(xiàn)，無論是管線周圍還是遠端各點，考慮累積孔壓與海床應力耦合效應時的累積孔壓發(fā)展速度均快于非耦合模型，且兩種模型的計算結(jié)果在管線周圍各點處的差異相對更為顯著，例如，在管底接觸點處，波浪作用100T時耦合模型的累積孔壓相對于非耦合模型增大了約19%，而在對應遠點處耦合模型的累積孔壓較非耦合模型增大約6%。

圖 6 不同位置處累積孔壓隨時間變化對比圖Fig. 6 Comparison of development of residual pore pressures at different locations

兩種計算模型下累積孔壓的水平分布和發(fā)展特征，可以在圖 7 中得到更為清晰的展示。如圖 7 所示，耦合模型下管線對海床累積孔壓的發(fā)展影響相對較大，一方面管線周圍累積孔壓的數(shù)值遠大于遠端各點，另一方面，水平方向上管線對累積孔壓的影響范圍隨著波浪作用時間的增加而有所增大，例如，當波浪作用時間為10T時，管底處的累積孔壓為13.6 kPa，遠點處為8.6 kPa，管底處的累積孔壓相對于遠點增大了58%，此時水平方向管線的影響范圍共約14D； 當波浪作用時間為40T時，管底處的累積孔壓為27 kPa，相應遠點的累積孔壓為16.6 kPa，管底處累積孔壓相對于遠點增大了63%，此時水平方向管線的影響范圍約32D。非耦合模型的計算結(jié)果表明，管線對累積孔壓的影響相對較小，首先，管線周圍累積孔壓的數(shù)值與遠端相比差異較小，其次，水平方向上管線的影響范圍也相對較小，且該影響范圍基本不受波浪作用時間的影響，例如，波浪作用時間為10T和40T時，管底處累積孔壓分別為9.4 kPa和19 kPa，相對于遠端累積孔壓均增加了約22%，且水平方向上管線的影響范圍約為9D。累積孔壓在管線頂面和中部水平截線上的分布特征與管線底部截面基本一致，只是在管線頂部端點處的累積孔壓相對有所減小，具體如圖 7 所示。

圖 7 管線附近海床的累積孔壓水平分布圖Fig. 7 Lateral distribution of residual pore pressure in the vicinity of the pipeline

耦合模型與非耦合模型下累積孔壓沿深度的分布特征如圖 8 所示，圖 8a和圖8b分別為耦合模型與非耦合模型下管線中心豎向截面與遠端海床豎向截面的累積孔壓分布曲線。從圖 8 中可以明顯看出，當考慮海床應力與累積孔壓的耦合效應時，管線周圍海床的累積孔壓明顯增大，且豎直方向上管線對累積孔壓的影響范圍隨著波浪作用時間的增加而有所增大，波浪作用為10T時管線對其下部海床累積孔壓的影響深度約為6D，當波浪作用時間為40T時，管線對其下部海床累積孔壓的影響深度擴大至18D。而當不考慮海床應力與累積孔壓的耦合效應時，管線對其下部海床累積孔壓的影響深度約為2D，其影響深度相對較小且基本不隨波浪作用時間而變化。

圖 8 管線中心豎向截面與海床遠端豎向截面的累積孔壓分布對比圖Fig. 8 Comparison of residual pore pressures along the vertical section across the central of the pipeline and that in the far field

沿管土交界面的海床累積孔壓分布如圖 9 所示，耦合模型與非耦合模型的計算結(jié)果均表明，管周下半部分海床的累積孔壓均較其上半部分發(fā)展更快，其中尤以管線底端點處孔壓的累積最為顯著，例如，考慮累積孔壓與海床應力的耦合效應， 30T時管線底端點累積孔壓相對于10T時累積孔壓的增加量為83%，而管線頂端點處的相應累積孔壓增加量為33%?？傮w而言，在波浪作用的前30個周期管周孔壓增加較快，之后隨著波浪作用時間的增加，海床累積孔壓的發(fā)展速度有所減緩，這從圖 6 中累積孔壓隨時間的發(fā)展變化曲線中也可以看出。

圖 9 累積孔壓沿管土交界面分布對比圖Fig. 9 Comparison of residual pore pressures along the interface of pipeline and soils

3.2.2 海床累積孔壓發(fā)展機制

由前述分析可知，當考慮累積孔壓與海床應力的耦合效應時，管線周圍和遠端海床的累積孔壓發(fā)展速度均有所增大，尤其是管線附近區(qū)域，累積孔壓的增大較為顯著。根據(jù)Liu et al. (2019)的分析可知，在耦合模型下，水平方向累積孔壓的不均勻分布會導致海床循環(huán)剪應力的不斷增加，從而會在較大程度上促進累積孔壓的快速發(fā)展。如前述圖 7 所示，在管線周圍海床的累積孔壓沿水平方向的差異相對較大，這會導致該區(qū)域剪應力的增加，如圖 10 所示，可以看出，隨著波浪作用周期的增加，管線周圍的海床剪應變(剪應力)也有所增大，而剪應變的增大會明顯促進累積孔壓的發(fā)展，因此，在管線附近區(qū)域海床的累積孔壓相較于遠端海床有明顯的增大。當不考慮累積孔壓與海床應力的耦合效應時，海床剪應力因不受累積孔壓的影響而保持不變，因此其累積孔壓的發(fā)展相對較慢，且管線對海床累積響應的影響范圍也相對較小。綜上所述，對于管線附近區(qū)域的海床，累積孔壓沿水平方向分布不均勻而導致的循環(huán)剪應力的增大，是耦合模型相對于非耦合模型累積孔壓明顯增大的主要原因。

圖 10 海床剪應變水平分布圖Fig. 10 Lateral distribution of soil shear strain

對于遠離管線的區(qū)域，管線對累積孔壓的影響可以忽略，這些區(qū)域的累積孔壓在水平方向上分布均勻，因而剪應力不隨波浪作用時間而變化； 這些區(qū)域內(nèi)耦合模型計算的海床累積孔壓大于非耦合模型計算得到的累積孔壓，其根本原因在于累積孔壓豎向分布的不均勻而導致的海床平均有效正應力的差異(Liu et al.，2019)，不同深度處耦合模型和非耦合模型平均有效正應力的差異曲線如圖 11 所示，深度越大，平均有效正應力的差異就越大，因而兩種模型計算所得到的累積孔壓的差異也就越大。

圖 11 遠點處海床平均有效正應力差值Fig. 11 Difference of the mean effective stress in the far field

3.3 含埋置管線海床漸進液化特征

當海床某深度位置的超靜孔隙水壓力大于其初始豎向有效應力時，該深度位置的土體會發(fā)生液化(Okusa， 1985)，液化準則的表達式為(有效應力以拉為正)：

pres+σ′z0>0

(16)

基于上述液化準則，對不同波浪作用時間下的海床液化區(qū)進行分析，結(jié)果如圖 12 所示。耦合模型下管線周圍海床的液化深度相對較大，隨著波浪作用時間的增加，液化區(qū)在豎向自上而下不斷發(fā)展，水平方向上液化區(qū)從管線附近逐漸向兩側(cè)不斷擴展，至20T時，管線頂部周圍海床發(fā)生液化； 至60T時，管線底部附近的海床開始液化，此時管線周圍海床的最大液化深度約為3.5 m。而對于非耦合模型，由于累積孔壓發(fā)展相對較慢，相應波浪作用時間海床尚未液化。因此，當不考慮累積孔壓與海床應力的耦合效應時，可能會一定程度地低估管線周圍海床的液化深度，不利于管線的安全設計。

圖 12 管線埋深2 m時海床的液化區(qū)分布圖Fig. 12 Liquefaction zone of the seabed around the pipeline at a buried depth of 2 m

管線的埋深對其周圍海床累積孔壓的發(fā)展和液化具有重要影響(Zhao et al.，2014)，一般而言，管線埋置深度越深，其周圍海床液化發(fā)展速度越慢。如圖 13 所示，為管線埋置深度為3 m時不同波浪作用時間下海床的液化區(qū)分布，在20T時，液化深度尚未達到管線頂端，而在60T時，管線底部區(qū)域的海床也尚未全部液化，與管線埋置深度為2 m時的液化區(qū)分布圖對比可知，管線埋置深度越深，越有利于管線的安全。

圖 13 管線埋深為3 m時海床的液化區(qū)分布圖Fig. 13 Liquefaction zone of the seabed around the pipeline at a buried depth of 3 m

需要注意的是，上述分析是基于滲透系數(shù)較小的相對松散的深厚均質(zhì)粉砂或細砂層地基，且波浪條件相對極端情況下的計算結(jié)果，由于海床滲透系數(shù)較小且相對松散，孔壓極易累積，因此，液化深度相對較大。在實際工程中，宜充分考慮管線周圍海床的實際地層條件，準確確定計算參數(shù)，在保證管線安全的前提下避免設計過度保守。

4 結(jié) 論

對非線性行進波作用下含埋置管線海床的累積響應特征和液化進行了數(shù)值模擬分析，探討了孔壓累積與海床應力耦合效應下管線周圍海床累積響應特征及其發(fā)展機制，并與不考慮孔壓累積與海床應力耦合效應的結(jié)果進行了對比分析，主要得到以下結(jié)論：

(1)行進波作用下，埋置管線附近海床的累積孔壓發(fā)展相對較快，其中管周下半部分海床的累積孔壓發(fā)展速度較其上半部分發(fā)展更快，尤以管線底端點的累積孔壓發(fā)展最為顯著。

(2)孔壓累積與海床應力的耦合效應會顯著增大管線對其周圍海床累積孔壓特征的影響，這種影響一方面體現(xiàn)在管線周圍海床累積孔壓的更快發(fā)展，另一方面體現(xiàn)在管線對孔壓影響范圍的擴大。結(jié)果分析表明，考慮累積孔壓與應力的耦合效應時，管線的影響范圍隨波浪作用時間會有所增大； 當忽略累積孔壓與應力的耦合效應時，管線的影響范圍基本不隨波浪作用時間而變化。

(3)由于累積孔壓與應力的耦合效應，管線附近累積孔壓的不均勻分布會引起海床循環(huán)剪應力的增大，從而導致管線周圍海床累積孔隙水壓力的快速發(fā)展，并促進海床的液化，忽略這種耦合效應會在一定程度上低估海床液化深度，不利于管線安全。實際工程設計中，應充分考慮管線周圍海床的實際地層條件，合理確定計算參數(shù)，保證管線的安全。