基于粘聚力單元的平整冰與海洋樁柱結(jié)構(gòu)相互作用機理研究

俞同強，劉 昆，劉俊杰，王自力

（1.江蘇科技大學(xué)船舶與海洋工程學(xué)院，江蘇鎮(zhèn)江 212003；2.中國船舶科學(xué)研究中心，江蘇無錫 214082）

0 引 言

北極地區(qū)是一個冰封已久的寶藏，不僅蘊藏著豐富的石油、天然氣等化石能源，還具有可觀的風(fēng)力、水、生物、土地、旅游等廣義資源。隨著全球氣候變暖的加劇以及技術(shù)的進步，北極資源開發(fā)利用的難度不斷降低，其巨大的經(jīng)濟和人文價值日益顯現(xiàn)。極地海域廣泛存在的漂浮平整冰是極地海洋裝備面臨的主要安全挑戰(zhàn)，其巨大的推動力和產(chǎn)生的冰激振動往往會給海洋裝備造成結(jié)構(gòu)性損害，危及平臺安全，與之相關(guān)的的事故屢見不鮮[1]。樁柱結(jié)構(gòu)作為目前海上平臺設(shè)施基座的基礎(chǔ)結(jié)構(gòu)，具有受波浪載荷作用小、結(jié)構(gòu)連續(xù)性好的特點，得到了廣泛應(yīng)用，其在冰載荷下的結(jié)構(gòu)性能也直接影響平臺整體結(jié)構(gòu)的安全性。因此，開展平整冰與海洋樁柱結(jié)構(gòu)相互作用研究，探究結(jié)構(gòu)的動態(tài)響應(yīng)和海冰的破壞機理，對于評估冰載荷的大小和變化特點、提高結(jié)構(gòu)性能、保障海上平臺的安全性具有十分重要的意義。

海洋結(jié)構(gòu)物與平整冰的相互作用機理一直是國內(nèi)外學(xué)者關(guān)注的主要問題，也是冰區(qū)結(jié)構(gòu)設(shè)計和疲勞壽命評估的重要依據(jù)。早期冰力的預(yù)報方法往往是基于試驗的回歸公式來計算最大靜冰力，很多學(xué)者和各國規(guī)范也都給出了相關(guān)的計算公式[2]。然而，事實上，海冰的運動、擠壓、破碎、重結(jié)晶等活動往往導(dǎo)致其動冰力顯著大于靜冰力，對冰激振動響應(yīng)和疲勞壽命的評估已成為目前海冰研究的主要工作[3-4]。動態(tài)冰載荷問題是一個時域性的問題，在理論方法上進行全面研究目前仍有較大的困難。十九世紀(jì)七十年代，有學(xué)者提出了接觸面冰材料非同時性破壞的特征[6]，即高壓力區(qū)和低壓力區(qū)交替出現(xiàn)，高壓力區(qū)冰材料發(fā)生破碎和脫落。在此基礎(chǔ)上，Jordaan 等[6-7]通過結(jié)構(gòu)與平整冰的壓縮試驗，分析了平整冰在平面擠壓作用下的破壞機理，建立了海冰與結(jié)構(gòu)相互作用面內(nèi)冰載荷高壓區(qū)概率數(shù)學(xué)模型。

現(xiàn)場測試[8-9]和模型試驗[10-11]被證明是研究海洋樁柱結(jié)構(gòu)動冰力的有效方法，但是近年來數(shù)值仿真逐漸成為海洋結(jié)構(gòu)物前期設(shè)計的重要輔助手段，其成本低、效率高、周期短的特點，在不同類型的結(jié)構(gòu)-平整冰作用研究中均得到了廣泛的應(yīng)用[12-13]，同時隨著計算機技術(shù)的迅猛發(fā)展，其準(zhǔn)確性和適用性也在不斷提高。但是，仿真分析中海冰材料的選擇和仿真方法的確定是影響計算結(jié)果準(zhǔn)確性和科學(xué)性的兩個最為重要的問題。粘塑性、粘彈性、彈塑性以及泡沫材料等均被學(xué)者嘗試應(yīng)用于結(jié)構(gòu)-冰碰撞的有限元仿真分析計算中[14-16]，在事故冰載荷的預(yù)報方面有較大的科學(xué)價值。但是傳統(tǒng)有限元方法基于連續(xù)介質(zhì)力學(xué)理論，使用單元強度模擬海冰強度，使用單元失效模擬海冰破碎的機制會導(dǎo)致失效單元的大量刪除，既違背了守恒定理，也難以真實模擬海冰的碎裂、堆積等物理現(xiàn)象。

