脆性聚甲基丙烯酸甲酯板動態裂紋傳播有限元模擬

張振亞，麻鴛鴛，周風華

(寧波大學沖擊與安全工程教育部重點實驗室，浙江寧波315211)

0 引言

20 世紀80 年代以來，脆性材料中動態裂紋傳播問題得到廣泛關注，取得了大量的實驗數據和理論研究成果［1-6］。但脆性材料在裂紋傳播過程中特有的物理現象，比如鏡面-拋物線-周期凹槽-分叉現象，以及特有的耗能模式，從數學上或經典斷裂力學上都還無法解釋，因此實驗手段和有限元分析是探索這些問題很好的手段。

對于動態裂紋傳播的行為，盡管有很多種研究方法，但數值模擬計算方法不失是一種很有效、合理的估計方法。目前有很多種模擬動態傳播過程，傳統方法使用動態應力強度或者J 積分作為斷裂準則。Nishioka 等把動態J 積分引入有限元網格成功模擬裂紋的傳播和扭結現象［7-8］，此方法能夠抓住裂紋尖端的奇異性，但要想獲得裂紋尖端的較好網格需要重新劃分網格技術。Belytschko 等［9］、Moes等［10］提出擴展的有限元法(XFEM)，這種方法的優點不需要重新劃分網格就能模擬任意裂紋，但它的算法過多依賴設定路徑的質量，對于比較復雜裂紋模式這個算法的精度難以保證。近年來，內聚力單元方法被廣泛用于模擬裂紋傳播行為，Camacho 等首先提出了一個簡單的剛性不可逆的內聚力準則［11-12］，pandolfi 等［13-14］運用此準則成功進行3D分析一些工程問題。基于材料的本身性質后來還出現了各種內聚力模型，如率無關和率相關等，這些數值方法都是基于自己的程序，對于動態裂紋傳播問題收斂較困難，需要調節參數。本文基于通用的大型軟件ABAQUS 處理器，該軟件算法比較成熟，容易收斂而且精度相對較高，結合開發的用戶子程序，這樣更加符合實際材料的物理背景。

本文是基于實驗的基礎上，基于內聚力單元方法，結合編制材料的用戶材料子程序，嵌入有限元軟件ABAQUS 求解器來模擬裂紋擴展行為。有限元模擬結果顯示:能夠抓住脆性材料聚甲基丙烯酸甲酯(PMMA)動態裂紋傳播中特有的狀態量(包括裂紋最小傳播速度、裂紋傳播速度、微分叉，局部分叉、裂紋傳播極限速度)。

1 動態裂紋傳播實驗

1.1 實驗條件

選用無色透明的PMMA 板材作為實驗材料，材料厚度D=3 mm，質地均勻，結構較為致密。材料的彈性模量E=3 000 MPa，泊松比υ =0.35，密度ρ =1.2 g/cm3. 將PMMA 板材裁為5 種不同尺寸的矩形，夾持在特制夾具中［15］。夾持后的試件實驗區域為長、寬分別是L×h 的5 種矩形，具體尺寸分別為:A 型(320 mm×240 mm)、B 型(320 mm ×200 mm)、C 型(320 mm ×160 mm)、D 型(240 mm ×240 mm)和E 型(240 mm ×200 mm). 利用絲網印刷技術在裂紋傳播區域順次布置等間距d =2.5 mm 的銀漆導電線路，每根線寬度為0.6 mm，導電線的厚度為10 μm. 實驗的其他方面在文獻［16］有詳細的描述。

1.2 實驗結果

實驗結果主要包括3 個部分:裂紋傳播速度，裂紋傳播速度與斷裂能關系和裂紋斷面留下的遺跡(局部分叉，宏觀分叉)。

如圖1 所示，裂紋速度經歷一個短暫加速階段進入傳播速度穩定段，且裂紋傳播速度最小值大約為215 m/s，這時的斷裂能為950 J/m2;裂紋傳播速度的最大值為650 m/s，對應的斷裂能大約為14 000 J/m2.

