破冰船阻力性能預報方法綜述

王 超，孫文林，康 瑞，王國亮

(哈爾濱工程大學 船舶工程學院，黑龍江 哈爾濱150001)

0 引 言

破冰船的阻力問題作為破冰船設計過程中的關鍵問題之一，對于破冰船的設計有著決定性的影響。雖然國外學者已經對其有了深入的研究，在國內卻尚屬起步階段，為此，本文將就實船試驗、模型試驗、經驗或半經驗公式以及理論預報等幾種破冰船阻力研究的主要方法進行系統地闡述和分析。

1 實船試驗

早期的破冰船阻力研究主要是在實船試驗的基礎上對試驗結果加以分析。時至今日，雖然模型試驗技術已經比較成熟，結果也比較可靠，卻仍然難以完全模擬自然條件下的破冰過程。所以，在破冰船的阻力研究中，實船試驗不可取代。

實船試驗前要進行環境條件的觀測、船舶參數的觀測和船舶狀態的觀測三方面工作。為了評估破冰船的阻力性能，通常選擇推力一定時在某一冰層厚度下的航速來進行試驗。試驗區域一般會選擇冰層厚度盡可能均勻的地點，以保證試驗盡可能接近于船舶在平整冰中的航行狀態，如圖1所示。

圖1 破冰船實船試驗Fig.1 Full scale test of icebreaker

破冰船的實船試驗過程與常規船舶類似，但為了研究海冰對船舶航行性能的影響，需要對海冰的相關參數及有關的氣象數據進行詳細記錄。

實際工作中冰的強度常用彎曲強度來衡量，但由于測量工作耗時耗力，因此多直接運用公式或圖譜由鹽度和溫度等參數得出，如Timco和O＇Brien［1］建議的海冰抗彎強度計算公式為:

式中:S 為海水鹽度;t 為海水溫度;Vb為鹽水體積。

傳統航道和冰脊的檔案通常由鉆探獲得，但隨著遙感衛星的普及和應用，如今的冰情檔案多是通過衛星數據得到，遙感衛星獲得的圖像如圖2所示。

為了實時測得冰層的厚度，試驗時還會用到特殊的電磁探測系統(見圖3和圖4)。該系統會向海水發出一個電磁場。由于冰的導電率相對于海水而言很低，因此可以忽略冰的導電率，通過測量電磁場強度得出設備距離海面的距離。再使用1 臺聲波測高儀或激光測高儀測出設備距離冰面的高度，就可以將兩者測得的距離相減獲得冰層的厚度。

圖2 遙感衛星可視波長圖像中的冰層Fig.2 Ice sheet inmODIS visual length image

圖3 首部的冰層厚度電磁探測裝置Fig.3 Electromagnetic system on the bow

圖4 冰厚測量裝置示意圖Fig.4 Sketch of ice thickness measuring system

實船試驗雖然在冰阻力值的測量和冰層破壞模式的觀測上有著不可取代的作用，但是試驗成本很高。因此，隨著模型試驗技術的發展，大量的試驗工作都已經改為采用模型試驗的方式進行。

2 模型試驗

隨著1955年在列寧格勒建立起世界上的第一座冰池，冰池模型試驗技術開始進入人們的視線。到20 世紀70年代人們開始重視模型試驗，并陸續在各地開始了冰池的建造［2］。

相比于實船試驗，模型試驗有準備時間短、花費低的優勢。隨著冰池試驗技術的提高和試驗設備的改進，目前的冰池試驗結果已經能夠滿足一般科研和設計工作的需要。圖5 為HSVA的冰池結構圖。

圖5 HSVA 冰池結構圖Fig.5 Structure drawings of HSVA＇s ice tank

冰池試驗過程中，由于涉及到海冰與船體間的相互作用，需要對海冰進行模擬并保證模型冰符合相似性準則的要求，其試驗參數的相似性關系如表1所示。

表1 冰池試驗相似性關系表Tab.1 Similarity relationship of ice tank test

