俞同強(qiáng)，劉 昆，劉俊杰，王自力

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

0 引 言

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

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

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

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

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

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

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

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

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

在基本的海冰物理/力學(xué)參數(shù)測(cè)量之外，海冰的單軸和多軸試驗(yàn)是海冰強(qiáng)度測(cè)量最為重要的試驗(yàn)之一，其中海冰的多軸試驗(yàn)對(duì)于建立海冰的本構(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)變?cè)隽?，Δ?為靜水應(yīng)力增量。在塑性階段，由于屈服面的外凸性，塑性應(yīng)變?cè)隽繛?/p>

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

牽引力法則分為三個(gè)階段，第一階段為初始損傷之前的彈性階段，此時(shí)應(yīng)力應(yīng)變線性增長(zhǎng)；第二階段為損傷之后的軟化階段，此時(shí)單元?jiǎng)偠入S損傷的發(fā)展逐漸降低；第三階段為到達(dá)最大分離值時(shí)，單元?jiǎng)h除?？梢姡尘哿δＰ陀袃蓚€(gè)重要材料參數(shù)，分別為最大粘聚力值與斷裂能，這兩個(gè)參數(shù)均為材料常數(shù)，其大小取決于材料的破壞強(qiáng)度及粘聚力區(qū)域的尺寸。

在幾何方面，理想粘聚力模型一般采用厚度為0的粘聚力單元來表示，但也可采用厚度值極小的單元來替代，如圖6所示。有限元軟件ABAQUS 提供了粘聚力單元功能，定義雙線性軟化粘聚力模型本構(gòu)，其牽引力法則曲線如圖7 所示，該本構(gòu)模型不可逆，可以看作是塑性區(qū)裂紋積累所導(dǎo)致的材料軟化現(xiàn)象。此外粘聚單元不代表任何物理材料，僅描述了發(fā)生斷裂時(shí)的粘聚力，在仿真中即使單元厚度不為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ù)開展基本模型測(cè)試驗(yàn)證，如表1所示。

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

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

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

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

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

2.1 試驗(yàn)及仿真模型

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

表2 實(shí)尺度下仿真模型相關(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é)果對(duì)比與分析

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

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

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

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

2.3 參數(shù)敏感性

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

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

圖13 速度為0.134 m/s和0.268 m/s時(shí)平整冰破碎情況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時(shí)接觸力時(shí)程曲線Fig.14 Time histories of contact forces at velocitise of 0.134 m/s and 0.268 m/s

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

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

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

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

圖17 不同傾角下接觸力時(shí)程曲線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

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

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

圖20 不同直徑下接觸力時(shí)程曲線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

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

3 結(jié) 論

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

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

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

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

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