汪春輝 王嘉安 王 超 郭春雨 朱廣元

(哈爾濱工程大學(xué)船舶工程學(xué)院，哈爾濱 150001)

引言

在極區(qū)，平整冰與海洋結(jié)構(gòu)物間的相互作用形式主要為擠壓破壞和彎曲破壞，冰層受結(jié)構(gòu)物作用產(chǎn)生斷裂破碎、裂紋擴(kuò)展和滑移清除過程[1-4].在某些特殊情況，結(jié)構(gòu)物需要上浮穿透冰層破出水面，對冰層施以自下而上的垂向載荷，而冰層垂直貫穿破裂過程十分復(fù)雜.

學(xué)界對冰與結(jié)構(gòu)物垂直相互作用的研究還不充分，國內(nèi)外相關(guān)研究主要以試驗(yàn)和理論方法為主.由于試驗(yàn)條件要求苛刻及數(shù)值技術(shù)的發(fā)展，一些具備完善冰作用機(jī)理的數(shù)值方法成為解決這些問題的好方法.并且海冰是一種復(fù)雜的材料，有許多因素會(huì)影響冰的物理特性，例如，溫度孔隙率和鹽度、速率依賴性.由于冰可發(fā)生大而復(fù)雜的變形，模擬冰與結(jié)構(gòu)的相互作用是一個(gè)巨大的挑戰(zhàn).

目前，有限元方法(FEM)、離散元方法(DEM)、光滑粒子法(SPH)和近場動(dòng)力學(xué)(PD)等數(shù)值模擬方法已被用于模擬冰與結(jié)構(gòu)物間的相互作用，本文采用有限元方法開展圓柱體垂直出水破冰研究.針對水下結(jié)構(gòu)物在冰層下方上浮至破出水面的過程的研究在世界范圍內(nèi)仍處于起步階段，國內(nèi)外公開資料極其稀少.因此利用有限元方法分析結(jié)構(gòu)物與冰層垂直相互作用進(jìn)展可以參考船與冰相互接觸的研究資料[5-7].

針對冰與結(jié)構(gòu)物耦合作用有限元數(shù)值分析方法，國際上已有大量的研究探索.Ehlers 和Kujala[8]在LS-DYNA 有限元軟件中使用彈塑性應(yīng)變率材料模擬冰的彎曲行為，并與4 點(diǎn)彎曲試驗(yàn)結(jié)果進(jìn)行對比.Polach 和Franzvon[9]應(yīng)用LS-DYNA 有限元軟件采用MAT_DAMAGE_3 材料模擬冰的力學(xué)特性，并基于損傷力學(xué)對數(shù)學(xué)模型進(jìn)行分析，將模擬結(jié)果與冰池試驗(yàn)結(jié)果對比.Gagnon[10]提出了一種彈性線性硬化材料，可模擬冰力學(xué)特征，但不能模擬冰的破壞斷裂行為.Wang 和Derradji-Aouat[11]基于罰函數(shù)法及拉格朗日歐拉法模擬碎冰對系泊結(jié)構(gòu)的作用過程，以評估系泊結(jié)構(gòu)的整體載荷.Ranta 等[12]通過FEM-DEM 耦合方法構(gòu)建冰與結(jié)構(gòu)物相互作用二維數(shù)值模型，利用多元回歸線性數(shù)學(xué)模型和變量消元法對冰載荷特征值分析，探究冰載荷影響因素.國內(nèi)方面，郭春雨等[13]基于LS-DYNA 有限元軟件模擬了極地船和浮冰的耦合作用過程，實(shí)現(xiàn)了冰區(qū)船舶在碎冰區(qū)域的航行阻力預(yù)報(bào).Wang 等[14]應(yīng)用FEM方法與DEM-CFD 方法模擬了冰區(qū)船舶在冰緣區(qū)的航行過程，并將數(shù)值模擬結(jié)果與模型試驗(yàn)對比，分析了碎冰載荷的主要影響因素.黃焱等[15]基于FEM方法構(gòu)建并分析了固定約束及摩擦約束邊界條件下的冰溫度膨脹力，發(fā)現(xiàn)固定約束會(huì)造成高度的應(yīng)力集中.任奕舟和鄒早建[16]、王峰[17]基于FEM 方法，通過粘聚力單元和彈塑性冰本構(gòu)模型構(gòu)建船舶與海冰碰撞場景的數(shù)值仿真，驗(yàn)證了粘聚元方法模擬海冰的準(zhǔn)確性與可行性.Chen 等[18]應(yīng)用FEM-SPH 方法建立了錐體和破冰船與平整冰相互作用的數(shù)值模型，預(yù)測了破冰阻力及碎冰堆積過程.

