航海模擬器中破冰船破冰阻力數(shù)值計(jì)算

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

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

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

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

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

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

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

1 船冰碰撞數(shù)值仿真理論

數(shù)值模擬法作為目前船冰碰撞問題的最主要研究方法,能夠有效地計(jì)算出船舶在發(fā)生碰撞過程中船舶的破壞變形、碰撞力的變化和能量變化等各種響應(yīng)。ANSYS/LS-DYNA是現(xiàn)階段最優(yōu)秀的顯性動力分析的有限元軟件,非常適合求解結(jié)構(gòu)的非線性高速碰撞等動態(tài)沖擊問題。

1.1 有限元分析方法

應(yīng)用有限元法求解船冰碰撞非線性動力學(xué)問題時,系統(tǒng)的一般方程描述為:

(1)

當(dāng)采用中心差分時間積分顯式求解加速度向量時,則計(jì)算結(jié)構(gòu)系統(tǒng)各節(jié)點(diǎn)在第n個時間步結(jié)束時刻tn的加速度矢量可描述為:

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

(2)

式中Fint為內(nèi)力矢量,其計(jì)算為:

(3)

式中:BTσdΩ、Fhg、Fcontact分別是當(dāng)前時刻單元應(yīng)力場等效節(jié)點(diǎn)力、沙漏阻力和接觸力矢量。

在系統(tǒng)中的節(jié)點(diǎn)速度和位移矢量為:

(4)

時間步和時間點(diǎn)定義為:

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

(5)

(6)

下一個模型構(gòu)型為:

xt+Δt=x0+ut+Δt

(7)

1.2 LS_DYNA時間步控制

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

Δt=0.9l/c

(8)

式中c為材料聲速,計(jì)算公式為:

(9)

l為實(shí)體單元特征長度,計(jì)算公式為:

l=V/Amax

(10)

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

2 船冰碰撞數(shù)值模擬

2.1 有限元模型建立

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

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

表1 楔形體與平整冰的主要參數(shù)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

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

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

2.2 模型材料選取

1)破冰船材料選取。

本文采用塑性動態(tài)模型模擬破冰船船體材料,具體的塑性動態(tài)模型材料的參數(shù)設(shè)置如表2所示。

表2 塑性動態(tài)模型材料的參數(shù)設(shè)置Table 2 Parameter setting of plastic dynamic model material

2)平整冰材料選取。

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

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

2.3 船-冰碰撞模型數(shù)值計(jì)算

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

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

圖7 彎曲破壞Fig.7 Bending failure

圖8 擠壓破壞Fig.8 Extrusion failure

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

(11)

(12)

(13)

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

總阻力系數(shù)的表達(dá)式為:

CIt=Cbr+Cc+Cb

(14)

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

圖9 破冰阻力系數(shù)與冰厚傅汝德數(shù)關(guān)系Fig.9 Ice breaking resistance coefficient and Fri

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

圖10 冰阻力Fig.10 Ice resistance

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

圖11 經(jīng)驗(yàn)公式計(jì)算的冰阻力Fig.11 Ice resistance calculated by empirical formula

3 船舶破冰阻力參數(shù)敏感性分析

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

3.1 不同平整冰厚度

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

圖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所示,數(shù)值模擬得到的平均破冰阻力與由經(jīng)驗(yàn)公式(見文獻(xiàn)[1-3])計(jì)算得到的平均阻力值進(jìn)行對比,數(shù)據(jù)結(jié)果如表4所示,與Lindqvist經(jīng)驗(yàn)公式計(jì)算所得平均值相對誤差約為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

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

3.2 不同碰撞速度

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

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

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

3.3 參數(shù)靈敏度

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

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

(15)

3.4 不同碰撞角度

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

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

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

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

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

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

3.5 不同旋回角度

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

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

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

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

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

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

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

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

4 結(jié)論

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

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

3)相比經(jīng)驗(yàn)公式方法,數(shù)值計(jì)算在可視化、船舶初始環(huán)境設(shè)定等問題上有很大優(yōu)勢,今后將利用破冰船實(shí)船試驗(yàn)數(shù)據(jù)進(jìn)行仿真驗(yàn)證,以提高數(shù)學(xué)模型的仿真精度。