早期的冰池試驗大多只考慮到傅汝德相似準則。直至1980年，Timco［3］指出柯西相似準則的滿足與否影響著模型冰的破壞形式、船體的縱搖運動和冰層的破碎與浸沒，所以適當的柯西數對于實驗結果的可靠性有著重要的影響。試驗中海冰彈性模量(E )與強度極限(σ)的比值如果與實尺度數據保持一致，根據柯西相似準則可以得到式(3):

模型試驗對于相似性的嚴格要求使模型冰制備的好壞成為試驗結果準確與否的一項重要影響因素。已有的模型冰可以分為低溫凍結模型冰和非凍結模型冰兩大類［4］。二者的區別就在于是否需要采用降溫過程來進行模型冰的制備。目前各地冰池中得到廣泛采用的是低溫凍結模型冰。這是因為低溫凍結模型冰不僅可以較好地滿足試驗對模型冰強度的要求，使用了氯化鈉溶液或其他溶液凍結而成的模型冰往往可以更好地還原海冰的內部晶體結構，得到的測量結果也更為準確。

最早得到應用的模型冰是鹽水冰，但是鹽水冰的使用受限于抗彎強度而無法滿足大多數冰池對于縮尺比的要求。1979年Timco［5］提出了尿素模型冰的制備方法。作者在對多種摻合物進行研究、比較后提出了尿素的水溶液制成的凍結模型冰并被各地的冰池所廣泛采用。

自然條件下的海冰隨著厚度的增加，晶體成長的方向會由隨機分布變為一致的垂向生長，使海冰下部區域的強度大為下降。實際工作中會利用這一特點，在室溫達到冰點以下后向冰池上空噴水，使液滴在空中凝結形成冰晶。這些冰晶會落到水面上成為晶核，加快模型冰的凍結;促使冰晶的垂向生長，使其強度降低;同時，減小晶體的直徑，防止因晶體直徑過大而造成冰層破壞形式的改變［6］。

由于冰層強度隨溫度升高而降低，制冰時也會在冰層的厚度達到目標厚度前關閉冰池的制冷系統，使整個冰池的溫度升高到冰點附近，從而達到降低冰層強度的目的。

隨著模型冰制作技術的發展，這些年也出現了一些新型的模型冰。如1986年Timco［7］制作出了新型的EG/AD/S 模型冰。EG/AD/S 即乙二醇、脂肪族洗滌劑和白糖。但是這種模型冰因為含有糖分，會導致細菌繁殖，使冰池難于維護。而且EG/AD/S冰制備成本高，因此沒有得到普及。

國內的模型冰研究方面大連理工大學于2000年制成了采用聚丙烯粉、水泥、塑料微珠等成分的DUT-1 型非凍結模型冰，并通過實驗證實該模型冰的物理和力學指標已能較好地滿足實驗的要求［4］。

另外，試驗相似性要求實船與試驗條件下船體與海冰間的摩擦系數相等。為此，船模的表面需要噴涂特制的涂料，并于試驗前測出二者間的滑動摩擦系數。進行摩擦試驗的裝置通常如圖6所示。

試驗過程中還要通過懸臂梁試驗或者平板彎曲試驗來測量模型冰的彈性模量和強度極限。進行懸臂梁實驗時一般在冰面上將冰層切割成如圖7所示的6Hi×2Hi×Hi的懸臂梁，并通過對懸臂梁的自由端加載來測得模型冰的彎曲強度σf［8］。

式中:F 為懸臂梁的斷裂載荷;L 為梁的長度(加載點至根部);B 為梁的寬度(根部處);Hi為冰的厚度。

圖6 摩擦試驗裝置圖Fig.6 Friction testing apparatus

圖7 懸臂梁概念圖Fig.7 Conceptual figure of cantilever beam test

根據梁上分布的位移探測裝置可以得到梁的撓曲線，并推出模型冰的彈性模量(圖8)。為了保證模型冰的質量，每次實驗后都需重新測量模型冰的彈性模量。

冰層的密度也是模型冰的重要性質之一。一般需要選取一塊矩形試樣，通過測量樣本干舷與厚度的比值可以根據式(5)得到冰的密度。

圖8 懸臂梁試驗裝置圖Fig.8 Cantilever beam test apparatus