本文應(yīng)用 LS-DYNA 有限元軟件采用結(jié)構(gòu)化-任意拉格朗日歐拉 (S-ALE)流固耦合方法及罰函數(shù)接觸算法，建立了冰、水、結(jié)構(gòu)物耦合作用數(shù)值模擬方法，自主搭建了圓柱體垂直貫穿冰層試驗(yàn)臺(tái)架.通過冰梁結(jié)構(gòu)4 點(diǎn)彎曲及楔形體入水?dāng)?shù)值模擬，驗(yàn)證了彈塑性應(yīng)變率本構(gòu)材料表征冰力學(xué)特性的準(zhǔn)確性和流固耦合方法的可行性，對比了無水環(huán)境下圓柱體垂直貫穿冰層和圓柱體垂直出水破冰過程中的冰載荷、冰層撓度變化、冰層應(yīng)力響應(yīng)及斷裂現(xiàn)象，以期為冰層與海洋結(jié)構(gòu)物垂向相互作用研究提供參考.

1 理論基礎(chǔ)

LS-DYNA 有限元軟件被廣泛應(yīng)用于航天、航海和武器裝備等工程制造領(lǐng)域.由于LS-DYNA 有限元軟件包在結(jié)構(gòu)計(jì)算、材料損傷破壞等領(lǐng)域的優(yōu)勢，越來越多的學(xué)者將其應(yīng)用于船冰碰撞、冰槳切削等冰工程方面的研究中，運(yùn)用Lagrange 格式的顯式中心差分法計(jì)算大變形的結(jié)構(gòu)物與冰層耦合作用問題.本文采取結(jié)構(gòu)化任意拉格朗日歐拉方法，基于罰函數(shù)流固耦合算法模擬圓柱體上浮垂直出水破冰中冰、水、結(jié)構(gòu)物耦合作用的強(qiáng)非線性動(dòng)態(tài)過程.

1.1 結(jié)構(gòu)化-任意拉格朗日歐拉流固耦合方法

結(jié)構(gòu)化任意拉格朗日歐拉方法原理上是一種基于統(tǒng)一連續(xù)介質(zhì)的動(dòng)量守恒方程的二階精度的平流算法[19].通過接觸算法獲取固體之間的接觸作用時(shí)的界面力，能夠有效解決物體在空間上存在大位移同時(shí)本身存在大變形問題.在S-ALE 計(jì)算中，定義了材料域和參考域，材料的速度為u，參考系的移動(dòng)速度為v，速度差u-v表示為w，由此可對S-ALE 計(jì)算方法的控制方程進(jìn)行推導(dǎo)，如下所示

式中J為雅可比矩陣，表示兩個(gè)域的相對體積微分.

S-ALE 公式中的材料時(shí)間導(dǎo)數(shù)寫成公式

式中，br表示b的參考域函數(shù).

S-ALE 方法的質(zhì)量、動(dòng)量以及能量守恒方程如下

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

罰函數(shù)接觸算法將可能發(fā)生接觸作用的兩個(gè)表面分別定義為主表面和從表面，發(fā)生穿透時(shí)，類似在主表面和從表面之間放置一系列法向彈簧，彈簧的作用是為了限制穿透，彈簧的力被稱為界面接觸力.

罰函數(shù)法的接觸力由下面計(jì)算公式得出

式中k為單元接觸面剛度，δ為穿透量.

單元接觸面剛度k由如下公式?jīng)Q定

式中pf為單元接觸面剛度的懲罰因子，K是接觸單元體積彈性模量，A是接觸段面積，V是主段體積.

在數(shù)學(xué)方法上，用歐拉變量描述流體的動(dòng)態(tài)，用拉格朗日變量描述結(jié)構(gòu)的運(yùn)動(dòng)邊界，通過分布節(jié)點(diǎn)力和插值速度來表示流場和結(jié)構(gòu)物的交互作用.該方法適用于廣泛的流固耦合應(yīng)用，界面接觸力為估算求得且可能存在數(shù)值振蕩[20].

1.3 冰本構(gòu)模型

冰的特性取決于溫度、應(yīng)變率和載荷方向，冰在受力變形破壞過程中，受加載速率的影響，會(huì)呈現(xiàn)韌脆轉(zhuǎn)換現(xiàn)象[21].本文選用冰本構(gòu)材料為彈塑性應(yīng)變率本構(gòu)材料*MAT _PLASTICITY_COMPRESSION_TENSION，該材料為一種各向同性的彈塑性材料，能夠模擬不同的壓縮和拉伸響應(yīng).此外，它還考慮了應(yīng)變率敏感性[22-24].

該材料模型需要定義兩條用于拉伸與壓縮行為的屈服應(yīng)力對有效塑性應(yīng)變的曲線，以及一條應(yīng)變率對壓縮屈服應(yīng)力比例系數(shù)(CYSF)的曲線，用于表征應(yīng)變率對壓縮強(qiáng)度的影響.假設(shè)彎曲強(qiáng)度與應(yīng)變率無關(guān)，拉伸和壓縮下的屈服應(yīng)力(σy)與有效塑性應(yīng)變(εp)的曲線均由以下方程推導(dǎo)得出

式中σ0為初始屈服應(yīng)力.

假定壓縮和拉伸過程中冰材料的塑性硬化模量(Ep)相同.通過Sánchez 等[25]給出的冰材料壓縮強(qiáng)度應(yīng)變率敏感性擬合函數(shù)可以計(jì)算得出應(yīng)變率()為10-3s-1時(shí)冰材料的初始壓縮屈服強(qiáng)度()

