胡 昕，詹成勝

(武漢理工大學，湖北武漢430063)

0 引言

隨著極地航線的開辟，對冰區(qū)船舶與海冰的相互作用進行研究變得十分必要。冰排是面積足夠大的漂浮冰體，通常是大面積的冰層，要想通過冰區(qū)水域必須與冰層發(fā)生擠撞［1］。國內(nèi)外的研究熱點主要在冰排與海工結構物如防波堤、閘墩、海洋平臺樁腿的相互作用等［2—3］方面。目前學者們對于船舶與冰山、塊狀冰體或者冰排碰撞［1，4—6］結構響應研究的比較多，直接考慮船冰水三者耦合作用模擬船舶破冰過程的研究很少。因此，研究船舶與冰、水耦合系統(tǒng)作用具有重要的理論意義與實際工程價值，可以為冰區(qū)航行船舶的抗冰沖擊強度分析提供依據(jù)。

由于傳統(tǒng)的有限元方法在涉及流固耦合問題時遇到很多的困難，而FEM-SPH(有限元法-光滑質(zhì)點流體動力學)耦合算法可以作為一種求解復雜流固耦合問題的新思路，它克服了因傳統(tǒng)有限元網(wǎng)格嚴重畸變而引起的數(shù)值計算困難。本文以在層冰中破冰航行的KVLCC2船型為數(shù)值模擬對象，用光滑質(zhì)點流體動力學這種無網(wǎng)格法模擬水，從而可以實現(xiàn)流體和固體系統(tǒng)的動態(tài)耦合分析。用合適的有限元網(wǎng)格模擬層冰，賦予冰適當?shù)牟牧蠈傩裕浖梢阅M層冰的具體破壞形式以及破壞的動態(tài)過程。

1 算法簡介

1.1 SPH 算法

1977年Lucyt和Gingold等分別提出了光滑質(zhì)點流體動力學(Smoothed Particle Hydrodynamics，SPH)方法，并且成功應用于天體物理領域中。近年來該方法被廣泛應用于處理結構動力學問題，主要是用于結構大變形解體、碎裂等分析(如高速碰撞、流固耦合碰撞等)［7］。

SPH基本思想是用一系列粒子表示問題域，用積分表示法來近似場函數(shù)，基于當前支持域分布的粒子，在每一時步通過粒子近似求解場函數(shù)及其導數(shù)，將粒子近似法應用于所有的場函數(shù)相關項中，得到只與時間相關的離散化形式，最后可以進行顯式積分求解。這些使得SPH方法具有無網(wǎng)格、自適應、穩(wěn)定以及拉格朗日性質(zhì)。

函數(shù)核近似式、函數(shù)導數(shù)核近似式:

在粒子i處函數(shù)粒子近似式、函數(shù)導數(shù)粒子近似式:

