基于碎點法的動態斷裂分析1)

沈寶瑩 王松 李明凈 董雷霆

(北京航空航天大學航空科學與工程學院,北京 100191)

引言

沖擊防護結構廣泛存在于航空航天、汽車、船舶、核能乃至國防等工程領域.這類結構在抵抗如撞擊、爆炸等極端的沖擊載荷時可能會發生動態斷裂直至最終破壞.準確預測結構在動態載荷下的斷裂破壞行為對防護結構的設計和改進具有重要意義.數值仿真方法是預測結構動態斷裂的重要手段,學者在這個方面已經進行了長期、大量的研究,然而由于動態斷裂問題的復雜性,這仍然是計算固體力學領域一個具有挑戰性的研究方向.

有限元法是當前工程領域應用最廣泛的數值仿真方法,有限元法基于伽遼金弱形式且使用可以使用簡單多項式的形函數,具有穩定性強、數值積分簡單等優點[1].然而有限元法在模擬動態斷裂導致的結構破壞時存在兩個方面的難題.第一個方面,有限元法本質上是基于位移連續性假設的,因此難以在模型中顯式引入裂紋.當前解決這個問題的一種方法是使用網格刪除技術,然而刪除網格會導致和物理過程不一致的材料損失[2].另一種方法是使用內聚力單元模擬內界面失效,但是固有型內聚力單元存在人工柔度問題,非固有型內聚力單元需要使用較大的內界面懲罰系數,在模擬復雜斷裂問題時存在困難[3-4].擴展有限元法能在有限元模型中顯式引入裂紋,是一種受到廣泛關注的計算方法,但該方法在引入裂紋時會引入額外的自由度,且通常需要對模型中的裂紋區域進行網格細化和重構,計算效率和穩定性有待進一步提高[5].第二個方面,使用有限元法模擬復雜動態斷裂問題時往往會伴隨裂紋附近局部區域的網格畸變,由此導致顯式動力學求解的臨界時間步長顯著減小,解的精度下降[6-9].過高的網格畸變甚至會因負體積等問題導致計算終止[9].當前通常的解決辦法是在計算過程中刪除模型中的畸變,但這同樣會導致與物理過程不符的材料損失.

以上所述的傳統有限元方法在模擬動態斷裂中面臨的諸多制約,其根源在于有限元方法使用單元離散求解域,并基于單元節點計算形函數,因此有限元對網格具有很強的依賴性.無網格方法是部分或完全擺脫網格依賴的一類數值方法,其核心思想是采用節點群對結構進行離散,并使用節點描述被離散結構的力學響應[10].因為無網格類方法能有效避免傳統有限元法普遍存在的網格畸變問題,因此這類方法一直是國內外計算力學界的研究熱點之一.從20世紀90年代起,學者們提出的無網格方法已達20余種,這些無網格法大體上可以分為兩種:粒子類方法和基于弱形式的無網格方法.

光滑粒子流體動力學(SPH)方法[11-12]是沖擊領域應用較為廣泛的一種粒子類無網格法,該方法利用光滑核函數近似模擬場函數,其離散控制方程是配點形式的[13].每個SPH 粒子具有質量、速度、體積等物理屬性,這種特性使得SPH 在顯式動力學仿真中具有很大優勢.但由于SPH 基于配點的強形式,其穩定性難以獲得嚴格的數學證明,因此SPH 往往需要較大的支持域,導致SPH 法的計算量通常高于有限元法[14].物質點法(MPM)是另一種代表性的粒子類方法.MPM 采用歐拉和拉格朗日雙重描述將物質離散為一組在空間網格上運動的質點,這些質點同樣攜帶質量、體積、速度等物理信息.動量方程在歐拉描述的空間網格中離散,避免了大變形造成拉格朗日網格畸變的問題[15].文獻[16-18]采用物質點法在物體超高速碰撞、沖擊侵徹、爆炸、動態斷裂等方面開展了一系列研究,取得了重要進展.

