修正內聚區長度計算公式在冰I 型裂紋擴展中的應用

倪寶玉，徐瑩，黃其，尤嘉，薛彥卓

1 哈爾濱工程大學 船舶工程學院，黑龍江 哈爾濱 150001

2 中國船舶及海洋工程設計研究院，上海 200011

3 中國船級社 江蘇分社，江蘇 南京 210011

0 引 言

結構物與冰的相互碰撞研究在冰區結構物的設計制造、極地地區資源的開發和利用，以及冰區船舶和結構物的安全航行及作業等領域均受到普遍重視。目前，在結構物?冰碰撞的研究中使用較廣泛的方法包括離散元法[1-2]和有限元法[3]。其中，離散元法用于冰裂紋擴展的研究具有一定的優勢，但對于海冰變形的研究不夠充分；有限元法用于海冰與結構變形的模擬較成熟，但對于冰裂紋的擴展模擬不準確。基于此，有研究者將內聚力單元法（cohesive element method, CEM）應用于冰裂紋的模擬。CEM 法是數值模型中模擬裂紋萌生和擴展的專用有限元法，最早可以追溯到研究者對裂縫尖端Dugdale-Barenblatt（D-B）[4-5]模型的開發和研究。關于內聚力模型的最初思想是由Barenblatt[5]提出的原子牽引分離定律，應用該定律可以避免裂紋尖端處不切實際的應力奇異性。后來，Dugdale[4]使用內聚力模型來描述受到法向應力影響的理想彈塑性材料在裂紋尖端附近的塑性變形。而后，Needleman[6]將D-B 模型中的裂紋區域稱為內聚力區（cohesive zone）。內聚力模型在許多研究者的研究中得以改進并應用于多個領域，與其他裂紋計算模型相比，內聚力模型具有能夠模擬多個裂紋路徑而無需采用裂紋路徑跟蹤的優點。此外，CEM 法無需預先知道裂紋擴展的方向，可通過任意放置內聚力單元來模擬裂紋路徑的傳播。本文將主要研究CEM 法在冰裂紋擴展方面的應用，這對于冰區結構物設計制造有著重要意義。

內聚力模型憑借在裂紋擴展和材料斷裂方面的優勢，近年來，被越來越多地應用于冰力學和結構物與冰碰撞的模擬中。例如，Mulmule[7]和Dempsey[8]采用I 型斷裂模式的冰材料對內聚定律進行了反算；Kuutti 等[9]通過將二維內聚力單元插入三角形冰單元周圍，模擬了垂直鋼板與厚冰塊碰撞的試驗，結果發現密集的有限元模型網格可更好地模擬試驗中冰塊的破壞過程；Lu 等[10]和Wang 等[11]均采用CEM 法模擬了錐形體與層冰的碰撞過程；王峰等[12]、劉路平等[13]和蔣昱妍等[14]則采用CEM 法模擬了平整冰與豎直圓柱體、海洋平臺等海洋結構物的碰撞過程，驗證了CEM 法在處理此類碰撞問題時的可靠性。

目前，盡管CEM 法在研究結構物?冰相互作用方面已得到一定的應用，但關于冰的內聚區長度公式的研究仍十分缺乏。有部分研究者直接采用了其他材料的公式，例如Lu 等[10]采用彈性材料的公式結合冰裂紋試驗，初步估算了內聚區長度。然而，大部分研究者并未關注內聚區長度的計算問題，例如Gürtner 等[15]、Kuutti 等[9]、Wang 等[11]在處理各種結構物?冰碰撞問題中均未提及內聚區長度與網格的關系，而在劃分網格時多采用試錯的方法。事實上，基于其他材料的研究成果在冰這種材料上不一定適用，且以往研究內聚區長度主要集中在低厚度模型[16-17]，例如復合材料的夾層，而對高厚度模型的研究很少。對于冰的斷裂模式而言，冰通常都是高厚度模型，現有的適用于薄板類型的內聚區長度不再適用。

對于上述研究中存在的不足，本文擬針對冰材料的內聚區長度計算和網格劃分問題進行新的適用性研究。首先，重現及推導現有的內聚區長度計算公式，并對公式進行修正；然后，基于有限元法建立雙懸臂梁數值模型進行模擬，以研究單個內聚區長度內的有限元模型網格密度；最后，基于LS-DYNA 軟件將相關結果應用于冰的三點彎曲試驗，以驗證CEM 法應用于模擬冰的I 型裂紋萌生和擴展時的有效性，期望能夠促進CEM法在冰材料力學特性研究領域中的應用。