國內(nèi)外學(xué)者也探索了很多具有創(chuàng)新性的仿真方法，其中以離散元方法、粘聚力單元、粘聚力接觸、擴展有限元（XFEM）、自由網(wǎng)格法及分支裂紋陣模型等（wing crack model）[17-19]為代表的新型仿真方法近年來也得到了一定的嘗試和應(yīng)用，但距離成熟仍有一定的距離?？梢?，目前對于動態(tài)冰載荷的研究還沒有一個達(dá)成共識的分析方法和理論，真正的動冰力計算還處于“初級階段”。

考慮到海冰材料特征的復(fù)雜性與仿真方法的有效性，本文采用基于自定義材料本構(gòu)和粘聚力單元相結(jié)合的仿真方法，將海冰的擠壓和彎曲剪切破壞結(jié)合，建立平整冰-海洋樁柱結(jié)構(gòu)的相互作用數(shù)值仿真模型并進行分析，重點關(guān)注相互作用過程中海冰的變形、失效以及冰載荷的時域特征，并對相關(guān)仿真參數(shù)影響進行一定的探討。

1 基于自定義材料本構(gòu)和粘聚力單元方法的平整冰仿真模型

如圖1 所示，當(dāng)平整冰遭遇直立/傾斜式海洋樁柱阻攔后，其主要發(fā)生的變形可分為三種：其一為平整冰自身發(fā)生擠壓破壞，此時接觸面邊緣破碎脫落，而接觸面中央的海冰在三維應(yīng)力作用下則會顯著強化；其二為平整冰的彎曲破壞，此時平整冰的破壞范圍通常較大，破壞的海冰堆積于表面或沉入平整冰之下，樁柱前端形成敞水區(qū)，冰載荷大幅卸載直至下一次接觸；其三為平整冰大范圍整體開裂，此種破壞往往涉及到平整冰內(nèi)部缺陷和大范圍剪切載荷的傳遞，在本文中不予考慮。

圖1 平整冰的幾種破壞形式Fig.1 Damage modes of level ice

1.1 考慮接觸界面擠壓失效的海冰材料本構(gòu)模型

在基本的海冰物理/力學(xué)參數(shù)測量之外，海冰的單軸和多軸試驗是海冰強度測量最為重要的試驗之一，其中海冰的多軸試驗對于建立海冰的本構(gòu)理論模型尤為重要，國內(nèi)外很多學(xué)者都做了相關(guān)方面的工作[20-21]，其中Derradji-Aoua[22]通過總結(jié)歸納得出的多曲面屈服方程等到了廣泛的接受和使用，該屈服方程可寫為

基于上述屈服方程，建立相應(yīng)彈塑性本構(gòu)模型，本構(gòu)模型在彈性階段，應(yīng)力應(yīng)變關(guān)系滿足各向同性廣義胡克定律

式中：G為剪切彈性模量，G=E/2( 1+μ)；E和μ分別為彈性模量和泊松比；γj為工程切應(yīng)變，γj=2εj(j=xy,yz,zx)；Δε0為靜水應(yīng)變增量，Δσ0為靜水應(yīng)力增量。在塑性階段，由于屈服面的外凸性，塑性應(yīng)變增量為

式中，dγ為塑性一致性參數(shù)增量，采用關(guān)聯(lián)流動法則，總應(yīng)變?yōu)閺椥詰?yīng)變與塑性應(yīng)變之和，可寫作