為了定義應(yīng)變率對壓縮屈服應(yīng)力比例系數(shù)的曲線，不同應(yīng)變率下的壓縮強(qiáng)度由應(yīng)變率為10-3s-1時(shí)的初始壓縮強(qiáng)度進(jìn)一步歸一化得出，本文采用的海水冰材料應(yīng)變率為10-3s-1時(shí)的初始壓縮強(qiáng)度為5.8 MPa.海水冰應(yīng)變率和相應(yīng)的壓縮屈服應(yīng)力比例系數(shù)如表1 所示.對于拉伸行為，實(shí)驗(yàn)結(jié)果顯示脆性破壞及應(yīng)變率對冰拉伸強(qiáng)度的影響可以忽略不計(jì)[26].

表1 海水冰應(yīng)變率和壓縮屈服應(yīng)力比例系數(shù)Table 1 Strain rates and compressive yield stress scale factors of sea ice

該材料模型采用等效塑性應(yīng)變破壞準(zhǔn)則，當(dāng)單元的等效塑性應(yīng)變達(dá)到定義的臨界值(εf)時(shí)，該單元將被刪除.可以通過調(diào)整εf值以獲得適當(dāng)?shù)膹澢憫?yīng).

2 數(shù)值計(jì)算模型

2.1 無水環(huán)境下圓柱體垂直貫穿冰層有限元模型

無水環(huán)境下圓柱體垂直貫穿冰層數(shù)值模型包含以下兩個(gè)部分:冰層和圓柱結(jié)構(gòu)體，均采用拉格朗日算法進(jìn)行描述.模型的幾何尺寸見下圖1，模型冰厚為3 cm，長寬均為50 cm;圓柱體直徑為3 cm，高為5 cm.冰層邊界條件采用四邊約束，圓柱體置于冰層正中心下1 cm 處.

圖1 數(shù)值模型尺寸Fig.1 Size of numerical model

本文將圓柱體定義為剛體，冰層本構(gòu)選用彈塑性應(yīng)變率材料模型，根據(jù)數(shù)值模擬需求將冰材料分為淡水冰及海水冰兩類.主要材料參數(shù)如表2所示[24-26].

表2 主要材料參數(shù)Table 2 Main material parameters

圓柱體與冰層的相互作用以及冰單元間相互作用通過侵蝕接觸算法實(shí)現(xiàn)，冰層單元不斷發(fā)生破壞并刪除失效單元形成裂紋，一旦網(wǎng)格單元失效且接觸表面發(fā)生變化，通過侵蝕接觸可重新建立表面形狀[21].

2.2 圓柱體垂直出水破冰有限元模型

圓柱體垂直出水破冰數(shù)值模型設(shè)置與文2.1 節(jié)開展的圓柱體垂直貫穿冰層試驗(yàn)工況保持一致.包含4 部分:空氣域、水域、圓柱結(jié)構(gòu)體和冰層，其中，空氣域和水域采用歐拉算法進(jìn)行描述，圓柱結(jié)構(gòu)體和冰層仍采用拉格朗日算法進(jìn)行描述.空氣域和水域的幾何尺寸及網(wǎng)格劃分如圖2，空氣域和水域的采用漸進(jìn)式網(wǎng)格劃分，即應(yīng)用*ALE_STRUCTURED_MESH 和*ALE_STRUCTURED_MESH_CONTROL_OINTS 關(guān)鍵字生成S-ALE 結(jié)構(gòu)化網(wǎng)格.水和空氣通過*MAT_NULL 材料來描述，水的狀態(tài)方程為OS-GRUNEISEN 來實(shí)現(xiàn)，空氣的控制方程為EOSINER-POLYNOMIAL，空氣與水的相關(guān)模型參數(shù)設(shè)置如表3，圓柱體結(jié)構(gòu)及冰層結(jié)構(gòu)參數(shù)設(shè)置與無水環(huán)境的數(shù)值模型參數(shù)設(shè)置保持一致.

表3 空氣、水相關(guān)參數(shù)設(shè)置Table 3 Setting of air and water related parameters

圖2 圓柱體垂直出水破冰數(shù)值模型Fig.2 Numerical model of vertical upward water of cylinder breaking through level ice

針對圓柱體垂直出水破冰流固耦合問題，通過INITIAL_HYDROSTATIC_ALE 關(guān)鍵字實(shí)現(xiàn)流體區(qū)域靜水壓力分層，使用*CONSTRAINED_AGRANGE_IN_SOLID 關(guān)鍵字實(shí)現(xiàn)冰層、圓柱體和流體區(qū)域間的流固耦合作用.

3 計(jì)算方法驗(yàn)證

3.1 冰本構(gòu)模型驗(yàn)證

Sodhi[27]指出冰層在垂向載荷作用下的破壞方式主要以彎曲破壞為主.通過對冰梁的4 點(diǎn)彎曲模型試驗(yàn)進(jìn)行數(shù)值模擬，驗(yàn)證本文所用應(yīng)變率彈塑性材料模擬冰彎曲破壞過程的可靠性.以Kujala 等[28]冰梁彎曲破壞模型試驗(yàn)為依據(jù)，模型幾何尺寸如圖3，上方的支撐架剛性固定，下方的支撐架以2.756 mm/s勻速向上移動(dòng).冰材料參數(shù)與表2 中的海水冰材料一致.

