袁洪魏,唐 維

(1.中國工程物理研究院 化工材料研究所,四川 綿陽 621900;2.中國工程物理研究院 研究生院,北京 100193)

引 言

高聚物黏結炸藥(Polymer Bonded Explosive, PBX)作為先進武器的重要炸藥,在武器殺傷、破壞與動力能源方面發(fā)揮著關鍵性作用。由于PBX材料強度較低[1],在制造、運輸、服役過程中易受載荷刺激而產(chǎn)生裂紋,嚴重影響武器系統(tǒng)的安全性和可靠性[2]。作為武器裝備安全性和可靠性評估的重要手段,數(shù)值預測PBX炸藥在典型載荷下的裂紋起裂擴展規(guī)律具有重要意義。

PBX材料是由體積占比95%的含能晶體及少于5%的黏結劑高溫壓制而成,其力學特性具有顯著的黏彈性、拉壓不對稱性、壓力相關性以及溫度相關性等[3]。其黏彈性表現(xiàn)在加載速率效應[4]、蠕變[5]及松弛等行為;拉壓不對稱性表現(xiàn)在其壓縮強度遠高于拉伸強度,且拉伸強度一般小于10MPa[1],因此在結構件中PBX材料常以拉伸失效為主。壓力相關性表現(xiàn)在材料應力—應變響應與靜水壓力(圍壓)相關,靜水壓越大,則材料強度越高[6-8]。因此,PBX材料往往能承受長時壓縮載荷,而在拉伸載荷作用下容易發(fā)生低應力短時失效。與此同時,由于PBX材料脆性特征明顯(破壞應變較低[1]),在環(huán)境溫度沖擊[9]以及動態(tài)載荷[10]作用下易發(fā)生多裂紋并生且分叉融合等現(xiàn)象。為準確仿真PBX結構斷裂行為,一方面需要充分考慮PBX材料的黏彈性、拉壓不對稱性及壓力相關性等特殊力學性能,另一方面需要計算方法能進行多裂紋萌生、擴展、分叉和融合等復雜斷裂行為的仿真計算。

針對PBX結構斷裂行為仿真,首先從損傷力學和斷裂力學兩個角度分析了目前斷裂計算方法的現(xiàn)狀與困境,然后針對新興的斷裂相場法闡述了其基本思想、優(yōu)勢及在PBX炸藥中應用的不足,接著針對相場法的不足提出了考慮拉壓不對稱性的黏彈性斷裂相場計算方法,并通過多個實際算例展示了其對于PBX材料力學行為及結構裂紋行為仿真分析的優(yōu)勢與潛力。

1 現(xiàn)有斷裂計算方法

從力學學科上區(qū)分,斷裂行為的仿真分析方法一般可分為損傷力學和斷裂力學兩個分支。

1.1 損傷力學

損傷力學方面常采用含損傷本構模型的方式進行斷裂行為仿真,即在本構模型中引入一個損傷因子,當損傷因子達到1時便認為失效,材料失效區(qū)域的連通便被認為是裂紋,所以該方法可以較為便利地進行裂紋的起裂、擴展、分叉和融合仿真分析。

羅景潤[11]針對PBX拉伸力學行為建立了指數(shù)型損傷演化模型;李丹[12]采用唯象分析和細觀統(tǒng)計相結合的方法提出了PBX材料的率型損傷本構模型;李英雷[13]、黃垂藝[14]基于損傷型Z-W-T模型研究了考慮PBX材料非線性黏彈性的本構關系,這些損傷本構模型大都致力于應變率效應及單軸應力—應變曲線的高精度描述上,尚未考慮PBX材料拉壓不對稱性、圍壓相關性等。