式中:▽為梯度算子;W(x-x'，h)為光滑核函數(shù)，h為光滑長度，用來定義光滑核函數(shù)的影響區(qū)域;mj、ρj(j=1，2，...，N)為粒子j的質(zhì)量和密度，N 為粒子j支持域內(nèi)的粒子總數(shù)。

利用上面的粒子近似式可以將流體質(zhì)量、動量和能量控制方程離散為下面的近似式:

式中:dρi/dt、dUi/dt、dEi/dt分別為粒子 i相應的隨體導數(shù)。

當粒子i周圍2倍光滑長度h內(nèi)的粒子分布比較均勻時，用粒子近似式計算可以得到相當好的結果［8］。

1.2FEM-SPH耦合算法

FEM-SPH耦合算法實質(zhì)是一種接觸(自我接觸、點面接觸、面面接觸等)算法。當SPH粒子與有限單元節(jié)點之間間距超過接觸厚度(人工變量)時，將會在接觸界面之間產(chǎn)生穿透，此時SPH粒子和有限元節(jié)點之間將會有接觸力作用。任何位于有限元節(jié)點支持域內(nèi)的SPH粒子都會對該節(jié)點產(chǎn)生接觸力，反之任何位于SPH粒子支持域內(nèi)的有限元節(jié)點也都會對該粒子產(chǎn)生接觸力。

界面穿透量為Penetration，接觸力FC大小可用式(8)、式(9)計算:

式中:SLFACM為接觸力懲罰系數(shù);Penetration為穿透厚度;Δt為時間步長;λ為接觸剛度。

圖1為接觸穿透示意圖，其中，Hcont為材料的接觸厚度，為人工變量。

圖1 接觸穿透示意圖

求出接觸力后，將接觸力以外力的形式加入到SPH動量方程和有限元動力學方程中［9］:

式中:M為結構質(zhì)量矩陣;C為結構阻尼矩陣;K為結構剛度矩陣;F為外載荷矩陣;FC為接觸力;ü、U、d分別為某時刻的加速度、速度和位移。

SPH數(shù)值積分過程中采用變光滑長度的方法消除光滑長度改變造成的影響［10］。根據(jù)粒子之間接觸力施加的原理(和光滑長度相關)，當界面之間沒有穿透時，粒子節(jié)點之間就沒有接觸力，所以在計算過程中，不需要定義接觸面以及接觸面的法線方向。該算法對復雜接觸面以及有限元網(wǎng)格變形較大的情況都適用。

2 數(shù)值模擬模型及參數(shù)

2.1 計算模型

本算例選用KVLCC2船型的縮尺(縮尺比為15.4)模型作為計算船體，模擬該船在層冰中連續(xù)破冰的工況。船模參數(shù)見表1。計算域模型和有限元網(wǎng)絡、SPH模型如圖2所示。

兼顧計算效率和計算結果，計算域范圍設定為12 m×5 m×0.896 5 m。

表1 KVLCC2船模主尺度

圖2 計算域模型和有限元網(wǎng)格、SPH模型

圖2 (a)、(b)分別示意了計算域模型和網(wǎng)格模型，其中圖2(a)示意了模型由船體、冰層、水(SPH粒子)、池壁邊界4部分組成，池壁邊界是一個沒有頂面的空心長方體，水可以在池中自由運動，冰層恰在自由液面上且冰層與池壁接觸的單元進行了全約束;圖2(b)展示了遮蔽池壁后模型局部網(wǎng)格模型，船體網(wǎng)格為shell殼單元，數(shù)量為29 930，單元平均大小為0.017 75 mm;冰網(wǎng)格為solid實體單元，單元平均大小為0.017 75 mm，冰厚為38.9 mm，網(wǎng)格數(shù)量為328 192;水為SPH單元，數(shù)量為135 806，粒子的生成可以由六面體網(wǎng)格轉(zhuǎn)化。

2.2 材料模型

冰力模型試驗既具有液流模型試驗的一切特點，同時也是一種材料力學特性的模型試驗，因此根據(jù)冰力模型試驗中的動力相似(重力相似、彈性力相似)可以導出試驗中一些主要物理量的比尺。冰力模型試驗物理量比尺見表2［11］。

表2 冰力模型試驗物理量比尺

根據(jù)以上比尺，可以計算得到縮尺模型所需要的冰材料屬性參數(shù)。模型冰材料屬性見表3。

表3 模型冰材料屬性

冰的破壞隨應變速率的不同表現(xiàn)為韌性破壞和脆性破壞，冰在高應變速率時的破壞表現(xiàn)為脆性破壞。本文計算中冰的應變速率遠大于冰的韌脆破壞分界點的應變速率，因此采用高應變速率時的材料模型［12］。

本文數(shù)值模擬中，船體看作剛體模型，材料采用一般結構用鋼;水體為SPH粒子，材料采用Murnaghan狀態(tài)方程表示:

式中:P0為參考壓力，P0=0;B為體積模量，B=27 300;ρ0為參照密度，ρ0=1 000 kg/m3;c為材料聲速;ρ為流體計算密度;γ為指數(shù)系數(shù)，γ=7。

2.3 接觸模型和模型約束

接觸模型用來模擬計算模型中不同部分之間的相互作用，相互作用的面有主、從之分，主、從面之間不能有初始穿透。從屬面的網(wǎng)格要比主面網(wǎng)格細密，防止主從穿透，因為穿透可能會引起沙漏現(xiàn)象的發(fā)生，后者導致計算出錯。

本文計算中的接觸模型定義有:船冰接觸、船水接觸、冰水接觸、水池壁接觸。上述接觸模型中船冰接觸，定義冰為從屬面，其余定義水為從屬面。接觸力比例因子為人工變量，取為默認值0.1。接觸深度為判斷主從面之間是否接觸的一個閥值，取值小于節(jié)點間距或者粒子間距，為人工變量，取值原則是要保證不能有初始穿透。

SPH單元屬性設置中，核函數(shù)取為默認4階貝塞爾核函數(shù)。軟件中核函數(shù)光滑長度是變化的，變化范圍一般取值為初始光滑長度的0.2～5倍。初始光滑長度和粒子半徑有關。

本文計算中，對池壁單元和冰排最外圍單元進行六自由度約束，對船舶進行垂向約束即不考慮船舶的升沉和縱傾。計算中也不考慮重力場的影響。剛性船舶經(jīng)過0.1 s加速至1.401 m/s，然后勻速破冰前行。

3 數(shù)值模擬結果分析

為了和文獻［13］中的數(shù)據(jù)作對比，除船舶主尺度外，文獻中的主要參數(shù)經(jīng)過換算后與本文計算一致。文獻［13］中的船舶主尺度換算后的尺寸為(長×寬×吃水):5.517 2 m×1.168 m×0.422 m。

經(jīng)過計算輸出狀態(tài)、求解器、程序版本、結束時間設置等，計算時間2 s，船舶受力穩(wěn)定，結果收斂。其中航行方向的冰力時程曲線結果見圖3，0.5 s后冰阻力穩(wěn)定后平均值為180 N，文獻［13］實驗測得的船舶冰阻力換算后為225 N。由于二者船型、船寬和吃水不同，且KVLCC2船型沒有做過冰區(qū)航行破冰試驗，并不能確定本文計算結果的準確性，但是在冰力的量級上經(jīng)過換算符合實測數(shù)據(jù)。本文計算船舶寬度小，計算得出的冰力極值小于實驗值，佐證了結論:接觸寬度的增大極大地影響著極值冰力的大?。?4］。

圖3 冰力時程曲線

從圖3可知，船冰剛開始接觸時冰力最大超過了600 N，隨著冰排的破裂，冰力迅速下降，最后在180 N上下波動，穩(wěn)定后的冰力峰值接近400 N。圖中冰力穩(wěn)定后第1次峰值大約在0.8 s，第2次峰值接近1.4 s，因為這時有新的橫向裂紋產(chǎn)生(見圖4中④、⑧)。隨著裂紋的擴展，處于裂紋線上的單元被刪除，冰力時程曲線呈現(xiàn)規(guī)則波動，一旦裂紋擴展結束，要產(chǎn)生新的裂紋時，冰排在受到擠壓過程中并沒有單元被刪除，此時冰力達到峰值。第3次峰值在1.6 s，由于前面初始裂紋的作用，此時的峰值明顯小于前2次。

圖4中各圖顯示了不同時刻冰排的動態(tài)斷裂過程。冰排斷裂過程中，在本文的船舶尺度和船速下冰排的斷裂主要以彎曲斷裂為主。通過刪除失效的冰單元實現(xiàn)冰排的斷裂，而冰單元的失效是由于裂紋尖端的應力集中，冰排局部Von Mises應力高達1.5×105Pa。

從圖4可以知道:當冰排和船首接觸時，船舶艏肩和艏部船舯最先出現(xiàn)裂紋并進行裂紋擴展，0.583 s船首兩側也出現(xiàn)了2道橫向裂紋并擴展，接著在1.155 s船首部從邊界斜(環(huán))向?qū)挾确较虍a(chǎn)生了2道裂紋把首部冰排分成了2個扇形的浮冰塊，再往后還是產(chǎn)生類似的斜(環(huán))向裂紋把扇形冰塊兒后面的左右2塊冰排切斷。

