基于參數(shù)化設(shè)計(jì)的浮冰區(qū)船舶冰阻力研究

童 波，涂勛程，谷家揚(yáng)，陶延武，張忠宇

（1.中國(guó)船舶工業(yè)集團(tuán)公司第七〇八研究所，上海200011；2.江蘇科技大學(xué) 船舶與海洋工程學(xué)院，江蘇 鎮(zhèn)江 212003；3.江蘇科技大學(xué) 海洋裝備研究院，江蘇 鎮(zhèn)江 212003）

0 引 言

隨著北極航道的初步通航，冰區(qū)船舶航行日益增加，但極地區(qū)域分布較廣的浮冰給船舶安全航行帶來(lái)一定威脅。浮冰的幾何形態(tài)及分布規(guī)律具有不確定性，船舶航行于浮冰區(qū)涉及多體碰撞和流固耦合等相關(guān)技術(shù)問(wèn)題，國(guó)內(nèi)外學(xué)者在該領(lǐng)域開(kāi)展了大量研究。

目前，浮冰阻力數(shù)值計(jì)算和試驗(yàn)研究主要關(guān)注浮冰厚度、船冰相對(duì)速度、浮冰密集度等因素對(duì)船舶冰阻力的影響[1-2]。國(guó)內(nèi)外研究人員發(fā)現(xiàn)，基于LS-Dyna軟件計(jì)算得到的船舶浮冰阻力，僅在較低速度時(shí)與試驗(yàn)測(cè)試結(jié)果比較符合，在高速度點(diǎn)時(shí)與實(shí)驗(yàn)值結(jié)果相差較大[3-4]。冰緣區(qū)可分為邊緣區(qū)、過(guò)渡區(qū)和內(nèi)陸區(qū)[5]，邊緣區(qū)到內(nèi)陸區(qū)的浮冰尺度從0.1~100 m不等，具有很強(qiáng)的空間離散性。Jeong等[6]在韓國(guó)KRISO海洋工程研究所冰水池中開(kāi)展了破冰通道寬度和冰塊尺度對(duì)冰阻力的影響研究。van den Berg等[7]采用數(shù)值計(jì)算和試驗(yàn)研究相結(jié)合的方法，對(duì)垂直形狀結(jié)構(gòu)與浮冰場(chǎng)的相互作用進(jìn)行了研究。研究發(fā)現(xiàn)，浮冰形狀對(duì)結(jié)構(gòu)運(yùn)動(dòng)方向上的冰荷載平均值和標(biāo)準(zhǔn)差影響較大。Liu[8]提出了一種采用擴(kuò)張多面體單元離散元模擬浮冰的方法，對(duì)北極浮冰區(qū)域的某浮式結(jié)構(gòu)進(jìn)行了冰載荷計(jì)算，分析了冰厚、冰濃度、浮冰尺度等冰荷載影響參數(shù)，但在計(jì)算過(guò)程中未使用可破碎的冰體材料，與船冰實(shí)際作用情況有所區(qū)別。盡管越來(lái)越多的學(xué)者[9]關(guān)注浮冰形狀效應(yīng)和尺度分布規(guī)律對(duì)冰阻力的影響，但浮冰的生成方法沒(méi)有通用性，無(wú)法滿足當(dāng)前浮冰建模的需求。

本文采用Grasshopper參數(shù)設(shè)計(jì)工具對(duì)浮冰區(qū)進(jìn)行了建模，參照真實(shí)冰區(qū)測(cè)量數(shù)據(jù)，基于遺傳算法對(duì)浮冰尺度概率分布開(kāi)展了優(yōu)化。采用可破碎的各向同性彈性斷裂失效模型模擬浮冰，基于LS-Dyna流固耦合方法對(duì)船舶在不同浮冰尺度范圍的冰區(qū)航行進(jìn)行了數(shù)值研究，探討了浮冰尺度效應(yīng)對(duì)船舶浮冰阻力的影響。

1 浮冰區(qū)參數(shù)化建模及浮冰尺度概率分布優(yōu)化