在所有的PBX材料損傷本構模型中,不可避免地需要討論viscoSCRAM模型,該模型考慮了微裂紋、微孔洞等細觀損傷機制,反映了炸藥從變形到破壞的物理過程,能較好地描述簡單應力狀態(tài)下材料斷裂破壞的特征。Dienes[15]首先提出了著名的統(tǒng)計微裂紋模型(SCRAM),后又進一步發(fā)展[16],SCRAM模型基于應變率可加性,包含一個描述基體行為的彈塑性部分,加上描述微裂紋網(wǎng)絡(形核、長大和聚并)的各向異性損傷。Addessio和Johnson[17]對SCRAM模型進行了簡化,提出了具有各向同性損傷和無微裂紋聚合的IsoSCRAM模型。Bennett[18-19]在IsoSCRAM模型的基礎上附加黏彈性部分,從而提出了viscoSCRAM模型,該模型宏觀上考慮了炸藥的黏彈性,細觀上以微裂紋作為主要損傷機制。柳明[20]對viscoSCRAM模型進行了修正,增加了Bodner-Partom黏塑性,引入壓力相關拉伸損傷增強因子,使得模型可以較為準確地描述PBX材料的拉壓不對稱性。其他損傷本構模型大都采用類似思路進行模型構建[21]。

可以看出,目前損傷本構模型重點在于描述材料從損傷到失效的行為描述而非斷裂行為本身的準確預測,同時現(xiàn)有的損傷本構模型尚不能完全同時考慮PBX材料的拉壓不對稱性和壓力相關性;除此之外,基于損傷力學的仿真方法無法準確考慮裂紋尖端的應力奇異性,因此無法準確預測PBX結構的斷裂過程。

1.2 斷裂力學

目前基于斷裂力學的PBX結構裂紋仿真方法主要分為基于有限元框架的數(shù)值計算方法以及非有限元框架的新型數(shù)值計算方法。

基于有限元方法的斷裂數(shù)值方法主要包括界面單元法(CFEM)及擴展有限元法(XFEM)。界面單元法是在單元之間插入內(nèi)聚力單元以模擬裂紋的起裂擴展,Tan H[22]和郭虎[23]均開展了界面單元法在準靜態(tài)加載下PBX結構力學響應的應用研究,界面單元法雖然方法簡單,但其只允許裂紋在網(wǎng)格邊界上擴展,具有強烈的網(wǎng)絡依耐性和計算不穩(wěn)定性[24]。擴展有限元方法通過Enrich函數(shù)補充自由度,裂紋可在單元內(nèi)部進行擴展,彌補了內(nèi)聚力有限元法的不足。黃西成[25]、戴開達[26]、李生濤[27]等均采用XFEM對PBX材料在準靜態(tài)或動態(tài)加載條件下的斷裂過程開展了仿真研究,但是XFEM需要追蹤裂紋面,處理過程復雜,且在模擬三維裂紋任意擴展時存在困難[28]。其他的單元刪除法、廣義有限元法、網(wǎng)格重構法等,大都存在需要追蹤裂紋面、判據(jù)眾多、難以處理三維問題等缺點。

在非有限元框架的新型裂紋計算方法方面,主要包括離散元法(DEM)、近場動力學方法(PD)、數(shù)值流形元法(NMM)。離散元法是一種無網(wǎng)格方法,其將固體離散為一系列帶質(zhì)量的物質(zhì)點,通過求解物質(zhì)點之間相互作用的動力學方程以開展裂紋仿真分析。傅華[29-30]采用DEM模擬了PBX細觀尺度的動態(tài)加載試驗,但是離散元方法假設較多,缺乏理論嚴密性,同時細觀參數(shù)難以確定,計算收斂性較差[31]。近場動力學法是另一種無網(wǎng)格方法,通過用物質(zhì)點之間的關系積分方程代替微分方程求解物質(zhì)點的運動,以此仿真分析裂紋問題。李潘[32-34]建立起可應用于PBX材料的近場動力學模擬方法并應用于PBX試樣的準靜態(tài)破壞行為計算,但是PD計算效率低下,存在數(shù)值不穩(wěn)定性,難以進行三維問題分析,且傳統(tǒng)物性參數(shù)與PD模型參數(shù)轉換困難[34]。數(shù)值流形元法是通過兩套獨立覆蓋系統(tǒng)和截斷不連續(xù)形函數(shù)將連續(xù)與不連續(xù)分析統(tǒng)一在一個框架內(nèi),其對于裂紋擴展模擬沒有網(wǎng)格依賴性,易于處理復雜裂紋。北京理工大學陳鵬萬團隊[35-37]系統(tǒng)研究了基于數(shù)值流形元法的PBX材料宏細觀破壞過程仿真計算方法,但是NMM存在計算效率低下、收斂性較差、難以處理三維情形的問題[35]。因此,新興的裂紋數(shù)值計算方法由于均存在計算效率低下、難以計算復雜情形、宏觀物性參數(shù)難以轉換為仿真所需參數(shù)等問題,難以進行實際工程應用。