Atluti等[19-20]提出的無網格局部彼得羅夫-伽遼金(MLPG)方法和Belytschko等[21-22]提出的無單元伽遼金(EFG)方法是兩種較為典型的弱形式無網格方法.這兩種無網格法基于弱形式方程,因此具有較好的數值穩定性[19,21,23].然而該類方法的形函數通常為復雜有理函數,這導致了在積分時無法籍由簡單的高斯積分直接計算,積分方法不當甚至會影響方法的穩定性.文獻[24-26]提出一系列積分方法,在降低計算量、提高積分精度及方法穩定性等方面具有明顯效果.

近年來,針對動態載荷作用下的復雜結構破壞問題,國內外學者在近場動力學、邊界元、間斷伽遼金有限元法等方面也開展了一定工作.劉立勝等采用近場動力學開展了濕熱環境下復合材料沖擊損傷及動態斷裂的研究[27],章青等進行了爆炸載荷作用下混凝土結構破壞過程的近場動力學模擬[28-29],Han等[30]研究了耦合損傷力學與近場動力學的材料失效仿真方法,高效偉等[31]采用邊界元法分析功能梯度材料動態斷裂力學問題,于福臨等開展了基于間斷伽遼金方法的船體板架爆炸沖擊響應數值模擬研究[32],相關的研究工作仍在進行當中.

綜上所述,國內外學者針對動態斷裂問題已經開展了多種數值仿真方法的研究,但由于該問題物理過程的復雜性,對動態斷裂過程的準確仿真預測仍然是一項具有挑戰性的工作.

近年來,作者所在研究團隊提出了一種不連續型無網格方法:碎點法(FPM)[33],該方法一方面參考弱形式無網格方法,采用弱形式方程,具有較好的數值穩定性,且基于支持域節點群定義形函數,使子域具有抵抗畸變的能力;另一方面參考間斷伽遼金有限元法,不對相鄰子域間的位移作連續性要求,因此該方法易于在模型中顯式引入裂紋.

碎點法使用空間中的節點群對問題域進行離散,基于節點群劃分具有明確幾何形狀的子域,如圖1(a)所示.碎點法模型中僅考慮節點的自由度,子域頂點信息由節點插值獲得.節點反映的是相應子域的物理狀態,因此節點具有體積、質量、速度等明確的物理含義.對于碎點法模型中的任意一個節點P0及其子域E0,其支持域為子域E0的所有相鄰子域E1～E6,則子域E0中的位移形函數是基于支持域節點群P1～P6定義的,如圖1(b)所示.由于碎點法使用支持域節點群而非子域頂點構建形函數,計算過程中子域形狀的畸變不會對形函數的計算造成影響,因此碎點法具有抵抗子域形狀畸變的能力.當然碎點法的支持域定義具有很大的靈活性,可以根據實際計算的需要進行修改,本文僅以圖1(b)所示的方法為例進行說明.

圖1 (a) 碎點法模型的離散和(b) 支持域及其節點群的定義Fig.1 (a) Discretization of FPM model and(b) definition of support range and its point cloud

碎點法不對相鄰子域間的位移形函數和試函數作連續性要求,即所使用的形函數和試函數可以是分片連續的,如圖2 所示.由于只要求函數在子域中局部連續,碎點法可以使用簡單多項式形式的位移形函數和試函數,因此子域內的積分計算可以使用簡單的高斯積分進行求解,從而避免了大多數弱形式無網格法因需要使用復雜有理數形函數而導致積分困難的問題.另一方面,由于碎點法使用的位移試函數是分片連續的,可以自然地在任意相鄰子域之間的內部界面處顯式引入裂紋,因此碎點法具備模擬大范圍復雜裂紋萌生和擴展問題的能力.在碎點法模型中顯式引入裂紋的具體方法將在下文中詳細介紹.

圖2 (a) 碎點法形函數示意圖(b) 位移場為的碎點法位移試函數[34]Fig.2(a) Shapefunctionof FPM(b) FPM’strialfunction of a field of displacement:

碎點法的方程是基于伽遼金弱形式構建的,因此該方法具有較好的數值穩定性.然而,碎點法的位移形函數和試函數是分片連續的,模型中相鄰子域間的位移連續性缺乏約束,因此方法的一致性無法得到保證.為了解決這個問題,參考間斷伽遼金有限元法,在弱形式方程中任意內界面處引入數值通量修正,以弱形式約束相鄰子域的位移連續性.在作者所在團隊的前期工作中,已經證明了在弱形式方程中引入內界面數值通量修正可以有效保證碎點法計算的一致性[33].

