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

張振亞，麻鴛鴛，周風(fēng)華

(寧波大學(xué)沖擊與安全工程教育部重點(diǎn)實(shí)驗(yàn)室，浙江寧波315211)

0 引言

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

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

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

1 動態(tài)裂紋傳播實(shí)驗(yàn)

1.1 實(shí)驗(yàn)條件

選用無色透明的PMMA 板材作為實(shí)驗(yàn)材料，材料厚度D=3 mm，質(zhì)地均勻，結(jié)構(gòu)較為致密。材料的彈性模量E=3 000 MPa，泊松比υ =0.35，密度ρ =1.2 g/cm3. 將PMMA 板材裁為5 種不同尺寸的矩形，夾持在特制夾具中［15］。夾持后的試件實(shí)驗(yàn)區(qū)域?yàn)殚L、寬分別是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). 利用絲網(wǎng)印刷技術(shù)在裂紋傳播區(qū)域順次布置等間距d =2.5 mm 的銀漆導(dǎo)電線路，每根線寬度為0.6 mm，導(dǎo)電線的厚度為10 μm. 實(shí)驗(yàn)的其他方面在文獻(xiàn)［16］有詳細(xì)的描述。

1.2 實(shí)驗(yàn)結(jié)果

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

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

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

如圖2 所示，裂紋穩(wěn)定傳播速度v0是預(yù)加能量的單調(diào)遞增函數(shù)，并且材料表現(xiàn)出速度增韌現(xiàn)象。

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

圖2 裂紋的動態(tài)傳播能與穩(wěn)定傳播速度之間的關(guān)系Fig.2 The relationship between the dynamic propagating energy and the crack propagation velocity

圖3 裂紋傳播的軌跡和分叉圖(圖片上導(dǎo)電線間隔包括線寬為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 動態(tài)裂紋傳播模擬

2.1 內(nèi)聚力準(zhǔn)則

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

式中:β 表示法向張開位移和剪切張開位移占有的權(quán)重。有效拉應(yīng)力σeff表示為

式中:σs表示剪切拉應(yīng)力;σn為法向拉應(yīng)力。內(nèi)聚力準(zhǔn)則是描述δeff與σeff之間的關(guān)系，普遍形式如下:

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

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

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

式中:δc0、和k 都是材料參數(shù);δc0是靜態(tài)的裂紋張開位移是裂紋張開位移變化率;k 是率相關(guān)系數(shù)。率相關(guān)的內(nèi)聚力模型如圖5 所示。

圖5 率相關(guān)的內(nèi)聚力準(zhǔn)則Fig.5 The rate-dependent cohesive law

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

式中:Γc表示產(chǎn)生裂紋面的表面能。

2.2 材料參數(shù)

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

利用(8)式可以計(jì)算出內(nèi)聚力區(qū)的大小為0.24 mm.

2.3 有限元模型

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

圖6 網(wǎng)格形成與有限元模型Fig.6 FEM model and mesh

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

3 計(jì)算結(jié)果

3.1 率無關(guān)模擬結(jié)果

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

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

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

3.2 率相關(guān)模擬結(jié)果

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

圖9 時(shí)間t=30 μs 時(shí)在相同載荷下的裂紋擴(kuò)展云圖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 實(shí)驗(yàn)結(jié)果與模擬結(jié)果對比圖Fig.11 The comparison of the experimental results and simulation results

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

4 結(jié)論

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

1)相同載荷下，采用率無關(guān)內(nèi)聚力準(zhǔn)則模擬的穩(wěn)定裂紋傳播遠(yuǎn)大于實(shí)驗(yàn)測出的穩(wěn)定裂紋傳播速度。

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

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

3)利用率相關(guān)的內(nèi)聚力模型能夠較好地模擬裂紋傳播行為，模擬結(jié)果與實(shí)驗(yàn)結(jié)果基本一致(見圖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］張振亞，周風(fēng)華. 脆性PMMA 帶板中的動態(tài)裂紋傳播實(shí)驗(yàn)［J］.實(shí)驗(yàn)力學(xué)，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］張振亞，段忠，周風(fēng)華. 脆性裂紋動態(tài)傳播過程中的速度振蕩現(xiàn)象和理論分析［J］. 力學(xué)學(xué)報(bào)，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.