因此,現(xiàn)有的PBX結構斷裂行為計算方法,其均存在固有的不足而不能既考慮PBX材料的力學特性又可以計算PBX結構的復雜斷裂行為。

2 斷裂相場法的優(yōu)勢與不足

斷裂相場法(Phase Field Method to Fracture)自本世紀初誕生,便由于其理論完備且相對簡單、數(shù)值實現(xiàn)方便、且易于計算復雜裂紋等優(yōu)點,在近二十年得到了迅猛發(fā)展。

2.1 相場法的發(fā)展現(xiàn)狀與優(yōu)勢

1998年,Francfort和Marigo[38]將裂紋表面能作為連續(xù)體勢能的一部分,通過最小勢能原理建立了斷裂變分理論。隨后Bourdin和Francfort[39]于2000年結合相場理論,通過引入序參量(后來稱為相場標量),給出了斷裂面的彌散表達式,為斷裂相場模型的發(fā)展奠定了堅實的理論基礎。斷裂相場法的基本思想(以線彈性斷裂系統(tǒng)為例)為:其自由能同時包含彈性能和裂紋表面能,裂紋的起裂、擴展、分叉和融合等行為受自由能最小化原理控制。因此該方法具有3方面的優(yōu)勢:其一,理論完備且簡單,斷裂相場法是Griffith斷裂理論[40]的延伸,滿足能量守恒定律,并且僅需要臨界能量釋放率GIC斷裂判據(jù),避免了傳統(tǒng)斷裂理論需要定義起裂、擴展、分叉等多個準則的缺陷;其二,數(shù)值實現(xiàn)方便,斷裂相場法獨立于數(shù)值處理方法,因此可以基于成熟的有限元框架進行程序開發(fā),減少了方法實現(xiàn)的成本,易于推廣應用[41];其三,易于處理復雜問題,斷裂相場法中裂紋的萌生、擴展、融合、分叉等演化行為是一組耦合偏微分方程數(shù)值解的自然結果,均不需要額外處理,能方便地仿真各種復雜裂紋問題,避免了跟蹤復雜的裂紋幾何。

由于斷裂相場法獨特的優(yōu)勢,其自誕生便迅速發(fā)展。Miehe[42]成功模擬了I、II型裂紋在不同加載強度下裂紋的擴展長度和發(fā)展方向;Borden[43-44]針對動態(tài)加載下脆性材料斷裂行為提出了相應的數(shù)值計算策略并實現(xiàn)了二維和三維動態(tài)裂紋擴展仿真模擬;Mesgarnejad[45]通過多種實驗對不同脆性斷裂模型進行了定性和定量對比。為了進一步滿足壓剪情況下的起裂擴展行為,清華大學柳占立[46]、同濟大學Zhou等[47]基于Morh-Column失效準則修正了相場驅(qū)動能量項采用顯式相場法實現(xiàn)了脆性材料動態(tài)及準靜態(tài)壓剪失效仿真模擬。目前斷裂相場法在線彈性斷裂方面已發(fā)展較為成熟,線彈性斷裂相場理論系統(tǒng)如圖1所示,同時在彈塑性分析[48]、力熱耦合[49]、流固耦合(水力壓裂方面)[50]、力電耦合等多場耦合領域[51],通過在自由能泛函中進一步引入塑性耗散能、熱能等并結合相場變分理論開展了相場的擴展應用研究。

圖1 線彈性斷裂相場法理論Fig.1 Theory of the linear elastic fracture phase field method

2.2 線彈性相場法中的拉壓不對稱性

2.3 黏彈性相場法的現(xiàn)狀與不足