綜上所述,碎點法具有如下優勢:(1)采用基于節位移形函數,因此具有抵抗畸變的能力;(2)放棄了相鄰子域間的形函數連續性要求,因此可以采用簡單多項式形式的位移試函數并用高斯積分求解弱形式方程,且易于在模型中顯式引入裂紋;(3)子域具有質量、體積、速度等物理屬性,因此便于顯式動力學的計算.基于上述特點,作者認為碎點法適合用于預測結構的動態斷裂行為.

碎點法已經在多個領域得到了檢驗和應用.Yang等[34]應用碎點法進行準靜態應力分析和斷裂分析,驗證了碎點法求解上述問題的有效性和準確性;Wang等[35]應用碎點法預測含U 型缺口脆性試件的斷裂強度和裂紋形貌,證明碎點法可以有效分析復雜的準靜態脆性斷裂問題;Guan等[36-37]證明了在分析撓曲電材料的斷裂時,碎點法相比于傳統有限元法具有計算簡單、易于顯式引入裂紋的特點;Guan等[38-39]應用碎點法求解非均質材料與復雜薄壁結構的瞬態熱傳導問題,展示了碎點法在解決復雜問題域及邊界的瞬態熱傳導問題方面的準確性、高效性及穩定性.

本文旨在提出基于碎點法核心思想的動力學碎點法理論并開發相應程序,同時驗證該方法分析應力波傳播和動態裂紋擴展的有效性.本文的理論推導和算例驗證都是基于二維空間小變形假設.本文的結構安排如下:第一章介紹碎點法的基本理論,包括碎點法基本理論、求解域離散方法、弱形式動量方程及其離散、顯式動力學求解格式;第二章通過算例驗證動力學碎點法在預測應力波傳播和動態裂紋擴展方面的能力;第三章對本文做出結論.

1 碎點法基本理論

1.1 求解域離散

碎點法采用布置節點的方式對求解域 Ω 進行離散,在建立碎點法方程時僅考慮節點處的自由度.基于節點進一步將求解域 Ω 劃分若干個子域 Ωi,每個子域內部包含一個節點,子域具有明確的幾何形狀,且相鄰子域之間互不相交重疊.碎點法的離散方式如圖1(a)所示.

碎點法的離散有兩種基本方式,其一是先在求解域中布置節點,然后基于節點定義Voronoi 多邊形作為子域,如圖1(a)所示.另一種是借助成熟的有限元前處理軟件,先在求解域中劃分有限元網格,以有限元單元作為碎點法的子域,然后在子域中定義節點.

對于碎點法模型中的任意子域,在子域節點處對其位移函數u進行泰勒級數展開,從而將該子域內的位移函數表示為節點位移和節點位移梯度的函數,由此便定義了子域內的多項式形式的位移試函數,如式(1)所示.以包含節點P0的子域E0為例,在P0處進行泰勒級數展開,則子域E0內的位移試函數如下

節點處的位移梯度可以采用多種方法進行求解,本文所采用的是廣義有限差分法.首先定義目標節點P0和子域E0的支持域:本文定義與子域E0相鄰且具有公共邊界的子域作為支持域,這些子域的節點P1,P2,···,Pm組成支持域節點群,如圖1(b) 所示.在定義了子域的支持域之后,可通過廣義有限差分法來求解位移梯度.

定義L2范數J如下式所示

其中a為位移梯度向量,u0和um分別為組裝后的P0和P1～Pm節點位移向量,A是包含P0和P1～Pm之間坐標差的張量,W是權函數張量,具體表達式如下

其中,uE為包含P0節點支持域節點群(P0～Pm)位移的組裝位移向量,I1和I2為兩個轉換張量,表達式如下

將式(4)代入式(3)可得節點P0處位移梯度向量a與支持域節點群位移向量uE的表達式如下