圖1 E 型試件在不同載荷作用下的裂紋傳播速度隨傳播距離的變化曲線Fig.1 Crack propagation velocities of E-type specimen under different loadings

如圖2 所示，裂紋穩定傳播速度v0是預加能量的單調遞增函數，并且材料表現出速度增韌現象。

如圖3 所示，這是裂紋傳播后留下的遺跡，從圖3(a)可以看出v0＜580 m/s，裂紋傳播路徑是一條直線，沒有出現局部分叉現象;圖3(b)描述了裂紋傳播速度v0在580 ～600 m/s 這個范圍，裂紋傳播路徑發生了彎曲，伴有對稱局部分叉產生，局部分叉的特征隨著裂紋傳播速度的增加，局部分叉越來越多，但裂紋的厚度沒有穿透試樣的厚度。這個階段從某種程度上反映了裂紋傳播開始出現失穩;圖3(c)代表了裂紋傳播速度達到610 m/s，裂紋傳播路徑開始出現宏觀一次分叉，隨著裂紋速度的提高，會出現多次分叉和次級分叉現象。由于提前裂紋出現分叉，無論施加的能量多大，單一主裂紋速度都不會增加下去，此時裂紋傳播速度的最大值為665 m/s，大約為材料的瑞麗波速cR的70%(PMMA材料的cR=930 m/s).

圖2 裂紋的動態傳播能與穩定傳播速度之間的關系Fig.2 The relationship between the dynamic propagating energy and the crack propagation velocity

圖3 裂紋傳播的軌跡和分叉圖(圖片上導電線間隔包括線寬為2.5 mm)Fig.3 Straight paths，local bifurcations and macro bifurcations of the propagating crack (the interval of conductive lines，including line width，is 2.5 mm)

2 動態裂紋傳播模擬

2.1 內聚力準則

混合型內聚力模型用于內聚力單元的本構關系，這種情況存在內聚力單元存在法向變形和剪切變形，即法向張開位移δn和剪切張開位移δs，有效張開位移δeff可以表示為

式中:β 表示法向張開位移和剪切張開位移占有的權重。有效拉應力σeff表示為

式中:σs表示剪切拉應力;σn為法向拉應力。內聚力準則是描述δeff與σeff之間的關系，普遍形式如下:

在本文有限元模擬中，利用一個簡單不可逆的線性衰減內聚力模型，在該模型中，首先忽略率相關的影響。內聚力的模型的形狀如圖4 所示。

圖4 初始剛度線性衰減不可逆內聚力準則Fig.4 The initial-rigid linear-decay irreversible cohesive law

假設脆性材料PMMA 的臨界應力σc是定值，而臨界張開位移δc隨著張開速度增加，這樣的率相關內聚力準則聯系著裂紋尖端斷裂過程，最終將導致伴隨微裂紋過程消耗更多的能量，這樣的內聚力模型能夠準確地描述脆性材料的斷裂過程。因此用下面的式子來描述δc和的關系

式中:δc0、和k 都是材料參數;δc0是靜態的裂紋張開位移是裂紋張開位移變化率;k 是率相關系數。率相關的內聚力模型如圖5 所示。

圖5 率相關的內聚力準則Fig.5 The rate-dependent cohesive law

σc是最大的內聚力，δc是臨界張開位移。圖4和圖5 中曲線所圍的面積就是裂紋面單位面積消耗的能量，即斷裂能Gc.

式中:Γc表示產生裂紋面的表面能。

2.2 材料參數

選用的材料是有機玻璃PMMA，材料性質:密度ρ=1 200 kg/m3，楊氏模量E =3 000 MPa，泊松比ν=0.35. 對于率無關的內聚力模型，最大內聚力σc和臨界張開位移δc是定值，在室溫狀態下材料PMMA表現出較小的塑性。而材料PMMA 的最大抗拉強度范圍是75 ～76 MPa，計算時選擇75 MPa. 實驗時發現材料PMMA 斷裂能Gc在420 ～480 N/m，因此選擇斷裂能Gc=450 N/m. 臨界張開位移δc根據(7)式計算出為12 μm. Gc也就是材料PMMA 靜態斷裂韌性。靜態裂紋前端內聚區的大小可以用Rice［17］公式估算如下:

利用(8)式可以計算出內聚力區的大小為0.24 mm.

2.3 有限元模型

計算模型是三維有限元模型。因為實驗中的幾何尺寸太大，選擇縮小的板來模擬。板的幾何尺寸:L =32 mm(x 軸)，高h =16 mm(y 軸)，厚度t =0.5 mm(z 軸)。有限元模型中所用的單元為:內聚力單元COH2D4，一般單元CPE3(三節點線性平面應變單元)，網格排列如圖6(a)所示，節點總數為51 884，單元總數為43 027. 板沿著垂直中心線有4 mm邊緣裂紋，板在y 軸施加一個相等邊界位移載荷Δ+和Δ-，試樣上的x 軸方向是自由的。在這種載荷下，裂紋通常是沿著徑向方向直線傳播。儲存在預加載荷板的單位面積的應變能W0計算如下:

圖6 網格形成與有限元模型Fig.6 FEM model and mesh

為了模擬實驗的實際工況，模擬工作分為兩個步驟:1)運行Plate_Preload.inp 文件，目的是計算出準靜態預加載荷Δ+和Δ-下算出板內應力應變(即板內儲存應變能);2)復制上步計算的結果文件Plate_Preload.Odb 和Plate_Preload. Prt 到當前文件夾，與Plate.inp 和裂紋起裂準則用戶子程序Vumat一塊運行。這樣做的好處是一方面邊界條件保持不變，同時把第1 步預加載荷算出的結果作為當前步已知條件(板內已經儲存了應變能)，啟動Vumat 即開動裂紋起裂準則，裂紋沿著設定的路徑就會依次釋放應變能。為了保證計算的精度和穩定性，模擬計算中設定時間步長為0.005 μs.

3 計算結果

3.1 率無關模擬結果

要模擬裂紋在不同載荷下的傳播速度情況，確定預加邊界位移δ0分別為0.045 mm、0.050 mm、0.055 mm、0.060 mm、0.065 mm、0.070 mm、0.080 mm，儲存在板內的能量分別為759 N/m、965 N/m、1 134 N/m、1 360 N/m、1 584 N/m、1 837 N/m、2 400 N/m. 裂紋傳播速度v0確定:在計算后處理中定義固定輸出的時間Δt，裂紋傳播的距離近似等于相鄰裂紋尖端輸出坐標差值Δl(裂紋直線傳播是精確滿足，傳播路徑發生彎曲近似滿足)，因此裂紋傳播速度為v0= Δl/Δt. 依此方法計算結果如圖7 所示。從圖7 可以看出:1)對于每一次實驗，動態裂紋傳播都是經歷了一個短暫的加速階段，裂紋傳播速度達到一個穩定速度v0，而且裂紋傳播路徑幾乎是直線傳播一直到斷裂;2)裂紋穩定傳播速度的最小值大約為235 m/s，正如圖中藍色虛線所示，也就是裂紋傳播速度小于這個值裂紋就會止裂不向前傳播，最大裂紋傳播速度為665 m/s，圖上黑色虛線表示的，無論外載荷再高，裂紋傳播都不會超過這個速度，并且裂紋傳播在低速時，即v0≤500 m/s，裂紋傳播速度變化較小，高速傳播時裂紋傳播速度變化較劇烈。只有當高載荷的情況下，即板內儲存較高的能量下，裂紋傳播速度達到580 m/s 左右，正如圖8(b)所示裂紋的傳播路徑會發生一點彎曲，并且沿著裂紋傳播路徑兩邊對稱出現微分叉。這些裂紋傳播速度以及裂紋傳播特有現象(微分叉、分叉)與實驗結果相契合，差別之處實驗中出現這些結果較模擬需要更高的能量，正如圖11 所示，在相同載荷下，率無關模擬結果遠遠大于實驗結果，出現的原因:裂紋傳播速度較大，模擬選擇Gc=450 N/m，對于PMMA這個材料太小，不能很好地反映材料的慣性效應。

