叢 濱，慕 鵬，宋星晨，李陳峰

(1. 中國人民解放軍92941 部隊，遼寧 葫蘆島 125000；2. 哈爾濱工程大學(xué)，黑龍江 哈爾濱 150001)

0 引 言

與開闊水域所遭受的風(fēng)浪流載荷作用不同，極地地區(qū)冰情復(fù)雜，存在許多大小形狀不定的浮冰，艦船在極地地區(qū)還要承受冰載荷的碰撞與擠壓作用，對船體結(jié)構(gòu)造成極大的威脅。因此對船舶在極地海域航行過程中可能遭受的載荷成分進(jìn)行分析，對船舶在冰區(qū)航行操縱性與安全性有重要意義。

目前國內(nèi)外學(xué)者開展了大量的船-冰碰撞數(shù)值模擬分析，隨著計算機(jī)技術(shù)的發(fā)展，形成了以網(wǎng)格單元建模的有限元方法(FEM)、無網(wǎng)格化建模的離散元方法(DEM)、光滑粒子流體動力學(xué)方法(SPH)及近場動力學(xué)方法(PD)[1]。蔣昱妍[2]基于粘聚單元法進(jìn)行了層冰-結(jié)構(gòu)物碰撞的數(shù)值仿真，對海冰失效模式和碰撞冰阻力進(jìn)行了研究，但忽略了海水流固耦合效應(yīng)的影響。Gagnon[3]基于流體動力學(xué)原理，建立了油船右舷舷側(cè)與冰山碰撞的數(shù)值模型，討論了碰撞過程中右舷結(jié)構(gòu)上的負(fù)載與壓力分布。Jou[4]對破冰船壓潰式破冰方式進(jìn)行了研究，通過DEM 方法對海冰及其裂縫進(jìn)行建模，預(yù)測了不同形狀破冰船船首與浮冰接觸時可能產(chǎn)生的浮冰破壞模式。Khayyer[5]提出了1 種增強(qiáng)的ISPH 流體模型和新的半隱式ISPH-SPH 耦合方法，并針對典型流固耦合問題進(jìn)行了模擬，做了嚴(yán)格的驗證和性能研究，是SPH 方法的巨大進(jìn)步。薛彥卓[6]將PD方法用于研究海洋結(jié)構(gòu)物與冰的相互作用問題，如基于近場動力學(xué)進(jìn)行了冰的3 點彎曲試驗?zāi)M、船舶與水平層冰相互作用的數(shù)值模擬等。上述文獻(xiàn)對于海冰-海洋結(jié)構(gòu)物相互作用的研究多集中分析船舶所受冰載荷作用，對浮冰的運(yùn)動狀態(tài)進(jìn)行了約束，部分忽略了水域的作用，未對水阻力及水域作用下浮冰-碰撞的運(yùn)動效果進(jìn)行分析。

為了研究船-冰-水三相耦合時的船體載荷響應(yīng)情況，本文根據(jù)船在冰緣區(qū)浮冰漂浮特性和船舶航行特性，在Ls-dyna 中通過S-ALE 流固耦合方法進(jìn)行船-冰-水相互作用機(jī)理研究。在此基礎(chǔ)上，研究不同尺寸浮冰與艦船相互作用的水阻力與冰載荷響應(yīng)規(guī)律。

1 數(shù)值方法原理

1.1 S-ALE 流固耦合方法

船舶與浮冰發(fā)生碰撞，由于體積與質(zhì)量差異，浮冰運(yùn)動姿態(tài)的改變更為明顯，因此船-冰-水的流固耦合模擬效果更多體現(xiàn)在浮冰運(yùn)動動態(tài)的真實性與船體受載的準(zhǔn)確性方面。本文在研究過程中對流體與結(jié)構(gòu)的耦合問題描述，采用結(jié)構(gòu)化-任意拉格朗日-歐拉(structured arbitrary Lagrangian-Euler algorithm，SALE)方法。