自2018年來多位研究人員開展了黏彈性相場法研究。Liu[59]提出了線性黏彈性斷裂相場方法,可以描述不同加載速率下的斷裂行為,該方法以彈性能作為裂紋驅(qū)動力。Loew[60]將廣義Maxwell黏彈性模型和Yeoh超彈性模型相結合,建立了橡膠材料的速率相關相場損傷模型。Yin和Kaliske[61]將橡膠材料的標準Maxwell黏彈性模型與Neo-Hooke超彈性模型相結合,建立了黏彈性斷裂相場模型,其采用Maxwell模型中平衡分支和非平衡分支的彈性能作為裂紋驅(qū)動力。Shen[62]提出了基于黏性耗散能驅(qū)動裂紋演化的黏彈性斷裂相場方法,不同于Loew[60],在該模型中,僅部分黏性耗散能驅(qū)動裂紋的演化,并研究了耗散能對蠕變、松弛、應變率效應和循環(huán)載荷作用下材料力學性能的影響。Huang[63]在Shen[62]的基礎上針對PBX材料建立了細觀斷裂相場仿真分析模型。Damma?[64]建立了統(tǒng)一的線性黏彈性斷裂相場模型,其在Shen[62]的基礎上考慮了體積變形和偏變形的黏彈性,并對彈性能和耗散能采用不同的退化函數(shù)。

值得注意的是,Liu[59]、Yin和Kaliske[61]以及Shen[62]均采用體積偏分分解法求得彈性能的拉伸分量和壓縮分量,其主要原因是黏彈性本構模型通常分解為體積部分和偏壓部分,因此用這種方法分解拉壓能量較為方便。然而,由于體積偏分分解中固有強度比范圍的限制,采用上述方法難以準確地描述PBX材料的拉壓不對稱性。

3 考慮拉壓不對稱性的黏彈性相場法理論修正

3.1 相場變分理論

基于斷裂相場法基本思想,本研究提出的黏彈性斷裂系統(tǒng)能量泛函為:

(1)

(2)

(3)

3.2 含損傷黏彈性本構模型

針對黏彈性本構模型,以廣義Maxwell模型為例進行了研究。同時考慮體積變形黏彈性和剪切變形黏彈性的廣義Maxwell模型如圖2所示。

圖2 廣義Maxwell模型示意圖Fig.2 Schematic diagram of the generalized Maxwell model

應變εij一般分解為體積應變volεij和偏應變eij:

(4)

應力σij也表述為體量volσij和偏量devσij之和:

(5)

根據(jù)元件定義,同時為提高數(shù)值計算穩(wěn)定性,采用Ambati[54]提出的混合相場模型,因此考慮相場引起損傷的本構模型為:

(6)

(7)

耗散能密度表示為:

(8)

以上即為本研究基于已有的黏彈性相場,考慮材料拉壓不對稱性的理論修正。

4 修正相場法的數(shù)值應用

針對上述建立的考慮拉壓不對稱性的黏彈性斷裂相場模型,本研究采用COMSOL進行了數(shù)值實現(xiàn),并通過多個實際算例來證明提出方法的有效性。同時在本研究中不研究耗散能貢獻,因此βv=0,且非平衡項體積模量取為1Kne=0。

4.1 拉壓過程仿真

4.1.1 單軸拉壓過程仿真

首先是針對單軸拉壓不對稱性的預測準確性驗證。針對PBX9502材料進行拉壓過程仿真分析,具體幾何尺寸、加載參數(shù)及試驗結果參考文獻[65-66]。仿真計算中采用中心對稱模型,約束及加載方式同試驗過程。采用最小二乘方法尋優(yōu)獲得的模型參數(shù)如下:Keq=1425MPa,μeq=380MPa,1μne=1130.5MPa,1τ=25.12s,Gc=13J/m2,l0=0.2mm,單元大小取為l0/he=2。計算得到的應力—應變曲線如圖3所示。

圖3 實驗與仿真得到的PBX9502應力—應變曲線對比Fig.3 Stress—strain curves for PBX9502 obtained from the experiment and simulation

