拖曳錨在海床中運動特性的大變形有限元分析

李培冬，劉海笑，趙燕兵

(天津大學 建筑工程學院，天津 300072)

拖曳錨在海床中運動特性的大變形有限元分析

李培冬，劉海笑，趙燕兵

(天津大學 建筑工程學院，天津 300072)

拖曳錨由于其承載性能和深水中便于安裝被廣泛應用于海洋工程系泊系統中，如：適用于懸鏈式系泊系統的傳統拖曳錨和適用于繃緊式系泊系統的法向承力錨。拖曳錨安裝過程中涉及諸多運動特性：錨板運動方向、系纜點處拖曳力和拖曳角及運動軌跡。基于大變形有限元分析技術耦合的歐拉-拉格朗日法，并引入纜繩方程，建立起錨-纜繩-海床土耦合作用的有限元分析模型；模擬了拖曳錨在均質和線性強度黏土中的嵌入安裝過程，研究了錨板運動方向、系纜點處拖曳力和拖曳角及運動軌跡等運動特性；通過與已有的有限元分析方法及理論方法進行對比，驗證了該分析模型的有效性；與已有的有限元分析方法相比，提出的分析模型有效地提高了計算效率。

拖曳錨；運動特性；耦合的歐拉-拉格朗日法；纜繩方程；大變形有限元

Abstract:Because of better performance both in pullout capacity and deepwater installation,drag anchors are widely utilized in mooring systems for offshore applications,such as the conventional drag anchor in the catenary mooring system and the vertically loaded plate anchor in the taut-wire mooring system.The penetration behaviors,such as the movement direction of the fluke,the drag force and drag angle at the shackle and the trajectory,are involved during the anchor installation.A large deformation finite element analysis using a coupled Eulerian-Lagrangian approach is performed to simulate the installation process of drag anchors in the clay with uniform and linear strengths,in which a chain equation is introduced to reflect the interaction between the installation line and the anchor.By comparing the results of the published FE method and theoretical method,including the movement direction of the fluke,the drag force and drag angle at the shackle and the trajectory,the efficiency of present work is well verified.Compared with the published FE method,the FE model of present work has a better performance in computational efficiency.

Keywords:drag anchor; kinematic behavior; coupled Eulerian-Lagrangian; chain equation; large deformation FE method

隨著海洋油氣開發不斷向深水及超深水發展，大型浮式海洋平臺及適用于深水油氣開發的繃緊式系泊方式得到了日趨廣泛的應用，同時也對深水系泊系統提出了更高的要求。法向承力錨具有承載性能優異和便于深水安裝的特點，因而被廣泛應用在深水系泊系統中[1]。

法向承力錨與傳統拖曳錨的安裝方式基本相同，拖曳錨與拖船通過纜繩連接，拖曳錨被拋錨至海床面后，拖船沿一定方向移動，拖曳錨逐漸嵌入海床土直至目標深度，同時錨板方位角逐漸減小[1]。在拖曳錨的整個安裝過程中，纜繩形態及通過纜繩傳遞的拖曳力也在不斷發生變化。拖曳錨安裝完成后，其承載性能高度依賴于最終嵌入深度和錨板方位角。嵌入深度和錨板方位角由拖曳錨的運動軌跡決定，而運動軌跡在安裝過程中是不可見的，因此研究錨的嵌入運動軌跡具有重要意義。

拖曳錨嵌入過程中，除運動軌跡[2-5]外，還涉及諸如錨板運動方向[6-8]、系纜點處拖曳力[9-11]和拖曳角[2,7,11]等運動特性，這些運動特性直接影響錨的嵌入運動軌跡。目前，針對拖曳錨運動特性的研究主要包括數值增量求解法、理論解析法和試驗研究，Zhang[1]對此作了詳細的總結。隨著有限元技術的發展，越來越多的研究者將有限元技術應用于錨的數值分析中，但對于拖曳錨運動特性的研究較少。有限元分析拖曳錨運動特性有兩個難點：其一，錨與土體的耦合作用會使土體產生極大變形，傳統的有限元方法在模擬大變形問題時會產生網格畸變，甚至出現計算不收斂的情況；其二，在安裝過程中根據纜繩形態可以把纜繩分成懸鏈纜、臥底纜及嵌入纜三部分[4]，嵌入纜形成的反懸鏈形態對拖曳錨的嵌入性能影響顯著，在有限元分析中必須采用恰當的方法來模擬嵌入纜的反懸鏈特性。