S-ALE 方法是在ALE 方法上發(fā)展的更為高階的流固耦合模擬方法。ALE 方法描述的流固耦合是指發(fā)生在流體與固體介質(zhì)交界處的協(xié)調(diào)關(guān)系和平衡，采用拉格朗日變量對固態(tài)結(jié)構(gòu)進(jìn)行運(yùn)動邊界的描述、采用歐拉變量對流體的運(yùn)動特性進(jìn)行描述，可以有效跟蹤物質(zhì)結(jié)構(gòu)邊界運(yùn)動和防止網(wǎng)格產(chǎn)生畸變。S-ALE 方法在原有算法和邏輯的基礎(chǔ)上，進(jìn)行了規(guī)律性簡化和自適應(yīng)網(wǎng)格生成，在計算時間上比ALE 方法減少20%～40%。用光滑delta 近似函數(shù)通過添加“源項”、分布節(jié)點力和插值速度來表示流場和結(jié)構(gòu)物的交互作用，實現(xiàn)結(jié)構(gòu)與流場的作用及反作用模擬，優(yōu)化了ALE 算法由于流體-結(jié)構(gòu)之間接觸定義產(chǎn)生的流體滲漏問題和壓力分布的連續(xù)性問題[7-8]。一般不可壓N-S 方程和力源項的分布函數(shù)為：

式中：ρ為流體密度；u為流場速度；p為流場壓力；t為時間；μ為動態(tài)粘度系數(shù)。F(x,t)為流場受到邊界的所有體積力，其中x為浸入邊界位移；f(s,t)為浸入邊界產(chǎn)生的單位力；X(s,t)為Dirac delta 函數(shù)；δ為光滑化的Dirac 函數(shù)。

1.2 罰函數(shù)接觸算法

船冰碰撞過程中存在大量接觸行為，如浮冰與船舶初始行為階段的碰撞，首柱舷側(cè)兩邊受到冰體的擠壓行為，因此需要通過接觸理論定義求解。船體與浮冰均為彈塑性結(jié)構(gòu)，接觸設(shè)置選擇為面對面接觸，接觸類型選擇為侵蝕接觸，接觸算法為罰函數(shù)法[9]。罰函數(shù)法從判斷物體接觸是否穿透的角度進(jìn)行分析，將要接觸的2 個面設(shè)置為主面與從面，在接觸進(jìn)行的各時間段內(nèi)，自動檢測主從面是否發(fā)生穿透。未穿透時，接觸繼續(xù)進(jìn)行；當(dāng)發(fā)生穿透時，在主面與從面之間引入垂直于主面法向的接觸力，作用方式類似于彈簧，如圖1 所示[10]。這一接觸力F可由下式得出：

式中：k為主面剛度；δ為穿透深度值。

主面剛度k與接觸面剛度的懲罰因子相關(guān)，在Lsdyna 中，懲罰因子的設(shè)置關(guān)乎數(shù)值模擬精度[11]。k的計算公式如下：

式中：pf為懲罰因子；V為主體積；K為接觸單元體積彈性模量；A為主從面接觸段面積。

圖1 罰函數(shù)法作用原理Fig. 1 Principle of penalty function method

2 船-冰-水耦合模型

以某船首部為研究對象，如圖2 所示。首部鋼材為907A 高強(qiáng)鋼，名義屈服強(qiáng)度為390 MPa，材料模型選擇各向同性塑性模型，材料參數(shù)如表1 所示。

圖2 船舶首部模型Fig. 2 Model of ship bow

表1 船體材料Tab. 1 Ship material parameters

船-冰-水耦合數(shù)值模型由船體、浮冰、海水、空氣四部分構(gòu)成。在Ls-dyna 中，S-ALE 方法對于流體域PART 的設(shè)置通過網(wǎng)格PART 和材料PART 的組合進(jìn)行定義。同時，水及空氣流體材料的物理特性通過材料與狀態(tài)方程結(jié)合來描述的，表達(dá)式如下：