1 內聚區長度計算公式及原理

CEM 法在提出后主要應用于混凝土結構和復合材料等領域，展現出了良好的裂紋模擬能力，可以形象地模擬材料的張開、滑開、撕開等開裂破壞過程，以及很好地模擬混凝土的碎裂和復合材料的分層失效過程。目前，控制內聚區的常見牽引力?分離關系包括4 種形式：雙線性型[18]、多項式型[19]、梯型[20]及指數型[21]，如圖1 所示。圖中，σ0為裂紋形成的最大應力， δn為損傷起始對應的上下面位移值， δf為“梯形”應力開始線性下降對應的上下面位移值， δ0為裂紋形成時上下面的最大位移值。

圖1 牽引力?分離曲線Fig. 1 Traction-separation curve

圖2 所示為不同材料的內聚力模型方向示意圖。圖中，定義裂紋擴展方向為長度方向、裂紋向外張開方向為厚度方向。由兩圖的對比可見，海冰屬于典型的高厚度模型。

圖2 不同材料內聚力模型方向定義Fig. 2 Defined dircections of cohesion model of different materials

內聚區長度可通過J積分推導計算。J積分是彈塑性斷裂力學中一個與路徑無關的積分[22]，可作為裂紋或缺口頂端應變場的平均度量。本文應用張開破壞的雙懸臂梁模型來簡述裂紋擴展區域長度的推導過程，如圖3 所示。圖中，Lcz為裂紋擴展區長度（也即內聚區長度），L為無初始裂紋的雙懸臂梁長度，Lpc為初始裂紋長度，L0為整個雙懸臂梁長度，h為單個懸臂梁的厚度。雙懸臂梁左端的初始裂紋用于約束裂紋的產生位置和擴展方向，外扭矩M施加于雙懸臂梁左端。

圖3 雙懸臂梁撕裂示意圖Fig. 3 Schematics of tearing of double cantilever beam

J積分的回路積分定義公式如下[22]：

為將所有常系數（含無量綱系數a的項及常數項）統一到一起并用字符代替，使公式具有一般性，式（3）可進一步變形為

式中：c為關于L和h的無量綱系數。將式（4）進一步簡化為

其中，

Hillerborg 等[23]基于D-B 模型得出原始的裂紋擴展區長度與裂紋位移和應力間的關系為

式中，f為裂縫形成時的最大應力，與式（6）中的σ0等價，即f=σ0；G=kδ0f，為斷裂能釋放率，其中k為常數，G適用于I 和II 型斷裂模式，可以理解為圖1 中各類型曲線與坐標系包圍的面積，該面積又可表達為最大應力、最大應變與一個常數的乘積。對比式（6）和式（7），可知文獻[23]的與式（6）中的c2是一致的，同樣也適用于撕裂（內聚力單元的I 型拉伸失效）破壞。為此，式（5）可改寫為

上述相關推導起初源于復合材料的裂紋擴展和計算裂縫區域長度，但因推導過程中并未限制材料的類型及屬性，所以上述描述裂紋的思想也適用于冰的I 型裂紋的萌生和擴展。文獻[23]提及的雙線性δ-σ 曲線適用于混凝土這類脆性或準脆性材料，本文應用CEM 法，在冰的力學模擬中也采用了雙線性δ-σ 曲線。

應用式（8）時，系數c1的選取是一大難點。以往相關研究者主要針對薄板模型，根據材料性質、內聚力模型本構關系等，對c1的取值給出了不同的推薦值，例如0.5[24]，π/8，0.731 或2.92[25]等。但是，這些參數均采用的是簡單的固定值，忽略了c1本身也會隨著模型尺寸變化的特征，且它們只對于低厚度的薄板材料（圖2（a））較為適用，而對于高厚度冰材料（圖2（b））的適用性尚未見諸于相關文獻報道。為此，本文將針對系數c1在冰材料中應用的取值進行研究。

2 應用于冰材料的內聚區長度