圖4 海冰斷裂動態(tài)過程

由于海冰斷裂過程中裂紋尖端某些單元被刪除，導致了冰排受力不均，使得裂紋現(xiàn)象并非完全對稱，船舶兩側冰排也不是同時發(fā)生裂紋，呈現(xiàn)了一定的隨機性。

一般來說，質(zhì)量較好、均勻強度較高的海冰破碎周期呈正態(tài)分布，常見的破碎周期為3～4 s［15］。本文破冰周期取為從橫向裂紋形成到扇形浮冰塊兒形成的時間間隔，經(jīng)過換算到實船尺度大約為3.2 s，與統(tǒng)計的最常見的海冰的破碎周期相符。

考慮船舶破冰工況下水的阻力變化這方面的資料較少，且缺乏實驗數(shù)據(jù)，本文計算雖然直接用SPH粒子模擬了水的作用，但是其準確性還有待驗證。圖5給出了計算得出的水的阻力時程曲線，穩(wěn)定后水阻力均值為80 N，根據(jù)經(jīng)驗此值偏大。由于本文對流固耦合計算采用的是接觸碰撞算法，和傳統(tǒng)CFD計算原理區(qū)別很大，所以水動力計算的準確性和專業(yè)CFD軟件相比還有不小的差距。

圖5 水阻力時程曲線

