盛澤強， 王毓杰， 張振南

(上海交通大學 船舶海洋與建筑工程學院，上海 200240)

1 引 言

脆性材料斷裂模擬一直深受關注，目前大部分模擬方法依賴于經(jīng)典連續(xù)介質斷裂力學理論，但其具有一定局限性。材料發(fā)生斷裂破壞，由連續(xù)介質轉變成非連續(xù)介質，原本構關系不再適用；傳統(tǒng)斷裂理論在斷裂模擬時需要考慮外部斷裂準則和網(wǎng)格重劃分等棘手問題。粘結單元法[1]CFEM(Cohesive finite element method)為解決這些問題提供了一種有效方法。該方法在連續(xù)介質單元邊界上插入界面單元，引入蘊含材料斷裂和強度準則的粘結法則(Cohesive law)來描述界面行為裂紋沿界面單元擴展，避免了斷裂準則選取與網(wǎng)格重構問題。粘結法得到了國內(nèi)外學者的廣泛關注和應用，如利用CFEM研究斷裂能隨著裂紋速度增加與裂紋分叉引起的比表面積增加之間的關系[2]及脆性材料裂縫擴展模擬等[3]。因為該方法考慮了界面特性，所以也用于模擬分層復合材料的裂紋起裂和擴展[4,5]。

雖然CFEM在模擬脆性材料斷裂方面表現(xiàn)出很大的優(yōu)勢，但仍存在一些不足，主要表現(xiàn)在隨著粘結單元(Cohesive interface)的數(shù)量增加，宏觀剛度折減，這種剛度折減受制于粘結單元尺寸和粘結剛度[6]。Tomar等[7]建議采用單元尺寸的上下界限來平衡剛度折減和準確性關系。也有學者通過提高粘結單元和體單元的剛度比來緩解剛度折減問題[8,9]。為了減弱全域范圍插入粘結單元造成的剛度折減和計算效率降低的問題，Camacho 等[10]提出了局部CFEM。規(guī)定只有應力達到材料強度時才在該應力區(qū)插入界面單元。這種方法使得剛度折減問題得到了解決，但同時也帶來了局部網(wǎng)格重劃分和結點重疊等一系列問題。Zhang等[6]提出了約束粘結單元法CFEM(constrained)，當粘結面狀態(tài)滿足特定條件時激活粘結單元，提高了計算效率，也緩解了剛度折減問題。

粘結單元法數(shù)值計算收斂問題主要由兩方面造成。(1) 隱式算法中無厚度的粘結單元容易在軟化行為中造成能量釋放問題，引起局部不穩(wěn)定[11]；(2) 粘結單元設置過大的罰剛度引起的數(shù)值振蕩[12]。有學者在粘結法則中引入粘性系數(shù)以消散材料軟化階段的能量釋放[11,13]，但這類方法涉及額外起裂準則和粘性系數(shù)的選取，增加了計算成本。Sarrado等[14]提出了將無厚度粘結面嵌入線彈性基質單元的方法，將材料的彈性模量和粘結面的初始罰剛度解耦，從而解決了罰剛度過大引起的收斂性問題。Manzoli等[15]將粘結面考慮成實體單元，即在體單元間插入由兩個大長寬比HAR(High-aspect-ratio)的三角單元組成的粘結單元對，通過損傷變量反映HAR單元的退化。HAR單元采用與體單元相同的強度，不需要單獨考慮粘結面的罰剛度，避免過大的初始罰剛度引起的收斂性問題。該方法應用于準脆性材料的斷裂[15]和水力壓裂問題的模擬[16,17]，然而只考慮了HAR單元法向變形的影響，對粘結面交叉點的空白處未作處理，因此仍有待完善。

