基于FEM-SPH 耦合算法的船舶破冰數(shù)值模擬

胡 昕，詹成勝

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

0 引言

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

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

1 算法簡介

1.1 SPH 算法

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

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

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

在粒子i處函數(shù)粒子近似式、函數(shù)導(dǎo)數(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相應(yīng)的隨體導(dǎo)數(shù)。

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

1.2FEM-SPH耦合算法

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

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

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

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

圖1 接觸穿透示意圖

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

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

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

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

2.1 計(jì)算模型

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

兼顧計(jì)算效率和計(jì)算結(jié)果，計(jì)算域范圍設(shè)定為12 m×5 m×0.896 5 m。

表1 KVLCC2船模主尺度

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

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

2.2 材料模型

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

表2 冰力模型試驗(yàn)物理量比尺

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

表3 模型冰材料屬性

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

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

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

2.3 接觸模型和模型約束

接觸模型用來模擬計(jì)算模型中不同部分之間的相互作用，相互作用的面有主、從之分，主、從面之間不能有初始穿透。從屬面的網(wǎng)格要比主面網(wǎng)格細(xì)密，防止主從穿透，因?yàn)榇┩缚赡軙鹕陈┈F(xiàn)象的發(fā)生，后者導(dǎo)致計(jì)算出錯(cuò)。

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

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

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

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

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

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

圖3 冰力時(shí)程曲線

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

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

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

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

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

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

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

圖5 水阻力時(shí)程曲線

4 結(jié)論

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

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

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

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

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

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

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

［3］ 單思鏑.河冰對橋墩撞擊作用的數(shù)值模擬分析［D］.哈爾濱:東北林業(yè)大學(xué)，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］.大連:大連海事大學(xué)，2013.

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

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

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

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

［11］ 孫金亮.海上結(jié)構(gòu)冰載荷的有限元分析和試驗(yàn)研究［D］.天津:天津大學(xué)，2007.

［12］ 張宿峰.流冰與橋墩的相互作用［D］.哈爾濱:東北林業(yè)大學(xué)，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ìn).冰-錐相互作用破壞全過程的有限元模擬［J］.中國海洋平臺，2007，22(2):22-27.

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

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