系泊系統耦合動力分析不同解法的比較研究

丁佐鵬，馬 山，段文洋，劉 昊

(哈爾濱工程大學船舶工程學院，黑龍江哈爾濱150001)

深海浮式平臺總體分析中，系泊纜索和立管的水動力阻尼效應對平臺低頻運動有重要影響，同時浮體搖蕩運動對系泊線誘導動態張力作用，需采用動力耦合分析方法同時求解浮體和系泊纜索的動力響應［1－2］.由于系泊浮體動力耦合分析涉及的學科領域較多，力學機理和分析技術復雜，尚有諸多問題需進一步探索和完善，其中耦合分析效率和算法的穩定性研究便是一個亟待解決的關鍵問題.為探索在浮體和系泊/立管耦合動力分析中提高計算效率，文中采用“弱耦合”模型，基于異步耦合動力分析［3］的求解思想，開發了深海浮體與系泊、立管的動力耦合程序，在單點系泊FPSO，SPAR平臺動力響應模擬上，其結果同商業軟件Aqwa獲得了很好驗證［4－5］.同時在數值模擬中發現針對浮體與系泊/立管動力耦合分析過程，通常柔性構件的穩定求解需要較小的時間步長，導致泊纜索/立管動力響應分析相對于浮體動力求解耗時最多.但是不能為了節約計算時間而隨意增大柔性桿件求解步長，這樣容易影響算法穩定性，因此有必要對系泊線/立管動力響應分析方法進行深入研究以提高耦合計算效率.

細長桿理論用于系泊線/立管動力響應分析，最早由 Garret［6］提出，然后 Webster等［7］對其改進，考慮桿件小拉伸變形效應.基于細長桿理論進行纜索動力分析時常用的有Adams-Moulton單步積分法［8］和 Newmark-β 迭代法［9］.文中針對這兩種方法，基于細長桿理論，編制了有限元數值求解的動力響應分析程序，分析了兩種算法的數值精度，在此基礎上開展了兩種算法計算穩定性和計算效率的對比分析.

1 細長桿理論模型簡介

1.1 纜索控制方程

圖1所示系泊線典型桿單元的運動控制方程可以寫成如下形式:

圖1 纜索分析空間坐標系示意Fig.1 Coordinates for cable analysis

式中:B為抗彎剛度;r(s，t)為纜索中心線位置向量，s為弧長變量，t為時間變量;λ為拉格朗日乘子，λ=T-Bκ2，κ 為曲率，T(s，t)=F·r′為軸向張力;ε為應變.借助Galerkin法對纜索控制方程中的高階偏導數降階，同時利用插值函數將沿桿單元的變量用節點變量近似表示，可得如下所示的控制方程:

式中:M為質量矩陣;U為節點位移坐標;q為外載荷;γikm，αikm，βikm，μim，τm，ηlm為積分常數;i=1，2，3，4;k=1，2，3，4;m=1，2，3;n=1，2，3.

1.2 水動力載荷和海底邊界條件的模擬

作用在纜索單元上的力主要有附加質量力qI，拖曳力qD和Froude-Krylov力qF-K.

上式中的前兩項可通過莫里森公式求得，在知道了水質點的加速度之后，Froude-Krylov力也可求得.

海底剛度對桿單元的剛度矩陣有影響，同時海底支持力也應該作為單元外載荷的一部分來考慮.文中通過在海底布置彈簧來模擬海底彈性基礎的作用，海底支持力為:

式中:c為海底彈簧剛度系數，文中取c=w/d2，w為纜索單位長度浸水重量，d為纜索直徑;r2為纜索單元垂向坐標，h為水深，e2為垂向單位向量.

2 兩種不同的動力學算法分析

2.1 Adams-Moulton 積分法

將式(2)中的二階微分方程用關于變量的兩個一階微分方程表示:

對上式在時間步內進行積分，并利用梯形公式進行積分運算:

通過聯立求解方程(7)和拉伸控制方程，即可求得節點位移和張力.由數值分析的知識可知，梯形公式的截斷誤差滿足以下關系:

由式(8)可知，Adams-Moulton法的截斷誤差為O((Δt)3)，具有二階精度，可以滿足工程計算需要.

2.2 Newmark-β 法

文中采用平均加速度方法，假設加速度在每個時間步內為常數，其值為步長始末兩點的平均值，在步長內通過積分獲得速度和位移的表達式［10］.在初始時刻，運動控制方程如下:

第k步的位移和速度表示為:

M(k)和q(k)的值可以通過u(k)，·u(k)和u¨(k)的值計算得到.

將式(10)代入式(9)，經整理可將第k步的增量方程變成如下形式:

將其與拉伸控制方程聯立求解可獲得第k步內的位移增量δU(knk)和拉格朗日乘子增量δλ(mk).取γ=1/2，β=1/4，并利用加速度為常數的假設可將式(10)中位移寫成如下形式:

可見，平均加速度假設下的位移具有二階泰勒展開形式，即具有二階精度.平均加速度法得到的加速度計算結果在節點處不連續，但是它是無條件穩定的.文中使用其迭代形式，即將前一步的結果作為初值代入控制方程(11)進行迭代求解直到δu(k)和δλ(k)足夠小，停止迭代，進行下一時間步的計算.初步比較 Adams-Moulton積分法和Newmark-β法，可知在進行纜索動力響應步進求解時，兩種方法具有相同的計算精度.下面將兩種方法分別應用到實際海洋工程系泊系統的動力響應求解中，對兩者結果進行比較.

2.3 兩種算法計算精度對比

本算例中分別利用Adams-Moulton法和Newmark-β迭代法研究了頂端受迫振動情況下系泊纜索頂端張力的時歷變化，并將結果與商用軟件Orcaflex［10］進行對比.系泊纜索為鋼鏈，工作水深223.5 m，鏈長762 m，名義直徑0.14 m，抗拉剛度2.21×109N，單位長度質量235.2 kg/m，法向附加質量系取2.8，法向拖曳力系數取3.2，頂端預張力1558.8 kN.系泊線頂端受迫運動形式如下:

式中:x0為纜索頂端初始橫坐標值，t為時間變量，G(t)為緩載函數.文中為如下形式:G(t)=3α2-2α3，α=t/Tload，Tload為緩載區間，a為受迫運動幅值(m)，ω為受迫運動頻率(rad/s)，T為受迫運動周期(s).圖2，3給出了文中采用的兩種算法計算得到的纜索頂端動態張力結果與商業軟件Orcaflex結果的比較.

圖2 振幅0.9144 m、周期4 s下頂端張力對比Fig.2 Comparison of top tension(a=0.9144 m，T=4 s)

圖3 振幅2.75 m、周期8 s下頂端張力驗證Fig.3 Comparison of top tensions(a=2.75 m，T=8 s)

從圖2，3中可以看出對不同振幅和頻率下的受迫振動Adams-Moulton法和Newmark-β迭代法所進行的系泊纜索張力響應分析結果與Orcaflex的結果具有很好的一致性，從而證明在實際的動力響應求解過程中，這兩種方法也具有相似的計算精度.

2.4 動力學算法穩定性和計算效率比較研究

在深海浮體動力耦合響應分析時，動力分析時間步長是一個關鍵參數，其對動力分析計算效率和穩定性有顯著影響.下面通過改變時間步長研究兩者收斂性和計算效率的比較規律.

選取某系泊線，其材料和工作水深見表1，EA為抗拉(壓)剛度，錨鏈和鋼纜的附加質量系數取1.0，拖曳力系數取2.0，按照式(13)分別選取了小幅高頻和大幅低頻兩種典型的受迫運動形式.

圖4，5分別給出了在幅值為1 m，振動周期為5 s以及幅值為4 m，振動周期為20 s兩種簡諧受迫運動情況下，兩種動力分析算法獲得的纜索頂端張力隨時間步長變化的結果.

表1 系泊線參數Table 1 Paramaters of the mooring line

圖4 小幅高頻運動時收斂性比較(振幅1 m，周期5 s)Fig.4 Convergence analysis of the two algorithms under different time step(a=1 m，T=5 s)

圖5 大幅低頻運動時收斂性比較(振幅4 m，周期20 s)Fig.5 Convergence analysis of the two algorithms under different time step(a=4 m，T=20 s)

可以看出，在保證計算結果收斂的前提下，改變動力分析的時間步長對結果的影響很小，說明這兩種算法都具有很好的收斂性.另外，比較兩個算例最大時間步長可以發現，Newmark-β法適用的步長范圍明顯比Adams-Moulton法大.通過分析可知，由于每次迭代過程相當于將結果重新代入平衡方程對加速度的結果進行誤差修正，減小了誤差累積，穩定性因而更高.以此為基礎，進一步研究了在不同的時間步長下系泊系統動力響應的數值求解效率.

文中算例通過模擬相同時間的纜索動力響應過程，對兩種算法計算效率進行比較(圖6，7).