在已有研究基礎上，本文提出了考慮厚度的局部粘結單元法TCFEM(Cohesive Finite Element Method with Thickness)。TCFEM在潛在裂紋擴展區(qū)內(nèi)插入粘結單元，將無厚度的粘結單元轉換成大長寬比的四邊形實體單元，采用拓展虛內(nèi)鍵(AVIB)本構[18]。當粘結單元完全開裂后，將自動轉變成只有抗壓剛度的Goodman節(jié)理單元。同時，將粘結面交叉處的空穴作為一個離散虛內(nèi)鍵(DVIB)模型[19]的多邊形鍵元胞來處理。相較于已有的考慮厚度的粘結面方法，TCFEM模型不需要單獨考慮粘結面的初始強度和損傷變量，且虛內(nèi)鍵模型考慮了法向和切向變形的貢獻，模型的幾何完整性也得到了完善。

2 TCFEM建模

TCFEM模型由體單元、考慮厚度的粘結單元和多邊形鍵元胞組成(圖1)。體單元在加載卸載中始終遵循線彈性本構關系；而材料的非線性行為由粘結單元和多邊形鍵元胞承擔?？紤]厚度的粘結單元采用AVIB本構[18]，多邊形鍵元胞滿足DVIB模型[19]。

圖1 TCFEM的單元組成

2.1 考慮厚度的粘結單元本構關系

AVIB源于虛內(nèi)鍵(VIB)，VIB是Gao等[20]提出的一種固體宏微觀本構模型。由于虛內(nèi)鍵勢函數(shù)蘊含了微觀斷裂機理，因此在模擬斷裂問題時不需要額外的斷裂準則。AVIB通過引入Xu-Needleman勢函數(shù)表征顆粒間法向及切向相互作用，克服了原VIB模型固定泊松比限制。AVIB考慮了脆性材料拉壓強度差異，適合準脆性材料的斷裂模擬。

根據(jù)文獻[18]，虛內(nèi)鍵的Xu-Needleman勢函數(shù)描述為

(1)

(2)

式中l(wèi)0為虛內(nèi)鍵的原長，ε為應變張量，ξ為變形前虛內(nèi)鍵的方向向量，上標T表示轉置。對于體積為V的代表單元體，其應變能密度為

(3)

根據(jù)超彈理論，應力張量和切線模量張量可表示為

(4)

在平面應力情況中，微觀虛內(nèi)鍵參數(shù)與宏觀材料常數(shù)之間關系為

(5)

式中E為楊氏模量，ν為泊松比。

為了避免AVIB在應變軟化的單元尺寸敏感性問題，Zhang等[18]提出在不改變材料強度條件下調整材料參數(shù)的方法，即軟化后的材料彈性模量和應變強度為

(6)

(7)

式中J為材料臨界J積分，h為粘結單元的厚度。

按式(6)調整材料參數(shù)主要是保證軟化過程中材料斷裂能守恒。對于裂尖單元不能一開始就應用調整后的材料參數(shù)，而應該是在進入峰后軟化階段才調整。在進入軟化階段前，材料還是由原來的材料參數(shù)控制。因此，裂尖材料開始調整參數(shù)的條件可表示為

(8)

當有厚度粘結單元的應變超過一定值時就意味著裂紋生成，此時粘結單元將自動轉化為接觸單元，以反映裂紋面之間的接觸和摩擦。粘結單元轉化為接觸單元的條件可表示為

(9)

式中λ為大于1的系數(shù)?？紤]計算效率問題，隱式計算的時間步長不能取得過小。在此情況下，如果λ取值剛好大于1，則在一個時間步內(nèi)，裂尖會有很多粘結單元應變滿足式(9)，而這些粘結單元并非真的已發(fā)生斷裂；如果λ取值太大，則會出現(xiàn)已發(fā)生斷裂的粘結單元沒有識別出來；同時也考慮到粘結面峰后能量釋放問題，本文取經(jīng)驗值λ≈3.0。

