靜態漸近解在動態斷裂問題中的適用性分析

諶赫， 鄒廣平

(哈爾濱工程大學 航天與建筑工程學院，黑龍江 哈爾濱 150001)

材料動態斷裂實驗的基本思路是測量載荷、加載點位移等物理量來計算動態應力強度因子(dynamic stress intensity factor，DSIF)。較早的研究工作主要是測量試樣在加載過程中的載荷、位移，代入相應的經驗公式計算動態應力強度因子[1-4]；隨著光學測量技術的發展，裂紋擴展過程也可以被觀察到[5-7]，甚至可以得到試樣表面任一時刻的應變分布[8-9]；應變片法在動態斷裂測試中也得到廣泛應用[10-11]，常用來標定其他方法；近年來，部分學者將實驗與數值模擬相結合，發展了實驗-數值法[12-14]。這種方法通過實驗測得載荷、起裂時間等參數，采用有限元軟件計算動態應力強度因子，進而確定材料動態斷裂韌性。相對于靜態斷裂力學，動態斷裂力學需要考慮材料的應變率效應與慣性效應。由于數學上的困難，人們對動態斷裂的研究尚不深入，動態斷裂力學的理論框架并不完備[15]。對于動態載荷作用下裂紋起裂的問題，主流觀點認為裂紋尖端應力應變場的形式與靜態問題相似，具有分離變量形式的解[16]。該解只有在接近裂紋尖端的區域才能成立，故稱漸近解。漸近解是聯系宏觀量與微觀量的橋梁，前面提到的各種動態應力強度因子測定與求解方法均以靜態漸近解為基礎。裂紋尖端塑性區的存在嚴重地限制了漸近解的應用。此外，對于不同構型的試樣，其載荷工況與邊界條件各不相同，邊界對裂尖應力應變場的影響也不能一概而論。

本文基于分離式Hopkinson拉桿(split-Hopkinson tension bar)實驗裝置，提出一種改進的緊湊拉伸剪切(compact tension shear specimen loaded,MCTS)試樣，用于II型動態斷裂測試[17]。本文采用有限元軟件ABAQUS對改進的緊湊拉伸剪切試樣在I型與II型載荷工況下進行了數值模擬，通過分析試樣裂尖應變場與漸近解的差異來討論靜態漸近解在動態問題中的適用性。

1 裂尖應力應變場漸近解

裂尖應力應變場的漸近解的表達式為：

(1)

平面應力狀態下，式(1)中fij(θ)、gij(θ)的表達式分別為：

(2)

(3)

漸近解滿足靜態問題的平衡方程σij,j=0，而對于動態問題慣性效應必須要考慮，平衡方程化為σij,j=ρui,tt，二者形式明顯不同，那么靜態問題的漸近解應用于動態問題的依據何在？

Sun[18]指出，裂紋尖端區域應力具有奇異性，應力對坐標的導數也是奇異的，而位移是有限的，因此平衡方程左邊遠遠大于右邊，基于上述假設，在裂尖附近將靜態問題的漸近解應用于動態問題是合理的；而在距離裂尖較遠處漸近解不適用。文獻[18]雖然給出定性的解釋，但并未提及裂尖“附近”的具體范圍。

2 改進的緊湊拉伸剪切試樣數值模擬

改進的緊湊拉伸剪切試樣的幾何構型與尺寸如圖1所示，通過特殊設計的夾具能夠實現I型與II型加載，試樣與分離式Hopkinson拉桿裝置的裝配圖見圖2。II型加載時需要在試樣側面施加垂直于側面的位移約束。

圖1 改進的緊湊拉伸剪切試樣幾何尺寸Fig.1 Geometry of MCTS specimen

圖2 改進的緊湊拉伸剪切試樣裝配圖Fig.2 Assembly of MCTS specimen

改進的緊湊拉伸剪切試樣與分離式Hopkinson拉桿裝置的材料參數如表1所示，其中采用雙線性隨動強化模型來表征試樣的塑性，分離式Hopkinson拉桿裝置的其他部分采用線彈性模型。

表1 有限元模型材料屬性Table 1 Material properties of finite element models

將沖擊拉伸試驗測得的波形進行濾波處理，作為入射桿端部施加的載荷[19]，如圖3所示；試樣裂尖網格劃分如圖4所示，在裂尖采用三棱柱奇異單元，半徑1 mm，圍繞裂尖32等分；其他部分采用六面體單元，寬度0.25 mm。

圖3 SHTB入射波Fig.3 Incident wave of SHTB

圖4 改進的緊湊拉伸剪切試樣裂紋尖端網格劃分Fig.4 Mesh at crack tip of MCTS specimen

3 動態應力強度因子計算

動態應力強度因子的定義為：

(4)

式(4)中坐標均為裂尖局部坐標。對于II型載荷工況，局部坐標與全局坐標相同，對于I型載荷工況，局部坐標相對于全局坐標順時針旋轉90°。

改進的緊湊拉伸剪切試樣裂尖塑性區如圖5所示。圖中標示出來的節點周圍單元塑性應變均為零，且距離裂尖最近。定義為塑性影響區的邊界，記作rp。由圖可見，相同的入射波作用下，I型裂尖塑性區rp的最大值是II型的2倍。

圖5 裂紋尖端塑性區Fig.5 Plastic zone at crack tip

根據定義式(4)，計算出裂紋延長線上各節點每一時刻的“應力強度”并線性擬合，其截距為這一時刻的動態應力強度因子數值解[20]。數據處理過程如下：從塑性影響區邊界開始，選取相鄰的7個節點進行線性擬合，依次重復此步驟，當截距變化量小于1%時記錄截距值，再重復3次，取截距的平均值作為動態應力強度因子值。計算出動態應力強度因子數值解如圖6所示。