式中：dσij為應(yīng)力增量；δij為Kronecker符號，其值為0（i≠j）或1（i=j）。

當(dāng)本構(gòu)模型進入塑性階段后，還需要對本構(gòu)模型定義合理的失效準(zhǔn)則。本文在冰材料本構(gòu)模型中僅考慮海冰的擠壓破壞，而不考慮海冰的彎曲或剪切破壞。這是由于當(dāng)多種破壞形式共同存在時，往往會導(dǎo)致單元的提前失效并刪除，其本質(zhì)上是過高估計了海冰材料的抗剪強度，使擠壓失效難以起到作用。本文基于Jordaan 的經(jīng)驗失效方程[6]，如式（7）所示，建立了以率無關(guān)和累積塑性應(yīng)變?yōu)楸碚鞯膭討B(tài)失效準(zhǔn)則來定義單元的失效，

積分率本構(gòu)方程的數(shù)值算法，即應(yīng)力更新算法，是通過本構(gòu)方程計算與反饋應(yīng)力狀態(tài)的主要途徑。本文采用Simo 等[23]提供的Cutting-Plane 算法來求解相關(guān)問題，首先計算測試應(yīng)力σ和屈服判定，若屈服方程fyield( )p,J2＜0，則單元處于彈性階段，測試應(yīng)力σ即為真實應(yīng)力，直接返回分應(yīng)力值；否則，單元進入塑性階段，進行塑性階段內(nèi)的迭代步計算，計算塑性參數(shù)Δγ并更新塑性迭代步內(nèi)的應(yīng)力值：

大型商用有限元軟件ABAQUS提供了完整的用戶子程序接口，依據(jù)上述分析選擇VUMAT子程序編寫海冰材料本構(gòu)子程序，建立單元測試模型，單元底部約束，頂部施加恒定的壓縮速度，單元四周施加恒定的圍壓，如圖2所示。

圖2 單元測試模型及載荷Fig.2 Single unit test model and load application

（1）輸出J2-p曲線，如圖3所示，可以看到，無論圍壓為0 MPa還是圍壓為0～100 MPa，在到達(dá)屈服點后（a、b 點），輸出曲線和理論曲線完全重合，頂點坐標(biāo)與理論解一致，有圍壓存在時，屈服點明顯后移，可見屈服階段的迭代計算收斂性較好。

圖3 屈服方程理論曲線及輸出值Fig.3 Theoretical and output value of yield equation

（2）圍壓為0 時，輸出εf-p曲線，如圖4 所示，可以看到，在單元到達(dá)失效點（c 點）之前，通過單元內(nèi)部計算得到的失效應(yīng)變曲線與理論失效應(yīng)變曲線完全重合，而單元自身塑性應(yīng)變在單元測試階段通過兩種方法求得，即通過塑性參數(shù)計算，以及通過總應(yīng)變減去彈性應(yīng)變的方法計算（εp=ε-εe），也取得了較好的一致性，可見在塑性階段的失效計算是準(zhǔn)確的。

圖4 單元失效應(yīng)變理論曲線及輸出值Fig.4 Theoretical and output of failure strain

（3）考察階梯變化圍壓下的失效應(yīng)變值，其中圍壓分別為10～80 MPa，每隔10 MPa 為一個圍壓工況，所得結(jié)果如圖5所示，可以看到，實際輸出值與所設(shè)圍壓值有一定的偏移，表明靜水應(yīng)力和圍壓共同作用在實際的屈服和失效準(zhǔn)則判定中。無論圍壓條件如何，失效時總應(yīng)變均稍大于塑性應(yīng)變，失效應(yīng)變隨圍壓的增大而減小，可見圍壓對于失效應(yīng)變的影響是顯著的。

圖5 不同圍壓下的塑性失效應(yīng)變值Fig.5 Plastic failure strain under different confining pressures

1.2 基于粘聚力單元的平整冰模型

