向量式有限元法在管土相互作用中的應用

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

（1.山東科技大學 土木工程與建筑學院，山東 青島266590；2.山東科技大學 山東省土木工程防災減災重點實驗室，山東 青島266590；3.中國海洋大學 工程學院，山東 青島266100）

0 引 言

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

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

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

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

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

1.1 彈簧模型

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

圖1 彈簧模型Fig.1 Spring model

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

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

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

1.2 向量式有限元分析模型

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

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

1.3 質(zhì)點控制方程

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

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

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

1.4 控制方程的求解

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

代入（4）式得：

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

取一個標準途徑單元tn-1≤t≤tn，桿件αβ 在tn時刻的兩端點位置向量是，方向向量為：

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

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

當質(zhì)點是αβ的起始點時τ=1，否則τ=2。

1.6 質(zhì)點所受外力

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

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

2 數(shù)值求解

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

圖3 單元虛擬運動示意圖Fig.3 Fictitious motion of the element

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

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

根據(jù)文獻[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，來代替懸垂段運動對立管觸地段的影響。因為主要變形發(fā)生在靠近自由端的部分，故圖像只取靠近自由端的30 m。

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

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

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

圖7 線性截斷彈簧模型下立管變形過程Fig.7 Riser deformation based on tension cut-off spring model

圖8 線性截斷彈簧模型下立管彎曲應力變化過程Fig.8 The variation of bending stresses based on tension cut-off spring model

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

3 參數(shù)分析

3.1 彈簧模型的影響

圖9 是采用向量式有限元法模擬三種不同彈簧模型得到的位形圖。從圖中可以看出線性彈簧模型溝槽最深處最接近自由端，因為線性彈簧模型是線性的，沒有考慮管土分離，所以立管一直承受海床土體吸力，溝槽最深處距離自由端最近，所取得的最大彎曲應力值也是三種模型中最大的，如圖10所示。相比之下，非線性彈簧模型雖然也未考慮管土分離，但假定了最大土體吸力以及最大抗力，所以溝槽最深處離自由端遠些；線性截斷彈簧模型考慮了土體最大抗力以及管土分離，當土體吸力達到最大值時，立管不再受吸力，所以溝槽最深處距自由端最遠，其最大彎曲應力也是三者中最小的。可以看出線性截斷彈簧模型更符合海床土體實際情況。

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

圖10 不同模型的立管彎曲應力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ù)的彎曲應力Fig.12 The bending stress of the riser with soil suction force coefficient

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

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

3.3 海床土體剛度的影響

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

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

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

4 結(jié) 論

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

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

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

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