式中:h0為樣本的干舷(見圖9)，水面以上任一平面距離試樣頂部的距離L0與其距離水面的距離L1作差即為試樣的干舷［9］;h 為模型冰試樣的總高度。

圖9 干舷法概念圖Fig.9 Conceptual figure of the free board method

冰池試驗的設備和數據記錄與常規拖曳水池試驗類似，但是冰池試驗中因為要進行冰層斷裂模式的觀察和海冰運動的研究，從而需要使用放置于多個方位的攝像機對試驗過程進行記錄。

目前的破冰船模型試驗技術雖然已經能夠滿足阻力測量的要求，現有試驗技術卻難以滿足對模型冰斷裂韌性的相似性要求。同時，試驗中需要用到大型的制冷設備，冰層的凍結也會花費大量時間。這不但使冰池試驗的成本相比普通的船模試驗大為提高，也增加了試驗的準備時間。另外，由于制冰過程中的環境條件難以控制，冰層的性質會因為溫度的擾動，以及表面的氣流等因素影響，與試驗的目標值存在微小的偏差。因此，當前的冰池試驗技術仍然有待改進。

3 理論預報

破冰船的理論預報方法分為分析方法和數值方法2 種［10］。其中分析方法采用經驗或半經驗公式對冰阻力數值進行預報，而數值方法則借助計算機程序對冰阻力數值進行求解。

3.1 分析方法

分析方法靠研究船體所受冰載荷與各變量之間的關系，對各個參數之間的關系做出一些基本假設，再以此為基礎通過對于實船試驗或模型試驗數據的回歸分析得到公式中的經驗系數，最終形成完整的阻力計算經驗公式或半經驗公式。分析方法的準確性與樣本數據有很大關系，船型與試驗采用的母型船相近時結果較為準確，船型相差較大時則結果偏差較大。另外早期的經驗公式大多存在不合理的假定［11］，因此多數經驗公式的可靠性較差。

隨著冰阻力研究的深入，后來出現了更為復雜的半經驗公式。這些半經驗公式結合了對于破冰過程的理論分析，并將大量相關的影響因素引入到了公式中。這些公式已經可以較為準確地反映出不同設計參數對于冰阻力的影響趨勢，對破冰船的設計選型具有一定的指導意義，也有助于減少模型試驗的工作量，進而節省設計成本。

分析方法中出現較早的是Jansson 于1956年給出的冰阻力估算的簡化公式。但是該公式并沒有為人們所廣泛接受。1970年Lewis和Edwards［12］在前人工作的基礎上推出了破冰阻力計算公式并得到較好的阻力計算結果。之后于1972年由Enkvist［13］提出的計算方法也得到了較好的應用。這一時期的經驗公式還有Lewis 公式和改進的Edwards 公式等。

這些早期公式的理論模型尚不完善，公式中的影響因素也相對較少。所以，這些公式往往不能正確地反映冰阻力同各項影響因素之間的關系。

1973年，Milano［14］嘗試將理論分析應用于冰阻力的計算上。根據結構理論和全尺度試驗的觀測結果將海冰假定為按照固定的模式彎曲，他的海冰破壞模型為后來的大量半經驗公式開辟了道路。

1979年，Enkvist［11］指出海冰的破碎是包含擠壓、剪切和撓曲等多種破壞的復雜破壞形式，同時為了簡化冰阻力的計算，他將破冰船的阻力劃分為破冰阻力、浸沒阻力和摩擦阻力3個部分。

1989年，Lindqvist［15］結合Enkvist的工作提出了Lindqvist 公式。公式針對破冰阻力基于結構力學模型進行了分析，給出了破冰阻力的計算公式。公式還考慮了首柱處擠壓破壞造成的阻力，并使用經驗系數給出了具體的計算公式。

破冰阻力Rc:

彎曲阻力RB:

式中:hice為冰厚;E 為楊氏模量;ν 為泊松比;σf為彎曲強度;B 為船寬;ρw為海水密度。

浸沒阻力Rs:

式中:ρi為冰密度;htot為冰與雪的厚度和;T 為溫度;B為船寬。

總的冰阻力Ri:

圖10 水線進角示意圖Fig.10 Water entrance angle