粘聚力模型理論最早于上世紀(jì)60 年代提出[24]，目前主要被應(yīng)用于混凝土及復(fù)合材料的研究中。本質(zhì)上，粘聚力是原子或者分子之間的相互作用力，這種力客觀存在，通過實際測量，定義合適的參數(shù)，可以有效分析界面損傷失效過程。在粘聚力區(qū)內(nèi)，通常將裂紋面上各向應(yīng)力定義為裂紋面上位移之間的關(guān)系，稱為牽引力法則（Traction Separation Law，TSL），可寫作

開裂過程中釋放的能量定義為斷裂能，可寫作

牽引力法則分為三個階段，第一階段為初始損傷之前的彈性階段，此時應(yīng)力應(yīng)變線性增長；第二階段為損傷之后的軟化階段，此時單元剛度隨損傷的發(fā)展逐漸降低；第三階段為到達(dá)最大分離值時，單元刪除?？梢?，粘聚力模型有兩個重要材料參數(shù)，分別為最大粘聚力值與斷裂能，這兩個參數(shù)均為材料常數(shù)，其大小取決于材料的破壞強度及粘聚力區(qū)域的尺寸。

在幾何方面，理想粘聚力模型一般采用厚度為0的粘聚力單元來表示，但也可采用厚度值極小的單元來替代，如圖6所示。有限元軟件ABAQUS 提供了粘聚力單元功能，定義雙線性軟化粘聚力模型本構(gòu)，其牽引力法則曲線如圖7 所示，該本構(gòu)模型不可逆，可以看作是塑性區(qū)裂紋積累所導(dǎo)致的材料軟化現(xiàn)象。此外粘聚單元不代表任何物理材料，僅描述了發(fā)生斷裂時的粘聚力，在仿真中即使單元厚度不為0，其質(zhì)量也應(yīng)當(dāng)設(shè)為0，因此粘聚力單元的侵蝕不違反守恒定律。

圖6 粘聚力單元圖Fig.6 Cohesive element

圖7 雙線性軟化粘聚力模型本構(gòu)Fig.7 Constitutive model of bilinear softening cohesive element

參考相關(guān)文獻(xiàn)資料[25]，確定輸入?yún)?shù)開展基本模型測試驗證，如表1所示。

表1 粘聚力單元材料參數(shù)設(shè)置Tab.1 Material parameters of cohesive element

建立基本三單元測試模型，其中中間是厚度為0 的粘聚力單元，分別施加恒定的壓縮與拉伸速度，考察單元失效形式，如圖8 所示?？梢钥吹剑趬嚎s載荷下，隨著擠壓作用的增強，冰體單元應(yīng)力不斷上升，冰體單元達(dá)到失效準(zhǔn)則后單元失效，此時單元刪除，單元應(yīng)力為10 MPa左右。當(dāng)載荷為拉伸載荷時，粘聚力單元達(dá)到失效準(zhǔn)則后刪除，有限元單元保留，單元應(yīng)力約為2 MPa左右，與設(shè)置的參數(shù)吻合，可見應(yīng)力狀態(tài)是合理的。

粘聚力單元的使用本質(zhì)上是一種唯像的方法，理論上粘聚力單元的分布形式很大程度上限制了裂紋的擴展路徑。Turon 等[26]指出，在粘聚區(qū)任意方向上至少需要三個網(wǎng)格才能滿足預(yù)測裂紋擴展的要求，對于各向同性多年冰以及橫觀各向同性一年冰，其結(jié)晶方向由表層向上進行，因此采用正交分布的粘聚力單元是合理的。同時，冰單元的應(yīng)力狀態(tài)也和網(wǎng)格尺寸有直接聯(lián)系，參照文獻(xiàn)[27]的方法進行網(wǎng)格敏感性研究，將1 m×1 m×1 m的立方體劃分為不同的網(wǎng)格大小并嵌入粘聚力單元，施加相同的壓縮載荷，觀察冰單元應(yīng)力狀態(tài)的變化。結(jié)果表明，單元尺寸越小，單元的應(yīng)力越大，當(dāng)單元尺寸小于0.14 時，單元應(yīng)力有明顯的收斂值（見圖9 所示）。在實際計算中，還應(yīng)根據(jù)計算時間的成本進行綜合考慮和試算，確定合理的網(wǎng)格大小。

