航海模擬器中破冰船破冰阻力數值計算

王曉雪, 張秀鳳, 劉兆春, 孟耀

(1.浙江師范大學 工學院,浙江 金華 321004;2.大連海事大學 航海學院,遼寧 大連 116026)

隨著北極航道的開通,極地船舶航行安全極為重要。由于北極海域環境惡劣,海冰成為影響極地船舶安全航行的最主要威脅。冰區航海模擬器在船舶駕駛員培訓中起到重要作用,其關鍵技術之一是冰區船舶運動數學模型。而冰區船舶運動數學模型的建立需要確定船舶與冰碰撞時的受力。因此,研究船舶與冰碰撞時的受力具有十分重要的意義。

對破冰阻力的預報常采用理論方法、數值方法和試驗方法。試驗方法最簡單直接,但是不足之處是海冰與海洋結構物之間的相互作用極其復雜,實船觀測的數據有限,且這些數據受到船型和海域的影響,會受到實驗條件和測試水平的限制,不能被廣泛應用。而經驗公式與數值模擬的應用范圍更為廣泛,實用性強,可在大量假設條件下以及近似情況的基礎上進行。所以學者們一般都采用理論與經驗公式相結合的方式來研究。

Lindqvist[1]將總阻力劃分為3個部分:破冰阻力、水下阻力和速度相關的阻力,根據試驗數據提出半經驗公式對冰阻力進行預報。Riska方法[2]是在一系列經驗數據基礎上,根據基礎理論相結合推導出來的公式。其中經驗系數是根據波羅的海的不同船型的多組實船試驗得到的。Jeong方法[3]是基于模型試驗數據,對Spencer所建立的近似估算方法進行了驗證并推廣。王超[4]利用經驗方法對連續破冰過程中船舶的冰阻力進行估算,對平整冰中船舶操縱性能進行了研究。

與理論方法和試驗方法相比,數值模擬方法作是目前最主要的研究方法,Luo等[5]考慮了船-水-冰的相互作用,利用CFD-DEM耦合數值方法對冰區航道航行的船舶阻力進行了研究,并與實船試驗數據作了對比。Zhang等[6]通過拖曳水池試驗和計算流體力學軟件,確定了船體斜拖MMG模型中的水動力和力矩,得到了船體4種前進速度下的水動力導數,最后總結了水動力導數與前進速度之間的關系。Su等[7]采用數值方法對船舶在平整冰中的操縱性作了相應預報,并針對在船舶與平整冰的接觸過程中分析了利用船體重力將海冰壓碎的過程。Su等[8]研究了船舶在冰區航行的推力,研究結果表明船舶的推力在冰區中相較于開闊水域中有一定的減額,并對船體將平整冰壓碎的過程進行了受力分析。Raed等[9]研究了船與冰相互作用過程的數值模擬,使用新的封閉式求解方式描述了海冰破裂過程。Wang等[10]采用有限元方法對航行于碎冰區的散貨船進行了數值分析,模擬了4種速度和3種碎冰密度下的冰區船舶受力,分析了航行速度和碎冰濃度對冰載荷的影響。張健等[11]并對封凍在層冰中的船舶受到層冰擠壓的損傷變形進行了分析,比較了不同結構形式的外板對船舶所受冰載荷的影響。劉為民等[12]應用LS-DYNA流固耦合方法對船舶在平整冰條件下航行的冰載荷進行了數值模擬,分析了速度和冰厚對破冰阻力的影響,發現在時域上各個方向的冰載荷都具有一定周期性,且破冰阻力隨航速和冰厚的增加基本呈線性增加的趨勢。