該公式的提出不但提高了經驗公式的準確程度，而且公式在對傳統楔形船首和匙型船首的計算結果對比中，很好地反映出了船首形狀對于冰阻力的影響，因此在之后得到了廣泛的應用。

1991年，Daley［16］提出了破冰活動周期的概念。他根據全尺度試驗和模型試驗中對于冰層破壞過程的觀測，將冰層的破壞簡化為擠壓破壞，彎曲破壞和冰層破碎3個不斷循環的過程。

Daley的接觸模型雖然對冰層的破壞過程進行了一定的簡化，但是這種破冰過程中的周期性規律卻與實驗結果基本吻合，而且將連續的破冰過程分割開來，使得冰阻力的計算模型得到了簡化。

1997年，由Riska 等［17］提出了Riska 公式，公式在原有的經驗公式基礎上做了改進，并根據波羅的海的大量實船數據分析得出了公式的經驗系數。實際應用證明，該公式對于波羅的海航行的破冰船能夠得到較為準確的結果。

式中:L，B，T 分別為船長、船寬和吃水;V 為航速;hi為冰厚;φ 為首柱傾角，(°);Lpar和Lbow分別為平行中體長和首部長度。經驗系數如表2所示。

表2 Riska 公式常數表Tab.2 Constants in Riska formula

雖然隨著冰阻力研究的深入，冰阻力的計算模型已經得到了長足的進步，但是經驗公式和半經驗公式普遍將冰阻力和敞水阻力區分為2個獨立的阻力成分而忽略了二者之間的耦合作用［11］。同時，在公式中普遍采用的冰阻力成分劃分方法雖然能幫助人們簡化冰阻力的計算，卻難以在試驗中測得各個成分的準確數值，這使得各阻力成分的經驗系數不夠準確。加之公式包含的參數有限，使得各參數的變化對公式準確性有著顯著的影響［18］。

也有部分公式為了計算的簡便，過分簡化了船體表面和冰層邊緣的幾何形狀和二者的接觸形式。此外，經驗公式雖然能夠較為準確地預報常規的破冰船在平整冰中航行時的阻力大小，對于破冰船經過冰脊或與冰山碰撞等特殊情況下受到的阻力卻無法進行準確的預報，因此適用范圍也較為有限。雖然其中的一部分公式得到了實驗數據的支持并且直到現在還被廣泛運用，卻只能夠用于設計階段的船型選擇，不足以替代模型試驗。

3.2 數值方法

由于以上原因，近幾年來人們建立了許多數值方法以提高冰阻力計算方法的準確性和可靠性。在這些數值模擬方法中應用最多的是有限元方法(FEM)和離散元方法(DEM)。根據這些數值模擬方法，計算機不但可以將破壞模式的影響考慮在內，還能比經驗公式更加明確地反映出船體首部形狀變化等因素所導致的阻力改變。為了發揮計算機強大計算性能的優勢，這些數值模型也比經驗公式所采用的模型復雜的多，因此模型中作的簡化相對較少，結果也比經驗公式更為準確。

這其中有限元方法的應用則稍早于離散元方法，是由M??tt?nen和Hoikkanen 在1990年最先應用于冰載荷模擬上的;近年來，Konuk 等 (2009)，Kolari 等(2009)和Ranta et 等(2010)也相繼對有限元方法在冰載荷模擬上的應用進行了研究。

1992年，Evgin 等最先將離散元方法用于冰載荷的模擬，其后的Hopkins (1997)和Paavilainen et al.(2006)等也對離散元方法進行了進一步的研究。

但是受制于早期計算機的計算能力，當時的數值模型相對于經驗公式而言并沒有顯現出特別的優勢，直到近來隨著人們對于破冰船破冰機理研究的深入和計算機性能的飛速發展，使得具有復雜計算模型的數值方法顯現出了前所未有的優勢，這才令冰阻力計算的數值模擬方法成為了人們研究的熱點。