基于Rhino3D平臺(tái)參數(shù)化設(shè)計(jì)工具Grasshopper生成具有隨機(jī)分布和非規(guī)則形態(tài)的浮冰區(qū)模型，運(yùn)用Grasshopper內(nèi)置的Voronoi運(yùn)算器生成浮冰二維模型。Voronoi圖是關(guān)于空間劃分的基礎(chǔ)數(shù)據(jù)結(jié)構(gòu)[10]，由一組連接相鄰兩點(diǎn)直線的垂直平分線所組成的連續(xù)多邊形。為實(shí)現(xiàn)初始冰區(qū)尺度、浮冰密集度、浮冰數(shù)量等參數(shù)控制和初始點(diǎn)集的隨機(jī)分布，建立相關(guān)參數(shù)和運(yùn)算器后，生成“Voronoi-浮冰區(qū)”腳本，該腳本共包含五個(gè)變量：浮冰區(qū)長(zhǎng)度、浮冰區(qū)寬度、浮冰數(shù)量、浮冰位置和Voronoi圖縮放因子。由該腳本生成的浮冰密集度為50%、浮冰數(shù)量為116塊的二維浮冰區(qū)模型（200 m×100 m）的可視化參數(shù)建模主要流程如圖1所示。

圖1 Grasshopper可視化參數(shù)建模主要流程Fig.1 Visualized parameterized modeling main process of Grasshopper

自然界的浮冰均為非規(guī)則多邊形，在研究浮冰尺度分布之前，必須選擇一個(gè)能表征浮冰尺度的特征量，Yulmetov等[11]定義浮冰等效直徑（Mean Calliper Diameter，簡(jiǎn)稱MCD）為：MCD=l/π（l為浮冰表面周長(zhǎng)，m）。

圖2浮冰MCD概率分布規(guī)律Fig.2 Distribution of probability of brash ice’s MCD

由“Voronoi-浮冰區(qū)”腳本生成的浮冰區(qū)二維模型是基于隨機(jī)點(diǎn)分布構(gòu)成的，隨機(jī)點(diǎn)的產(chǎn)生屬于泊松過(guò)程，是一種統(tǒng)計(jì)與概率學(xué)里常見(jiàn)的離散概率分布類(lèi)型，因此浮冰MCD概率分布的擬合曲線也符合泊松分布規(guī)律，浮冰MCD概率分布直方圖及擬合曲線如圖2所示。

浮冰模型的建立對(duì)于冰阻力預(yù)報(bào)至關(guān)重要?；赩oronoi圖生成的浮冰依賴于隨機(jī)點(diǎn)的生成方式，浮冰MCD概率分布服從泊松分布特性有別于真實(shí)浮冰的分布規(guī)律。自然環(huán)境下浮冰區(qū)MCD概率分布基本符合負(fù)指數(shù)冪函數(shù)，為研究不同的浮冰MCD概率分布對(duì)船舶冰阻力的影響，本文采用遺傳算法對(duì)Voronoi圖進(jìn)行了優(yōu)化（概率分布由正態(tài)分布轉(zhuǎn)化為負(fù)指數(shù)冪函數(shù)分布）。通過(guò)分析極地浮冰尺度的測(cè)量數(shù)據(jù)，浮冰MCD概率分布的負(fù)指數(shù)冪函數(shù)表達(dá)式[11]為：

式中：β為受冰區(qū)地理位置影響的參數(shù)；Dmin為浮冰最小MCD；Dmax為浮冰最大MCD。

浮冰MCD概率分布優(yōu)化借助于Grasshopper中的Galapagos運(yùn)算器，遺傳算法是計(jì)算數(shù)學(xué)中常用的搜索算法，采用該算法建立“遺傳算法優(yōu)化-Voronoi-浮冰區(qū)”腳本，優(yōu)化流程主要概括如下：在已知浮冰概率分布函數(shù)和浮冰總數(shù)的前提下，生成數(shù)量確定同時(shí)具有一定尺度范圍的浮冰，運(yùn)用遺傳算法運(yùn)算器產(chǎn)生2倍于浮冰數(shù)量的“浮點(diǎn)”，將“浮點(diǎn)”累加作用在隨機(jī)點(diǎn)x、y坐標(biāo)值上，隨機(jī)點(diǎn)相對(duì)位置的改變影響到各Voronoi圖的相對(duì)尺度，將基于目標(biāo)概率分布函數(shù)獲得的各浮冰面積與優(yōu)化后的各Voronoi圖面積作“差值”計(jì)算，當(dāng)差值之和取最小極限時(shí)，則確定為全局最優(yōu)解。

基于測(cè)量數(shù)據(jù)[11]（Okhotsk海域，2000年2月，Dmin=7 m，β=1.8），對(duì)浮冰密集度為50%的目標(biāo)浮冰區(qū)（200 m×100 m）進(jìn)行優(yōu)化，優(yōu)化前后的浮冰區(qū)模型如圖3所示，浮冰MCD概率分布直方圖及擬合曲線如圖4所示，可以看出優(yōu)化后的浮冰MCD概率分布（7 m≤MCD≤20 m）與目標(biāo)冪函數(shù)基本重合，“遺傳算法優(yōu)化-Voronoi-浮冰區(qū)”腳本生成的浮冰模型基本符合自然冰區(qū)的浮冰分布方式。