Zhao和Liu[12]基于大變形有限元分析技術耦合的歐拉-拉格朗日法(coupled Eulerian-Lagrangian,CEL)，建立了能夠模擬拖曳錨在海床土中嵌入安裝過程的有限元分析模型。CEL法克服了傳統有限元方法在模擬大變形問題時常見的接觸問題和網格畸變問題。通過與已有理論預測結果對比，驗證了CEL法分析拖曳錨運動特性的有效性。為模擬嵌入纜的反懸鏈特性，Zhao和Liu通過LINK連接單元連接柱形塊體來構建安裝纜繩，安裝纜繩與錨之間同樣通過LINK連接單元連接。對比結果驗證了這種方法在模擬嵌入纜反懸鏈特性時的有效性[12,13]，但缺點在于：纜繩的存在使得模型的網格規模擴大，由此導致模型的計算效率較低。因此，提出一種能夠有效提高計算效率的分析模型勢在必行。

本文基于大變形有限元分析技術耦合的歐拉-拉格朗日法，引入纜繩方程，建立錨-纜繩-海床土耦合作用的有限元分析模型。其中，系纜點處纜繩與錨的相互作用由纜繩方程來描述，并通過ABAQUS用戶子程序VUAMP來實現。分析模型模擬了拖曳錨在均質和線性強度黏土中的嵌入安裝過程，研究了錨板運動方向、系纜點處拖曳力和拖曳角及運動軌跡等運動特性。通過與已有的有限元分析方法進行對比，來驗證該分析模型的有效性和在計算效率方面的優勢。

1 數值方法

1.1耦合的歐拉-拉格朗日法(CEL)

有限元分析中，描述物體運動和變形的方法有兩種：拉格朗日法和歐拉法。拉格朗日法(Lagrangian)是一種網格依附于材料的計算方法，網格隨材料變形而變形。網格內充滿單一材料，能夠很好地追蹤材料邊界的變化。但對于大變形問題，網格會發生嚴重變形，將導致計算不收斂。歐拉法(Eulerian)采用網格固定而材料可在網格中自由流動的方式，網格不隨材料變形而變形，克服了拉格朗日法網格畸變所帶來的數值奇異，但歐拉法不能很好地追蹤材料邊界的變化。

為了克服傳統有限元法的缺陷，大變形有限元法CEL被提出。該方法最先由Noh[14]提出，Benson和Okazawa[15,16]又進一步將其完善。CEL法集合了歐拉算法和拉格朗日算法的優點。在分析土和結構物耦合作用時，土體采用歐拉算法，即網格固定，土體在網格中流動，不存在網格畸變問題。通過計算每個單元的歐拉材料體積分數(EVF)來追蹤土體在網格中的流動變形，EVF=1表示網格被土體填滿，EVF=0表示網格中沒有土體。結構采用拉格朗日算法，通過結構的邊界來追蹤土與結構接觸面的變化。目前，CEL法已經被廣泛應用于模擬結構物與土體耦合作用的大變形問題[12,17-20]。

1.2子程序VUAMP

有限元軟件ABAQUS提供了FORTRAN語言編寫的子程序接口，供用戶根據需求進行二次開發。為了克服直接構造安裝纜繩造成的計算耗時，本文通過纜繩方程來描述系纜點處纜繩與錨的相互作用，即在拖曳錨系纜點處施加不斷變化的拖曳力。VUAMP子程序接口使得用戶可以定義相關變量的幅值曲線，如載荷條件、邊界條件和預定義場等，描述這些變量隨時間或位移等參數的變化情況[21]。

為了實現在拖曳錨系纜點處施加不斷變化的拖曳力，需要定義集中載荷隨時間變化的幅值曲線。現以本文所研究問題為例，介紹VUAMP使用過程中的關鍵技術：

1)定義與集中載荷(拖曳力)相關的歷史輸出變量，并通過傳感器(Sensor)將相應的輸出變量傳遞給子程序VUAMP，用于計算下一增量步的載荷值。拖曳錨系纜點處拖曳力與系纜點的水平和豎向位移相關，因此，需將系纜點的水平和豎向位移作為輸出變量；

2)創建集中載荷(拖曳力)，并定義相應的用戶幅值(User amplitude)，在定義用戶幅值時指定所需狀態變量(Svars)的數目，狀態變量用于在兩次調用子程序間傳遞變量；

3)ABAQUS在每個增量步(Increment)計算開始時，讀取子程序VUAMP設置的幅值，進而更新系纜點處的集中載荷(拖曳力)。

2 數值方法流程

2.1纜繩方程