本文在給出系數c1的修正公式之前，首先探討內聚區長度Lcz隨模型尺寸的變化規律。這里，將關注兩個值：厚度值h和長厚比L/h。在有限元模型中，可以通過單元的應力或應變是否超過模型中的設定值來判斷單元是否失效[26]。本文內聚力單元的變形位移量是通過追蹤每個時刻單元節點坐標計算得到的?；谖灰品▌t，通過節點坐標的變化來判斷內聚力單元是否處于不可恢復的變形狀態，結合牽引力?分離曲線（圖1），來判定內聚力單元是否失效。為降低誤差，本文盡可能減小了坐標數據輸出的時間間隔，將輸出間隔時間設置為5×10?5s。

構建雙懸臂梁數值模型。雙懸臂整個長度L0=150 mm，單個懸臂梁厚度h=1.55 mm，模型受力端有一個初始裂紋長度Lpc=35 mm 的切口，無初始裂紋的雙懸臂梁長度L=115 mm。將裂紋擴展長度設為內聚區長度Lcz，雙懸臂梁受方向相反、大小相等的拉力F作用。模型中的內聚力單元受力變形如圖4 所示。由圖可見，當加載時間由0 s 增加至0.35 s 時，內聚力單元隨著力的加載逐漸被拉伸變形；當加載時間達到0.5 s 時，內聚力單元因受力達到極限并失效。

圖4 內聚力單元拉伸變形Fig. 4 Tensile deformation of cohesive element

在保證初始裂紋長度Lpc=35 mm 不變，而改變無初始裂紋的雙懸臂梁長度L（115，100 和85 mm）的情況下，繪制了如圖5 所示雙懸臂梁破壞過程中內聚區長度Lcz隨單個懸臂梁厚度h的變化曲線，并進一步繪制出如圖6 所示Lcz隨Ln(L/h)的變化曲線。

圖5 Lcz -h 曲線Fig. 5 Lcz-h curve