4 結論

本文用SPH方法模擬水對船舶和冰排的作用，驗證了SPH方法用于解決流固耦合問題是可行的。粒子之間的接觸判斷特別費時，因此在計算效率上還需要提升。

本文模擬了模型尺度下的船舶破冰過程，可以為模型試驗提供參考。

文中數(shù)值模擬了海冰彎曲斷裂過程，其冰力時程曲線與文獻［13］中船舶冰力數(shù)據(jù)相符。在本文破冰速率下，冰排的破壞呈現(xiàn)脆性破壞。冰力時程曲線峰值對應于特定形式的裂紋(如本文中的橫向裂紋)產(chǎn)生之前的一個過渡狀態(tài)。

船型會影響流場及船艏興波，進而對冰的破壞產(chǎn)生影響。文中船艏的形式并非破冰型艏，導致冰排向上彎折斷裂。由于KVLCC2船型艏部干舷很小，部分冰塊甚至沖上了甲板，因此合適的破冰型艏有利于破冰航行，更加有利于航行安全。

本文對KVLCC2船型破冰過程進行了數(shù)值模擬，對船舶的抗冰載荷結構設計和船舶冰區(qū)航行安全有一定的參考價值。

［1］ 張健，陳聰，張淼溶，等.船舶與冰排碰撞結構響應研究［J］.船舶工程，2014，36(6):24-26.

［2］ 路衛(wèi)衛(wèi).冰排與建筑物擠壓破壞有限元模擬分析研究［D］.天津:天津大學，2007.

［3］ 單思鏑.河冰對橋墩撞擊作用的數(shù)值模擬分析［D］.哈爾濱:東北林業(yè)大學，2011.

［4］ Daley C，Kim H.Ice collision forces considering structural deformation［C］//ASME 2010 29th International Conference on Ocean，Offshore and Arctic Engineering.American Society of Mechanical Engineers，2010:817-825.

［5］ Liu Zhenhui，Amdahl J，Loset S.Integrated Numerical Analysis of An Iceberg Collision With a Foreship Structure［J］.Marine Structures，2010(1):1-19.

［6］ 何偉.基于ANSYS/LS-DYNA的船舶與冰碰撞分析研究［D］.大連:大連海事大學，2013.

［7］ 張洪濤，趙美英，等.SPH和FEM耦合方法分析機翼前緣鳥撞的響應問題［J］.科學技術與工程，2009，9(7):1802-1806.

［8］ 徐志宏，湯文輝，羅永.光滑粒子模擬方法在超高速碰撞現(xiàn)象中的應用［J］.爆炸與沖擊，2006，26(1):53-58.

［9］ 張志春，強洪夫，高巍然.SPH-FEM接觸算法在沖擊動力學數(shù)值計算中的應用［J］.固體力學學報，2011，32(3):319-324.

［10］ 強洪夫，高巍然.完全變光滑長度SPH法及其實現(xiàn)［J］.計算物理，2008，25(5):569-575.

［11］ 孫金亮.海上結構冰載荷的有限元分析和試驗研究［D］.天津:天津大學，2007.

［12］ 張宿峰.流冰與橋墩的相互作用［D］.哈爾濱:東北林業(yè)大學，2014.

［13］ Su B，Riska K，Moan T.A numerical method for the prediction of ship performance in level ice［J］.Cold Regions Science and Technology，2010，60(3):177-188.

［14］ 韓雷，李鋒，岳前進.冰-錐相互作用破壞全過程的有限元模擬［J］.中國海洋平臺，2007，22(2):22-27.

［15］ 武文華，于佰杰，許寧，等.海冰與錐體抗冰結構動力作用的數(shù)值模擬［J］.工程力學，2008，25(11):192-196.

［16］ 張雄，劉巖.無網(wǎng)格法［M］.北京:清華大學出版社，2004.