系泊纜索動力分析中Newmark-β迭代解法的計算效率

， ， ，

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

0 引 言

在深海系泊浮體分析中，系泊纜索和立管的動態(tài)響應對浮體有很大影響，需采用動力耦合方法求解浮體與系泊纜索及立管的相互作用。系泊浮體耦合運動包括3個方面：浮體水動力建模、系泊纜索力學分析建模和系泊纜索的耦合算法。按照JACOB等[1]的討論，耦合算法有“強耦合”方法和“弱耦合”方法?！皬婑詈稀狈椒▽⒏◇w和系泊線動力學方程通過彈簧反力(矩)和阻尼力(矩)耦合在一起，形成一個浮體和細長桿件動力學系統(tǒng)同時求解的多自由度方程組，在這種方法中浮體和系泊纜索是整體求解的。“弱耦合”方法將浮體動力學模塊和柔性桿件動力學模塊獨立求解，兩者通過位移和受力的傳遞進行耦合。RAN[2]和GARRET[3]分別開發(fā)了時域全耦合動態(tài)分析程序。KIM等[4]將聚酯纜材料的非線性數(shù)學模型融合進系泊纜索的運動和控制方程，并與深海平臺進行耦合分析。ZHANG等[5]分別采用準靜態(tài)耦合、半耦合和全耦合3種方法對某Spar平臺進行數(shù)值模擬，并將結(jié)果與試驗進行對比。袁夢[6]研究了聚酯纜的動力特性，并對某張緊式Spar平臺進行耦合動力分析。楊敏冬[7]對風、浪、流共同作用下的Truss Spar平臺分別進行了靜態(tài)耦合模擬和動態(tài)模擬?；贘ING等[8]提出的異步耦合方法，馬山等[9-11]開發(fā)了深海浮體與系泊、立管的動力耦合分析程序，在對單點系泊FPSO、Spar平臺的動力響應分析中得到了很好的驗證。該方法屬于“弱耦合”方法，由于系泊系統(tǒng)作為獨立模塊求解，適合采用并行計算。

在浮體動力耦合分析數(shù)值模擬中系泊纜索、立管的動力分析時間步長較小，耗時較長，因此，提高系泊纜索動力分析的計算效率對于減小系泊浮體耦合分析的計算耗時至關重要。JACOB等分別采用“強耦合”和“弱耦合”方法對某半潛式平臺進行動力分析對比，發(fā)現(xiàn)“弱耦合”方法結(jié)合并行計算能大幅提高計算效率。丁佐鵬等[12]基于細長桿理論將Adams-Moulton法和Newmark-β迭代法兩種算法用于系泊纜索動力分析，并對這兩種算法進行對比研究工作，發(fā)現(xiàn)在相同時間步長的情況下Adams-Moulton法計算效率更高，但是其只適用于小步長，而Newmark-β迭代法穩(wěn)定性更好，適用于大步長。

在采用系泊浮體與動力分析異步耦合算法時，針對系泊纜索、立管動力分析常采用細長桿理論。該方法最早由GARRET[13]提出，適用于不可拉伸細長桿件，后經(jīng)MA等[14]改進，考慮了桿件的拉伸變形。唐友剛等[15]基于細長桿理論建立了三維非線性細長桿單元剛度矩陣，實現(xiàn)有限元軟件ABAQUS對纜索單元的調(diào)用，該方法非常適用于工程實際。馬剛[16]結(jié)合細長桿理論分析了某單點系泊FPSO在波浪中的運動和受力。袁夢等[17]對控制方程進行修改，將細長桿理論應用于深水輕質(zhì)纜索的動力分析。細長桿理論動力求解中經(jīng)常采用的一種方法是Newmark-β迭代法[18]，本文采用該方法，以提高系泊纜索動力分析計算效率為目標進行以下工作：(1)研究單元數(shù)、時間步長、控制精度和每個時間步的迭代次數(shù)等要素對計算效率的影響。(2)將OpenMP并行計算應用中系泊纜索的動力分析中，研究其對計算效率的影響。發(fā)現(xiàn)適當減小單元數(shù)、增加時間步長、降低動力分析迭代求解控制精度以及減小每個時間步的迭代次數(shù)有利于提高計算效率，采用并行計算對提高計算效率也是有利的。

圖1 纜索分析空間坐標示例

1 基于細長桿理論的Newmark-β迭代法

1.1 基于細長桿理論的動力學響應控制方程簡介

圖1給出了纜索動力分析的弧長坐標系，系泊纜索的運動及控制方程如下：

r′·r′=(1+ε)2(2)

引入插值函數(shù)，利用Galerkin法離散后的有限元運動及拉伸控制方程如下：

μimqmn+fin(3)

1.2 Newmark-β迭代法求解格式

系泊纜索有限元離散動力學方程組采用Newmark-β迭代法求解。具體求解格式如下：

在t=0時，或者時間步k=0時，有

(5)

假設初始時刻結(jié)構(gòu)靜止，初始條件為

(6)

(7)

時間步k時各變量的值用前一時間步k-1步的值通過式(9)進行預估：