圖3浮冰區(qū)模型優(yōu)化前后對(duì)比 Fig.3 Comparison of brash ice zone models before and after optimization

圖4浮冰MCD概率分布規(guī)律Fig.4 Distribution of probability of brash ice’s MCD

2 浮冰MCD概率分布對(duì)浮冰阻力的影響

2.1 數(shù)值方法

利用LS-Dyna軟件研究“船-浮冰-水”流固耦合問(wèn)題時(shí)，通常采用ALE方法（即Arbitrary Lagrange-Euler）。該方法將Lagrange和Euler方法進(jìn)行結(jié)合，ALE方法突破了固體大變形數(shù)值計(jì)算的難題，目前已經(jīng)成為分析大變形問(wèn)題的重要數(shù)值方法。

任意物理量f在ALE描述下對(duì)時(shí)間的導(dǎo)數(shù)由兩部分組成：

式中：Xi和xi分別為拉格朗日坐標(biāo)和歐拉坐標(biāo)，Δvi為相對(duì)速度，ui和wi分別為流體質(zhì)點(diǎn)速度和參考坐標(biāo)系下的網(wǎng)格速度。通過(guò)上式又可導(dǎo)出ALE描述下的流體連續(xù)性方程及動(dòng)量方程，由質(zhì)量守恒原理用流體密度ρ描述（2）式中的物理量f可得：

（3）式即為流體連續(xù)性方程，對(duì)于不可壓縮流體，可簡(jiǎn)化為：

結(jié)合牛頓第二定律對(duì)流體微元運(yùn)動(dòng)情況進(jìn)行分析，可導(dǎo)出流體單元本構(gòu)方程：

式中：p為流場(chǎng)壓強(qiáng)，μ 為粘性系數(shù)，τij為應(yīng)力張量，δij為 Kronecker δ函數(shù)。

根據(jù)流體單元上壓力、質(zhì)量力、粘性力以及慣性力相平衡的條件可以推出：

將（5）式中的σij代入（6）式中可得ALE描述下的流體動(dòng)量方程：

LS-Dyna采用罰函數(shù)約束的方法來(lái)實(shí)現(xiàn)流固耦合，程序?qū)⒆詣?dòng)追蹤拉格朗日節(jié)點(diǎn)（船體結(jié)構(gòu)和浮冰）和歐拉流體（水和空氣）位置間的相對(duì)位移，檢查每個(gè)從節(jié)點(diǎn)對(duì)主物質(zhì)表面的貫穿。若發(fā)生從節(jié)點(diǎn)對(duì)主物質(zhì)表面的貫穿，界面力fs就會(huì)分布到歐拉流體的節(jié)點(diǎn)上，耦合界面會(huì)在ALE的每個(gè)輸運(yùn)載荷步中進(jìn)行重構(gòu)，對(duì)未出現(xiàn)該情況的則不進(jìn)行任何操作。

界面力大小受貫穿點(diǎn)的相對(duì)位置影響，滿足如下關(guān)系式：

式中：δi為穿透量；ki為接觸剛度因子。

2.2 計(jì)算模型和參數(shù)

計(jì)算模型的水域和空氣域長(zhǎng)寬為280 m×110 m，水深12 m，空氣域高度8 m，建模時(shí)將二者處理為共面，劃分網(wǎng)格時(shí)將兩種材料處理為共節(jié)點(diǎn)，水和空氣均采用Null材料模型，狀態(tài)方程分別采用GRUNISEN和LINEAR_POLYNOMIAL描述，流場(chǎng)外邊界全部采用無(wú)反射邊界條件*BOUNDARY_NON_REFLECTING。流體域?qū)嶓w單元采用變間距網(wǎng)格劃分，網(wǎng)格在z方向上以自由液面處為起始面向兩端漸變式增大，網(wǎng)格數(shù)量為400 400個(gè)。

通過(guò)Grasshopper“遺傳算法優(yōu)化-Voronoi-浮冰區(qū)”腳本得到二維浮冰模型后，建立厚度為1 m的浮冰區(qū)有限元模型，浮冰密集度為50%，設(shè)置優(yōu)化前和優(yōu)化后工況：A1~A5、B1~B5，共計(jì)10個(gè)工況。優(yōu)化后工況除浮冰尺度外的其余變量與優(yōu)化前保持一致，優(yōu)化工況參照冰區(qū)信息如表1所示。