將式(5)代入式(1)中,可以得到子域E0中的位移試函數uh與目標子域的支持域節點群自由度uE之間的關系如下

矩陣N為子域E0內任意一點x=[x,y]T處的位移形函數,其表達式如下

1.2 碎點法的弱形式動量方程

在求解域 Ω 中,強形式動量方程和邊界條件為:

式中 ρ為密度,是加速度向量,是施加的體力向量,σij是應力張量.邊界條件中的 Γt和Γu分別代表施加應力邊界條件和位移邊界條件的外邊界,ni是邊界上的外法線單位向量,和分別為在外邊界Γt和Γu上施加的面力和位移.

以虛位移 δui作為檢驗函數,采用伽遼金法推導弱形式動量方程如下

對式(8) 括號內第一項使用分部積分和高斯公式可得

式中 ?E表示子域E的邊界,nj為邊界?E上單位外法向量.

在碎點法模型中,任意子域E都需滿足式(9),對所有子域上的方程進行累加可得

定義 Γ 表示所有子域邊界的集合,則Γh=ΓΓt-Γu表示所有相鄰子域間內邊界的集合.假定任意相鄰子域E1和E2的公共邊界為e,對于任意物理量w,其在邊界e上的跳躍算子[]和平均算子{ }定義為

將跳躍算子和平均算子以及面力邊界條件引入式(10)的第四項中可得

式中 ηh內邊界 Γh處的懲罰系數,用于施加相鄰子域間的弱形式位移連續性約束,ηu為外邊界 Γu處的懲罰系數,用于施加弱形式位移邊界條件,lh和lu分別為 Γh和Γu處子域界面的特征長度.

合并式(14)中等號右邊的第三和第四項,并引入內部界面數值通量的定義如下

最終推導得到包含數值通量的碎點法弱形式動量方程如下

式中等號右邊第三項成為數值通量修正項,可以看出在碎點法弱形式方程中,相鄰子域之間的相互作用僅由數值通量修正項決定,因此式(15)中定義的數值通量包含相鄰子域間相互作用面力的物理含義.

1.3 碎點法弱形式方程的離散

1.3.1 離散格式的弱形式方程離散格式的弱形式方程

為了便于表達,以下推導過程均采用Voigt 標記法的矩陣表達式.

由位移試函數與支持域節點自由度之間的關系式(7)可得速度和加速度的試探函數表達式如下

由于試探解與檢驗函數采用相同形函數,對虛位移 δu有

為子域E0的應變形函數.

將式(7)及式(17)～式(20)代入式(16)中,對碎點法求解域中所有子域進行組裝,并在等號兩邊同時消去全局虛位移向量 δuE,可得到離散形式的碎點法動量方程

其中,矩陣M是全局質量矩陣,是全局加速度矩陣,fext是全局等效節點外力矩陣,fint是全局等效節點內力矩陣,其表達式如下

式中,為子域的體力矩陣,為基于應力邊界條件施加的面力,為子域的應力,是內邊界上的數值通量,是基于內邊界單位法向量將子域應力轉換成內邊界面力的投影矩陣,為子域節點處的位移,為基于位移邊界條件施加的位移.

1.3.2 引入裂紋對弱形式方程的處理

碎點法使用的是分片連續的弱形式方程,因此易于在任意內部界面處顯式引入裂紋.在離散格式的碎點法弱形式方程式(22)中,任意相鄰子域之間的聯系包含兩個部分,一部分是數值通量上的聯系,如上所述,數值通量t*具有相鄰子域相互所用面力的物理含義;另一部分是位移形函數上的聯系,因為碎點法的形函數是基于支持域節點群計算的,而相鄰子域互為對方的支持域.

由于數值通量的物理含義是相鄰子域的相互作用面力,可以自然地基于數值通量定義內界面斷裂準則.本文使用最大拉應力準則作為斷裂準則,即當碎點法模型中任意內界面處的法向數值通量達到材料的拉伸強度時,則判斷該內界面發生斷裂.當然,本文僅以最大拉應力準則為例,其他斷裂準則也可以用作碎點法的內界面斷裂判據.