圖3 冰梁四點(diǎn)彎曲數(shù)值模擬幾何模型Fig.3 Geometric model for numerical simulation of four-point bending of ice beam

冰梁受到彎曲載荷后在加載支撐處形成斷裂裂紋如圖4(a)，通過對比文獻(xiàn)中基于SPH 方法的斷裂現(xiàn)象即圖5(b)發(fā)現(xiàn)，本文模擬的冰梁4 點(diǎn)彎曲斷裂現(xiàn)象與試驗(yàn)[28]及SPH 方法[23]計(jì)算結(jié)果吻合度極高.

冰梁4 點(diǎn)彎曲過程的彎曲載荷-時(shí)間曲線如圖5，試驗(yàn)及數(shù)值計(jì)算結(jié)果的曲線趨勢以及峰值點(diǎn)一致性良好.結(jié)合圖4 冰梁4 點(diǎn)彎曲破壞斷裂現(xiàn)象可以說明本文所用應(yīng)變率彈塑性材料可以準(zhǔn)確模擬冰層彎曲破裂問題.

圖5 冰梁4 點(diǎn)彎曲數(shù)值驗(yàn)證模擬結(jié)果Fig.5 Numerical verification of simulation results of four-point bending of ice beam

3.2 圓柱體垂直貫穿冰層試驗(yàn)

基于自主開發(fā)的試驗(yàn)臺(tái)架、以凍結(jié)模型冰為試驗(yàn)對象、開展無水環(huán)境下圓柱體垂直貫穿冰層試驗(yàn)，探究該動(dòng)態(tài)破冰過程的現(xiàn)象及破冰載荷.試驗(yàn)臺(tái)架如圖6，主要儀器包括伺服電動(dòng)缸、數(shù)據(jù)采集與分析處理系統(tǒng)和接觸式壓力傳感器等，試驗(yàn)臺(tái)架主要用來固定模型冰以及搭載伺服電動(dòng)缸

圖6 試驗(yàn)用主要儀器Fig.6 Main instruments for test

試驗(yàn)使用凍結(jié)型模型冰作為模型材料，將煮沸后的純凈水倒入模具中冷卻凝固至定型，凍結(jié)24 h定型后取出進(jìn)行尺寸修正，再將修正過的模型冰放于低溫試驗(yàn)箱中儲(chǔ)存24 h 后進(jìn)行試驗(yàn)，低溫試驗(yàn)箱環(huán)境溫度恒定為-25 °C.為測量制備模型冰的力學(xué)性能參數(shù)，分別進(jìn)行了單軸壓縮試驗(yàn)和巴西圓盤試驗(yàn)測量其抗壓強(qiáng)度及拉伸強(qiáng)度.單軸壓縮試驗(yàn)和巴西圓盤試驗(yàn)壓頭加載速率均為4 mm/s，與圓柱體垂直貫穿冰層試驗(yàn)加載速率保持一致.

單軸壓縮試驗(yàn)試驗(yàn)選用長方體模型冰試樣，冰試樣尺寸為40 mm×40 mm×100 mm，盡可能減少冰體試樣本身對試驗(yàn)結(jié)果的影響，測得模型冰試樣平均抗壓強(qiáng)度為2.41 MPa.單軸壓縮試驗(yàn)典型破壞過程如圖7.

巴西圓盤法是一種間接測量抗拉強(qiáng)度的方法，能夠有效避開直拉法試驗(yàn)過程中存在的冰試件夾持和加載困難的缺點(diǎn)[29].試驗(yàn)用圓盤狀冰體試樣直徑為10 cm 厚度為5 cm，測得模型冰試樣平均拉伸強(qiáng)度為0.39 MPa.巴西圓盤試驗(yàn)典型破壞過程如圖8.

Sodhi[27]發(fā)現(xiàn)冰層的垂向貫穿破裂過程中的大撓度變形僅發(fā)生于施加載荷附近區(qū)域，通常小于冰厚度的10～ 20 倍，且固定約束(四邊約束)下冰載荷穩(wěn)定性優(yōu)于自由約束下冰載荷穩(wěn)定性，但固定約束(四邊約束)下冰載荷均值要大于自由約束下冰載荷均值.為獲得可靠的試驗(yàn)結(jié)果及規(guī)律性結(jié)論，選取符合尺度效應(yīng)的冰層尺寸(50 cm×50 cm×3 cm)的模型冰板并采用固定約束(四邊約束)進(jìn)行多組試驗(yàn)，如圖9.將模型冰板置于試驗(yàn)臺(tái)架上，鋁型材扣置于模型冰上方約束冰層運(yùn)動(dòng).試驗(yàn)用圓柱體試樣直徑3 cm，高5 cm，壓載速率為4 mm/s，圓柱體試樣勻速上升與冰層接觸，通過數(shù)據(jù)采集系統(tǒng)和高速攝像機(jī)，記錄模型冰板斷裂破壞過程中冰載荷的時(shí)歷曲線以及圓柱體試樣與模型冰接觸作用時(shí)冰層破裂現(xiàn)象.