圖7 裂紋速度隨試樣長度的變化曲線(利用率無關內聚力準則)Fig.7 Crack propagation velocity vs. specimen length. Simulations using rate-independent cohesive law

圖8 時間t=40 μs 時在不同載荷下的裂紋擴展云圖Fig.8 Crack propagation paths under different pre-loadings at 40 μs

3.2 率相關模擬結果

在這個率相關內聚力模擬計算中有4 個參數σc、δc0、k 需要確定。與前面計算相似，指定臨界應力σc=75 MPa，靜態臨界張開位移δc0=0.012 mm.可見靜態斷裂能Gc=450 N/m. 為了避免復雜的參數調試，假設k=1. 裂紋張開變化率決定率相關的程度，一個非常小的就會體現很明顯的率相關性。在計算中選定10 m/s，圖9 給定了t=30 μs 時刻時相同載荷、不同應變率的裂紋擴展云圖。從圖9 可以看出，值越大，裂紋傳播距離就越遠，意味著裂紋傳播速度就越大。

圖9 時間t=30 μs 時在相同載荷下的裂紋擴展云圖Fig.9 Crack propagation paths under W0 =1 837 N/m preloading at 30 μs. The rate-dependent cohesive law is used

圖10 不同載荷下裂紋傳播速度隨裂紋長度的變化曲線=5 m/s)Fig.10 The crack propagation velocity vs. specimen length.Simulations using rate-dependent cohesive law =5 m/s)

3.3 裂紋分叉

圖11 實驗結果與模擬結果對比圖Fig.11 The comparison of the experimental results and simulation results

采用率相關的內聚力模型來模擬裂紋分叉，在該模擬計算中，內聚力參數如前面一樣:σc=75 MPa，δc0=0.012 mm，k =15 m/s. 計算結果如圖12. 從計算結果可以看出，在W0≈7 350 N/m，裂紋傳播速度大約為v0=580 m/s 時，裂紋開始出現對稱的微分叉，微分叉幾乎貫穿了試樣的厚度，隨著裂紋速度進一步加大，正如圖12(b)，裂紋傳播路徑發生彎曲，有的微分叉開始發展成為一次分叉，有的保持微分叉，微分叉的數量明顯減少，可見微分叉發展為一次分叉消耗了大量的能量，與圖3(b)的結果基本是一致的。裂紋傳播速度繼續加大，即W0≈12 180 N/m，v0≈645 m/s 時，裂紋傳播路徑發生彎曲，接著出現一次分叉、二次分叉，速度達到660 m/s，裂紋還出現次級分叉，由于分叉過早的出現，消耗了一部分能量，導致了材料PMMA 的上極限速度永遠達不到瑞利波速cR.

4 結論

利用Cohesive 單元法，開發用戶材料子程序，結合ABAQUS 求解器模擬脆性材料PMMA 動態裂紋傳播行為。模擬計算過程中考慮采用率無關和率相關兩種情況，其模擬結果如下:

1)相同載荷下，采用率無關內聚力準則模擬的穩定裂紋傳播遠大于實驗測出的穩定裂紋傳播速度。

2)無論是采用率無關還是率相關的內聚力準則，模擬出的裂紋傳播速度的上限速度都是665 m/s，下限速度為235 m/s，與實驗測試結果是吻合的，這也說明了這是材料PMMA 固有的材料屬性，與其他無關。

圖12 =5 m/s 時在不同載荷下裂紋分叉圖Fig.12 The bifurcations of cracks under different loads