圖6 動態應力強度因子數值解Fig.6 Numerical solution of DSIF

4 裂尖應變場分析

根據定義已經求得動態應力強度因子的數值解，將其代入靜態漸近解的表達式，任一節點應變都可求出。將式(2)寫成矩陣形式：

(5)

Rittel[21]提出一種雙應變片求解動態應力強度因子的實驗方法，其中r1、θ1、r2、θ2分別為應變片1、2在裂尖局部坐標系內的坐標，反解式(5)即可求出動態應力強度因子。本文將通過式(5)計算的值稱為漸近解，通過有限元分析得到的解稱為數值解，二者之間的關系如圖7所示。可以看出，若節點應變漸近解與數值解不相同，則動態應力強度因子漸近解與數值解也會存在誤差。本節在不同方向上選取節點，分析其應變漸近解與數值解之間的差異。

圖7 應變分析流程Fig.7 Flow graph of strain analysis

實驗結果表明，裂紋往往在動態應力強度因子達到極大值之前開始擴展[12]，在這個時間段內討論應變漸近解與數值解的誤差是有意義的。本文討論的范圍從動態應力強度因子達到0.1倍最大值開始，達到最大值為止。限于篇幅，文中只列出±90°、±67.5°、±45°、±22.5°和0°方向上節點應變漸近解與數值解在討論范圍內絕對誤差的最大值、平均值以及絕對誤差最大值與應變分量最大值之比(記為相對誤差，在括號內標注)。I型與II型2種載荷工況分別如表2、表3所示。

表2 節點應變誤差(I型載荷工況)Table 2 Error of nodal strain (mode I loading)

表3 節點應變誤差(II型載荷工況)Table 3 Error of nodal strain (mode II loading)

ε11、ε22在θ=0°及θ=±22.5°方向上相對誤差平均值小于20%。ε11誤差最小的節點位于θ=0°，r=2.52 mm處，編號為1 023；除該節點外，ε22誤差最小的節點位于θ=-22.5°，r=3.02 mm處，編號為1 051。其應變曲線分別如圖8(a)、(b)所示。將節點1 023的應變分量ε11，節點1 051的應變分量ε22代入式(5)反解出動態應力強度因子，并與數值解對比。如圖8(c)所示。可見KΙ漸近解與數值解誤差很小，而KΙΙ誤差很大。由于KΙΙ數值解約等于零，且系數矩陣存在零元素，故KΙΙ漸近解明顯偏離數值解。

圖8 I型載荷工況節點應變及動態應力強度因子漸近解與數值解對比Fig.8 Nodal strain and comparison of asymptotic solution and numerical solution of DSIF for mode I loading

對于II型載荷工況，應變誤差相對I型載荷工況更小。ε11在θ=-45°、θ=-67.5°方向上相對誤差小于10%；ε22在θ=45°方向上相對誤差小于10%。ε11相對誤差最小的節點位于θ=-67.5°，r=3.52 mm處，編號為1109；ε22誤差最小的節點位于θ=45°，r=2.77 mm處，編號為962。其應變曲線分別如圖9(a)、(b)所示。將節點1109的應變分量ε11，節點962的應變分量ε22代入式(5)，解出動態應力強度因子如圖9(c)所示。可見動態應力強度因子漸近解與數值解幾乎完全重合。

圖9 II型載荷工況節點應變及DSIF漸近解與數值解對比Fig.9 Nodal strain and comparison of asymptotic solution and numerical solution of DSIF for mode II loading

考察各方向上各節點應變絕對誤差最大值與應變分量最大值之比，取其平均值如表4。可見多數方向上節點應變分量漸近解與數值解相差較大，一些方向上相差極大，其原因是這些方向上相應的節點應變分量漸近解或數值解約等于零。

表4 各方向節點應變相對誤差平均值Table 4 Average error of nodal strain in different orientations

對于I型載荷工況，由于塑性影響區尺度較大，所選節點距離裂尖較遠，受邊界影響更大，因此各方向上節點應變誤差相對II型載荷工況較大。對于II型載荷工況，裂尖塑性影響區尺度較小，應變誤差相對較小。

圖10顯示了2種載荷工況下節點應變漸近解與數值解相對誤差小于15%的區域，稱為應變理想區。應變理想區位于rp～2rp。對于I型載荷工況，ε11與ε22的理想區重合，位于θ=0°附近；而對于II型載荷工況，ε11與ε22的理想區分別位于θ=-45°與θ=45°附近。可見靜態漸近解適用的區域非常有限。

圖10 應變理想區Fig.10 Suitable zone of strain solution

在應變理想區內選取4組節點計算動態應力強度因子并分析漸近解與數值解的誤差，2種載荷工況分別如表5所示。可見在節點應變相對誤差不大的條件下，動態應力強度因子相對誤差與節點應變相對誤差量級相同；限于篇幅文中不再列出動態應力強度因子曲線。

表5 動態應力強度因子漸近解與數值解的誤差Table 5 Error of DSIFs between asymptotic and numerical solutions

5 結論

1)裂尖應變場的解析解在特定方向上與數值解較為接近(Δε≤10%)，而在其他方向上有較大差異。對于I型載荷工況，誤差較小的方向是0°附近；對于II型載荷工況則是±45°附近。

2)DSIF解析解與數值解的相對誤差取決于節點應變的相對誤差，2者在相同數量級。

3)靜態漸近解在動態斷裂問題中的適用范圍非常有限，且與試樣載荷工況有關。應針對具體問題分別討論。