基于參數化設計的浮冰區船舶冰阻力研究

童 波，涂勛程，谷家揚，陶延武，張忠宇

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

0 引 言

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

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

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

1 浮冰區參數化建模及浮冰尺度概率分布優化

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

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

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

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

由“Voronoi-浮冰區”腳本生成的浮冰區二維模型是基于隨機點分布構成的，隨機點的產生屬于泊松過程，是一種統計與概率學里常見的離散概率分布類型，因此浮冰MCD概率分布的擬合曲線也符合泊松分布規律，浮冰MCD概率分布直方圖及擬合曲線如圖2所示。

浮冰模型的建立對于冰阻力預報至關重要。基于Voronoi圖生成的浮冰依賴于隨機點的生成方式，浮冰MCD概率分布服從泊松分布特性有別于真實浮冰的分布規律。自然環境下浮冰區MCD概率分布基本符合負指數冪函數，為研究不同的浮冰MCD概率分布對船舶冰阻力的影響，本文采用遺傳算法對Voronoi圖進行了優化（概率分布由正態分布轉化為負指數冪函數分布）。通過分析極地浮冰尺度的測量數據，浮冰MCD概率分布的負指數冪函數表達式[11]為：

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

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

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

圖3浮冰區模型優化前后對比 Fig.3 Comparison of brash ice zone models before and after optimization

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

2 浮冰MCD概率分布對浮冰阻力的影響

2.1 數值方法

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

任意物理量f在ALE描述下對時間的導數由兩部分組成：

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

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

結合牛頓第二定律對流體微元運動情況進行分析，可導出流體單元本構方程：

式中：p為流場壓強，μ 為粘性系數，τij為應力張量，δij為 Kronecker δ函數。

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

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

LS-Dyna采用罰函數約束的方法來實現流固耦合，程序將自動追蹤拉格朗日節點（船體結構和浮冰）和歐拉流體（水和空氣）位置間的相對位移，檢查每個從節點對主物質表面的貫穿。若發生從節點對主物質表面的貫穿，界面力fs就會分布到歐拉流體的節點上，耦合界面會在ALE的每個輸運載荷步中進行重構，對未出現該情況的則不進行任何操作。

界面力大小受貫穿點的相對位置影響，滿足如下關系式：

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

2.2 計算模型和參數

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

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

表1優化工況參照冰區信息Tab.1 Referenced ice area information of optimized operating conditions

冰體材料采用*MAT_ISOTROPIC_ELASTIC_FAILURE（MAT_013）。該材料是帶有塑性應變失效準則的各向同性彈性斷裂失效模型，當有效塑性應變達到失效應變或當壓力達到失效截斷壓力時，單元失去承載應力的能力，偏應力變為零，即材料表現為流體狀態，冰材料參數[12]如表2所示。

表2冰體材料參數Tab.2 The material parameters of ice

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

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

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

2.3 結果分析

選取優化后浮冰尺度最大的B1工況和浮冰尺度最小的B5工況為例，船舶在浮冰區航行如圖6所示。從模擬結果來看，船舶在浮冰區航行時，浮冰均有不同程度破碎，破碎主要發生在船艏附近區域。浮冰在船艏有明顯的堆積現象，浮冰發生破碎的同時伴隨著翻轉，隨后貼合船體表面下沉并沿船體滑動，最終到船體艉部被尾流場清除。總體而言，大尺度浮冰的破碎更為劇烈，但翻轉和堆積現象相對小尺度浮冰則較弱。

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

利用LS-prepost對計算結果進行后處理可得到船舶與浮冰作用的受力，提取x方向分量，即為船舶在浮冰區航行時的浮冰阻力，船舶進入浮冰區后的浮冰阻力-時間歷程曲線（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工況，更大尺度的浮冰意味著具有更大的質量，這是優化后各工況的浮冰阻力峰值整體大于優化前的主要原因。

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

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

圖9浮冰阻力統計值與MCD平均值關系Fig.9 Relationship between statistical value of brash ice resistance and MCD average value

從圖9（b）可以看出優化后的浮冰阻力平均值均小于優化前，優化前后的浮冰阻力平均值隨MCD增大而呈現負指數冪函數的減小趨勢，減小趨勢在6 m

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

3 結 論

本文采用Voronoi圖對不規則幾何形態的浮冰進行了參數化設計，并參照真實冰區環境下浮冰MCD概率分布規律，基于遺傳算法對其進行優化，以獲得接近自然環境條件下的浮冰區模型，將浮冰尺度作為研究變量，采用有限元方法對船舶在優化前后的浮冰區航行開展了數值計算研究，通過數值計算結果結合經驗值對比分析得出如下結論：

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

（2）優化后的浮冰阻力峰值整體大于優化前，浮冰阻力峰值受到浮冰形態不規則性及姿態隨機性等因素的影響，并未隨浮冰MCD增大而呈現出明顯趨勢。

（3）數值模擬結果在較大MCD范圍內與經驗值比較符合，浮冰阻力平均值隨MCD增大而呈現負指數冪函數的減小趨勢，相比大尺度浮冰而言，小尺度浮冰較高頻率的沖擊載荷對浮冰阻力平均值起到了明顯的加強作用。

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