圖9 試驗(yàn)用凍結(jié)模型冰Fig.9 Frozen model ice used in the test

圓柱體試樣與冰層接觸后，冰層試樣不會(huì)瞬間破裂，如圖10(a)，隨著圓柱體壓頭的上壓，冰層中間部分受壓發(fā)生向上彎曲，冰層四周向上翹曲，如圖10(b).隨著圓柱體試樣的進(jìn)一步上壓，裂紋從冰層中部產(chǎn)生，瞬間擴(kuò)散至冰層四周，并在四邊約束的作用下，冰層先會(huì)沿約束內(nèi)側(cè)產(chǎn)生周向斷裂，最終被頂起形成碎冰堆，如圖10(c).

圖10 試驗(yàn)現(xiàn)象Fig.10 Experimental phenomenon

圖10 試驗(yàn)現(xiàn)象(續(xù))Fig.10 Experimental phenomenon (continued)

由于在冰層的制取過程中，冰的形狀、厚度存在一定不可避免的差異，無法保證每次試驗(yàn)?zāi)Ｐ屯耆嗤瑢?dǎo)致冰載荷極值與作用時(shí)間存在一定的差異.選取4 組3 cm 厚度冰層進(jìn)行試驗(yàn)重復(fù)性驗(yàn)證，其試驗(yàn)結(jié)果如圖11，結(jié)果表明整體試驗(yàn)數(shù)據(jù)一致性良好.

圖11 冰層垂直貫穿破裂試驗(yàn)重復(fù)性驗(yàn)證Fig.11 Repeatability verification of vertical penetr-ation test of level ice

3.3 圓柱體垂直貫穿冰層數(shù)值模型驗(yàn)證

網(wǎng)格質(zhì)量影響數(shù)值計(jì)算精度，過多的網(wǎng)格往往會(huì)極大的延長數(shù)值計(jì)算時(shí)間，而數(shù)值計(jì)算采用的模型計(jì)算單元需要保證結(jié)果準(zhǔn)確性的同時(shí)提高計(jì)算效率.由于圓柱體模型為剛體材料，網(wǎng)格質(zhì)量影響相對較小，所以為驗(yàn)證網(wǎng)格合理性，分別對冰層數(shù)值模型劃分大小為2.5 mm，3.5 mm 以及4 mm 的正方形網(wǎng)格，進(jìn)行網(wǎng)格無關(guān)性分析.

如圖12 所示，圖中3 條曲線總體來看趨勢一致且峰值點(diǎn)近似相同，呈振蕩上升趨勢且冰載荷特征穩(wěn)定，2.5 mm 網(wǎng)格單元曲線與3.5 mm 網(wǎng)格單元曲線擬合度良好，4 mm 網(wǎng)格單元曲線載荷卸載時(shí)間相對提前.所以選取冰層網(wǎng)格單元大小為3.5 mm，既可以保證數(shù)值精度又可以提高計(jì)算效率.

圖12 不同網(wǎng)格尺寸下的冰載荷時(shí)歷曲線Fig.12 Time history curve of ice load under different element sizes

開展無水環(huán)境下圓柱體垂直破冰數(shù)值仿真以驗(yàn)證數(shù)值方法準(zhǔn)確性，數(shù)值計(jì)算模型見1.2 節(jié)，冰層材料選用表2 中淡水冰材料，計(jì)算工況與上述試驗(yàn)一致，圖13 為數(shù)值模擬結(jié)果.從圖13(a)應(yīng)力云圖中可以看出，當(dāng)圓柱體與冰層接觸而冰層并未發(fā)生斷裂現(xiàn)象時(shí)，冰層發(fā)生彈性變形，應(yīng)力集中在冰層中心，由于冰層采用四邊約束應(yīng)力區(qū)域呈環(huán)狀.隨著圓柱體與冰層的持續(xù)作用，部分冰層單元達(dá)到失效極限，出現(xiàn)徑向裂紋向四周擴(kuò)展，裂紋相互作用進(jìn)一步形成楔子，如圖13(b)所示.當(dāng)圓柱體繼續(xù)上升，由于應(yīng)力集中于冰層中央，導(dǎo)致冰層中部破碎形成碎冰堆積，冰層產(chǎn)生周向裂紋，隨后冰層向中心塌陷，如圖13(c)所示，對比圖10(c)可以看出冰層向中心塌陷整體現(xiàn)象以及破壞范圍均呈現(xiàn)良好的一致性，且均存在約束下的徑向裂縫在其形成后不再進(jìn)一步擴(kuò)展，并且由徑向裂縫形成的大塊楔形物隨圓柱體進(jìn)一步加載而隨之推動(dòng).

圖13 圓柱體垂直破冰數(shù)值模擬冰層破裂過程Fig.13 Numerical simulation of level ice rupture process by cylinder vertical ice breaking