式(9)實際上為裂紋識別準則，其與一般材料破壞準則(如王來貴等[21])不同。本文采用AVIB本構，單元是否破壞由AVIB本構來控制。式(9)是將已斷裂的粘結單元識別出來。

2.2 粘結面交叉多邊形空穴本構關系

在體單元邊界插入有厚度的粘結單元后，在粘結面交叉點會形成不規(guī)則的多邊形空穴，如圖2所示。這些空穴的存在使得計算模型與原完整模型不匹配，人為地在原材料中引入了缺陷，空穴周圍引起應力集中，削弱了原材料的強度。為了彌補這一模型缺陷，本文將此多邊形空穴當作一個鍵元胞補上，采用DVIB[19]方法對其進行描述，即鍵元胞的節(jié)點力向量F為

(10)

鍵元胞的剛度矩陣Κ為

(11)

式中ui為鍵元胞中結點的位移向量。

鍵元胞中的鍵能參數(shù)與宏觀材料常數(shù)的關系為

圖2 粘結面交叉點處的多邊形鍵元胞

(12)

式中Ω為鍵元胞的虛內(nèi)鍵數(shù)量。

考慮到裂紋只能沿著粘結單元發(fā)生，粘結單元端點間通過多邊形鍵元胞銜接，故也需要引入鍵元胞破壞準則。為了簡化處理，認為結點鍵元胞為彈脆性。因而可采用鍵的應變準則來判斷鍵是否發(fā)生斷裂，即

(13)

當鍵元胞的虛內(nèi)鍵斷裂總量達到一定值時[19]，可認為鍵元胞發(fā)生破裂。

3 數(shù)值實現(xiàn)

3.1 TCFEM網(wǎng)格生成算法

本文開發(fā)了在指定區(qū)域內(nèi)生成有厚度粘結單元的算法。如圖3(a)所示，在常規(guī)網(wǎng)格的紅圈范圍內(nèi)生成無厚度的粘結單元，將這些單元按比例縮小，從而在單元之間形成了一定厚度間隙，如圖3(b)所示。在單元之間插入有厚度粘結單元，對于粘結面交叉處形成的空穴區(qū)域，引入多邊形鍵元胞進行填充，最終形成一個完整的網(wǎng)格結構，如圖3(c)所示。

圖3 生成有厚度粘結單元的主要步驟

3.2 數(shù)值實現(xiàn)過程

在數(shù)值模擬過程中采用隱式算法，通過Newton-Raphson法進行迭代計算。在計算過程中分別得到體單元、粘結單元和多邊形鍵元胞的剛度矩陣，整合成總體剛度矩陣，采用有限元方法進行處理和求解，具體計算流程如圖4所示。

圖4 TCFEM模型計算流程

4 數(shù)值模擬與驗證

4.1 考慮厚度的粘結單元長寬比影響

本文引入的粘結單元對比常規(guī)四邊形單元具有很大的長寬比，達到了1∶10及以上。漆文邦[22]指出，四邊形單元的計算精度會隨著長寬比的增大而降低。本文的粘結單元主要是為了表現(xiàn)體單元間的粘結作用，以抗拉特性為主，因此只需保證粘結單元的應變計算精度在可接受范圍內(nèi)即可。

為了驗證單元長寬比對模型精度的影響，采用圖5(a)的懸臂梁進行模擬。梁長L=8 m，高D=1 m，懸臂端承受分布剪力P=1 kN。材料參數(shù)為E=30.0×109Pa，v=0.25，梁的厚度為一個單元長度。按照平面應變問題求解。模擬懸臂梁的單元網(wǎng)格如圖5(b)所示，在網(wǎng)格中間層插入單層粘結單元，自由端承受分布剪力。

該梁自由端中間點撓度的解析解[22]為

(14)