由圖3可知,計算得到的拉伸壓縮應力—應變曲線與試驗數(shù)據(jù)吻合良好。計算得到的拉伸強度和壓縮強度分別為4.65MPa和12.28MPa,壓縮拉伸強度比為2.64;實驗所得拉伸及壓縮強度分別為4.53MPa和12.52MPa,壓縮拉伸強度為2.76;拉伸和壓縮強度誤差分別為1%和-2.2%。在本算例中僅采用了一階廣義Maxwell模型即可較為準確地描述該材料的拉伸和壓縮曲線。

啞鈴拉伸過程裂紋萌生擴展至斷裂過程如圖4所示。通過單軸壓拉斷裂行為的仿真分析,采用應變張量譜分解的相場法結合合適的本構模型可以精確地描述PBX材料的拉壓不對稱性特征。

圖4 PBX啞鈴拉伸過程裂紋萌生、擴展及斷裂過程(云圖為相場自由度d)Fig.4 Crack initiation, extension and fracture processes during PBX dumbbell stretching (legend for phase field d)

4.1.2 圍壓壓縮過程仿真

進一步開展該模型對于圍壓相關性的預測精度驗證。針對名為EDC37的PBX材料,試驗獲得的不同圍壓下的應力—應變曲線如圖5所示[8],圍壓越大,壓縮強度越高。計算模型為一個三維單元,x-、y-和z-面采用簡支約束,x+、y+面自由,加載過程中首先對x+、y+和z+面施加等靜壓,在此基礎上,進一步對z+面通過位移進行等應變率加載。模型參數(shù)如下:Keq=444.6MPa,μeq=263.6MPa,1μne=167.7MPa,1τ=25.12s,Gc=5J/m2,l0=0.2mm,單元大小取為l0/he=2。

圖5 不同圍壓下的壓縮應力—應變曲線Fig.5 Compressive stress—strain curves at different isostatic pressures

不同圍壓下仿真與試驗所得的應力—應變曲線如圖6所示,二者規(guī)律一致。

圖6 實驗與仿真獲得的不同圍壓下壓縮強度對比Fig.6 Comparison of the compression strengths obtained experimentally and simulatively at different confining pressures

從圖6中進一步可以看出,壓縮強度隨著圍壓壓力的增大呈線性增加;實驗與仿真獲得的壓縮強度隨圍壓變化的斜率分別為2.11和2.19,誤差為2.8%,說明本研究提出的相場法可以較為準確描述不同圍壓下的壓縮應力—應變曲線及壓縮強度。

4.2 帶孔板壓縮過程仿真

為進一步驗證提出的模型對于結構的預測精度,仿真分析了心帶孔平板試樣的壓縮破壞過程。材料為PBX9502。Liu[67]通過實驗研究了該構件在50℃下的壓縮破壞試驗,加載速率為1.27mm/min,并通過DIC獲得了全場變形及裂紋演化過程。圖7(a)為DIC獲得的裂紋形貌,其為沿著上下圓弧處的豎直方向裂紋。

圖7 裂紋形貌:(a)試驗結果[67];(b)仿真結果-應變張量譜分解;(c)仿真結果-體偏分解[62]Fig.7 Crack morphology: (a) experimental results[67]; (b) simulation results of the proposed method; (c) simulation results of Shen et al[62]

仿真模型采用二維平面應力模型,約束同實驗過程,底部固支、頂部約束水平方向位移同時按照1.27mm/min速率施加壓縮位移,參數(shù)同4.1節(jié)。

首先從裂紋形貌上分析,本方法獲得的結果與實驗裂紋形貌較如圖7所示,二者較為一致。由圖7可知,其裂紋由上下圓弧中心開始萌生,然后沿著豎直方向開始擴展,形成上下對稱的宏觀裂紋,同時計算得到的裂紋上下、左右完全對稱,證明本方法具有良好的計算魯棒性。而采用體偏分解的方法[62]計算得到裂紋由左右圓弧處萌生并沿水平方向擴展,與試驗結果不一致。進一步對比裂紋擴展長度均為16mm時的壓縮載荷,試驗與仿真獲得的此刻載荷分別為8.62和9.62MPa,相對誤差為11.6%。所以此方法不僅能計算準確的裂紋形貌,結合合適的本構模型后可以準確預測失效載荷。