(9)

式中：Δt為時間步長。

M(k)，q(k)可以通過以上各值求得。

式中：ρ為水的密度;A為材料截面積;FEA為纜索抗拉剛度。

求解出以上方程后，通過式(12)修正各變量的值：

(12)

式中：γ、β為常數(shù)，γ=1/2，β=1/4。

1.3 海底邊界條件

海底對系泊纜索同時有支持力和摩擦力的作用，本文忽略摩擦力的作用，只考慮支持力的作用。將海底對纜索支撐力作用看成線彈性彈簧，參考文獻[18]，支持力的公式為

(13)

參數(shù)運動及控制方程，對支持力進行離散積分后得到其表達式：

(14)

對于完全觸底的單元，式(14)中的l1即為單元長度，而μim、γikm仍為積分常數(shù)；但是對于一部分觸底、一部分懸空的單元，式(14)中的l1不是單元長度，μim、γikm也不再是常數(shù)，積分只在觸底的那一部分區(qū)間上進行。

對于部分觸底的單元需先求解出觸地點，即求解式(15)中的零點：

(15)

2 各要素對計算效率的影響

基于以上理論，本文編制了系泊纜索動力求解程序?，F(xiàn)以單一成分系泊線的強迫運動動力分析為例，研究有限單元數(shù)、時間步長、控制精度和每個時間步的迭代次數(shù)對計算效率的影響。系泊線的具體參數(shù)見表1。

表1 單一成分系泊線參數(shù)

圖2 單一成分系泊線強迫運動示例

系泊線作業(yè)水深1 000 m，在系泊線頂端施加一水平強迫正弦運動，如圖2所示，運動位移關系：

x(t)=x0+a·sin(ωt)(16)

式中：x0為系泊線頂端初始時刻坐標；a為施加的強迫運動的幅值；ω=2π/T為振蕩圓頻率，T為振蕩周期。

坐標系原點在靜水面上，取x0=0，a=10 m，T=15 s，系泊線頂端位于水平面下8 m處。

2.1 單元數(shù)和時間步長的影響

現(xiàn)將整個系泊線分別劃分為10、20、40個單元，時間步長分別取0.1 s、0.2 s、0.3 s，共9種工況進行數(shù)值模擬，控制精度取10-2，并將系泊線頂端張力的時歷曲線與商業(yè)軟件OrcaFlex的計算結(jié)果進行比較，如圖3～圖5所示。

圖3 單元數(shù)為10的計算結(jié)果 圖4 單元數(shù)為20的計算結(jié)果

由圖3～圖5可以看出:劃分不同的單元數(shù)，取不同的時間步長，采用Newmark-β迭代法進行數(shù)值計算得出的結(jié)果總體與OrcaFlex模擬結(jié)果是一致的，說明本文采用的數(shù)值方法是可靠的。上述比較也說明：在相同單元數(shù)目情況下，不同時間步長下計算結(jié)果較為一致，本文采用的Newmark-β迭代法數(shù)值穩(wěn)定性較好，在較大的時間步長下纜索動力分析也能夠較為穩(wěn)定地進行下去。

圖5 單元數(shù)為40計算結(jié)果

現(xiàn)按以上各工況進行數(shù)值模擬，數(shù)值模擬時長均3 000 s，得到計算耗時見表2。

表2 各工況模擬3 000 s的計算耗時對比 s

由表2可以發(fā)現(xiàn)：在固定控制精度為10-2，單元數(shù)相同時，隨著時間步長的增加，計算耗時逐漸減小，隨著單元數(shù)的增加，這一趨勢越明顯。在時間步長相同時，單元數(shù)越多，計算耗時越長，而且隨著單元數(shù)的增加，耗費時間的增加速度越快。從這2點可知，單元數(shù)的劃分對模擬的耗時非常敏感，雖然增加單元數(shù)可以提高計算精度，但是單元數(shù)增加到一定程度時，數(shù)值計算的精度已不能提高很多，如在單元數(shù)取20時計算精度已經(jīng)很高，增加到40并沒有帶來太大的精度提升，反而由于分析的單元數(shù)增加了，擴大了整個系泊纜索求解的矩陣，導致計算耗時的增加。因此在采用Newmark-β迭代法進行數(shù)值計算時，在滿足計算精度的要求前提下盡量減少單元數(shù)，對于提高計算效率有很大的幫助。同時可以看到，提高時間步長也減小了計算耗時，這一現(xiàn)象與丁佐鵬等所述情況有所不同，主要是因為此處控制精度取的是10-2，比較容易達到要求，增加時間步長對迭代次數(shù)的影響并不是很大。而丁佐鵬等在進行纜索動力分析計算效率研究時，由于針對系泊纜索動力分析迭代控制精度比較高，不容易達到，增加時間步長導致迭代次數(shù)明顯增加，導致大時間步長反而比小步長更耗時。

2.2 控制精度的影響