綜上,目前大多學者利用數值模擬方法對船-冰碰撞過程開展研究,主要分析了破冰過程中不同船速、船艏形狀、冰厚等對破冰能力的影響,對船舶不同角度與平整冰的碰撞以及船舶在冰區旋回航行的破冰能力研究的文獻很少。本文基于ANSYS/LS-DYNA軟件,對破冰船與平整冰的碰撞進行數值求解,并研究了船舶與平整冰在不同速度、不同碰撞角度和旋回角度以及不同的平整冰厚度發生碰撞時的動態響應特性,定性分析了不同工況下船冰碰撞的船舶的碰撞力、能量損耗等規律,并且對船速與冰厚2個因素進行了定量分析,擬合出公式,進一步為破冰船在冰區的破冰作業提供一定的理論參考價值。

1 船冰碰撞數值仿真理論

數值模擬法作為目前船冰碰撞問題的最主要研究方法,能夠有效地計算出船舶在發生碰撞過程中船舶的破壞變形、碰撞力的變化和能量變化等各種響應。ANSYS/LS-DYNA是現階段最優秀的顯性動力分析的有限元軟件,非常適合求解結構的非線性高速碰撞等動態沖擊問題。

1.1 有限元分析方法

應用有限元法求解船冰碰撞非線性動力學問題時,系統的一般方程描述為:

(1)

當采用中心差分時間積分顯式求解加速度向量時,則計算結構系統各節點在第n個時間步結束時刻tn的加速度矢量可描述為:

a(tn)=M-1[P(tn)-Fint(tn)]

(2)

式中Fint為內力矢量,其計算為:

(3)

式中:BTσdΩ、Fhg、Fcontact分別是當前時刻單元應力場等效節點力、沙漏阻力和接觸力矢量。

在系統中的節點速度和位移矢量為:

(4)

時間步和時間點定義為:

Δtn-1=(tn-tn-1),Δtn=(tn+1-tn)

(5)

(6)

下一個模型構型為:

xt+Δt=x0+ut+Δt

(7)

1.2 LS_DYNA時間步控制

LS_DYNA軟件顯式有限元方法采用中心差分法離散化時間來求解運動微分方程。LS_DYNA軟件會通過檢查所有單元來自動計算所需時間步長:

Δt=0.9l/c

(8)

式中c為材料聲速,計算公式為:

(9)

l為實體單元特征長度,計算公式為:

l=V/Amax

(10)

式中:Amax是單元的最大面積;V為單元的體積。

2 船冰碰撞數值模擬

2.1 有限元模型建立

“雪龍2號”極地考察船是我國第1艘自主建造的極地科學考察破冰船,是全球第1艘采用船艏、船艉雙向破冰技術的極地考察船。本文以“雪龍2號”主尺度為原型,繪制一艘用于破冰作業的仿“雪龍2號”的船舶模型。只考慮船艏破冰[13],因此對船艉做流線型簡化,如圖1所示。

為了降低對計算機性能的要求,在保證精度的同時,縮減仿真計算運行時間,驗證試驗方法可行性。本文參照“雪龍2號”船艏幾何形狀,建立了簡化的楔形體模型,進行船舶破冰的仿真試驗。利用楔形體簡化模型,通過改變多項影響因子,對船冰碰撞進行多次模擬,得到大量數據。根據數據對各項參數進行靈敏度分析,總結各項參數影響下的冰阻力變化規律。圖2、3為船-冰的有限元模型。楔形體模型與平整冰的主要幾何參數如表1所示。

表1 楔形體與平整冰的主要參數Table 1 Main dimension of wedge and level ice m

圖2 破冰船-平整冰有限元模型Fig.2 Finite element model of icebreaker and level ice

圖3 楔形體與平整冰有限元模型Fig.3 Finite element model of wedge and level ice

在船體破冰模擬中,船體有限元模型采用殼單元,網格采用三角形網格,尺寸長度50 cm(為避免網格畸形,防止在仿真計算過程中產生負體積,定義最小網格尺寸不小于10 cm);平整冰設定為無限冰區,為了提高計算效率,故平整冰采用中心加密方式,即平整冰采用沿船舶中縱剖面方向為中心線,依次向兩側稀疏的方式。其中最小網格尺寸為10 cm×10 cm×20 cm,如圖4所示。

圖4 “雪龍2號”與平整冰的網格劃分Fig.4 The meshing of “Xuelong 2” and level ice