粘結單元的長寬比分別取10∶1，20∶1和 50∶1，將模型計算結果與解析解對比，誤差分別為 1.0853%，1.1131%和1.1340%。表明單元長寬比越大，引起的誤差也越大，與文獻[22]結論一致。但模擬結果誤差都在可接受范圍內(nèi)，因此單元的長寬比不是影響本方法準確性的重要因素?？紤]到單元過大的長寬比會影響收斂性，故本文模擬算例中，粘結單元長寬比的取值均為10∶1。

圖5 粘結單元長寬比對計算精度的影響算例

4.2 算例驗證

4.2.1 拉伸型裂紋擴展模擬

圖6 三點彎(TPB)實驗試件尺寸(D =75 mm)與加載方案

圖8為三點彎梁的網(wǎng)格圖，在試件槽口附近潛在的裂紋擴展路徑區(qū)域插入有厚度的粘結單元(圖8 詳圖)。粘結單元的具體厚度根據(jù)單元尺寸決定，約為單元邊長的1/10。同時為了驗證網(wǎng)格尺寸是否對模擬結果有影響，采用了兩種尺寸的網(wǎng)格，即粗網(wǎng)格(單元總數(shù)為7369，結點總數(shù)為3828(size1))和細網(wǎng)格(單元總數(shù)為13126，結點總數(shù)為6742(size2))。

圖7 通過AVIB本構模擬的單軸拉伸應力-應變曲線

圖8 三點彎梁(TPB)實驗試件網(wǎng)格及局部網(wǎng)格詳圖

模擬結果如圖9所示。與文獻[23,24]相比，模擬的力位移曲線(圖9(a))與實驗結果吻合較好；同時,試件的變形和裂紋擴展路徑也與試驗基本一致(圖9(b))，說明該方法可以模擬材料拉伸斷裂行為。而且粗細網(wǎng)格模擬的結果也基本一致，說明模擬過程中斷裂能是守恒的，在很大程度上消除了單元尺寸敏感性問題。圖9(a)中曲線的峰值存在微小下降，主要是因為插入粘結單元后材料會發(fā)生一定程度的軟化。與傳統(tǒng)CFEM相比，由于本文方法采用有厚度粘結單元，相當于實體單元，因而不存在剛度折減問題。

4.2.2 復合型裂紋擴展模擬

圖9 三點彎梁(TPB)試驗與模擬結果

圖10 復合型試驗試件尺寸與加載方案

梁變形如圖11(b)所示。

圖11 復合型裂紋試驗結果(方案1)

圖12為第二種加載方式下的兩種試件模擬結果。圖12(a)顯示，P -CMOD 曲線的模擬結果與實驗結果基本吻合。試件加載中，裂縫從預設裂縫尖端起裂并發(fā)生擴展。圖12(b)顯示裂紋擴展路徑向加載點位置偏轉延伸，與實驗結果基本一致。這說明TCFEM方法可以有效再現(xiàn)復合型裂紋擴展行為。

圖12 復合型裂紋試驗結果(方案2)

5 結 論

本文發(fā)展了一種考慮厚度的局部粘結單元法，用于脆性材料的斷裂模擬。數(shù)值實現(xiàn)過程中，將指定區(qū)域內(nèi)的體單元按比例縮小，使單元之間形成一定間隙，從而生成有厚度的粘結單元。采用AVIB本構模擬粘結單元。由于粘結單元有了厚度，在粘結面交叉處形成多邊形空穴，采用DVIB鍵元胞將其補上，以保持計算模型的幾何完整性。由于插入的粘結單元是有厚度的實體單元，在軟化行為前采用與體單元相同的剛度，避免了傳統(tǒng)CFEM方法的剛度折減問題以及罰剛度過大而引起的收斂性問題。采用DVIB模型彌補了由于粘結單元厚度引起的粘結面交叉處多邊形空穴問題，保證了計算模型的幾何完整性。通過懸臂梁撓曲變形、混凝土梁三點彎和四點彎試驗模擬及對比，表明本文方法是有效的。本方法將在準脆性材料斷裂問題研究中具有廣闊的應用前景。