圖6 Lcz-Ln(L/h)曲線Fig. 6 Lcz-Ln((L/h) curve

由圖5 可見，整體內聚力區長度Lcz會隨無初始裂紋雙懸臂梁長度L的增加而逐漸增加。在增加單個懸臂梁厚度h的情況下，Lcz相對于h增加的變化率逐漸減小。當h增加到一定數值后，Lcz基本不再發生變化，但對于低厚度條件下的Lcz變化影響較大。

由圖6 可見，長厚比L/h的增加對Lcz的影響也呈現很明顯的下降趨勢。當L/h>30（Ln(30)≈3.40）時，Lcz基本不再發生變化；當L/h減小時，Lcz受L/h變化的影響逐漸敏感；當L/h進一步減小，達到L/h<1.25（Ln（1.25）≈0.223）時，Lcz再次在某個穩定值附近波動，本文將出現該現象的點稱為拐點。根據3 種不同長度L的無初始裂紋雙懸臂梁模型，拐點出現時，L/h的比值均在1.25 附近。由圖5 和圖6 可見，當L/h在1.25～30 的范圍內取值時，Lcz受到L/h值變化的影響比較明顯；當L/h不在該范圍時，其對Lcz的影響比較小。

3 修正內聚區長度計算公式

目前，CEM 法主要應用于較薄，即L/h值偏大[16]的復合材料的力學特性研究中，其內聚區長度公式的修正系數以常數項修正為主，未考慮長厚比L/h值由1.25 趨近1 時內聚區長度不再隨厚度變化的現象。而在實際的模擬中，應考慮存在內聚區長度變化的拐點，即在L/h到達一定數值后，內聚區長度變化不再明顯的現象。本文擬在前人研究的基礎上對系數c1的選取進行修正，將其與L/h值的變化聯系起來，以提高內聚區長度公式的精確性。

根據不同長度L的無初始裂紋雙懸臂梁模型，本文選取L/h值大于1.25 的部分，并假設c1為關于L與h的函數，再根據式（8）對模擬中的Lcz–h散點進行多次數據擬合，以建立c1與L/h的關系式。通過對數據的擬合整理，可得：

代入式（8），可得到新的內聚區長度計算公式為

上式中：x為當前L/h的取值；hmin為模擬的最小模型厚度。當hmin遠小于模型長度時，hmin對計算結果影響很小，對于冰材料目前尚未見相關計算，本文參照文獻[5] 和[16] 中對于復合材料的取值，選擇hmin=1.5 mm。對于冰的I 型裂紋的萌生和擴展， σ0取為冰的抗拉強度，其值的選取參考了冰的抗拉試驗[27]。此外，根據上節所述，拐點出現后內聚區長度基本不變，故該公式僅適用于L/h大于1.25 的情況，對于L/h小于1.25 的情況，均取L/h=1.25 對應的內聚區長度。

下文將驗證系數c1對內聚區計算結果修正的有效性。重新構建一個高厚度的雙懸臂梁數值模型，如圖7 所示，其寬度B為20 mm，長度L分別為70 和210 mm，Lpc仍為35 mm，厚度2h不斷變化。在1～150 的L/h范圍內適當選取一定的散點繪制相應曲線，將模型計算結果與公式修正后的結果進行比較，如圖8 所示。

圖7 高厚度的雙懸臂梁模型Fig. 7 Model of high-thickness double cantilever beam

由圖8 可見，在2 種不同尺寸模型下，修正模型模擬結果與直接模擬結果在L/h的所有取值范圍內都有較好一致性。由分析可知，除驗證公式的有效性外，拐點處（L/h=1.25，即圖8（a）中的虛線位置）的內聚力區長度與雙臂梁模型整個長度也有一定的關系。隨著整個長度和內聚區長度的增加，公式的計算結果與模擬結果間的誤差逐漸變小，這點在L/h取值較小的范圍內尤其明顯。

圖8 模擬結果和修正內聚區長度對比Fig. 8 Comparison of simulation results and modified cohesive zone length

4 修正公式在冰材料上的應用

網格密度決定了計算量大小。在內聚區長度Lcz中內聚力單元的密度越大，裂紋形成和力的曲線越光滑，精確性越高，反之，力的曲線更容易出現鋸齒形波動。在單個內聚區范圍內設置合適的內聚力單元個數可以很好地模擬裂紋形成和力的傳遞過程?；谏衔腖cz計算方法，本節以內聚力單元的I 型拉伸失效為研究對象，研究一個Lcz中應布置的內聚力單元個數，即Lcz=mSc（其中，Sc為內聚力單元尺寸，m為單元數）。

本節狹長雙懸臂梁模型的主尺度（L×B×2h）為150 mm×20 mm×3.1 mm，如圖9 所示。冰參數如下：彈性模量5.72 GPa，剪切模量2.2 GPa，硬化模量4.26 GPa，體積模量5.26 GPa，上下牽引的運動速度為0.004 m/s。網格尺寸依次從2.0 mm 細化到0.5 mm，間隔0.05 mm。模擬在拉力帶動下狹長雙懸臂梁間的內聚力單元的I 型拉伸失效過程，并測試不同網格密度下網格密度對失效拉力和位移曲線的影響，如圖10 所示。

圖9 狹長雙懸臂梁模型主尺度示意圖Fig. 9 Schematics of main dimensions of narrow double cantilever beam model

由圖10 可見，隨著網格的細化，內聚力單元需傳遞的牽引力?時間曲線逐漸趨于光滑，整體波動趨于緩和，數值逐漸向中間靠攏。根據式（10）計算的內聚區長度Lcz=2 mm，結合上述模型數據，實際模擬得到的內聚區長度也在2 mm 附近。

根據圖10 所示不同網格密度下牽引力?時間曲線分布情況，當網格尺寸小于等于0.5 mm 時，曲線的波動基本上可以接受，且0.5 mm 的網格密度所需計算量也可以接受。綜上所述，內聚力單元在Ⅰ型拉伸失效模式下一個內聚區長度中至少存在4 個內聚力單元才能較為精確地描述材料的斷裂過程，并可較好地保證載荷曲線的精確性和穩定性。

圖10 不同網格密度下內聚力單元牽引力?時間曲線Fig. 10 Traction-time curves of cohesive element with different grid densities

在理論上推斷了CEM 法網格長度與模型尺寸關系的精確性并進行數值驗證后，本節將采用CEM 法對冰的三點彎曲試驗進行數值模擬，以呈現冰由變形到破壞的過程。數值模型[28]主尺度（L×B×2h）為70 mm×70 mm×650 mm，支點間距離為600 mm，如圖11 所示。相關參數如表1 所示。數值模型按照實際的試驗模型建立，冰參數的取值根據試驗結果推導得出。

表1 試驗參數Table 1 Experimental parameters

圖11 冰試樣模型主尺度示意圖Fig. 11 Schematics of main dimensions of ice specimen

為保證良好的網格穩定性和計算結果的精確性，根據式（10）計算所需內聚區長度Lcz，按照裂縫發展方向，以及上文中對于內聚區長度公式中厚度方向的定義，確定了三點彎曲試驗模型的長厚比L/h約為0.215，這明顯小于拐點出現時的L/h值1.25，故取拐點處的Lcz值。經計算，得到的拐點處Lcz=15.4 mm，并根據m=4，得到最大內聚力單元網格尺寸為3.85 mm。

為保證網格劃分的取整性，計算中裂紋發展方向的內聚力單元長度取值為3.5 mm，垂直于裂紋發展方向的內聚力單元長度為4.5 mm，這樣可以節省一定的計算量。在上述網格密度條件下，計算步長為3.68×10?7s，計算時長約4 h。

按照上述試驗數據和冰參數構建有限元模型，在可能斷裂的試樣中間部分添加內聚力單元，添加方式如圖12 所示。圖中，紅色為冰單元，灰白色為內聚力單元。冰力學模型僅提供變形和力的傳遞，冰的破壞和分離采用內聚力單元代替，故不考慮冰單元的破壞，而受破壞的僅有內聚力單元。冰力學模型為各向同性的線彈塑性模型，數值模擬過程如圖13 所示，試驗冰試樣斷裂形式如圖14 所示。

圖12 冰試樣截面間內聚力單元Fig. 12 Cohesive element between ice sample sections

圖13 有限元模擬三點彎曲過程Fig. 13 Finite element simulations of three-point bending process

圖14 試驗冰試樣斷裂[28]Fig. 14 Fracture of experimental ice specimen[28]

本節數值模擬所用冰模型是均勻的，內聚力單元牽引力法則選用的是雙線性型，且在冰力學模型建立過程中不考慮冰的原始缺陷，因此，圖13所示冰試樣斷裂面較光順，導致整體模擬過程中冰的破壞較為平滑。根據數值模擬中物體的撓度與壓力變化，繪制出如圖15 所示壓力?撓度曲線。

圖15 壓力?撓度曲線Fig. 15 Pressure-deflection curves

試驗中，冰試樣被破壞前的最大受力值為1 127.9 N，基于CEM 法模擬三點彎曲的極限載荷為1 161 N，斷裂點的誤差為2.9%，與試驗結果較為接近。從圖15 所示壓力?撓度曲線的不同階段來看，CEM 法模擬的曲線趨勢和試驗曲線基本一致。從最終斷裂撓度來看，CEM 法對應的極限撓度值相較于試驗數據有一定的誤差，但較小。此外，由于該數值模擬選擇的冰力學模型為線性模型，且內聚力單元建立過程中未考慮原始缺陷，使得模擬得到的曲線相較于試驗數據得到的曲線更光滑。從兩個曲線的對比可知，數值模擬過程中冰試樣的彎曲變形、破壞過程與試驗的基本一致，且峰值過后的載荷卸載過程很好地體現了冰的脆性和裂紋成型的連貫性?？傊敬文M數據與試驗數據的對比結果很好地驗證了CEM 法及相關網格劃分方式在模擬冰的力學特性方面的有效性。

5 結 語

本文從J積分出發，基于部分假設并結合前人的相關研究結果，重現了內聚區長度計算公式的推導過程，在原有的內聚區長度公式基礎上，增加了關于L/h值的修正函數。通過建立新的高厚度雙懸臂梁數值模型檢驗了內聚區長度公式的精確性，解決了該公式是否適用于冰力學模型的問題。通過構建雙懸臂梁有限元模型進行數值模擬，結果發現，在一個內聚區長度中至少存在4 個內聚力單元才能夠較為精確地描述材料的斷裂過程。最后，將研究結果應用于三點彎曲試驗模型的模擬中，結果顯示斷裂點的極限載荷誤差為2.9%且在合理范圍內，由此驗證了修正后的內聚區長度公式在模擬冰的I 型裂紋萌生和擴展時的有效性。