2001年，Valanto［19-20］提出了一個3D 有限元模型。在該模型中考慮了冰面下水流的流體動力學作用對于冰層破壞的影響。此外Valanto 還指出可以通過歐拉- 拉格朗日公式對冰層破壞產生的碎冰進行建模，這進一步促進了浸沒阻力和摩擦阻力的直接計算方法的發展。2001年Wang［21］提出了一種求得時域解的計算方法，計算模型中通過幾何網格來模擬固定的圓錐與移動的冰層之間的連續接觸過程。2011年Vegard Aksnes［22］提出了一種用于計算平整冰對系泊船只作用的半經驗數值計算方法。該方法將冰的運動分為水線附近的冰層破壞、碎冰翻轉以及水線以下由于冰與船體摩擦產生的運動，船體則被離散為平板。且最終計算得到的結果與模型試驗數據比較后證明了該方法具有較好的準確性。2013年Xiang Tan 等［23］提出了一個描述冰與船體之間相互作用的六自由度模型，從而將船體的升沉、橫搖和縱搖運動加入了模型的考慮范圍。

時至今日，人們已經提出了多種用于計算破冰船阻力的數值方法，上述的眾多方法中以有限元方法較為成熟。通常用于計算平整冰條件下冰阻力的有限元方法會在每一個時間步長內依照船體受力計算、運動方程求解、更新水線位置、接觸判定、受力計算、冰層破壞的循環進行。

不同的有限元程序對于冰阻力問題的求解流程雖然相似，但其受力計算公式和冰層破壞的判定卻有很大的差異。至于離散元方法，其主要應用于碎冰條件下的冰阻力計算，因此在完成初始設置和船體與碎冰間的接觸判定后需要對水域中碎冰的受力進行計算，據此得到各個離散元的運動方向和速度，確定其下一時間步長內的位置。

然而數值方法的出現并沒有完全解決分析方法中所存在的問題。盡管通過計算機的輔助，人們得以采用更為復雜的計算模型，其計算過程仍然涉及到冰層的破碎以及冰與海水間的流固耦合問題。為了避免復雜的問題，目前仍有大量的數值方法采用經驗公式的理論模型來對問題進行簡化，致使這些數值方法難以擺脫傳統分析方法的限制。

除了傳統的有限元方法和離散元方法，隨著近年來數值計算方法的迅速發展也有許多嶄新的數值方法被用于冰阻力的計算。2004年，Munjiza［24］結合了有限元和離散元2 種方法的優點，提出了FEM– DEM 方法。通過FEM – DEM 方法可以運用有限元方法來模擬冰層的碎裂并使用離散元方法來模擬碎冰的堆積過程，這樣就可以結合2 種方法的優勢對破冰過程進行計算，得到更為準確的結果。2009年，Poloj?rvi和Tuhkuri［25］首先將該方法應用于冰脊龍骨的沖撞模擬中。這一嘗試不但發展了FEM-DEM 方法在研究破冰過程方面的應用，而且推進了數值模型在模擬冰脊等特殊條件下的應用。

盡管目前的數值仿真方法獲得的冰阻力計算結果還不盡如人意，但是相比經驗公式而言，數值仿真方法便于對船體和冰層形狀進行細致的描繪并一定程度上計入流體動力學作用對冰阻力的影響。而隨著人們對于冰阻力研究的深入，數值方法也在不斷完善。因此，數值方法的研究有著廣闊的前景。

4 結 語

破冰船的阻力問題是影響破冰船能否滿足其設計破冰能力的關鍵問題之一。國外在冰阻力的研究上已經形成了成熟的理論體系，但是由于國內對破冰船領域的關注較晚，這一領域的研究尚在起步階段。

本文就破冰船阻力研究的幾種主要方法進行了簡要的介紹，對不同方法的發展現狀、自身優勢、現存問題及發展方向進行了簡要的分析與論述，對相關的研究者有一定的借鑒意義。

［1］TIMCO G W，O.BRIEN S.Flexural strength equation for sea ice［J］.Cold Regions Science and Technology，1994，22:285-298.

［2］孟廣琳.國外模型冰槽及冰中海上結構物的物理模擬［J］.海洋環境科學，1986(1):85-95.MENG Guang-lin.Physical simulation of foreignmodel ice tank and offshore structures［J］.Marine Environmental Science，1986(1):85-95.