式中：C0～C6為多項式擬合參數(shù)，?為體積參數(shù)，E0為初始內(nèi)能。

浮冰采用SOLID164 實體單元進(jìn)行建模，模型尺寸為10m×10m×1m，網(wǎng)格尺寸為0.2m×0.2m×0.2m。冰材料模型選擇遵循Mohr-Coulomb 屈服準(zhǔn)則[12]的各向同性彈塑性本構(gòu)模型。根據(jù)國內(nèi)外學(xué)者相關(guān)試驗研究和海冰數(shù)據(jù)統(tǒng)計，本文海冰材料參數(shù)如表2 所示。

表2 空氣、水狀態(tài)方程參數(shù)Tab. 2 The state equation parameters of air and water

船-冰-水耦合數(shù)值模型如圖3 所示

圖3 船-冰-水耦合模型Fig. 3 Ship-ice-water coupling model

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

3.1 船-冰-水耦合結(jié)果分析

船-冰水相互作用過程，水阻力與冰載荷是船體載荷響應(yīng)的主要成分。首先對敞水環(huán)境下水阻力進(jìn)行研究，對船體施加縱蕩方向的速度模擬船舶在恒定功率下勻速前進(jìn)。航速控制分別為2 kn，5 kn，8 kn，10 kn，15 kn。在Ls-dyna 可以實現(xiàn)船舶前進(jìn)過程中的興波現(xiàn)象模擬，如圖4 所示。船舶首部型破開水域前進(jìn)，將水向兩邊推開，自由液面發(fā)生改變，波高出現(xiàn)區(qū)域船身壓力增大，波谷區(qū)域船身壓力降低，形成興波現(xiàn)象。隨著船舶繼續(xù)前行，在船體水線下部分出現(xiàn)擾流，最終在尾部形成尾波現(xiàn)象，從首波、船身兩側(cè)的V-形擴(kuò)散波和尾波軌跡與開爾文船行波系一致[13]。

圖4 艦船航行數(shù)值模擬結(jié)果Fig. 4 Numerical simulation results of ship navigation

圖5 開爾文波波峰軌跡Fig. 5 Kelvin wave crest trajectory

表3 海冰力學(xué)參數(shù)Tab. 3 Mechanical parameters of sea ice

圖6 為15 kn 航速下船體外板的應(yīng)力云圖，水線下部分應(yīng)力隨著水深而逐漸變大，船身表面最大應(yīng)力出現(xiàn)在水線以下首柱破波區(qū)域和球鼻艏處。

圖6 水阻力作用下外板應(yīng)力云圖Fig. 6 Hull stress diagram under water resistance

圖7 為不同航速下敞水阻力時歷曲線，對數(shù)據(jù)進(jìn)行濾波處理來減小振蕩幅度?？梢钥闯?，在0～2 s內(nèi)，船體在推力作用下，導(dǎo)致靜水?dāng)_動，水阻力值周期振蕩提升，3 s 后當(dāng)水域趨于穩(wěn)定后，船體表面受力趨于固定值，逐漸收斂。對收斂后水阻力值與航速關(guān)系進(jìn)行分析，如圖8 所示，當(dāng)航速從2 kn 逐漸增加至15 kn，水阻力值從3.5×105N 增加至7.5×105N，本船未進(jìn)行船體模型水阻力試驗，將結(jié)果與S60 系列船型通過三因次法換算得到的船體總阻力對比，阻力值量級相同。并且在Fr值相同的情況下，水總阻力預(yù)報值與船體航速存在正相關(guān)的趨勢。

圖7 不同航速下水阻力時歷曲線Fig. 7 Time history curves of water resistance

圖8 水阻力穩(wěn)定值隨航速變化Fig. 8 Variation of water resistance with speed