2.2 模型材料選取

1)破冰船材料選取。

本文采用塑性動態模型模擬破冰船船體材料,具體的塑性動態模型材料的參數設置如表2所示。

表2 塑性動態模型材料的參數設置Table 2 Parameter setting of plastic dynamic model material

2)平整冰材料選取。

平整冰材料模型的選擇一直是船-冰碰撞仿真的重點與難點。如表3所示,本文選取各向同性的彈塑性斷裂模型[14]模擬冰材料,失效準則為Von Mises屈服準則,破壞模式為最大塑性應變模式,分離模式為恒定最小壓力模式。

表3 冰的材料模型參數Table 3 Material model parameters of level ice

2.3 船-冰碰撞模型數值計算

破冰過程主要分為破碎、翻轉、滑行。設定船舶初始速度為1.5 m/s,與1 m厚的平整冰發生接觸作用。如圖5、6所示,平整冰首先被船艏撞出一個缺口,隨著破冰船的持續破冰,缺口逐漸增大,直至達到冰層的失效應變或截斷壓力,平整冰破碎斷裂?！把?號”船舶破冰時,海冰的失效主要表現為彎曲破壞,如圖7所示;楔形體破冰時,海冰破損主要為擠壓破環,如圖8所示。

圖6 楔形體破冰應力應變效果Fig.6 Stress-strain effect of wedge ice breaking

圖7 彎曲破壞Fig.7 Bending failure

圖8 擠壓破壞Fig.8 Extrusion failure

為檢驗仿真結果和數據的可靠性,本文采用破冰船Terry Fox號的船模水池試驗[15-16]數據進行對比。由于Terry Fox船模與本文數值模擬所采用船模尺寸不同,需對兩者所得出的結果進行無因次化處理。根據相似原理和量綱分析法,定義阻力的無量綱系數為:

(11)

(12)

(13)

式中:ρice、hice分別是冰的密度和厚度;Δρ是冰和水的密度差;BWL是水線寬;T是吃水;Rb、Rc、Rbr分別是冰浮阻力、冰滑阻力和破冰阻力;Cb、Cc、Cbr分別是冰浮阻力系數、冰滑阻力系數和破冰阻力系數。

總阻力系數的表達式為:

CIt=Cbr+Cc+Cb

(14)

繪制破冰船Terry Fox號的破冰阻力系數與冰厚傅汝德數Fri的關系如圖9所示。由于本文在數值計算中沒有考慮冰清阻力和冰滑阻力,在Fri=0.5時,誤差最大為22%,具有較好的精度,說明了本文數值計算的可靠性。

圖9 破冰阻力系數與冰厚傅汝德數關系Fig.9 Ice breaking resistance coefficient and Fri

如圖10所示,破冰船與平整冰碰撞過程中,平整冰處于“接觸-擠壓-變形-失效-應力釋放-再接觸”狀態。船艏與冰層首先單點接觸,由于船艏接觸點突然受到冰阻礙,造成暫時停頓,船體受到的冰阻力此時產生較大的峰值,時歷第1個峰值后,船舶繼續前行破冰,破冰力逐漸穩定在一個范圍內;與此同時,冰層斷裂過程中,冰層因受到船艏沖擊在接觸處產生微裂紋,船艏對冰層施力時,冰層內部的受力并不均勻,大部分會集中在微裂紋處。冰阻力到達峰值,造成進一步開裂,裂紋擴展是一個正反饋過程,即裂紋越長,應力集中效應越明顯,進一步開裂也越容易。

圖10 冰阻力Fig.10 Ice resistance

對該時間段內的冰載荷進行數據統計得到船體受到的平均阻力為6.358×106N。將數值模擬得到的破冰阻力與Lindqvist[1]經驗公式及Jeong[2]經驗公式進行對比,如圖11所示,2種經驗公式計算得到的冰阻力值隨時間變化的總體趨勢與數值模擬結果相符,Lindqvist經驗公式的阻力平均值4.489×106N及Jeong經驗公式的阻力平均值1.230×106N,與數值模擬得到的平均破冰阻力6.358×106N為同一個數量級,且誤差值在20%以內。故本文對船舶破冰阻力預報數值模擬研究具有一定可靠性。