如圖1所示，纜繩嵌入海床土中的部分，即嵌入纜，由于土抗力和摩擦力的共同作用形成反懸鏈形態。嵌入纜與錨在系纜點處發生相互作用，通過纜繩傳遞的拖曳力Ta以及拖曳力與水平方向的夾角θah直接影響錨在海床土中的運動特性。拖曳錨安裝過程中，θah和Ta逐漸增大，使得作用在錨上的豎向力逐漸增大，促使錨板逐漸抬平。本文通過引入纜繩方程來模擬纜繩與土、纜繩與錨的耦合作用。

圖1 安裝纜繩示意Fig.1 Definition of installation line

Neubecker和Randolph[10]、Zhang等[11]提出了系纜點處拖曳力Ta與θah的關系式：

式中：μ為纜繩-土摩擦系數，在飽和黏土中，Neubecker和Randolph[10]、Zhang等[11]分別建議μ取0.4～0.6和0.1～0.6，DNV[22]建議鋼纜取0.1～0.3，鐵鏈取0.6～0.8；θ為某一深度z處纜繩與水平面的夾角；za為系纜點深度；En為有效承載系數，對于拖纜，Degenkamp等[23]基于模型試驗，建議鋼纜取1，鐵鏈取2.5；d為纜繩直徑；Nc為纜繩端阻力系數，取值范圍為7.6～14[22-25]；su為土體不排水抗剪強度，表示為：

式中：su0為海床面土體強度;k為土體強度梯度，對于均勻土，k=0。

對式(1)由海床面到系纜點深度za積分得：

式中：θe為海床面處纜繩與水平面的夾角，本文假設拖曳錨安裝過程中海床面處始終存在臥底纜，即θe=0°。隨著拖曳錨的不斷嵌入，系纜點處拖曳力Ta逐漸增大，通過式(3)求解相應的θah，然后在分析模型中將Ta以與水平面成一定角度θah的方式施加在系纜點處，從而直接影響拖曳錨的運動行為。

2.2實施流程

由于ABAQUS中只能加載水平和豎向載荷，因此，將作用在系纜點的拖曳力Ta沿水平方向和豎直方向分解(如圖2)，載荷值由式(4)決定。Tx和Tz在嵌入過程中不斷變化，需要通過子程序VUAMP定義Tx和Tz隨嵌入過程變化的幅值曲線。式(4)中的θah為拖曳力與水平方向的夾角，在子程序中通過式(3)更新。

有限元模型的分析流程如下：

1)給定初始時刻拖曳錨系纜點處的拖曳力Ta0和拖曳角θah0；

2)由式(4)計算作用在系纜點處的水平和豎向集中載荷及載荷幅值；

3)執行一個增量步計算，在下一個增量步計算開始之前調用子程序VUAMP，提取系纜點深度z，并利用式(3)計算新的θahi；

4)每隔一定時間tc，計算該段時間內系纜點水平運動的平均速度vd，比較控制速度vc和vd的相對大小，用于判斷是否更新系纜力。若vd≥vc，保持系纜力不變；若vd

5)判斷錨是否達到極限嵌入條件，是則停止計算，否則重復執行2)～5)。

整個分析流程在ABAQUS中通過調用子程序實現，不需要人為干預。在下文分析模型的驗證中，tc取0.01 s，vc取0.5 m/s，初始時刻系纜力Ta0取1kN，θah0取0°，每次更新系纜力時ΔT取0.1kN。

圖2 系纜點處拖曳力Fig.2 Drag force at the shackle

圖3 安裝纜繩模型[12]Fig.3 The installation line in the model[12]

3 數值方法驗證

為檢驗本文提出的有限元分析模型，將其應用于模擬拖曳錨在均質和線性強度黏土中的嵌入安裝過程，分析嵌入過程中錨板運動方向、系纜點處拖曳力和拖曳角及運動軌跡等運動特性，并與Zhao和Liu[12]的有限元分析結果對比。Zhao和Liu通過LINK連接單元連接柱形塊體來構建安裝纜繩，如圖3所示，由此來模擬嵌入纜在海床土中的反懸鏈特性，纜繩的存在嚴重降低了模型計算效率。本文通過引入式(3)來模擬纜繩與土、纜繩與錨的耦合作用，從而提高模型計算效率。

3.1數值模型

圖4 有限元模型Fig.4 FE model

本文通過纜繩方程來模擬纜繩與土、纜繩與錨的耦合作用，見式(3)。式(3)中參數的取值均與Zhao和Liu一致,具體取值見表3。由于不對纜繩建模，本文分析模型的網格規模較Zhao和Liu的小，有效提高了模型的計算效率。

表1 拖曳錨模型參數Tab.1 Parameters for anchor in the model

表2 土體模型參數Tab.2 Parameters for soil in the model

表3 方程(3)所涉及參數Tab.3 Parameters in Eq.3