圖9 網(wǎng)格收斂性驗證Fig.9 Mesh convergence verification

2 平整冰與海洋樁柱結(jié)構(gòu)相互作用機理研究

2.1 試驗及仿真模型

基于上述自定義本構(gòu)和粘聚力單元相結(jié)合的方法，建立平整冰與樁柱結(jié)構(gòu)相互作用的仿真分析模型，仿真模擬工況參考武海斌等[28]在天津大學(xué)冰水池開展的單樁風(fēng)機基礎(chǔ)與冰層作用的模型試驗，針對3 MW 單樁風(fēng)機的試驗結(jié)果進行仿真驗證。模型試驗所模擬的冰況為冰厚40 cm，擠壓強度2.06 MPa，冰速范圍0.05～1.2 m/s，縮尺比λ=20。由于粘聚力單元中有關(guān)模型冰的參數(shù)，尤其是斷裂能等難以獲得，而實尺度下有國外相關(guān)的測量結(jié)果，因此本文建立實尺度工況下的數(shù)值模型以保證輸入?yún)?shù)的科學(xué)性。綜合考慮計算時間等因素，實尺度下仿真模型的相關(guān)參數(shù)如表2所示。仿真分析模型如圖10 所示，其中樁柱結(jié)構(gòu)采用殼單元建模（S4R），網(wǎng)格大小為0.1 m，材料為剛體，不考慮其變形，且約束除運動方向外的其他自由度。平整冰采用實體六面體單元建模（C3D8R），網(wǎng)格大小為0.1 m，材料屬性為用戶自定義材料，基本材料參數(shù)如表3 所示，平整冰除接觸面外其他三個面邊緣節(jié)點約束所有自由度。粘聚力單元為0厚度體單元（COH3D8），材料參數(shù)仍如表1所列，接觸摩擦系數(shù)為0.15，粘聚力單元粘性系數(shù)設(shè)置為0.0001，以增強收斂性，采用36核主頻2.24 GHz工作機計算200 s耗時約320小時左右。

表2 實尺度下仿真模型相關(guān)參數(shù)Tab.2 Parameters of simulation model in full scale

表3 平整冰單元基本材料參數(shù)Tab.3 Material parameters of leveling ice

圖10 仿真分析模型與網(wǎng)格特征Fig.10 Simulation analysis model and mesh characteristics

2.2 仿真結(jié)果對比與分析

針對3 MW 風(fēng)機樁柱在平均水位破冰速度為0.05 m/s工況進行了仿真對比研究，該工況對應(yīng)于低速工況。

當(dāng)實際相對速度為0.05 m/s，即試驗相對速度為0.011 m/s時，圖11給出了仿真和試驗中海冰破壞的對比圖?？梢钥吹?，海冰的局部擠壓和小范圍壓屈破壞在仿真和試驗中均有體現(xiàn)，試驗中小范圍壓屈以環(huán)狀裂紋的形式擴展，但在仿真中以規(guī)整的單元斷裂為特征，這主要受到有限元網(wǎng)格劃分形式的影響。圖12給出了仿真和試驗中的時程接觸力圖，可以看到仿真和試驗中的接觸力曲線總體特征基本吻合，峰值力均在400 N 左右。兩者的接觸力平均值分別為86 N和91 N，仿真值略小于試驗值。在小范圍壓屈階段（圖中A 段）接觸力呈現(xiàn)整體加載—快速波動—整體卸載的變化軌跡，并在其后有一個較長時間的準(zhǔn)空載階段，可見在發(fā)生小范圍壓屈時，存在塊狀海冰非同時破壞的現(xiàn)象。在局部擠壓階段，接觸力總體呈現(xiàn)震蕩上升到達(dá)峰值力后震蕩下降的趨勢，總體來說，小范圍壓屈破壞所占比例相對較小，平整冰主要發(fā)生的還是擠壓破壞，這點在仿真和試驗中是一致的。但是，由于平整冰碎裂和破壞的隨機性，時程接觸力的局部特征難以吻合，試驗中除初始階段外沒有整體屈曲破壞發(fā)生，但是在仿真中這種破壞形式在整個作用過程中均有出現(xiàn)（圖中B段），且試驗中測量值的震蕩幅度較大，這在仿真后半段時間內(nèi)才有一定的體現(xiàn)。

