王 飛，李效民，馬芳俊，郭海燕

（1.山東科技大學(xué) 土木工程與建筑學(xué)院，山東 青島266590；2.山東科技大學(xué) 山東省土木工程防災(zāi)減災(zāi)重點(diǎn)實(shí)驗(yàn)室，山東 青島266590；3.中國海洋大學(xué) 工程學(xué)院，山東 青島266100）

0 引 言

隨著深海油氣資源的開發(fā)，鋼懸鏈線立管（SCR）已經(jīng)成為開發(fā)深海油氣資源的主要立管形式。鋼懸鏈線立管觸地點(diǎn)承受較大的彎曲應(yīng)力，與海床接觸區(qū)域易發(fā)生疲勞破壞。因此對于SCR 觸地段與海床間的相互作用分析越來越受到人們的重視。

目前，國內(nèi)外專家對立管與海床間的相互作用進(jìn)行了一些研究，主要有有限元法、有限差分法或通過實(shí)驗(yàn)研究的方法。Pesce 等[1]將海床等效為線性彈簧模型，分析了觸地點(diǎn)處海床硬度對立管運(yùn)動的影響。2H 海洋工程公司發(fā)起了工業(yè)聯(lián)合開發(fā)計(jì)劃（STRIDE JIP）[2]，對管土相互作用進(jìn)行了全尺寸模型試驗(yàn)。You[3]通過有限差分法對管土相互作用進(jìn)行了分析。Aubeny 等[4]根據(jù)小尺度模型實(shí)驗(yàn)提出了基于P-y 曲線的管土作用模型以及各個階段的經(jīng)驗(yàn)公式。Hu 等[5]通過試驗(yàn)表明立管反復(fù)作用在海床上會導(dǎo)致土體剛度衰減。郭海燕等[6]應(yīng)用ANSYS 中的接觸單元建立模型對立管-海床的作用進(jìn)行了分析。王坤鵬等[7]用ABAQUS 創(chuàng)建觸地單元分析管土相互作用。Wang 等[8]通過室內(nèi)大比尺模型試驗(yàn)研究觸地區(qū)的埋置深度對管內(nèi)應(yīng)力的影響。白興蘭等[9]基于改進(jìn)的程序CABLE3D RSI 對非線性管土作用下鋼懸鏈線立管觸地區(qū)疲勞進(jìn)行了分析。李敢等[10]用ANSYS 對管土相互作用引起的疲勞進(jìn)行了分析。

盡管人們對于管土相互作用已經(jīng)做了大量工作，但為保證立管在位運(yùn)行的安全性，在設(shè)計(jì)分析過程中運(yùn)用多種數(shù)值計(jì)算方法和模型進(jìn)行驗(yàn)證是極其必要的。由丁承先教授[11]提出的向量式有限元法，又稱為有限質(zhì)點(diǎn)法，以向量力學(xué)作為分析的理論基礎(chǔ)，并與數(shù)值計(jì)算相結(jié)合，是求解結(jié)構(gòu)的大變形、大變位、彈塑性、碰撞、倒塌等非線性或不連續(xù)性力學(xué)行為的新方法。本文將其應(yīng)用于鋼懸鏈線立管觸地段與海床土體相互作用的分析研究中。相較于傳統(tǒng)分析方法，向量式有限元是對立管的力學(xué)行為建立計(jì)算理論，計(jì)算中無需集成結(jié)構(gòu)的剛度矩陣，可方便更改單元和邊界條件，基于中央差分的顯式積分求解方案無需迭代求解非線性方程組，從而能夠降低立管分析的難度，同時分析初始時刻便可以將荷載直接加載到立管上，即使是對靜力問題的分析，亦是通過逐步施加外部載荷并不斷改變立管的位置和形態(tài)直至達(dá)到最終結(jié)果，因此可以更好地模擬真實(shí)海洋環(huán)境下立管的運(yùn)動變形過程，使分析相對準(zhǔn)確。

本文基于向量式有限元法采用彈簧模型模擬海床，將懸鏈線立管觸地段離散成一系列質(zhì)點(diǎn)的集合，質(zhì)點(diǎn)間用只承受內(nèi)力而無質(zhì)量的平面彎曲桿件元連接，采用一系列途徑單元描述立管的受力變形過程，通過虛擬的逆向運(yùn)動計(jì)算立管的純變形和內(nèi)力。采用牛頓第二定律來表達(dá)質(zhì)點(diǎn)的運(yùn)動行為，并用中央差分的顯式積分法求解質(zhì)點(diǎn)各時刻的控制方程，編制相應(yīng)Matlab 求解程序，驗(yàn)證該方法的可行性并分析管土相互作用。