4.3 PBX細觀壓縮過程仿真

本節(jié)通過PBX細觀模型的斷裂過程仿真來驗證方法對于復雜裂紋的計算能力。首先基于PBX細觀形貌分布特征[如圖8(a)所示]建立了球形填充的PBX炸藥細觀結構模型,模型大小為1mm×1mm;然后參考文獻[68]將界面進行彌散[見圖8(b)],彌散寬度取為0.002mm;然后將球形晶體、黏結劑及界面區(qū)域設置不同材料參數(shù);最后通過本方法進行豎直方向的單軸壓縮仿真,加載速率為6.67×10-3mm/s。

圖8 PBX壓縮斷裂演化圖:(a)試驗測得細觀形貌;(b)建立并界面彌散后的細觀仿真模型;(c)第一條裂紋萌生;(d)多裂紋萌生;(e)最終裂紋形貌;(f)試驗測得最終形貌Fig.8 PBX compression fracture evolution diagrams: (a) experimentally measured microscopic morphology; (b) microscopic simulation model after building and interface dispersion; (c) first crack initiation; (d) multi-crack initiation; (e) final crack morphology; (f) experimentally measured final morphology

模型及參數(shù)如下:晶體采用線彈性模型,彈性模量為13.375GPa,泊松比為0.32,臨界能量釋放率為0.25N/mm。黏結劑采用廣義Maxwell黏彈性本構模型,其彈性模量、剪切模量和松弛時間取值參考文獻[69],泊松比為0.495,臨界能量釋放率為2.5N/mm。界面采用線彈性模型,彈性模量為20MPa,臨界能量釋放率為0.12N/mm;晶體、黏結劑及界面的彌散裂紋寬度均取為0.002mm。

圖8(c)～(d)為壓縮過程微觀裂紋萌生、擴展-融合直至失效裂紋形貌圖。圖8(c)為第一條裂紋萌生,圖8(d)為多條裂紋萌生并擴展,圖8(e)[圖9為圖8(e)放大圖]為多裂紋分叉融合后形成的最終斷裂形貌,且與試驗結果[圖8(f)]較為吻合。通過圖8(c)～(d)可知,該方法可自然地進行多裂紋萌生、擴展、分叉和融合等復雜斷裂過程的仿真計算擴展,相對于擴展有限元等方法具有明顯優(yōu)勢。

圖9 最終時刻裂紋形貌(黃色虛線框為裂紋分叉/融合區(qū)域)Fig.9 Final crack morphology (yellow dashed box shows crack bifurcation/fusion area)

5 結 論

(1)理論方面,結合廣義Maxwell黏彈性本構模型,采用應變張量譜分解的方法對彈性勢能進行分解,認為只有拉伸彈性能和耗散能驅(qū)動損傷演化和裂紋形成,構建了考慮拉壓不對稱性的黏彈性斷裂相場模型。

(2)材料力學行為驗證方面,采用提出的模型對典型PBX材料的單軸拉伸壓縮過程進行仿真,拉伸壓縮強度描述誤差不大于2.2%,證明了提出的方法可精確描述PBX材料的拉壓不對稱性;同時對典型PBX材料不同圍壓壓縮過程進行模擬,強度隨圍壓變化斜率誤差不大于3%,證明了提出的模型可不經(jīng)額外處理而自然地描述不同圍壓下的力學性能。

(3)結構斷裂行為驗證方面,采用提出的模型計算了典型PBX帶孔板壓縮過程,從裂紋形貌上與試驗結果完全一致,從相同裂紋長度時載荷誤差小于12%,且計算具有較好的魯棒性;采用提出的模型仿真計算了PBX細觀結構模型的壓縮過程,證明了提出的模型依然可以自然地進行多裂紋萌生、擴展、分叉及融合復雜斷裂行為的仿真預測。

(4)黏彈性斷裂相場模型可同時考慮PBX材料復雜的力學特征以及PBX結構復雜的斷裂行為;針對PBX炸藥斷裂行為的仿真計算,相場法具有顯著的優(yōu)勢及潛力。