對于碎點法模型中的任意內部邊界,當邊界的數值通量滿足斷裂準則時,就在該邊界處引入裂紋.如上所述,任意相鄰子域之間的聯系包括兩個部分,因此需分別進行修正.首先,對于斷裂的內邊界,移除弱形式方程中該邊界處的數值通量修正項,即將數值通量置零t*=0,就相當于消除了數值通量部分的相互聯系,如圖3(b)所示.然后,對該內界面兩端子域各自的支持域進行修正,分別從各自的支持域中將對方移除,就相當于消除了形函數部分的相互聯系.

圖3 碎點法裂紋的引入Fig.3 Steps of introducing crack in FPM

通過以上步驟,就可以在碎點法模型中任意內部邊界處引入顯式裂紋.值得一提的是,采用上述方法在碎點法模型中引入裂紋時,只修正了裂紋兩端子域各自的支持域和移除了該內界面處的數值通量,這些修正只影響裂紋附近局部區域,對整體弱形式方程的影響較小,且不需要進行網格重構,也不會額外增加整體模型的自由度.由此可見,碎點法只需要通過簡單的修正就可以在模型中的任意內邊界處顯式引入裂紋,該操作僅對弱形式方程造成局部影響,且不增加自由度,因此碎點法可以簡單、高效地模擬結構中任意位置的裂紋萌生和任意形貌的裂紋擴展,在模擬斷裂、破碎等極端問題時具有一定的優勢.

1.4 碎點法的顯式動力學求解

本文采用中心差分法來對碎點法弱形式動量方程進行求解,中心差分法的時間軸如圖4 所示.

圖4 中心差分法時間軸示意圖Fig.4 Time axis of the Verlet method

將式(24)和式(25)整理可得

這種形式的中心差分法也通常被稱作蛙跳格式.基于該方法,定義動力學碎點法的求解過程如下:

(1) 由式(23)計算tn時刻的質量和節點力矩陣

(7) 重復上述(1)～(6)步,直至最后一個時間步.

由于中心差分法是條件穩定的,因此為了保證計算的穩定性,每一個時間步的時間增量Δt必須小于臨界時間步長Δtcr.

下面討論中心差分法的穩定性問題.為了維持算法的穩定性,時間步長不得超過臨界時間步長Δtcr,一般情況下取

式中 α 是安全系數,碎點法參考有限元的安全系數取值,取0.8 ≤α ≤0.98,臨界時間步長的計算公式如下

式中le代表碎點法模型中子域的特征長度,c是線彈性材料的絕熱聲速.

由于碎點法的弱形式方程引入了數值通量修正,用于保證方法的一致性,而數值通量修正項中包含懲罰系數 ηh,如式(15)所示.數值通量會使得整個系統產生一定的數值振蕩,而數值通量產生的作用力的大小取決于懲罰系數 ηh,因而懲罰系數會對臨界時間步長產生影響,較小的時間步長可以逐漸消除數值通量帶來的振蕩的影響.

對碎點法中中心差分法的臨界時間步長計算公式進行修正得到

使用上式對碎點法模型中的每一個子域進行計算,使用最小的臨界時間步長作為整體模型的臨界時間步長.

2 碎點法動力學算例

2.1 拉力作用下的平板應力波傳播

為了驗證動力學碎點法模擬應力波傳播的能力,首先應用動力學碎點法程序求解如圖5 所示平板的應力波傳播問題.參考文獻中的設置[40],該算例中使用矩形薄板作為研究對象,矩形板的面內尺寸為長L=4 m,寬D=2 m.平板材料屬性為,楊氏模量E=80kPa,泊松比 μ=0,而材料密度 ρ=1 kg/m3.邊界條件如圖5 所示,矩形板左端固支,右端受均布拉伸載荷P(t),載荷大小隨試件的變化關系如圖6所示,仿真的總時間設置為0.3 s.由于分析對象為薄板,因此可將該問題等效為平面應力問題,碎點法模型中使用40×20個均勻分布的節點對矩形板進行離散.由于該問題中材料的泊松比設置為零,因此本質上這是一個一維應力波傳播問題.

圖5 拉力作用下的矩形平板模型Fig.5 The model for rectangular plate under tension