1 鋼懸鏈線立管觸地段管土相互作用分析

1.1 彈簧模型

鋼懸鏈線立管觸地段與海床相互作用可用一系列彈簧模擬海床，從而可以將P-y 曲線[4]簡化為三種彈簧模型，如圖1 所示分別為線性彈簧模型、非線性彈簧模型和線性折斷彈簧模型。線性彈簧模型假設(shè)海床土體是線性的，立管觸地段所受土體抗力與土體剛度和入土深度成正比；非線性彈簧模型考慮了土體的最大土體抗力與土體最大吸引力；線性截?cái)鄰椈赡Ｐ涂紤]到管土分離后，不存在土體吸力。

圖1 彈簧模型Fig.1 Spring model

海床初始嵌入骨架曲線為：

其中：Su0為海床表面抗剪強(qiáng)度，Sg為抗剪強(qiáng)度變化梯度，a 和b 隨著溝槽形狀和管道的粗糙度變化，y為溝槽深度，D 為管道直徑。

定義最大土體吸力與極限承載力的比值：

1.2 向量式有限元分析模型

如圖2 所示，鋼懸鏈立管觸地段管土模型可以簡化為立管放置在一系列彈簧上，采用向量式有限元法進(jìn)行分析。立管左端自由，右端鉸接，立管觸地段整體坐標(biāo)系定義如下：取自由端為坐標(biāo)原點(diǎn)，x 軸沿立管觸地段軸線，y 軸為水深方向。采用N+1（從左向右編號依次為1，2，…，N，N+1）個等間距的質(zhì)點(diǎn)模擬立管的直線狀態(tài)，質(zhì)點(diǎn)之間利用平面彎曲桿單元相連接，單元的質(zhì)量平均分配到兩端的質(zhì)點(diǎn)上，立管的初始張力通過施加單元的初始軸力來實(shí)現(xiàn)。對觸地段與海床土體的靜力相互作用展開研究，通過對觸地段左端逐步施加豎向位移代表鋼懸鏈立管懸垂段的變形運(yùn)動過程，可以使1 號節(jié)點(diǎn)在R 時間內(nèi)運(yùn)動到指定點(diǎn)，并在R 時間內(nèi)將重力、浮力及土支撐力以斜坡函數(shù)加載的方式作用在各質(zhì)點(diǎn)上，總共作用T 時間，使管道達(dá)到平衡狀態(tài)。

圖2 管土相互作用模型及坐標(biāo)示意圖Fig.2 Diagram of the riser deformation and coordinate

1.3 質(zhì)點(diǎn)控制方程

對于平面問題，質(zhì)點(diǎn)有三個位移分量，亦即兩個平移量和一個轉(zhuǎn)動量，考慮阻尼的影響，質(zhì)點(diǎn)控制方程式為：

其中：Mi為i 質(zhì)點(diǎn)的質(zhì)量，為i 質(zhì)點(diǎn)加速度，分別為i 質(zhì)點(diǎn)所受的外力、內(nèi)力和阻尼力。其展開為：

其中：mi是質(zhì)點(diǎn)i 的質(zhì)量；Ii是其對z 軸的質(zhì)量慣性矩；βi是轉(zhuǎn)角；ζi為立管的阻尼系數(shù)；分別為質(zhì)點(diǎn)所受到的沿x、y 軸的外力及繞z 軸的外力偶矩，分別為質(zhì)點(diǎn)所受到的沿x、y 軸的內(nèi)力及繞z 軸的內(nèi)力偶矩。

1.4 控制方程的求解

采用中央差分的顯式時間積分法，質(zhì)點(diǎn)的速度和加速度可表示為：

代入（4）式得：

1.5 單元內(nèi)力計(jì)算

取一個標(biāo)準(zhǔn)途徑單元tn-1≤t≤tn，桿件αβ 在tn時刻的兩端點(diǎn)位置向量是，方向向量為：

以tn-1時刻單元內(nèi)力和變形作為基礎(chǔ)架構(gòu)，通過虛擬的逆向運(yùn)動計(jì)算立管下一時刻內(nèi)力和變形純增量，如圖3（a）所示，然后根據(jù)材料力學(xué)的理論求得桿件元αβ 在tn時刻單元內(nèi)力為：

最后通過虛擬的正向運(yùn)動如圖3（b）回到原先tn時刻位置。兩個節(jié)點(diǎn)內(nèi)力分量分別轉(zhuǎn)換成域坐標(biāo)分量和，彎矩不需轉(zhuǎn)換，將各單元節(jié)點(diǎn)力作用在各質(zhì)點(diǎn)，再對質(zhì)點(diǎn)內(nèi)力進(jìn)行集成，得到：