圖11 經驗公式計算的冰阻力Fig.11 Ice resistance calculated by empirical formula

3 船舶破冰阻力參數敏感性分析

破冰船進行破冰作業時,影響船-冰碰撞的因素有很多,本文選取不同平整冰厚度、不同船舶初速度以及不同船冰碰撞角度、旋回角度等易于量化的參數來分析對船舶破冰阻力的影響。

3.1 不同平整冰厚度

本文選取了平整冰厚度為0.5、0.8、1.0、1.2、1.5 m的5個不同的工況,破冰船以初速度為1.5 m/s向平整冰域航行破冰。破開冰的長度、動能、破冰阻力變化曲線分別如圖12～14所示,破冰船進行破冰作業的同一時間內,隨著平整冰厚度的增加,破冰船破開平整冰的長度漸漸變小,動能迅速下降,隨著破冰船與平整冰的不斷接觸,冰力呈現出鋸齒狀的非線性特征波動,并伴隨有多個峰值出現。

圖12 不同平整冰厚度下破冰長度Fig.12 Ice breaking length under different ice thickness

圖13 不同冰層厚度下動能Fig.13 Ship kinetic energy under different ice thickness

圖14 不同平整冰厚度下碰撞力Fig.14 Force under different ice thickness

如圖15所示,數值模擬得到的平均破冰阻力與由經驗公式(見文獻[1-3])計算得到的平均阻力值進行對比,數據結果如表4所示,與Lindqvist經驗公式計算所得平均值相對誤差約為10.85%,與Riska[3]方法的相對誤差約為28.21%,與Jeong公式的相對誤差約為8.15%。

表4 不同冰厚的平均破冰阻力值Table 4 Average ice breaking resistance under different ice thicknesses 106 N

圖15 不同冰厚的平均破冰阻力Fig.15 Average ice breaking resistance under different ice thickness

由于本文的仿真試驗忽略了破碎浮冰在滑行階段對船體造成的摩擦阻力,僅考慮了海冰環境干擾力,所以數值模擬得到的值比經驗公式計算所得值偏小,但誤差仍在一定范圍內,故本文對船舶破冰阻力數值模擬預報具有一定可靠性。

3.2 不同碰撞速度

破冰船以不同初速度1、1.5、2、2.5 m/s沖撞冰層時,破冰阻力和動能時歷曲線如圖16、 17所示,由于在船-冰碰撞過程中,不斷發生“碰撞-平整冰破碎-碰撞-平整冰破碎”這一現象,碰撞力曲線呈現出了一定程度內的波動。

圖16 不同速度下冰阻力Fig.16 Ice force at different speeds

圖17 不同速度下動能Fig.17 Kinetic energy at different speeds

3.3 參數靈敏度

將上述得到的有關不同冰厚與不同速度得到的冰阻力如表5所示的數據,為利于冰區航海模擬器數學模型的實時解算,采用矩形域的最小二乘法曲面擬合方法,得到適用于“雪龍2號”船舶運動預報的阻力表達式為:

表5 不同工況下船-冰碰撞計算冰載荷Table 5 Ice load in ship ice collision under different working conditions 106 N

(15)

3.4 不同碰撞角度

船舶在無限冰區航行時,不管船艏對著哪個方向都是直航破冰,但在要進入冰區航行時,破冰船船艏與所要接觸的平整冰區所呈角度不同,對碰撞的影響也不同。如圖18所示,破冰船船身與平整冰邊緣成θ角,船速V為2 m/s破冰。

圖18 船舶斜向破冰Fig.18 Oblique ship ice breaking

