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

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

引言

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

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

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

針對冰與結(jié)構(gòu)物耦合作用有限元數(shù)值分析方法，國際上已有大量的研究探索.Ehlers 和Kujala[8]在LS-DYNA 有限元軟件中使用彈塑性應(yīng)變率材料模擬冰的彎曲行為，并與4 點彎曲試驗結(jié)果進行對比.Polach 和Franzvon[9]應(yīng)用LS-DYNA 有限元軟件采用MAT_DAMAGE_3 材料模擬冰的力學(xué)特性，并基于損傷力學(xué)對數(shù)學(xué)模型進行分析，將模擬結(jié)果與冰池試驗結(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 有限元軟件模擬了極地船和浮冰的耦合作用過程，實現(xiàn)了冰區(qū)船舶在碎冰區(qū)域的航行阻力預(yù)報.Wang 等[14]應(yīng)用FEM方法與DEM-CFD 方法模擬了冰區(qū)船舶在冰緣區(qū)的航行過程，并將數(shù)值模擬結(jié)果與模型試驗對比，分析了碎冰載荷的主要影響因素.黃焱等[15]基于FEM方法構(gòu)建并分析了固定約束及摩擦約束邊界條件下的冰溫度膨脹力，發(fā)現(xiàn)固定約束會造成高度的應(yīng)力集中.任奕舟和鄒早建[16]、王峰[17]基于FEM 方法，通過粘聚力單元和彈塑性冰本構(gòu)模型構(gòu)建船舶與海冰碰撞場景的數(shù)值仿真，驗證了粘聚元方法模擬海冰的準確性與可行性.Chen 等[18]應(yīng)用FEM-SPH 方法建立了錐體和破冰船與平整冰相互作用的數(shù)值模型，預(yù)測了破冰阻力及碎冰堆積過程.

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

1 理論基礎(chǔ)

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

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

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

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

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

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

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

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

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

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

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

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

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

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

1.3 冰本構(gòu)模型

冰的特性取決于溫度、應(yīng)變率和載荷方向，冰在受力變形破壞過程中，受加載速率的影響，會呈現(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)變率對壓縮強度的影響.假設(shè)彎曲強度與應(yīng)變率無關(guān)，拉伸和壓縮下的屈服應(yīng)力(σy)與有效塑性應(yīng)變(εp)的曲線均由以下方程推導(dǎo)得出

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

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

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

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

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

2 數(shù)值計算模型

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

無水環(huán)境下圓柱體垂直貫穿冰層數(shù)值模型包含以下兩個部分:冰層和圓柱結(jié)構(gòu)體，均采用拉格朗日算法進行描述.模型的幾何尺寸見下圖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

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

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

圓柱體垂直出水破冰數(shù)值模型設(shè)置與文2.1 節(jié)開展的圓柱體垂直貫穿冰層試驗工況保持一致.包含4 部分:空氣域、水域、圓柱結(jié)構(gòu)體和冰層，其中，空氣域和水域采用歐拉算法進行描述，圓柱結(jié)構(gòu)體和冰層仍采用拉格朗日算法進行描述.空氣域和水域的幾何尺寸及網(wǎng)格劃分如圖2，空氣域和水域的采用漸進式網(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 來實現(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)鍵字實現(xiàn)流體區(qū)域靜水壓力分層，使用*CONSTRAINED_AGRANGE_IN_SOLID 關(guān)鍵字實現(xiàn)冰層、圓柱體和流體區(qū)域間的流固耦合作用.

3 計算方法驗證

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

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

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

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

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

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

3.2 圓柱體垂直貫穿冰層試驗

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

圖6 試驗用主要儀器Fig.6 Main instruments for test

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

4.1 冰載荷分析

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

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

4.2 冰層撓度變化分析

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

圖19 冰層撓度測量點Fig.19 Measuring point of level ice deflection

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

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

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

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

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

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

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

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

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

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

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

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

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

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