當(dāng)質(zhì)點(diǎn)是αβ的起始點(diǎn)時τ=1，否則τ=2。

1.6 質(zhì)點(diǎn)所受外力

立管質(zhì)點(diǎn)所受外力包括重力、浮力以及海床土體支撐力。表達(dá)式可以寫為：

其中：ω 為立管濕重，k 為土體剛度。

2 數(shù)值求解

基于上文所述理論及計(jì)算公式，采用Matlab編寫求解程序， 分析管土相互作用， 并與用ABAQUS[12]軟件分析的結(jié)果和You[3]基于有限差分法編制的程序進(jìn)行對比分析，計(jì)算流程如圖4所示。本文模型參數(shù)采用STRIDE JIP 實(shí)驗(yàn)[2]中的模型數(shù)據(jù)，如表1 所示。

圖3 單元虛擬運(yùn)動示意圖Fig.3 Fictitious motion of the element

圖4 向量式有限元法計(jì)算流程Fig.4 Flow chart of pipe/seafloor interaction analysis using the VFIFE

表1 模型參數(shù)Tab.1 Model parameters

根據(jù)文獻(xiàn)[5]取a=6.73，b=0.15，溝槽深度y=D（D 為立管外徑），則土體極限承載力為Pmax=3.96 kN/m，取m=0.5，則最大吸力為：Pt=0.5×Pmax=1.98 kN/m。在鋼懸鏈立管觸地段左端施加豎向位移D，來代替懸垂段運(yùn)動對立管觸地段的影響。因?yàn)橹饕冃伟l(fā)生在靠近自由端的部分，故圖像只取靠近自由端的30 m。

圖5 是分別采用線性彈簧模型、非線性彈簧模型以及線性截?cái)鄰椈赡Ｐ屯ㄟ^向量式有限元法(VIFIFE)、有限差分法(FDM)和ABAQUS 得到的位形變化對比圖。從圖中可以看出三種方法得到的三種彈簧模型下的位形結(jié)果高度吻合，管道在靠近自由端的位置處向上位移，沒有滲入海床，管土發(fā)生分離，隨著遠(yuǎn)離自由端管道入土深度逐漸加大，直到達(dá)到最深處，然后逐漸變小趨于穩(wěn)定。圖6 是分別在三種模型下用三種方法模擬得到的管道彎曲應(yīng)力分布圖，彎曲應(yīng)力在管道觸地點(diǎn)急劇增大到最大值，隨后迅速變小，在遠(yuǎn)離自由端后趨于平穩(wěn)，數(shù)值很小，幾乎可以忽略。從圖5-6 兩對比圖可以看出，三種方法吻合良好，證明向量式有限元法在管土相互作用模擬分析中是可行的。

圖5 三種彈簧模型位形對比圖Fig.5 Comparison of configurations for three spring model

圖6 三種彈簧模型彎曲應(yīng)力對比圖Fig.6 Comparison of bending stresses for three spring model

圖7 線性截?cái)鄰椈赡Ｐ拖铝⒐茏冃芜^程Fig.7 Riser deformation based on tension cut-off spring model

圖8 線性截?cái)鄰椈赡Ｐ拖铝⒐軓澢鷳?yīng)力變化過程Fig.8 The variation of bending stresses based on tension cut-off spring model

圖7 是采用線性截?cái)鄰椈赡Ｐ湍M海床時立管觸地段位形變化過程圖。觸地段初始狀態(tài)為靜態(tài)直線狀態(tài)，然后在觸地段左端逐步施加豎向位移使管道從水平位置移動到D 位置，最終管道達(dá)到靜止?fàn)顟B(tài)。因此向量式有限元法分析過程是一個動態(tài)的實(shí)現(xiàn)過程，是對立管力學(xué)行為過程的模擬，更加符合立管的真實(shí)變化過程；圖8 是與之對應(yīng)的彎曲應(yīng)力變化過程。相比于基于數(shù)學(xué)模型的有限差分法及其他有限元軟件，向量式有限元法基于物理模型，計(jì)算中無需集成結(jié)構(gòu)的剛度矩陣，也無需迭代求解，可方便更改單元和邊界條件，從而簡化了立管分析的難度。即使是對于靜力位形及內(nèi)力的計(jì)算亦是通過逐步動態(tài)實(shí)現(xiàn)的過程，變形及受力符合立管實(shí)際受力情況，所得結(jié)果相對精確。

3 參數(shù)分析

3.1 彈簧模型的影響