圖11 速度為0.011 m/s時平整冰破碎情況Fig.11 Breaking condition of level ice at a velocity of 0.011 m/s

圖12 速度為0.011 m/s時接觸力時程曲線Fig.12 Time history of contact force at a velocity of 0.011 m/s

2.3 參數(shù)敏感性

針對3 MW 風(fēng)機樁柱在平均水位破冰速度為0.6 m/s 和1.2 m/s（即試驗相對速度為0.134 m/s 和0.268 m/s）兩個工況進行仿真對比研究，對應(yīng)于中高冰速工況。

圖13 和圖14 分別為兩種工況下仿真和試驗中平整冰破壞的對比圖和接觸力時程曲線?？梢钥吹较鄬λ俣容^大時，平整冰局部擠壓和小范圍壓屈破壞兩種破壞形式同時存在的特征并沒有發(fā)生變化，單元最大應(yīng)力為8～10 MPa 左右，應(yīng)力水平隨速度的提升有一定的增強。此時試驗中平整冰屈曲破壞的范圍大幅減小，同時有小尺寸碎塊產(chǎn)生，這種趨勢隨著速度的增大而越發(fā)明顯；仿真中則表現(xiàn)為平整冰的破壞更加碎片化，碎片的飛濺也更加嚴(yán)重，海冰的彎曲變形更為明顯，但海冰的屈曲破壞仍然存在并占據(jù)一定的比例。從時程接觸力的對比來看，兩種形式破壞的接觸力特征區(qū)別已不明顯，除初始階段外接觸力曲線沒有顯著的卸載段出現(xiàn)，總體上接觸力表現(xiàn)為在一定幅值內(nèi)反復(fù)震蕩的特征，伴隨有規(guī)律性的個別較大峰值的出現(xiàn)。試驗和仿真結(jié)果都表明，速度的變化不是影響峰值接觸力大小的主要原因，無論是仿真還是試驗，速度增大24倍，但是峰值力的大小沒有顯著的變化，均維持在300～400 N 左右，但是速度的變化顯著改變了接觸力的平均值和時域變化特征，速度越大，平均力越大，同時規(guī)律性的振蕩也越明顯，接觸力振蕩的頻率以及峰值碰撞力出現(xiàn)的頻率均顯著提高，其往往是引起冰激共振、造成結(jié)構(gòu)損壞的主要原因。

圖13 速度為0.134 m/s和0.268 m/s時平整冰破碎情況Fig.13 Breaking condition of level ice at velocities of 0.134 m/s and 0.268 m/s

圖14 速度為0.134 m/s和0.268 m/s時接觸力時程曲線Fig.14 Time histories of contact forces at velocitise of 0.134 m/s and 0.268 m/s

標(biāo)記樁柱運動過程中失效的平整冰和粘聚力單元，記錄樁柱位移分別為2 m、4 m、6 m、8 m時各類型單元的失效數(shù)量和粘聚力失效單元所占的百分比，如圖15 所示，其中粘聚力單元和冰單元原始數(shù)量分別為131 120 和48 000 個?？梢钥吹?，失效單元的總數(shù)隨樁柱運動距離和速度的增加而呈近似線性增加的趨勢。從各類型單元失效占比來看，隨著運動的進行，粘聚力失效單元占比均表現(xiàn)為逐漸下降的趨勢，同時隨速度的增大而逐漸減小。這主要是由于速度較大時，冰單元往往難以形成穩(wěn)定的圍壓狀態(tài)，導(dǎo)致其更易失效。不同工況下，粘聚力失效單元數(shù)量占比均在95%以上，可見由粘聚力單元引起的彎曲剪切失效仍是主要的失效模式。同時需要指出的是，粘聚力單元的基數(shù)約三倍于冰單元數(shù)量，這也是造成粘聚力單元失效數(shù)量占比較大的一個重要原因。