改變角度θ,取值為60°、65°、70°、75°、80°、85°。如圖19所示,破冰船1/4船身進入平整冰區域的碰撞力時間歷程曲線。船舶沖撞式破冰,撞擊力過大,導致平整冰破損缺口過大,船舶行進時暫時與冰層無接觸,致使有冰力為0的狀態,而且船舶由于船體左右舷兩側受到的冰載荷大小差值過大,碰撞力大幅度波動。

圖19 船舶斜向破冰的碰撞力Fig.19 Force of ship oblique ice breaking

在破冰船以θ角在冰區航行時,將船冰碰撞過程中最大的碰撞力作為破冰船以θ角航行時的碰撞力。如表6所示,船冰碰撞角度為90°時,碰撞力最大,破冰船與平整冰成直角撞擊為最危險工況,由于船舶設計時,一般考慮最危險工況,故本文其他工況船舶與平整冰碰撞的角度均為90°。

表6 不同破冰角度的碰撞力Table 6 Force of different ice breaking angles

3.5 不同旋回角度

如圖20所示,本文設計破冰船船身與平整冰成直角,給予船體前進速度Vy和橫向速度Vx,使破冰船以合速度2 m/s旋回破冰。

圖20 船舶旋回破冰Fig.20 Ship turning ice breaking

但由于破冰船前期以合速度∑V斜向上與平整冰碰撞,如圖20的虛線船身所示,導致船艏右側受到的冰力遠大于左側,致使破冰船向左轉向旋回,直至船舶停止。這一過程不易界定船舶的破冰阻力,故將該過程的平均阻力作為船舶以α角在冰區航行時的破冰阻力。

船舶旋回破冰作業時,改變角度α,取值為5°、10°、15°、20°、25°、30°、10 s內的碰撞力時間歷程曲線如圖21所示,可見碰撞力隨著船舶逐漸進入冰區,船舶與平整冰接觸碰撞的面積逐漸增大而增大。在前期船舶向右前方破冰,碰撞力緩慢增加;在約5～6 s之后,船舶開始在冰區旋回破冰,這個階段的碰撞力時間歷程曲線與圖19的碰撞力狀況相似,運動趨勢一致。

圖21 不同破冰角度的碰撞力Fig.21 Force at different ice breaking angles

船舶在冰區旋回破冰的這7個工況的碰撞力值如表7所示,將船舶旋回破冰的碰撞力與直航連續式破冰的碰撞力對比,發現船舶在冰區旋回破冰比直航連續式破冰受到的碰撞力更大。碰撞力隨著旋回角度的增大,也小幅度的增大。

表7 不同破冰角度的碰撞力Table 7 Force at different ice breaking angles

由于冰區船舶操縱模擬器行為真實感的實時性的要求,依據上述有限元軟件數值計算出的破冰阻力值,進行了適合于航海模擬器實時性要求的處理:1)擬合出適合模擬器用的較高精度的實時估算的阻力公式;2)建立了適用于運動數學模型實時解算的破冰阻力離線數據庫,利用插值方法,實時獲得冰阻力,并在冰區船舶操縱模擬器中得到了應用。

4 結論

1)利用Terry Fox號的船模水池阻力試驗結果與“雪龍2號”船舶的數值計算的無量綱化結果比較,驗證了數值計算的可靠性,并將仿真得到的阻力平均值與3種經驗公式計算的平均冰阻力值進行了驗證對比,進一步說明了數值計算的可信度;

2)利用“雪龍2號”簡化后的楔形體模型,通過LS_DYNA軟件對楔形體模型以不同的碰撞速度、碰撞角度、旋回角度以及平整冰厚度,進行多次船-冰碰撞數值計算,利用矩形域的最小二乘法定量擬合了航速、冰厚的冰阻力計算公式,適用于船舶低速破冰的阻力值計算,本文擬合出的估算公式可應用到航海模擬器中破冰船開辟航道的作業模擬中。

3)相比經驗公式方法,數值計算在可視化、船舶初始環境設定等問題上有很大優勢,今后將利用破冰船實船試驗數據進行仿真驗證,以提高數學模型的仿真精度。