在確定船-水耦合下的水阻力規(guī)律后，對冰水耦合下的水阻力、冰阻力作用機(jī)理進(jìn)行研究。在多相耦合時，船體速度不再以恒定功率下航行，給船舶施加10 kn 的初速度，擁有初始動能的船首與浮冰發(fā)生沖撞，浮冰的存在導(dǎo)致水阻力與冰載荷作用機(jī)理之間存在相互影響。

圖9 為船-冰-水三相耦合下的作用過程，船與冰接觸發(fā)生碰撞。在碰撞力作用下，產(chǎn)生排水效應(yīng)，在浮冰周圍形成環(huán)形波浪，同時浮冰浮態(tài)發(fā)生改變，在波浪中呈現(xiàn)動態(tài)平衡效果。

圖9 船-冰-水相互作用數(shù)值模擬Fig. 9 Numerical simulation of ship-ice-water interaction

對碰撞時間內(nèi)的冰載荷時歷曲線進(jìn)行分析，浮冰與船體發(fā)生碰撞時，首柱處受到冰載荷的碰撞與擠壓作用，并分為與船體運(yùn)動方向相反的縱向力、沿船寬方向?qū)ΨQ分布的橫向力及沿Z軸向上的垂向力。時歷曲線如圖10 所示，而冰單元在碰撞后發(fā)生局部破碎失效導(dǎo)致峰值出現(xiàn)多次卸載現(xiàn)象。從時歷曲線可發(fā)現(xiàn)，縱向力是冰阻力的主要組成部分，并且在時間上先于橫向力出現(xiàn)。因為浮冰擠壓作用是關(guān)于船線對稱，橫向破冰力表現(xiàn)為沿0 附近的振蕩狀態(tài)。船體在與浮冰接觸過程中，船體自身重力對浮冰存在壓潰作用，在冰阻力上表現(xiàn)為持續(xù)存在的垂向冰載荷。

圖10 冰載荷時歷曲線Fig. 10 Time curve of ice load

將冰-水耦合作用下的冰載荷與無水域狀態(tài)時約束Z向自由度的冰載荷進(jìn)行對比，如圖11 所示。水域浮力的作用下弱化了船-冰碰撞時的接觸，使得峰值數(shù)據(jù)較無水狀態(tài)下，減小了32%，整體載荷趨勢也有所下降。因此計算船冰碰撞時，水域的影響效果是不能忽視的。

圖11 冰阻力時歷曲線對比Fig. 11 Comparison of ice load time curve

船-冰-水耦合作用下的水阻力時歷曲線如圖12 所示，呈逐漸衰減直至消失趨勢。這是因為水阻力在船舶與浮冰發(fā)生碰撞后，船體動能一部分轉(zhuǎn)化為浮冰內(nèi)能與動能，轉(zhuǎn)化為水域內(nèi)能部分減少。同時隨著船體速度下降，擊水作用力與排水速度降低，導(dǎo)致船體表面水阻力逐漸衰減。

圖12 水阻力時歷曲線Fig. 12 Time curve of water resistance

3.2 冰型尺寸影響分析

浮冰厚度與尺寸是規(guī)則浮冰與船舶首部碰撞的重要冰型參數(shù)，在船-冰-水相互作用機(jī)理研究的基礎(chǔ)上，改變浮冰厚度或浮冰尺寸，展開初速度為10 kn 下的碰撞仿真分析。工況設(shè)置如表4 所示，水域及碰撞參數(shù)設(shè)置不作改變。

表4 船-冰碰撞計算工況Tab. 4 Ship-ice collision calculation conditions

分別以冰厚和浮冰尺寸為變形分析case1～case7 工況船體所受冰阻力與水阻力，各冰厚和浮冰尺寸冰阻力時歷曲線對比如圖13 所示。

圖13 不同冰型參數(shù)下冰阻力時歷曲線Fig. 13 Ice load time history curves under different ice type parameters