表1優(yōu)化工況參照冰區(qū)信息Tab.1 Referenced ice area information of optimized operating conditions

冰體材料采用*MAT_ISOTROPIC_ELASTIC_FAILURE（MAT_013）。該材料是帶有塑性應(yīng)變失效準(zhǔn)則的各向同性彈性斷裂失效模型，當(dāng)有效塑性應(yīng)變達(dá)到失效應(yīng)變或當(dāng)壓力達(dá)到失效截?cái)鄩毫r(shí)，單元失去承載應(yīng)力的能力，偏應(yīng)力變?yōu)榱?，即材料表現(xiàn)為流體狀態(tài)，冰材料參數(shù)[12]如表2所示。

表2冰體材料參數(shù)Tab.2 The material parameters of ice

船舶主尺度如表3所示，船體材料為剛體，約束船體所有節(jié)點(diǎn)在z方向上的位移，設(shè)定船速16 kns，模擬時(shí)間長(zhǎng)度為25 s。船體-浮冰接觸定義為*CONTACT_ERODING_SURFACE_TO_SURFACE，浮冰自身接觸定義為*CONTACT_ERODING_SINGLE_SURFACE；定義流固耦合關(guān)鍵字為*CONTROL_ALE、*CONSTRAINED_LAGRANGE_IN_SOLID。在浮冰區(qū)前設(shè)置了長(zhǎng)75 m的緩沖區(qū)以消除波浪對(duì)浮冰姿態(tài)的影響，浮冰區(qū)其余三個(gè)方向距離水域邊界均為5 m。以A1、B1工況為例，優(yōu)化前后的船舶-浮冰區(qū)有限元模型（空氣域已隱去）如圖5所示。

表3船舶主尺度Tab.3 Main parameters of ship

圖5船舶-浮冰區(qū)有限元模型（隱去空氣域）Fig.5 Finite element model of ship-brash ice(Air domain is hidden)

2.3 結(jié)果分析

選取優(yōu)化后浮冰尺度最大的B1工況和浮冰尺度最小的B5工況為例，船舶在浮冰區(qū)航行如圖6所示。從模擬結(jié)果來(lái)看，船舶在浮冰區(qū)航行時(shí)，浮冰均有不同程度破碎，破碎主要發(fā)生在船艏附近區(qū)域。浮冰在船艏有明顯的堆積現(xiàn)象，浮冰發(fā)生破碎的同時(shí)伴隨著翻轉(zhuǎn)，隨后貼合船體表面下沉并沿船體滑動(dòng)，最終到船體艉部被尾流場(chǎng)清除。總體而言，大尺度浮冰的破碎更為劇烈，但翻轉(zhuǎn)和堆積現(xiàn)象相對(duì)小尺度浮冰則較弱。

圖6數(shù)值模擬船舶在浮冰區(qū)航行情境Fig.6 Numerical simulation of ship’s navigation in brash ice zone

利用LS-prepost對(duì)計(jì)算結(jié)果進(jìn)行后處理可得到船舶與浮冰作用的受力，提取x方向分量，即為船舶在浮冰區(qū)航行時(shí)的浮冰阻力，船舶進(jìn)入浮冰區(qū)后的浮冰阻力-時(shí)間歷程曲線（10~25 s）如圖7所示，分別對(duì)比A1~A5和B1~B5工況可以看出，船體受浮冰的沖擊頻率隨浮冰尺度的減小均呈明顯增大趨勢(shì)。

圖7浮冰阻力-時(shí)間歷程曲線Fig.7 Brash ice resistance-time history curve

從圖8所示的浮冰MCD范圍在各工況下的分布情況來(lái)看，B1~B5各工況的MCD分布范圍及MCD峰值均大于A1~A5工況，更大尺度的浮冰意味著具有更大的質(zhì)量，這是優(yōu)化后各工況的浮冰阻力峰值整體大于優(yōu)化前的主要原因。