圖6 不同時間步長下兩種算法計算效率比較Fig.6 Comparison of computational efficiency between the algorithms under different time steps

圖7 不同時間步長下兩種算法計算效率比較Fig.7 Comparison of computational efficiency between the algorithms under different time steps

從圖6，7中可看出隨著時間步長的增大，Adams-Moulton法的數值計算時間顯著降低，而Newmark-β法則未能表現出與Adams-Moulton法相似的規律，甚至在一些大步長下計算效率反而下降了.通過仔細研究發現，迭代過程其實對Newmarkβ法的計算時間影響很大.尤其是在振動初始階段，由于結構的響應比較大，為了滿足精度要求往往在一個時間步需要迭代幾十次，即使到了穩定計算階段，每步也需要數次的迭代，從而影響了整體的計算效率.

3 結論

文中利用細長桿理論建立了系泊纜索的有限元模型，分別根據Adams-Moulton法和Newmark-β法這兩種動力學算法編寫了相應的纜索動力學求解程序.通過與商業軟件的結果進行對比，表明在考慮了實際環境下的各種非線性因素(系泊線幾何參數，水動力載荷和海底接觸條件)之后，這兩種算法依然能夠得到滿足工程應用的結果，具有較好的數值精度.

通過算例表明，Newmark-β法能適應較大時間步長，算法穩定性更好;Adams-Moulton單步法的計算效率隨時間步長的增大提高很快，而增大步長對Newmark-β迭代方法的效率影響不大，步長太大了效率甚至會下降.分析發現，增大Newmark-β迭代法時間步長之后，為獲得穩定的結果每一步所需的迭代次數增多，每一步的計算時間就增加了，從而影響了整體的計算效率.總體來講，Newmark-β迭代方法算法穩定性好，更適合環境載荷瞬息萬變的海洋工程柔性桿件分析，但需要進一步研究如何在保證計算精度的前提下減小每一時間步的迭代次數，以提高求解效率.

References)

［1］ Ormberg H，Larsen K.Coupled analysis of floater motion and mooring dynamics for a turret moored ship［J］.Applied Ocean Research，1998，20:55 －67.

［2］ 竇培林，朱建，陳剛，等.風浪流成長全過程單點系泊FPSO運動低頻響應極值研究［J］.江蘇科技大學學報:自然科學版，2014，28(2):114 －119.Dou Peilin，Zhu Jian，Chen Gang，et al.Low frequency extreme response of single-point moored FPSO during the whole process of growing of wind/wave/current［J］.Journal of Jiangsu University of Science and Technology:Natural Science Edition，2014，28(2):114 －110.(in Chinese)

［3］ Jing X，Webster W C，Xu Q，et al.Coupled dynamic modeling of a moored floating platform with Risers［C］∥Proceedings of the ASME30th Inter:Conference on Ocean，Offshore and Arctic Engineering.OMAE2011 －49553.Rotterdam，Netherlands:ASME，2011.

［4］ Ma Shan，Duan Wenyang.Asynchronous coupling analysis of dynamic response for deepwater moored floating platform［C］∥Proceedings of the Eighth International Workshop on Ship Hydrodynamics.IWSH2013 －P3B －3.Seoul，Korea:September，2013:23 －25.

［5］ 馬山，段文洋.深海移動式平臺與系泊/立管系統動力響應異步耦合分析研究［C］∥第二十五屆全國水動力學研討會暨第十二屆全國水動力學學術會議論文集北京，2013:786－792.

［6］ Garret D L.Dynamic analysis of slender rods［J］∥Journal of Energy Resource Technology，1982，104(4)，302－306.

［7］ Pauling J R，Webster W C.A consistent larger amplitude analysis of the coupled response of a TLP and tendon system ［C］∥5th OMAE Conference，Tokyo，Japan，1986，3:126 －133.

［8］ Ran Z.Coupled dynamic analysis of floating structures in waves and current［D］.Texas:Texas A＆M University，2000.

［9］ Chen Xiaohong.Studies on dynamic integration between deep-water floating structures and their mooring/tendon systems［D］.Texas:Texas A＆M University，2002.

［10］ 李鴻晶，王通，廖旭.關于 Newmark-β法機理的一種解釋［J］.地震工程與工程振動，2011，31(2):55－62.Li Hongjing，Wang Tong，Liao Xu.An interpretation on Newmark beta methods in mechanism of numerical analysis［J］.Journal of Earthquake Engineering and Engineering Vibration，2011，31(2):55 －62.(in Chinese)