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

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

0 引 言

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

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

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

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

基于Rhino3D平臺參數(shù)化設(shè)計工具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ū)”腳本，該腳本共包含五個變量：浮冰區(qū)長度、浮冰區(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ī)則多邊形，在研究浮冰尺度分布之前，必須選擇一個能表征浮冰尺度的特征量，Yulmetov等[11]定義浮冰等效直徑（Mean Calliper Diameter，簡稱MCD）為：MCD=l/π（l為浮冰表面周長，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)生屬于泊松過程，是一種統(tǒng)計與概率學(xué)里常見的離散概率分布類型，因此浮冰MCD概率分布的擬合曲線也符合泊松分布規(guī)律，浮冰MCD概率分布直方圖及擬合曲線如圖2所示。

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

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

浮冰MCD概率分布優(yōu)化借助于Grasshopper中的Galapagos運(yùn)算器，遺傳算法是計算數(shù)學(xué)中常用的搜索算法，采用該算法建立“遺傳算法優(yōu)化-Voronoi-浮冰區(qū)”腳本，優(yōu)化流程主要概括如下：在已知浮冰概率分布函數(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)相對位置的改變影響到各Voronoi圖的相對尺度，將基于目標(biāo)概率分布函數(shù)獲得的各浮冰面積與優(yōu)化后的各Voronoi圖面積作“差值”計算，當(dāng)差值之和取最小極限時，則確定為全局最優(yōu)解。

基于測量數(shù)據(jù)[11]（Okhotsk海域，2000年2月，Dmin=7 m，β=1.8），對浮冰密集度為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)化前后對比 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概率分布對浮冰阻力的影響

2.1 數(shù)值方法

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

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

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

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

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

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

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

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

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

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

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

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

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

通過Grasshopper“遺傳算法優(yōu)化-Voronoi-浮冰區(qū)”腳本得到二維浮冰模型后，建立厚度為1 m的浮冰區(qū)有限元模型，浮冰密集度為50%，設(shè)置優(yōu)化前和優(yōu)化后工況：A1~A5、B1~B5，共計10個工況。優(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á)到失效截斷壓力時，單元失去承載應(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，模擬時間長度為25 s。船體-浮冰接觸定義為*CONTACT_ERODING_SURFACE_TO_SURFACE，浮冰自身接觸定義為*CONTACT_ERODING_SINGLE_SURFACE；定義流固耦合關(guān)鍵字為*CONTROL_ALE、*CONSTRAINED_LAGRANGE_IN_SOLID。在浮冰區(qū)前設(shè)置了長75 m的緩沖區(qū)以消除波浪對浮冰姿態(tài)的影響，浮冰區(qū)其余三個方向距離水域邊界均為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é)果來看，船舶在浮冰區(qū)航行時，浮冰均有不同程度破碎，破碎主要發(fā)生在船艏附近區(qū)域。浮冰在船艏有明顯的堆積現(xiàn)象，浮冰發(fā)生破碎的同時伴隨著翻轉(zhuǎn)，隨后貼合船體表面下沉并沿船體滑動，最終到船體艉部被尾流場清除。總體而言，大尺度浮冰的破碎更為劇烈，但翻轉(zhuǎn)和堆積現(xiàn)象相對小尺度浮冰則較弱。

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

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

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

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

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

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

圖9浮冰阻力統(tǒng)計值與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ù)的減小趨勢，減小趨勢在6 m

將數(shù)值模擬結(jié)果與Vance浮冰阻力經(jīng)驗(yàn)公式[13]計算值對比可以發(fā)現(xiàn)，當(dāng)浮冰MCD平均值>5 m左右時，優(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é)果明顯偏高，若不對原本服從正態(tài)分布的MCD概率分布向指數(shù)分布進(jìn)行優(yōu)化，將會得到過于保守的計算結(jié)果。

3 結(jié) 論

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

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

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

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

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