圖13 圓柱體垂直破冰數(shù)值模擬冰層破裂過程(續(xù))Fig.13 Numerical simulation of level ice rupture process by cylinder vertical ice breaking (continued)

圖14 為圓柱體垂直破冰時(shí)數(shù)值模擬與試驗(yàn)結(jié)果的冰載荷時(shí)歷曲線.在圓柱體接觸冰層且勻速上升過程中，數(shù)值模擬冰載荷數(shù)據(jù)呈震蕩上升的趨勢，伴隨著冰層裂紋產(chǎn)生和擴(kuò)展，整體趨勢、載荷峰值且突破時(shí)間均與試驗(yàn)值一致.

圖14 圓柱體垂直破冰載荷數(shù)值與試驗(yàn)對比曲線Fig.14 Numerical and experimental comparison curve of vertical ice breaking load of cylinder

3.4 基于S-ALE 方法的楔形體入水?dāng)?shù)值模擬

為驗(yàn)證使用S-ALE 方法分析流固耦合問題的準(zhǔn)確性與可行性，應(yīng)用此方法模擬了二維楔形體入水砰擊過程，并與陳光茂等[30]基于SPH 方法模擬二維楔形體入水問題的數(shù)值結(jié)果及Zhao 和Faltinsen[31]的實(shí)驗(yàn)結(jié)果進(jìn)行對比.

依據(jù)文獻(xiàn)[30]建立二維楔形體入水?dāng)?shù)值模型，水池長2.4 m，高1.3 m，以此保證入水砰擊的過程中壁面反射不會(huì)對楔形體受力產(chǎn)生干擾.楔形體形狀為等腰三角形，頂邊長度為0.5 m，斜升角(即斜邊與水平方向的夾角)為30°.在楔形體入水前的瞬間，它的垂向速度為6.15 m/s，隨后自由入水，計(jì)算使用的二維楔形體為實(shí)心剛體，重量為 241 kg.

圖15 為S-ALE 方法模擬的楔形體入水過程自由表面變形，與文獻(xiàn)[30]試驗(yàn)現(xiàn)象一致，通過圖16入水時(shí)間0.015 8 s 時(shí)壓力分布對比圖可知S-ALE 方法數(shù)值模擬結(jié)果與文獻(xiàn)[30]中應(yīng)用SPH 方法模擬結(jié)果相同，且圖17 中入水垂向速度隨時(shí)間變化對比曲線吻合性良好，證明了S-ALE 方法處理流固耦合問題的準(zhǔn)確性與合理性.

圖15 楔形體入水過程自由表面變形Fig.15 Deformation of free surface of wedge during water entry

圖16 入水時(shí)間0.015 8 s 時(shí)壓力分布對比圖Fig.16 Comparison diagram of pressure distribution when water entry time is 0.015 8 s

圖16 入水時(shí)間0.015 8 s 時(shí)壓力分布對比圖(續(xù))Fig.16 Comparison diagram of pressure distribution when water entry time is 0.015 8 s (continued)

圖17 入水垂向速度隨時(shí)間變化圖Fig.17 Variation diagram of vertical velocity of water inflow with time

4 圓柱體垂直出水破冰過程計(jì)算結(jié)果分析

圓柱體垂直出水破冰數(shù)值模擬有限元模型如2.2 節(jié)所示，模型冰材料選取及模擬工況均與無水環(huán)境下圓柱體垂直貫穿冰層一致.

4.1 冰載荷分析

圓柱體垂直出水破冰和無水環(huán)境下圓柱體垂直貫穿冰層冰載荷時(shí)歷對比曲線如圖18，出水破冰時(shí)圓柱體所受冰載荷整體趨勢與無水破冰時(shí)類似，冰載荷峰值基本一致;出水破冰時(shí)冰層與圓柱體相互作用時(shí)間比無水破冰時(shí)的作用時(shí)間長，原因是在圓柱體與冰層接觸之前的近場逼近過程中，結(jié)構(gòu)體與冰層之間的水由于受到結(jié)構(gòu)物與被約束的冰層的雙重?cái)D壓而預(yù)先產(chǎn)生的一個(gè)高壓力場，該壓力場使船-冰之間產(chǎn)生一個(gè)對冰層存在支撐作用的阻力，即水墊效應(yīng)[32].水墊效應(yīng)作用在冰層，使的冰層彈性變形階段更長，且冰層失效后未塌陷，仍然與圓柱體接觸產(chǎn)生作用力;出水破冰時(shí)圓柱體上浮過程中水會(huì)上升擠壓冰面，高壓力場導(dǎo)致圓柱體尚未與冰層接觸前發(fā)生彎曲變形從而產(chǎn)生撓度，所以有水環(huán)境下圓柱體與冰層的冰載荷要略晚于無水環(huán)境出現(xiàn)(0.08 s)，且有水情況下的冰載荷變化趨勢相對平緩.

圖18 有水和無水環(huán)境下圓柱體上浮破冰冰載荷時(shí)歷曲線Fig.18 Time history curve of ice breaking load on floating cylinder in water and waterless environment

4.2 冰層撓度變化分析