3.2均質土工況對比

圖5 均質土工況對比Fig.5 Comparison of results in uniform clay

系纜點處拖曳力在工程上同樣值得關注。圖5(d)表明，本文求得的拖曳力與Zhao和Liu的計算結果在整體趨勢上一致，但本文計算結果較Zhao和Liu的小。為了說明兩者的差異，將采用理論解(即式(3))與兩種有限元解進行對比，Neubecker和Randolph[10]及Zhang等[11]已分別通過試驗驗證了式(3)的有效性。式(3)中參數(θah,za)分別取自Zhao和Liu及本文的有限元解，得到兩種有限元解對應的理論解，如圖5(d)所示。由于本文求得的za小于Zhao和Liu的(見圖5(c))，使得本文對應式(3)的解小于Zhao和Liu的；本文求得的拖曳力與式(3)完全吻合，Zhao和Liu的有限元解較式(3)大，增大了兩種有限元解的差距。

Zhao和Liu有限元模型中引入的纜繩長37.8 m，土體網格單元數為875 420，本文模型通過引入纜繩方程減小了土體尺寸，土體網格單元數為575 276，較前者減少了34.3%。所使用的計算機配有2顆Intel E5-2690 v2處理器，32G內存。前者的CPU計算時間為50.2 h，本文模型計算時間較前者減少36.7%，有效提高了模型的計算效率。對于前者，如果引入更長的纜繩，土體網格單元數會相應增加，使得模型計算時間更長，此時本文有限元模型在計算效率方面的優勢將更明顯。

3.3線性土工況對比

本文線性土土體強度su=5+1.5zkPa，其余條件和均質土工況保持一致，計算結果如圖6所示。

圖6 線性土工況對比Fig.6 Comparison of results in the clay with linear strength

圖6(d)表明，本文求得的拖曳力與Zhao和Liu的計算結果在整體趨勢上一致，本文有限元解較小，造成兩者差異的原因與均質土中的一致。對比圖5(d)和圖6(d)可以發現，線性土中拖曳力隨著嵌入深度增加而線性增加，而均質土中拖曳力卻逐漸趨于穩定，這是因為作用在錨上的土抗力和土體強度密切相關。

本文有限元模型的計算時間比Zhao和Liu的減少了41.6%，計算效率較均質土工況提高更多。

4 結 語

1)本文基于大變形有限元分析技術耦合的歐拉-拉格朗日法，引入纜繩方程，建立了錨-纜繩-海床土耦合作用的有限元分析模型，其中，系纜點處纜繩與錨的相互作用由纜繩方程來描述，并通過ABAQUS用戶子程序VUAMP來實現。

2)分析模型模擬了拖曳錨在均質和線性強度黏土中的嵌入安裝過程，研究了錨板運動方向、系纜點處拖曳力和拖曳角及運動軌跡等運動特性。通過與Zhao和Liu的有限元分析方法對比表明：本文分析模型求得的錨板運動方向與Zhao和Liu的相同；拖曳角稍小于Zhao和Liu的計算結果，均勻土和線性土中分別相差2.4°和1.7°；拖曳角越小錨的嵌入深度越淺，因此導致錨的嵌入深度也稍淺于Zhao和Liu的計算結果，均勻土和線性土中分別相差6.5%和3.4%；由于嵌入深度較小，使得本文求得的拖曳力小于Zhao和Liu的，而Zhao和Liu的有限元解大于式(3)的理論解，增大了兩種有限元解的差距。

3)由于不對纜繩建模，本文有限元模型的網格規模比Zhao和Liu建立的模型要小，使得均質土中本文有限元模型的計算時間比Zhao和Liu的減少36.7%，線性土中減少41.6%，有效提高了計算效率，且拖曳纜繩越長本文有限元模型在計算效率方面的優勢將更明顯。

通過引入纜繩方程來描述系纜點處纜繩與錨的相互作用，模型的簡化在提高計算效率的同時使得本文有限元解與Zhao和Liu存在一定差異，通過合理分析產生差異的原因，證明了本文所提分析模型的可行性。

[1] ZHANG W.Penetration mechanism and kinematic behavior of drag anchors[D].Tianjin:Tianjin University,2011.

[2] AUBENY C P,MURFF J D,KIM B M.Prediction of anchor trajectory during drag embedment in soft clay[J].International Journal of Offshore Polar Engineering,2008,18(4):314-319.

[3] MURFF J D,RANDOLPH M F,ELKHATIB S,et al.Vertically loaded plate anchors for deepwater applications[C]//Proceedings of International Symposium on Frontiers in Offshore Geotechnics.2005:31-48.