3)利用率相關的內聚力模型能夠較好地模擬裂紋傳播行為，模擬結果與實驗結果基本一致(見圖11)，說明對于脆性材料PMMA，必須要考慮裂紋尖端張開位移變化率，這樣能夠消耗很多的能量。

References)

［1］Ravi-Chandar K，Knauss W G. An experimental investigation into dynamic fracture—Ⅳ. On the interaction of stress waves with propagation cracks［J］. International Journal of Fracture，1984，26(3):189 -200.

［2］Fineberg J，Marder M. Instability in dynamic fracture［J］. Physical Reports，1999，313:1 -108.

［3］Scheibert J，Guerra C，Dalmas D. Brittle-quasibrittle transition in dynamic fracture:an energetic signature［J］. Physical Review Letters，2010，104:045501.

［4］Xia K W，Chalivendra V B，Rosakis A J. Observing ideal“selfsimilar”crack growth in experiments［J］. Engineering Fracture Mechanics，2006，73:2748 -2755.

［5］Fineberg J，Gross S P，Marder M. Instability in dynamic fracture［J］. Physical Review Letters，1991，67:457 -460.

［6］Fineberg J，Gross S P，Marder M. Instability in the propagation of fast cracks［J］. Physical Review B，1992，45:5146 -5154.

［7］Nishioka T. Recent developments in computational dynamic fracture mechanics［C］∥Aliabadi M H. Dynamic Fracture Mechanics Int Series Comput Engng. Southampton，UK，and Boston，US:Computational Mechanics Publications，1995:1 -60.

［8］Nishioka T，Tokudome H，Kinoshita M. Dynamic fracture-path prediction in impact fracture phenomena using moving finite element method based on Delaunay automatic mesh generation［J］.International Journal of Solids and Structures，2001，38 (5):273 -301.

［9］Belytschko T，Black T. Elastic crack growth in finite elements with minimal remeshing［J］. International Journal for Numerical Methods in Engineering，1999，45(6):1 -20.

［10］Moes N，Dolbow J，Belytschko T. A finite element method for crack growth without remeshing［J］. International Journal for Numerical Methods in Engineering，2000，46(1):31 -50.

［11］Camacho G T，Ortiz M. Computational modelling of impact damage in brittle materials［J］. International Journal of Solids and Structures，1996，33(2):899 -938.

［12］Camacho G T，Ortiz M. Adaptive Lagrangian modelling of ballistic penetration of metallic targets［J］. Computer Methods in Applied Mechanics and Engineering，1997，142:269 -301.

［13］Pandolfi A，Krysl P，Ortiz M. Finite element simulation of ring expansion and fragmentation:the capturing of length and time scales through cohesive models of fracture［J］. International Journal of Fracture，1999，95:279 -297.

［14］Pandolfi A，Guduru P R，Ortiz M，et al. Three dimensional cohesive-elements of dynamic fracture in C300 steel［J］. International Journal of Solids and Structures，2000，37:3733 -3760.

［15］張振亞，周風華. 脆性PMMA 帶板中的動態裂紋傳播實驗［J］.實驗力學，2012，27(4):401 -407.ZHANG Zhen-ya，ZHOU Feng-hua. Experiment of dynamic crack propagation in brittle PMMA strip［J］. Journal of Experimental Mechanics，2012，27(4):401 -407.(in Chinese)

［16］張振亞，段忠，周風華. 脆性裂紋動態傳播過程中的速度振蕩現象和理論分析［J］. 力學學報，2013，45(5):729 -738.ZHANG Zhen-ya，DUAN Zhong，ZHOU Feng-hua. Experimental and theoretical investigations on the velocity oscillations of dynamic crack propagating in brittle material［J］. China Journal of Theoretical and Applied mechanics，2013，45(5):729 -738.(in Chinese)

［17］Rice J R. The mechanics of earthquake rupture［C］∥Dziewonski A M，Boschi E. Proceedings of International School of Physics Enrico Fermi. North Holland:Physics of the Earth’s Interior，1980:555 -649.