圖6 載荷 P(t) 隨時間的關系Fig.6 The curve for the load-time correlation

對于此問題,文獻[40]給出了自由端A點的水平位移uA、中點B點的水平位移位移uB、中點B點的水平正應力、以及固定端C點的水平正應力的精確解.本文以上述數據為參考對碎點法的計算結果進行驗證,同時與有限元結果進行對照.如圖7 所示,在碎點法模型中的相同位置讀取位移和應力信息,峰值處的相對誤差如表1和表2 所示.碎點法的仿真結果與有限元及理論解吻合良好,從而驗證了本文所提出的動力學碎點法在模擬一維應力波傳播問題方面的有效性.

圖7 拉力作用下平板應力波傳播問題精確解及FPM 結果Fig.7 Exact solutions and FPM’s results of the displacement and stress at different points

圖7 拉力作用下平板應力波傳播問題精確解及FPM 結果(續)Fig.7 Exact solutions and FPM’s results of the displacement and stress at different points(continued)

表1 有限元、碎點法的位移峰值與理論解的比較Table 1 Relative error of maximum displacement obtained from FEM and FPM

表2 有限元、碎點法的應力峰值與理論解的比較Table 2 Relative error of maximum stress obtained from FEM and FPM

2.2 自由端剪切載荷作用下懸臂梁動態響應

接著使用動力學碎點法程序對自由端剪切載荷作用下的懸臂梁動態響應進行仿真,通過這個算例來驗證所提出的方法在模擬二維應力波傳播問題方面的能力.該算例中使用圖8(a)所示的懸臂梁作為研究對象,梁的尺寸為,長度L=48 m,高度H=12 m.材料的彈性模量E=30MPa,泊松比為0.3,密度ρ=1.0kg/m3.邊界條件設置如圖8(a)所示,懸臂梁左端面固支,右端面受向下的瞬態剪切載荷P=100Pa,載荷作用時間為0～0.5 s,總仿真時間為 2 s.本算例同樣使用平面應力假設,碎點法模型使用 4 8×12 個均勻分布的節點對懸臂梁進行離散,如圖8(b).從模型中隨機選取一個節點A點,用于讀取位移和應力仿真結果.本算例中所設置的幾何、材料參數僅用于算法驗證,不代表任何具體的物理對象.

圖8 (a)懸臂梁示意圖和(b)計算模型節點分布Fig.8 (a) Configuration of the beam and(b) the scatter of subdomains

同時,還在商業有限元軟件ABAQUS 中建立了這個算例的有限元模型,用有限元模擬結果對碎點法程序進行驗證.有限元模型的網格與碎點法模型的子域完全一致.

在圖9～ 圖11 中分別展示了懸臂梁加載端(右端)中點的豎直位移、固定端(左端)中點各應力分量和A點水平正應力的有限元方法FEM和碎點法計算結果.圖11 中僅繪制A點的水平正應力曲線是因為該點處的其他應力分量基本為零.可以看出,碎點法預測的應力和位移狀態與有限元結果吻合良好,兩種方法得到的應力、位移響應頻率特征一致.各圖中,兩種方法得到的位移曲線幾乎完全重合,應力曲線峰值的相對誤差在 2.5% 以內,其余部分幾乎完全重合.由此可見,本文所提出的動力學碎點法程序可以有效模擬二維應力波傳播問題.

圖9 自由端中點豎直方向位移響應曲線Fig.9 Vertical displacement curve of the mid-point at the free end of the beam

圖10 固定端中點應力響應曲線Fig.10 Curves of stress components of the mid-point at the fixed end of the beam

圖11 懸臂梁A 點處應力瞬態響應曲線Fig.11 Stress-time curve of point A

2.3 含裂紋金屬板低速沖擊破壞行為仿真

最后,使用動力學碎點法程序對一個經典動態斷裂問題進行模擬仿真,以驗證該方法模擬裂紋動態擴展的能力.該算例基于 Kalthoff等[41-42]在1987 年和2000年所做的一系列試驗工作,試驗試件如圖12(a)所示,為長和寬分別為100mm和200mm的矩形金屬板,板上包含著兩條對稱的水平初始裂紋,裂紋間距為50mm,初始裂紋長度為50mm.試驗使用一塊底面半徑為25 mm 的圓柱形子彈以33 m/s的速度沖擊兩條裂紋之間的區域[42].試驗觀察到含裂紋金屬板在子彈沖擊下的主要失效模式為脆性斷裂,裂紋擴展路徑與初始裂紋的夾角約為70°,如圖12(b)所示.