圖15 失效單元數(shù)量與粘聚力失效單元占比Fig.15 Number of failure elements and proportion of cohesive failure elements

考察圓臺樁柱傾角對于平整冰失效特性的影響，在原型樁柱尺寸的基礎(chǔ)上，改變圓臺上底直徑至200 mm、150 mm、100 mm、50 mm，下底直徑不變，即傾角分別為86°、84°、82°、80°，平整冰的破碎情況及接觸力時程曲線如圖16 和圖17 所示，失效單元數(shù)量及計算得到的粘聚力失效單元占比如圖18 所示?？梢钥吹剑捎诮嵌茸兓^小，不同傾角下平整冰的破碎特征沒有顯著的差異，傾角較大時，整體破碎的現(xiàn)象較為明顯，冰單元應(yīng)力大小與典型工況基本相同，均維持在8 MPa左右。時程接觸力非線性、多峰值的特征沒有發(fā)生改變，總體上峰值接觸力隨傾角的增大而呈非線性減小的趨勢，第一個峰值力大小分別為343.73 N、258.71 N、207.60 N和167.34 N。穩(wěn)定破冰后前三種工況接觸力的峰值差異不明顯，傾角為80°時峰值接觸力明顯減小，同時，接觸力波動頻率和峰值接觸力出現(xiàn)的頻率也基本相同。平均接觸力隨傾角的增大變化明顯，平均力從傾角為86°時的76 N下降到傾角為80°時的54 N左右。從失效單元數(shù)量和比例上來看，總體上，失效單元數(shù)量隨傾角的增大而逐漸減小，在單個工況中失效單元數(shù)量隨移動距離呈現(xiàn)線性特征，粘聚力單元失效比例不斷降低，總體占比在85%～90%左右，由于各工況傾角僅有2°的變化值，可見粘聚力單元的失效對于傾角的變化是十分敏感的。

圖16 不同傾角下平整冰破碎情況Fig.16 Breaking condition of level ice at different inclination angles

圖17 不同傾角下接觸力時程曲線Fig.17 Time history of contact force at different inclination angles

圖18 不同傾角下失效單元數(shù)量與粘聚力失效單元占比Fig.18 Number of failure elements and proportion of cohesive failure elements at different inclination angles

考察圓臺樁柱接觸部分直徑對于平整冰失效特性的影響，在原型樁柱尺寸的基礎(chǔ)上，改變圓臺下底直徑至250 mm、200 mm、150 mm、100 mm，上底直徑按照比例改變，此時傾角不變，為87°，平整冰的破碎情況和接觸力時程曲線如圖19和圖20所示，失效單元數(shù)量及計算得到的粘聚力失效單元占比如圖21 所示。可以看到，平整冰的破壞形式隨樁柱直徑的減小顯著不同，但是整體上屈曲失效和破碎失效同時存在的特征沒有發(fā)生改變，單元應(yīng)力有小幅減小，下底直徑為100 mm 時單元最大應(yīng)力減小為8.26 MPa左右，但是總體應(yīng)力水平仍在8 MPa左右。在時程接觸力上，第一峰值力分別為376.12 N、352 N、342.92 N和246.0 1N，呈非線性下降的趨勢，這種趨勢在第二個峰值力體現(xiàn)得更為明顯，第二峰值力在下底直徑為250 mm 時為327.21 N，而在下底直徑為200 mm時大幅降低至202.66 N。接觸力的平均值也呈非線性減小的趨勢，四種工況下接觸力平均值分別為76.43 N、59.39 N、51.91 N、40.57 N。接觸力峰值出現(xiàn)的頻率隨直徑的減小有一定程度的提高，這主要是由于直徑較小時，圓柱面弧度增大，冰單元更易向兩側(cè)擠出，增加了撞擊幾率。由圖21 可以看到，隨著樁柱直徑的減小，無論是冰體單元還是粘聚力單元，失效單元的數(shù)量均顯著減少，粘聚力失效單元占比也逐漸降低，但總體占比仍在85%～90%左右，可見粘聚力單元的失效在平整冰整體失效破碎中起著十分重要的作用。