浮冰厚度越小，冰單元破壞程度越大，失效單元數(shù)越多，峰值卸載頻率越高，峰值與冰厚變化成正相關(guān)關(guān)系；冰阻力在浮冰邊長達(dá)到20 m 與50 m 時出現(xiàn)畸變，船-冰作用時間隨浮冰尺寸的增加呈現(xiàn)增加趨勢。當(dāng)浮冰尺寸為50 m 時，冰阻力持續(xù)存在，這是因為船舶本身動能不足以使巨型浮冰產(chǎn)生明顯速度變化，反之船舶受到浮冰質(zhì)量的阻力影響。其初速度下降更為明顯，二者之間一直處于接觸狀態(tài)，擠壓與碰撞冰載荷也一直存在。不同冰厚下的冰載荷峰值如表5 所示，不同浮冰尺寸下的載荷峰值如表6 所示。

表5 不同冰厚下的冰載荷峰值Tab. 5 The peak ice load under different ice thickness

表6 不同浮冰尺寸下的載荷峰值Tab. 6 The peak ice load under different floating ice sizes

水阻力與船舶航速關(guān)系密切，水阻力值幾乎不隨浮冰參數(shù)發(fā)生改變，趨勢與圖12 基本相同，僅隨船舶速度衰減而逐漸減小。從冰載荷在總載荷值中占比發(fā)現(xiàn)，浮冰厚度為0.5m 時，冰載荷峰值為2.11MN。當(dāng)浮冰厚度增加1 倍、1.4 倍和2 倍時，峰值分別增加29%、65%與138%。浮冰邊長增加給浮冰質(zhì)量的增加是二次方級的，因此，冰載荷增長速度存在躍遷趨勢，且當(dāng)浮冰尺寸大于型寬時，增長趨勢減緩，達(dá)到冰載荷閾值。假設(shè)浮冰尺寸無限增加，即可視為覆蓋在海面的層冰，不產(chǎn)生移動速度，船體初始動能對浮冰的作用效果幾乎完全轉(zhuǎn)化為浮冰破碎失效的變形能。而浮冰越厚，尺寸越大，水域?qū)Ρd荷的弱化效果越弱，冰載荷在總載荷中占比越高。

圖14 冰阻力峰值隨冰型參數(shù)變化Fig. 14 The peak value of ice load varies with ice parameters.

4 結(jié) 語

為了研究艦船在航行時遭遇浮冰，載荷響應(yīng)情況，本文基于S-ALE 流固耦合方法對冰區(qū)海域浮冰漂浮和船-冰-水相互作用進(jìn)行數(shù)值模擬，并對不同冰型參數(shù)下載荷影響進(jìn)行了分析。主要結(jié)論歸納如下：

1）通過對船-冰-水相互作用過程中船行波以及船-冰碰撞現(xiàn)象的準(zhǔn)確模擬，驗證了S-ALE 流固耦合方法處理船-冰問題的可行性。

2）船舶敞水域航行過程中，當(dāng)航速達(dá)到穩(wěn)定后，會產(chǎn)生開爾文船行波，并且水阻力趨于穩(wěn)定值，與航速存在正相關(guān)關(guān)系。同時將船-冰-水相互作用的冰載荷結(jié)果與無水狀態(tài)下約束浮冰Z向自由度的浮冰相比，冰載荷峰值減弱了32%，證明了在考慮船-冰作用時水域的重要性。

3）水阻力僅與航速相關(guān)，而冰載荷受到水域外部因素及冰厚、浮冰尺寸等自身因素的影響。水域在船冰碰撞過程對冰載荷起到了削弱作用，而這一削弱效果隨浮冰尺寸、冰厚的增加而減弱。浮冰尺寸、冰厚越大，冰載荷在艦船航行總載荷中占比越大，且尺寸因素影響大于浮冰厚度。