圖12 (a)含裂紋金屬板沖擊試驗裝置示意圖和(b)試驗觀測到的裂紋擴展路徑[43]Fig.12 (a) Experimental setup for the Kalthoff plate impact test and(b) the experimentally-observed crack path[43]

根據該問題的對稱性,建立二分之一模型以減少計算量,計算模型如圖13(a)所示,圖中紅色線條代表初始裂紋.文獻[33]研究了碎點法節點分布方式給計算結果帶來的影響,由于數值通量修正的引入,碎點法使用均勻分布節點或隨機分布節點均可以得到準確的結果,驗證了碎點法的魯棒性.基于該特點及試驗中觀測到的裂紋主要擴展區域,在布置節點時,矩形板右上區域布置的子域節點更加密集,在其他區域子域節點稍微稀疏.這種方式一方面可以更精細地捕捉裂紋的擴展路徑,另一方面可以節省計算資源.本文的工作中,鄰近子域節點群的選取為所有與該子域共邊界的子域節點,并未因節點分布密度不同而進行調整.

試驗中金屬板由型號18 Ni1900的鋼材制成,碎點法仿真中采用文獻[44]中的材料參數,楊氏模量E=190GPa,泊松比 μ=0.3,密度 ρ=8 t/m3.根據文獻[45],取內界面破壞的拉伸強度為1773 MPa.同時,還建立了使用更加致密節點的碎點法模型,如圖13(b)所示,用于分析網格密度對仿真結果的影響,并驗證圖13(a)所示的有效性.

圖13 計算模型:(a) 6516 個子域節點和(b) 11 995 個子域節點Fig.13 Domain discretization with(a) 6516 FPM subdomains and(b)11 995 FPM subdomains

由于沖擊載荷作用時間短,可認為子彈速度在沖擊過程中幾乎沒有降低,因此可以在沖擊區域用速度邊界條件等效代替子彈的沖擊載荷.根據文獻[45-46],試驗中所使用的子彈和試件由相同材料制備而成,二者擁有相同的彈性阻抗,因此加載速度取試驗中子彈速度的一半 1 6.5 m/s.

使用碎點法仿真計算得到的靜水壓力云圖如圖14 所示,可以看到圖14(a)中,在 8 μs 時,由于矩形板左側載荷的作用,應力波先是沿著矩形板下方初始裂紋與對稱面之間的區域,而矩形板其他區域暫時未受到擾動,相應的靜水壓力均為零.到12 μs時,裂紋尖端處靜水壓力較大,矩形板其他位置的靜水壓力幾乎為零,如圖14(b)所示,這是因為初始裂紋的存在導致了應力波逐漸在初始裂紋尖端處聚集,產生應力集中的現象,仿真的結果與實際的物理機制吻合.

由于裂紋尖端應力集中,裂紋尖端區域的內邊界達到強度而開始萌生裂紋,并沿著與初始裂紋夾角約70°角的方向擴展,如圖14(c)所示,此時應力集中的位置也隨之發生變化,應力集中現象發生在新的裂紋尖端處,而已生成裂紋的區域應力逐漸減小,仿真的結果吻合物理機制.

圖14 靜水壓云圖Fig.14 Configuration with hydrostatic pressure

圖14 靜水壓云圖(續)Fig.14 Configuration with hydrostatic pressure(continued)