圖19 不同直徑下平整冰破碎情況Fig.19 Breaking condition of level ice under different diameters

圖20 不同直徑下接觸力時程曲線Fig.20 Time histories of contact forces under different diameters

圖21 不同直徑下失效單元數(shù)量與粘聚力失效單元占比Fig.21 Number of failure elements and proportion of cohesive failure elements under different diameters

綜上可見，樁柱的傾角和直徑顯著影響平整冰的破壞和失效特征，但是屈曲破壞和擠壓破碎同時存在，且非同時性破壞特征明顯的情況沒有改變。應(yīng)力水平在不同工況下有規(guī)律性地變化，總體在8～10 MPa左右，可見在冰體受到擠壓時考慮圍壓的影響是十分必要的。失效單元隨速度、直徑的增加和傾角的減小而增加，粘聚力和平整冰單元同時存在失效，無論何種工況，粘聚力失效單元比例隨樁柱的運動距離逐漸降低，但是粘聚力失效單元占主要多數(shù)的情況沒有發(fā)生顯著變化，總的比例均維持在85%以上。

3 結(jié) 論

本文基于自定義材料本構(gòu)和粘聚力單元相結(jié)合的方法，建立了平整冰與海洋樁柱相互作用的數(shù)值仿真模型，深入分析了平整冰在與樁柱作用過程中平整冰的破壞特點和時程接觸力的變化特征，探究了速度、傾角和直徑的影響規(guī)律，得到如下結(jié)論：

（1）通過自定義材料本構(gòu)和粘聚力單元相結(jié)合的方法，將平整冰的擠壓破壞和彎曲破壞分開對待，解決了傳統(tǒng)有限元軟件難以模擬碎冰的問題。通過單元驗證工作，驗證了子程序的執(zhí)行正確性，同時通過簡單的粘聚力測試模型，闡述了粘聚力的作用機理和網(wǎng)格的收斂性特征。

（2）仿真結(jié)果表明，本文所提出的仿真方法能夠有效地模擬平整冰與樁柱相互作用時平整冰的破壞特征及接觸力的變化情況，尤其是局部擠壓和小范圍壓屈破壞的破壞特征在仿真中均得到了合理的反映。

（3）不同工況下平整冰破碎的特征基本相同，屈曲破壞和擠壓破碎同時存在，且非同時性特征明顯。粘聚力和平整冰單元同時存在失效，失效單元隨速度、直徑的增加和傾角的減小而增加。無論是仿真還是所參考的試驗結(jié)果均表明，速度的變化對于接觸力峰值的大小沒有明顯影響，但會顯著改變接觸力的平均值和時域變化特征，主要表現(xiàn)為規(guī)律性的振蕩和峰值，其往往會引起冰激共振，造成結(jié)構(gòu)疲勞損傷，帶來安全風(fēng)險。

需要指出的是，由于海冰的擠壓破碎現(xiàn)象和粘聚力客觀存在，本文提出的仿真方法具有較強的實用性和科學(xué)性，同時考慮到一年生平整冰表現(xiàn)為各向異性，通過給橫向和豎向粘聚力單元施加不同的材料屬性能夠模擬這一特征，結(jié)合自定義材料本構(gòu)的可拓展性優(yōu)勢，這一仿真方法在未來的科學(xué)研究中也具有較廣的應(yīng)用前景。