將纜索劃分為20個單元，時間步長取0.2 s，控制精度分別取10-2、10-3、10-4、10-5、10-6進行數(shù)值模擬，計算得到纜索頂端張力時歷曲線，如圖6所示。數(shù)值模擬時間取3 000 s，將各控制精度下每個時間步的平均迭代次數(shù)及計算耗時進行對比，見表3。

圖6 各控制精度下頂端張力時歷曲線對比

表3 各控制精度下每個時間步平均迭代次數(shù)及計算耗時

從圖6中可以看出，各控制精度下頂端張力的時歷曲線基本一致，表3顯示每個時間步的平均迭代次數(shù)隨著控制精度的提高而迅速增加，迭代次數(shù)的增加必然導致計算耗時的增加。因此，一味地提高控制精度并沒有太大的意義，反而會增加耗時，得不償失。在計算精度得到保證的前提下，降低控制精度有利于減少每個時間步的迭代次數(shù)，從而減小計算耗時，提高計算效率。

2.3 每個時間步迭代次數(shù)的影響

圖7 各迭代次數(shù)下位移的精度 圖8 各迭代次數(shù)下張力的精度

圖9 每個時間步迭代次數(shù)為3次和20次時 頂端張力時歷曲線(1 000 s)

圖10 每個時間步不同迭代次數(shù)計算耗時

從圖7和圖8可以看出：各迭代次數(shù)下位移的精度遠高于張力的精度，可以推斷控制迭代精度時其實主要控制的是張力的迭代精度。但無論是位移還是張力，隨著每個時間步迭代次數(shù)的增加，計算的精度均在不斷提高。在迭代次數(shù)較少的時候，精度提高的速度很快，但是當?shù)螖?shù)增多之后，精度的提高越來越慢，迭代次數(shù)達到10次之后精度的提高已不很明顯。即使每個時間步的迭代次數(shù)僅為3次，張力的精度相較于張力的幅值已經(jīng)很高了，且模擬時長達到1 000 s計算的結(jié)果如圖9所示，也能保持良好的穩(wěn)定性，并與迭代20次的結(jié)果吻合得很好。

從圖10可以看出：計算耗時隨著每個時間步的迭代次數(shù)的增加，基本呈線性增加趨勢，因此控制每個時間步的迭代次數(shù)是提高計算效率的一種有效方法。

圖11 并行計算

圖12 并行計算耗時

3 并行計算

系泊浮體的系泊系統(tǒng)通常由多根系泊線組成。在對系泊浮體動力耦合分析開展研究的“弱耦合”方法中，在系泊纜索求解時各根系泊纜索動力分析是獨立進行的，允許對系泊纜索并行求解提高計算效率，如圖11所示。本文應用Intel Fortran自帶的并行運算方法OpenMP實現(xiàn)多根系泊纜索并行計算，該方法相對簡單，容易實現(xiàn)。

為研究并行計算的效率，取12根系泊線，簡便起見，其參數(shù)均與第2節(jié)中的一致，且均做與上節(jié)中同樣的強迫運動。將泊線劃分為20個單元，時間步長取0.2 s，模擬時長取500 s。運行用計算機參數(shù)為Intel Core i5-4460，四核。并行計算時分別采用2線程、3線程和4線程進行計算，將并行計算與串行計算(即線程數(shù)為1)的耗時進行對比，如圖12所示。

從圖12可以看出：隨著線程數(shù)的增加并行計算的耗時顯著減小，說明將并行計算運用于系泊纜索的動力分析是十分有效的。

4 結(jié) 語

本文在系泊纜索動力分析細長桿理論的基礎上，建立系泊纜索有限元模型，采用Newmark-β迭代法求解，編制動力分析計算程序，并利用商業(yè)軟件OrcaFlex驗證程序的正確性。以此為前提，本文著重討論了各要素對計算效率的影響。其中系泊線單元數(shù)的劃分對計算耗時的影響比較敏感，隨著單元數(shù)的增加，計算耗時呈非線性增加趨勢，適當增加時間步長有利于減小計算耗時。隨著控制精度的提高，每個時間步的迭代次數(shù)增加速度越來越快，一味提高控制精度對計算結(jié)果精度的提高并沒有顯著的作用，反而造成計算耗時的大幅增加。控制每個時間步的迭代次數(shù)，能有效減小計算耗時，而且計算結(jié)果的精度也并未有太大損失。對于實際應用中深水懸鏈式系泊，系泊纜索的長度常在千米以上，根據(jù)前文算例分析及筆者經(jīng)驗，動力分析的時間步長可取0.2 s，單元數(shù)可取20，每個時間步的迭代次數(shù)可取3，動力分析結(jié)果精度能夠滿足實際工程的要求。另外，針對系泊浮體動力耦合分析的“弱耦合”方法，在對多根系泊線進行動力分析時可采用并行計算。本文針對系泊線動力分析開展的并行計算表明：隨著線程數(shù)的增加，計算耗時呈線性減少趨勢，采用并行計算能夠有效降低系泊纜索動力求解的耗時。下一步將集成系泊纜索與浮體動力耦合分析程序，在此基礎上開展系泊纜索并行計算下系泊浮體動力耦合算法計算效率研究，以便進一步提高計算效率。