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

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

0 引 言

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

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

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

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

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

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

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

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

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

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

基于上述屈服方程，建立相應彈塑性本構模型，本構模型在彈性階段，應力應變關系滿足各向同性廣義胡克定律

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

式中，dγ為塑性一致性參數增量，采用關聯流動法則，總應變?yōu)閺椥詰兣c塑性應變之和，可寫作

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

參考相關文獻資料[25]，確定輸入參數開展基本模型測試驗證，如表1所示。

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

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

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

圖9 網格收斂性驗證Fig.9 Mesh convergence verification

2 平整冰與海洋樁柱結構相互作用機理研究

2.1 試驗及仿真模型

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

表2 實尺度下仿真模型相關參數Tab.2 Parameters of simulation model in full scale

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

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

2.2 仿真結果對比與分析

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

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

圖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 參數敏感性

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

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

圖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

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

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

考察圓臺樁柱傾角對于平整冰失效特性的影響，在原型樁柱尺寸的基礎上，改變圓臺上底直徑至200 mm、150 mm、100 mm、50 mm，下底直徑不變，即傾角分別為86°、84°、82°、80°，平整冰的破碎情況及接觸力時程曲線如圖16 和圖17 所示，失效單元數量及計算得到的粘聚力失效單元占比如圖18 所示?？梢钥吹剑捎诮嵌茸兓^小，不同傾角下平整冰的破碎特征沒有顯著的差異，傾角較大時，整體破碎的現象較為明顯，冰單元應力大小與典型工況基本相同，均維持在8 MPa左右。時程接觸力非線性、多峰值的特征沒有發(fā)生改變，總體上峰值接觸力隨傾角的增大而呈非線性減小的趨勢，第一個峰值力大小分別為343.73 N、258.71 N、207.60 N和167.34 N。穩(wěn)定破冰后前三種工況接觸力的峰值差異不明顯，傾角為80°時峰值接觸力明顯減小，同時，接觸力波動頻率和峰值接觸力出現的頻率也基本相同。平均接觸力隨傾角的增大變化明顯，平均力從傾角為86°時的76 N下降到傾角為80°時的54 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 不同傾角下失效單元數量與粘聚力失效單元占比Fig.18 Number of failure elements and proportion of cohesive failure elements at different inclination angles

考察圓臺樁柱接觸部分直徑對于平整冰失效特性的影響，在原型樁柱尺寸的基礎上，改變圓臺下底直徑至250 mm、200 mm、150 mm、100 mm，上底直徑按照比例改變，此時傾角不變，為87°，平整冰的破碎情況和接觸力時程曲線如圖19和圖20所示，失效單元數量及計算得到的粘聚力失效單元占比如圖21 所示?？梢钥吹?，平整冰的破壞形式隨樁柱直徑的減小顯著不同，但是整體上屈曲失效和破碎失效同時存在的特征沒有發(fā)生改變，單元應力有小幅減小，下底直徑為100 mm 時單元最大應力減小為8.26 MPa左右，但是總體應力水平仍在8 MPa左右。在時程接觸力上，第一峰值力分別為376.12 N、352 N、342.92 N和246.0 1N，呈非線性下降的趨勢，這種趨勢在第二個峰值力體現得更為明顯，第二峰值力在下底直徑為250 mm 時為327.21 N，而在下底直徑為200 mm時大幅降低至202.66 N。接觸力的平均值也呈非線性減小的趨勢，四種工況下接觸力平均值分別為76.43 N、59.39 N、51.91 N、40.57 N。接觸力峰值出現的頻率隨直徑的減小有一定程度的提高，這主要是由于直徑較小時，圓柱面弧度增大，冰單元更易向兩側擠出，增加了撞擊幾率。由圖21 可以看到，隨著樁柱直徑的減小，無論是冰體單元還是粘聚力單元，失效單元的數量均顯著減少，粘聚力失效單元占比也逐漸降低，但總體占比仍在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 不同直徑下失效單元數量與粘聚力失效單元占比Fig.21 Number of failure elements and proportion of cohesive failure elements under different diameters

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

3 結 論

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

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

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

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

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