不同節點分布的仿真的裂紋擴展路徑在圖15中畫出.對于仿真結果1,在最初時刻,裂紋以大約70°角開始擴展.隨著裂紋擴展到水平位置60mm左右,夾角出現偏轉,此時角度約為64°;在水平位置80mm 處左右,裂紋近似沿垂直方向擴展,直至上邊緣附近,又重新以約70°角擴展;從初始裂紋尖端到裂紋最終位置的平均角度約為66°.而對于仿真結果2,裂紋幾乎一致沿著約70°的角度擴展,僅在接近板上邊緣附近出現了角度增大的情況.總體而言,兩種不同節點密度的模型預測的裂紋擴展形貌都和試驗觀測基本一致.其中模型2 使用更致密的節點,因此仿真結果與試驗結果吻合更加良好;模型1 的自由度顯著少于模型2,計算效率更高,且能夠給出關鍵性的應力傳播和裂紋擴展特性,也滿足動態裂紋擴展預測和分析的需求.

圖15 FPM 仿真裂紋擴展路徑Fig.15 Crack propagation from FPM simulation

下面繪制裂紋擴展的速度并與文獻[45]的仿真結果對比,討論碎點法仿真得到的裂紋擴展速度的響應特征.根據瑞利波速計算公式

式中cs為剪切波傳播速度,可以計算得到瑞利波速為cR=2799 m/s.

對裂紋長度隨時間的變化數據進行差分求導處理可求得裂紋擴展速度,如圖16 所示,裂紋擴展速度在25 μs 左右時開始增長,之后趨于穩定,在56 μs左右開始出現下降.碎點法得到的裂紋擴展速度略低于文獻仿真結果,但總體趨勢與文獻[45]中數據吻合良好.

圖16 裂紋擴展速度-時間曲線Fig.16 Curve of crack propagation speed versus time

以上算例表明,本文所提出的動力學碎點法能夠模擬典型的動態裂紋擴展行為,仿真得到的裂紋形貌和裂紋擴展速度與試驗觀測結果吻合良好.考慮到碎點法可以在任意內界面處判斷界面狀態,當界面破壞時只需要簡單的修正就可以顯式引入裂紋,且引入裂紋不會增加模型整體的自由度,因此動力學碎點法具備模擬大規模復雜裂紋萌生和動態擴展問題的能力.

3 結論及展望

本文根據碎點法的基本思想和原理,推導了碎點法弱形式動量方程、建立了顯式動力學碎點法求解格式并編制了相應的計算程序,還使用經典算例驗證了動力學碎點法模擬應力波傳播和動態裂紋擴展問題的能力.

碎點法具有如下特色:(1)碎點法是一種基于伽遼金弱形式的無網格法,具有較好的算法穩定性;(2)采用基于節點的位移形函數,因此具有抵抗子域畸變的能力;(3)碎點法的形函數是分片連續的,因此可以使用多項式形式的形函數,并使用簡單的高斯積分進行求解;(4)在弱形式方程中引入數值通量修正項,從而保證放棄了內界面位移連續性的碎點法的一致性和穩定性;(5)內界面數值通量具有相鄰子域間相互作用面力的物理含義,因此可用于定義內界面的斷裂準則;(6)只需要簡單的操作就可以在碎點法模型中顯式引入裂紋,即移除對應內界面的數值通量修正,并修改界面兩端子域的支持域,從而分別相鄰子域間數值通量和形函數兩方面的相互聯系;(7)子域節點具有質量、體積、速度等明確的物理屬性,便于顯式動力學計算.

隨后通過三個典型算例驗證了動力學碎點法對應力波傳播及動態裂紋預測的預測能力.前兩個算例驗證了碎點法模擬拉伸和剪切應力波傳播的有效性和準確性,所得結果與精確解和有限元結果吻合良好.第三個算例采用經典動態斷裂算例來驗證動力學碎點法模擬裂紋動態擴展過程的能力,計算結果與文獻中的試驗結果吻合良好,體現了碎點法模擬動態斷裂的問題的有效性.本文所提出的動力學碎點法理論及程序為動態斷裂問題的研究提供了一種簡單、高效的數值仿真方法.

最后需要指出的是,本文并未涉及如裂紋分叉、交匯等問題,但碎點法易于顯式引入裂紋的特點使得其在這類問題的仿真中具有很大潛力.這些問題背后的復雜物理機制意味著需要做進一步工作,如開發更加細化的界面模型,確定合適的斷裂準則,以及建立這些準則與界面數值通量的聯系等,這也將是作者團隊下一步的研究方向.