系泊系統(tǒng)耦合動力分析不同解法的比較研究

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

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

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

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

1 細(xì)長桿理論模型簡介

1.1 纜索控制方程

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

圖1 纜索分析空間坐標(biāo)系示意Fig.1 Coordinates for cable analysis

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

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

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

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

上式中的前兩項可通過莫里森公式求得，在知道了水質(zhì)點的加速度之后，F(xiàn)roude-Krylov力也可求得.

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

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

2 兩種不同的動力學(xué)算法分析

2.1 Adams-Moulton 積分法

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

對上式在時間步內(nèi)進(jìn)行積分，并利用梯形公式進(jìn)行積分運算:

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

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

2.2 Newmark-β 法

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

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

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

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

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

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

2.3 兩種算法計算精度對比

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

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

圖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-β迭代法所進(jìn)行的系泊纜索張力響應(yīng)分析結(jié)果與Orcaflex的結(jié)果具有很好的一致性，從而證明在實際的動力響應(yīng)求解過程中，這兩種方法也具有相似的計算精度.

2.4 動力學(xué)算法穩(wěn)定性和計算效率比較研究

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

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

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

表1 系泊線參數(shù)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)

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

文中算例通過模擬相同時間的纜索動力響應(yīng)過程，對兩種算法計算效率進(jìn)行比較(圖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法的數(shù)值計算時間顯著降低，而Newmark-β法則未能表現(xiàn)出與Adams-Moulton法相似的規(guī)律，甚至在一些大步長下計算效率反而下降了.通過仔細(xì)研究發(fā)現(xiàn)，迭代過程其實對Newmarkβ法的計算時間影響很大.尤其是在振動初始階段，由于結(jié)構(gòu)的響應(yīng)比較大，為了滿足精度要求往往在一個時間步需要迭代幾十次，即使到了穩(wěn)定計算階段，每步也需要數(shù)次的迭代，從而影響了整體的計算效率.

3 結(jié)論

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

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

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］ 竇培林，朱建，陳剛，等.風(fēng)浪流成長全過程單點系泊FPSO運動低頻響應(yīng)極值研究［J］.江蘇科技大學(xué)學(xué)報:自然科學(xué)版，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］ 馬山，段文洋.深海移動式平臺與系泊/立管系統(tǒng)動力響應(yīng)異步耦合分析研究［C］∥第二十五屆全國水動力學(xué)研討會暨第十二屆全國水動力學(xué)學(xué)術(shù)會議論文集北京，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］ 李鴻晶，王通，廖旭.關(guān)于 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)