［3］TIMCO G W.The mechanical properties of saline-dopedand carbamide (urea)-doped model ice［J］.Cold Regions Science and Technology，1980，3(1):45-56.

［4］王永學，李志軍，李廣偉.DUT-1 非凍結合成模型冰物模技術及應用［J］.大連理工大學學報，2001(1):94-99.WANG Yong-xue，LI Zhi-jun，LI Guang-wei.Technique of DUT- 1 non-refrigerated breakable ice materials and its applications ［J］.Journal of Dalian University of Technology，2001(1):94-99.

［5］TIMCO G.The mechanical and morphological properties of doped ice:A search for a better structurally simulated ice for model test basins［C］//POAC 79，Proceedings，1979:719-739.

［6］史慶增，徐繼祖，宋安.冰力模型實驗［J］.冰川凍土，1990(2):117-123.SHI Qing-zeng，XU Ji-zu，SONG An.The model test of ice forces［J］.Journal of Glaciology and Geocryology，1990(2):117-123.

［7］TIMCO G W.EG/AD/S:A new type of model ice for refrigerated towing tanks［J］.Cold Regions Science and Technology，1986，12(2):175-195.

［8］SCHWARZ J，FREDERKING R，GAVRILLO V，et al.Standardized testing methods for measuring mechanical properties of ice ［J］.Cold Regions Science and Technology，1981，4(3):245-253.

［9］RISKA K，KATO K M，BELJASHOV M V，et al.Report of the performance in ice-covered waters committee 21st［R］.1996.

［10］TIMCO G，CROASDALE K，WRIGHT B.An overview of first-year sea ice ridges［R］.PERD/CHC report，2000.

［11］ENKVIST E，VARSTA P，RISKA K.The ship-ice interaction ［C］//POAC 79， Proceedings of the 5th International Conference on Port and Ocean Engineering under Arctic Conditions.1980，2.

［12］LEWIS J W，EDWARDS Jr R Y.Methods for predicting icebreaking and ice resistance characteristics of icebreakers［R］.1970.

［13］ENKVIST E.On the ice resistance encountered by ships operating in the continuous mode of icebreaking［R］.1972.

［14］MILANO V R.Ship resistance to continuous motion in ice［D］.Stevens Institute of Technology，1972.

［15］LINDQVIST G.A straightforward method for calculation of ice resistance of ships［C］//Proceedings of POAC.1989:722-735.

［16］DALEYC G.Ice edge contact- a brittle failure process model［R］.1991.

［17］RISKA K.Performance of merchant vessels in ice in the Baltic［M］.Sj?fartsverket，1997.

［18］THORSEN I B.Estimation and computation of ice-resistance for ship hulls［D］.Norwegian University of Science and Technology，2012.

［19］VALANTO P，JONES S J，ENKVIST E，et al.The resistance of ships in level ice.Discussion.Author ＇s closure［J］.Transactions-Society of Naval Architects and Marine Engineers，2001，109:53-83.

［20］VALANTO P.On the cause and distribution of resistance forces on ship hulls moving in level ice［C］//Proceedings of the 18th International Conference on Port and Ocean Engineering Under Arctic Conditions，Ottawa，Ontario，Canada.2001:803-813.

［21］WANG S.A dynamic model for breaking pattern of level ice by conical structures ［J］.Acta Polytechnica Scandinavica Mechanical Engineering Me，2001，156.

［22］AKSNES V.A panel method for modelling level ice actions onmoored ships.Part 1:Local ice force formulation［J］.Cold Regions Science and Technology，2011，65(2):128-136.

［23］TAN X，SU B，RISKA K，et al.A six-degrees-of-freedom numerical model for level ice-ship interaction［J］.Cold Regions Science and Technology，2013，92:1-16.

［24］XIANG J，MUNJIZA A，LATHAM J P，et al.On the validation of DEM and FEM/DEM models in 2D and 3D［J］.Engineering Computations，2009，26(6):673-687.

［25］POLOJARVI A，TUHKURI J.3D discrete numerical modelling of ridge keel punch through tests［J］.Cold Regions Science and Technology，2009，56(1):18-29.