圓柱體垂直出水破冰過程可分為水下上浮階段和出水破冰階段兩部分.為研究圓柱體出水破冰過程中冰層撓度變化趨勢，分別選取冰層上表面中心點(diǎn)A和Y軸四等分點(diǎn)B和X軸四等分點(diǎn)C觀測撓度變化曲線，冰層撓度測量點(diǎn)A，B及C如圖19 所示.

圖19 冰層撓度測量點(diǎn)Fig.19 Measuring point of level ice deflection

圓柱體水下上浮階段未與冰層接觸時(shí)A，B，C點(diǎn)的撓度變化曲線如圖20.可以看出A，B，C，3 點(diǎn)撓度變化趨勢及斜率基本一致，撓度變化曲線均出現(xiàn)負(fù)值，呈現(xiàn)震蕩后趨于穩(wěn)定;A點(diǎn)撓度變化隨圓柱體上浮過程整體大于B和C兩點(diǎn)的撓度變化，說明水下上浮階段圓柱體運(yùn)動(dòng)過程中引起水上升擠壓冰面，導(dǎo)致圓柱體尚未與冰層接觸前發(fā)生變形從而產(chǎn)生撓度變化，且應(yīng)力集中在冰面中心點(diǎn);B和C兩點(diǎn)撓度變化曲線幾近吻合，說明水下上浮階段圓柱體對冰層的影響范圍以冰面中心點(diǎn)向四周擴(kuò)散且逐漸衰減，與冰層放置于水面之上且四周約束的工況設(shè)定一致.

圖20 圓柱體水下上浮階段冰層各點(diǎn)撓度Fig.20 Deflection of each point of level ice during underwater floating of cylinder

由圖20 知，圓柱體出水破冰過程中冰層A點(diǎn)撓度隨上浮時(shí)間變化最為明顯，因此選取A點(diǎn)對比無水破冰和出水破冰時(shí)冰層的撓度變化.

圖21 為無水環(huán)境和有水環(huán)境下圓柱體上浮破冰過程中冰層A點(diǎn)處的撓度變化曲線.由圖21 可知兩條撓度變化曲線整體趨勢一致，隨著圓柱體逐漸上升擠壓冰層，冰層撓度逐漸增大;出水破冰過程冰層撓度變化峰值為1.48 mm，高于無水破冰時(shí)冰層撓度變化峰值0.81 mm;圓柱體突破冰層的時(shí)間也有明顯差異，有水環(huán)境下的冰層突破時(shí)間遠(yuǎn)大于無水環(huán)境下的冰層突破時(shí)間;在2.5～ 2.58 s 間有一段平滑階段，由于結(jié)構(gòu)物與冰層擠壓水形成的高壓力場，使得在圓柱體與冰層接觸之前冰層就存在初始撓度，即圓柱體在2.5～ 2.58 s 間未接觸到冰層.圖21 與圖18 曲線表現(xiàn)一致，均證明了有水環(huán)境下結(jié)構(gòu)物-冰層耦合作用時(shí)存在水墊效應(yīng).

圖21 有水環(huán)境和無水環(huán)境下圓柱體上浮破冰過程中冰層A 點(diǎn)撓度變化Fig.21 Deflection change of point A of level ice in the process of cylinder floating and breaking ice in water environment and waterless environment

4.3 冰層應(yīng)力響應(yīng)及斷裂現(xiàn)象分析

圓柱體出水破冰過程中冰層未發(fā)生斷裂時(shí)von Mises 等效應(yīng)力云圖如圖22.在圓柱體水下上浮階段，圓柱體運(yùn)動(dòng)排開水?dāng)D壓冰層，對冰層產(chǎn)生影響，應(yīng)力主要集中于冰層中心即結(jié)構(gòu)物正上方，但也存在應(yīng)力并不集中的情況，由于結(jié)構(gòu)物處于短暫加速狀態(tài)會(huì)使應(yīng)力值不穩(wěn)定，如圖22(a)到圖22(b)過程;在圓柱體處于勻速上升階段，應(yīng)力完全集中在冰層中心，且冰層中心應(yīng)力隨著結(jié)構(gòu)物運(yùn)動(dòng)時(shí)間越來越大，如圖22(c)到圖22(d)過程.在結(jié)構(gòu)物出水破冰階段，冰層應(yīng)力極值點(diǎn)隨著圓柱體的運(yùn)動(dòng)迅速增大，當(dāng)t=2.58 s 時(shí)結(jié)構(gòu)物接觸冰層，應(yīng)力極值為542.9 kPa，如圖22(e);當(dāng)t=2.94 s 時(shí)結(jié)構(gòu)物即將突破冰層，應(yīng)力極值達(dá)到最大值，冰層突破應(yīng)力為2509 kPa，如圖22(f).圓柱體出水破冰過程中的應(yīng)力響應(yīng)云圖與冰層撓度變化趨勢表現(xiàn)一致，無水環(huán)境下圓柱體突破冰層應(yīng)力極值為2381 kPa，小于有水環(huán)境下冰層的突破應(yīng)力極值.

圖22 冰層未發(fā)生斷裂時(shí)von Mises 等效應(yīng)力云圖Fig.22 von Mises equivalent stress diagram when the level ice does not break