[4] LIU H X,LIU C L,YANG H T,et al.A novel kinematic model for drag anchors in seabed soils[J].Ocean Engineering,2012,49:33-42.

[5] LIU H X,LIU C L,ZHAO Y B,et al.Reverse catenary equation of the embedded installation line and application to the kinematic model for drag anchors[J].Applied Ocean Research,2013,43:80-87.

[6] DUNNAVANT T W,KWAN C T T.Centrifuge modeling and parametric analyses of drag anchor behavior[C]//Proceedings of the 25th Offshore Technology Conference.1993:OTC7202.

[7] O’NEILL M P,BRANSBY M F,RANDOLPH M F.Drag anchor fluke-soil interaction in clays[J].Canadian Geotechnical Journal,2003,40(1):78-94.

[8] LIU H X,ZHANG W,LIU C L,et al.Movement direction of drag anchors in seabed soils[J].Applied Ocean Research,2012,34:78-95.

[9] REESE L C.A design method for an anchor pile in a mooring system[C]//Proceedings of the 5th Annual Offshore Technology Conference.1973:OTC1745.

[10] NEUBECKER S R,RANDOLPH M F.Profile and frictional capacity of embedded anchor chains[J].Journal of Geotechnical Engineering,1995,121:797-803.

[11] ZHANG W,LIU H X,ZHAO Y B,et al.Interactional properties between drag anchor and installation line[J].Journal of Geotechnical and Geoenvironmental Engineering,2014,140(2):04013018.

[12] ZHAO Y B,LIU H X.Numerical simulation of drag anchor installation by a large deformation finite element technique[C]//Proceedings of the 33rd International Conference on Ocean,Offshore and Arctic Engineering.2014:OMAE2014-23476.

[13] ZHAO Y B,LIU H X.Large deformation finite element analysis of the anchor line embedded in seabed soils[C]//Proceedings of the 32nd International Conference on Ocean,Offshore and Arctic Engineering.2013:OMAE2013-10586.

[14] NOH W F.A time-dependent two-space-dimensional coupled Eulerian-Lagrangian code[J].Methods in Computational Physics,1964,3:117-179.

[15] BENSON D J.Computational methods in Lagrangian and Eulerian hydrocodes[J].Computer Methods in Applied Mechanics and Engineering,1992,99(2-3):235-394.

[16] BENSON D J,OKAZAVA S.Contact in a multi-material Eulerian finite element formulation[J].Computer Methods in Applied Mechanics and Engineering,2004,193(39-41):4277-4298.

[17] LIU H X,ZHAO Y B.Numerical study of the penetration mechanism and kinematic behavior of drag anchors using a coupled Eulerian-Lagrangian approach[J].Geotechnical Engineering,2014,45(4):29-39.

[18] QIU G,HENKE S.Controlled installation of spudcan foundation on loose sand overlying weak clay[J].Marine Structures,2011,24(4):528-550.

[19] QIU G,HENKE S,GRABE J.Application of a coupled Eulerian-Lagrangian approach on geomechanical problems involving large deformations[J].Computers and Geotechnics,2011,38(1):30-39.

[20] PUCKER T,GRABE J.Numerical simulation of the installation process of full displacement piles[J].Computers and Geotechnics,2012,45:93-106.

[21] Dassault Systemes.ABAQUS,Version 6.10 documentation[S].2010.

[22] DNV.Design and installation of fluke anchors in clay[S].2000.

[23] DEGENKAMP G,DUTTA A.Soil resistances to embedded anchor chain in soft clay[J].Journal of Geotechnical Engineering,1989,115(10):1420-1438.

[24] VIVATRAT V,VALENT P J,PONTERIO A A.The influence of chain friction on anchor pile design[C]//Proceedings of the 14th Offshore Technology Conference.1982:OTC4178.

[25] YEN B C,TOFANI G D.Soil resistence to stud link chain[C]//Proceedings of the 16th Offshore Technology Conference.1984:OTC4769.

Large deformation finite element analysis on the kinematic behavior of drag anchors in the seabed

LI Peidong,LIU Haixiao,ZHAO Yanbing

(School of Civil Engineering,Tianjin University,Tianjin 300072,China)

P751

A

10.16483/j.issn.1005-9865.2016.02.008

1005-9865(2016)02-0056-08

2015-01-19

天津市應用基礎與前沿技術研究計劃重點資助項目(14JCZDJC39900)；國家自然科學基金資助項目(51179124)

李培冬(1990-)，男，福建漳平人，碩士研究生，從事海洋工程系泊基礎數值研究。E-mail:zplpd@163.com

劉海笑。E-mail:liuhx@tju.edu.cn