船艏完全進(jìn)入浮冰區(qū)在15 s時(shí)刻，對(duì)15~25 s內(nèi)的阻力值進(jìn)行統(tǒng)計(jì)，統(tǒng)計(jì)值與MCD平均值關(guān)系如圖9所示，從圖9（a）可以看出各工況的阻力峰值在優(yōu)化前后的變化，B3工況與上述規(guī)律有一定出入，其浮冰阻力峰值相對(duì)A3工況更小。原因主要有以下兩點(diǎn)：（1）優(yōu)化浮冰尺度屬于全局優(yōu)化，針對(duì)單個(gè)浮冰的優(yōu)化而言仍具有較強(qiáng)隨機(jī)性，即區(qū)域中的某一塊浮冰被放大或縮小是隨機(jī)的，在該工況的船舶航線上可能恰好有較多的浮冰尺度被縮??；（2）因計(jì)算資源的限制，模擬船舶在浮冰區(qū)航行的時(shí)長(zhǎng)有限，浮冰運(yùn)動(dòng)和姿態(tài)變化在一定周期內(nèi)具有累加效應(yīng)，B3工況阻力峰值可能并未出現(xiàn)在計(jì)算時(shí)間內(nèi)。總體來(lái)看，浮冰阻力最大值的影響因素較多，并未隨MCD平均值變化呈明顯趨勢(shì)。

圖8浮冰MCD范圍分布情況Fig.8 Distribution of brash ice’s MCD range

圖9浮冰阻力統(tǒng)計(jì)值與MCD平均值關(guān)系Fig.9 Relationship between statistical value of brash ice resistance and MCD average value

從圖9（b）可以看出優(yōu)化后的浮冰阻力平均值均小于優(yōu)化前，優(yōu)化前后的浮冰阻力平均值隨MCD增大而呈現(xiàn)負(fù)指數(shù)冪函數(shù)的減小趨勢(shì)，減小趨勢(shì)在6 m

將數(shù)值模擬結(jié)果與Vance浮冰阻力經(jīng)驗(yàn)公式[13]計(jì)算值對(duì)比可以發(fā)現(xiàn)，當(dāng)浮冰MCD平均值>5 m左右時(shí)，優(yōu)化后數(shù)值結(jié)果的浮冰阻力平均值與經(jīng)驗(yàn)值相差較小。需指出的是，Vance提出的基于全尺度試驗(yàn)得出的經(jīng)驗(yàn)公式不含有浮冰尺度及浮冰密集度的變量，但該經(jīng)驗(yàn)值反映出的平均冰阻力在一定程度上也驗(yàn)證了本文數(shù)值模擬結(jié)果的正確性；另外，浮冰在優(yōu)化前的模擬結(jié)果明顯偏高，若不對(duì)原本服從正態(tài)分布的MCD概率分布向指數(shù)分布進(jìn)行優(yōu)化，將會(huì)得到過(guò)于保守的計(jì)算結(jié)果。

3 結(jié) 論

本文采用Voronoi圖對(duì)不規(guī)則幾何形態(tài)的浮冰進(jìn)行了參數(shù)化設(shè)計(jì)，并參照真實(shí)冰區(qū)環(huán)境下浮冰MCD概率分布規(guī)律，基于遺傳算法對(duì)其進(jìn)行優(yōu)化，以獲得接近自然環(huán)境條件下的浮冰區(qū)模型，將浮冰尺度作為研究變量，采用有限元方法對(duì)船舶在優(yōu)化前后的浮冰區(qū)航行開(kāi)展了數(shù)值計(jì)算研究，通過(guò)數(shù)值計(jì)算結(jié)果結(jié)合經(jīng)驗(yàn)值對(duì)比分析得出如下結(jié)論：

（1）大尺度浮冰的破碎更為劇烈，但翻轉(zhuǎn)和堆積現(xiàn)象相對(duì)小尺度浮冰則較弱，船體受浮冰的沖擊頻率隨浮冰尺度的減小呈明顯增大的趨勢(shì)。

（2）優(yōu)化后的浮冰阻力峰值整體大于優(yōu)化前，浮冰阻力峰值受到浮冰形態(tài)不規(guī)則性及姿態(tài)隨機(jī)性等因素的影響，并未隨浮冰MCD增大而呈現(xiàn)出明顯趨勢(shì)。

（3）數(shù)值模擬結(jié)果在較大MCD范圍內(nèi)與經(jīng)驗(yàn)值比較符合，浮冰阻力平均值隨MCD增大而呈現(xiàn)負(fù)指數(shù)冪函數(shù)的減小趨勢(shì)，相比大尺度浮冰而言，小尺度浮冰較高頻率的沖擊載荷對(duì)浮冰阻力平均值起到了明顯的加強(qiáng)作用。

（4）采用直接基于Voronoi圖生成的浮冰區(qū)計(jì)算的浮冰阻力結(jié)果過(guò)于保守，對(duì)浮冰MCD概率分布進(jìn)行優(yōu)化，可以有效降低船舶在小尺度浮冰區(qū)受到的平均冰阻力，提高浮冰阻力預(yù)報(bào)精確性。