圖22 冰層未發(fā)生斷裂時(shí)von Mises 等效應(yīng)力云圖 (續(xù))Fig.22 von Mises equivalent stress diagram when the level ice does not break (continued)

圖23 為圓柱體出水破冰過程中冰層發(fā)生斷裂時(shí)的von Mises 等效應(yīng)力云圖及破裂現(xiàn)象.圓柱體出水破冰過程中冰層斷裂現(xiàn)象大體上整體上與無水破冰相似，可以分為3 個(gè)階段.第1 階段，冰層受到圓柱體擠壓達(dá)到突破載荷超過冰的抗拉強(qiáng)度，浮冰板生成十字徑向裂紋，如圖23(a)所示，此時(shí)發(fā)生輕微卸載，徑向裂紋周圍的應(yīng)力小于圓柱體的突破應(yīng)力，隨著圓柱體運(yùn)動(dòng)十字徑向裂紋持續(xù)擴(kuò)展，徑向裂紋周圍的應(yīng)力也逐漸增大，如圖23(b)所示;第2 階段，冰層在出現(xiàn)十字徑向裂紋之后，隨著圓柱體繼續(xù)上浮，冰層上萌生更多徑向裂紋形成楔子，此時(shí)冰層仍未到達(dá)失效極限，如圖23(c)所示;第3 階段，隨著圓柱體繼續(xù)上升，達(dá)到冰層的失效極限，由于在冰層的上、下表面已經(jīng)存在大量裂紋，圓柱體突破冰層造成冰層失效，部分碎冰在水浮力的作用下沿圓柱體堆積形成碎冰堆，如圖23(d)所示.

圖23 冰層發(fā)生斷裂時(shí)von Mises 等效應(yīng)力云圖及現(xiàn)象Fig.23 von Mises equivalent stress diagram and phenomenon when level ice breaks

圖23 冰層發(fā)生斷裂時(shí)von Mises 等效應(yīng)力云圖及現(xiàn)象(續(xù))Fig.23 von Mises equivalent stress diagram and phenomenon when level ice breaks (continued)

5 結(jié)論

本文針對結(jié)構(gòu)物與冰層相互作用強(qiáng)非線性問題，應(yīng)用LS-DYNA 有限元軟件建立了基于S-ALE流固耦合方法及罰函數(shù)接觸算法的冰-水-結(jié)構(gòu)物耦合作用數(shù)值模擬方法，研究了圓柱體垂直出水破冰過程，并對比分析了無水環(huán)境及有水環(huán)境下圓柱體垂直貫穿冰層破裂的共同點(diǎn)與區(qū)別，得到如下結(jié)論.

(1)通過冰梁結(jié)構(gòu)4 點(diǎn)彎曲數(shù)值模擬，對比文獻(xiàn)中的試驗(yàn)結(jié)果及數(shù)值模擬現(xiàn)象，驗(yàn)證了彈塑性應(yīng)變率材料的冰力學(xué)特性.進(jìn)行了無水環(huán)境下圓柱體垂直貫穿冰層模型試驗(yàn)及數(shù)值模擬，對比破冰載荷和冰層破裂現(xiàn)象，驗(yàn)證了有限元方法研究結(jié)構(gòu)物-冰層耦合作用問題的準(zhǔn)確性.基于S-ALE 方法模擬了楔形體入水流固耦合過程，驗(yàn)證了流固耦合算法的可行性.

(2)開展了圓柱體垂直出水破冰數(shù)值模擬，通過與無水破冰對比可知:有水環(huán)境下結(jié)構(gòu)物-冰層間作用存在水墊效應(yīng).冰層突破載荷極值大小基本不會(huì)因有水環(huán)境或無水環(huán)境而產(chǎn)生明顯變化.由于水對冰層的支撐作用，有水環(huán)境下的結(jié)構(gòu)物突破冰層的冰載荷持續(xù)時(shí)間顯著長于無水環(huán)境下冰載荷持續(xù)時(shí)間.

(3)描述了圓柱體垂直出水過程中的冰層斷裂現(xiàn)象，與無水破冰相比裂紋萌生到擴(kuò)展現(xiàn)象相似，最終碎冰失效形式不同.無水環(huán)境下碎冰會(huì)迅速塌陷，而有水環(huán)境下的碎冰會(huì)沿圓柱體周圍堆積.有水環(huán)境下冰層彈性變形階段更長，撓度變化大于無水環(huán)境下的撓度變化.

綜上所述，由于在有水、無水環(huán)境下的結(jié)構(gòu)物突破冰層冰載荷極值大小基本一致，研究冰層突破載荷極值及結(jié)構(gòu)物強(qiáng)度相關(guān)問題時(shí)均可在無水環(huán)境下研究提高計(jì)算效率.而由于水墊效應(yīng)的存在有水環(huán)境下結(jié)構(gòu)物垂直出水貫穿冰層破裂過程與無水環(huán)境下結(jié)構(gòu)物垂直貫穿冰層破裂過程存在明顯區(qū)別，需要根據(jù)研究目標(biāo)選擇合適的數(shù)值模擬方法.