圖9 是采用向量式有限元法模擬三種不同彈簧模型得到的位形圖。從圖中可以看出線性彈簧模型溝槽最深處最接近自由端，因?yàn)榫€性彈簧模型是線性的，沒有考慮管土分離，所以立管一直承受海床土體吸力，溝槽最深處距離自由端最近，所取得的最大彎曲應(yīng)力值也是三種模型中最大的，如圖10所示。相比之下，非線性彈簧模型雖然也未考慮管土分離，但假定了最大土體吸力以及最大抗力，所以溝槽最深處離自由端遠(yuǎn)些；線性截?cái)鄰椈赡Ｐ涂紤]了土體最大抗力以及管土分離，當(dāng)土體吸力達(dá)到最大值時，立管不再受吸力，所以溝槽最深處距自由端最遠(yuǎn)，其最大彎曲應(yīng)力也是三者中最小的。可以看出線性截?cái)鄰椈赡Ｐ透虾４餐馏w實(shí)際情況。

圖9 不同模型的立管位形Fig.9 The configuration of the riser with different spring model

圖10 不同模型的立管彎曲應(yīng)力Fig.10 The bending stress of the riser with different spring model

圖11 不同吸力系數(shù)的立管位形Fig.11 The configuration of the riser with different soil suction force coefficient

圖12 不同吸力系數(shù)的彎曲應(yīng)力Fig.12 The bending stress of the riser with soil suction force coefficient

3.2 海床土體吸力系數(shù)的影響

采用更符合實(shí)際情況的線性截?cái)鄰椈赡Ｐ湍M海床，在立管觸地段自由端施加向上位移D，改變吸力系數(shù)m，分別取0、0.25、0.5、0.75 和1，得到不同立管位形和應(yīng)力，如圖11-12 所示。從圖中可以看出隨著吸力系數(shù)m 的增加，觸地段所形成的溝槽最深處越靠近自由端，同時彎曲應(yīng)力也相應(yīng)變大。說明海床土體吸力增大了立管觸地段的彎曲應(yīng)力。

3.3 海床土體剛度的影響

在立管觸地段自由端施加向上位移D，改變海床土體剛度，分別取100、200、300、400 和500 kN/m/m，得到不同剛度下立管位形和應(yīng)力，如圖13-14 所示。隨著海床土體剛度的增大，溝槽最深處距離自由端越遠(yuǎn)，相應(yīng)的彎曲應(yīng)力也逐漸變小。當(dāng)土體剛度小時，其對應(yīng)力的影響很明顯，但當(dāng)剛度變大時，對最大應(yīng)力影響變小。說明土體剛度變大會使溝槽最深處遠(yuǎn)離自由端，應(yīng)力變小，土體剛度變大達(dá)到一定值后，應(yīng)力值應(yīng)該趨于穩(wěn)定，接近剛性海床。

圖13 不同海床土體剛度的立管位形Fig.13 The configuration of the riser with different soils stiffness

圖14 不同海床土體剛度的彎曲應(yīng)力Fig.14 The bending stress of the riser with different soils stiffness

4 結(jié) 論

本文采用彈簧模型模擬海床，將向量式有限元法用于鋼懸鏈線立管觸地段與海床土體相互作用的模擬中，并與ABAQUS 及有限元差分法模擬的結(jié)果作對比，并分析了不同參數(shù)對立管觸地段與海床相互作用的靜力影響。通過本文研究可以得出以下結(jié)論：

（1）三種方法模擬的結(jié)果吻合，表明向量式有限元法在管土相互作用的模擬中是可行的；

（2）即使是對于管土相互作用的靜力問題分析，向量式有限元亦是動態(tài)實(shí)現(xiàn)的過程，即由靜止開始到運(yùn)動再到靜止結(jié)束的過程。分析過程中可以任意改變荷載，可以更好地模擬真實(shí)海洋環(huán)境下立管的運(yùn)動行為，采用中央差分的顯式時間積分求解，不必形成剛度矩陣，避免了迭代求解和收斂問題，計(jì)算中容易改變單元與邊界條件，計(jì)算程序中只存在質(zhì)點(diǎn)位置計(jì)算和單元內(nèi)力計(jì)算兩個循環(huán)，整體步驟呈系統(tǒng)化；能夠直接求解立管的變形和內(nèi)力，避免了由變形計(jì)算內(nèi)力的復(fù)雜過程和可能的誤差；

（3）向量式有限元法可以有效地分析不同吸力系數(shù)和海床剛度時的管土相互作用，能夠很好地契合三種模型，但還需進(jìn)一步驗(yàn)證該方法在循環(huán)作用下土體剛度的衰減P-y 曲線模型中的應(yīng)用。