三維裂紋擴展時動態應力強度因子的有限元研究

王光軍 熊峰

（中國船舶重工集團公司712研究所，武漢 430064）

工程中很多構件受動載荷作用，例如：柴油機氣缸內氣體燃燒時受動態內壓作用，氣缸內腔經氣體燃燒膨脹將出現龜裂，繼續燃燒膨脹裂紋將不斷擴展，可能產生疲勞破壞等等。實際上，三維裂紋擴展時動態應力強度因子的求解是一個較為困難的問題。由于增加了一個時間變量，不僅在數學處理上困難，在物理上也很復雜。對于三維問題應力強度因子，目前尚不能獲得精確的解析解。而三維問題的數值分析，由于裂紋尖端奇異性的復雜性，以及計算工作量大等因素，使研究受到一定的限制。本文主要利用有限元軟件ANSYS計算三維裂紋擴展時動態應力強度因子[1]，能為動態應力強度因子的計算、分析提供幫助和指導。

1 三維裂紋擴展時動態應力強度因子有限元分析

1.1 分析思路

對于線性彈性均勻介質，當載荷隨時間變化，穩定裂紋頂端的漸近位移場與靜態情況完全類似。對于平面I型問題

對于動態載荷下靜止裂紋，裂紋尖端應力場和位移場公式與靜態載荷下的類似。通過裂紋面上各點處垂直于裂紋面的位移，在平面應變情況下，由式(1)，式(2)，式(3)可得

(4)式中位移μy(t)為相應時刻裂紋面上各點處垂直于裂紋面的位移，可由上述的有限元法求得，外推到裂紋尖端求出r=0處的K1(t)值，即為所求動態應力強度因子。

1.2 三維裂紋擴展基本方程

動態隨機載荷導致裂紋的隨機擴展，如果將構件的動態干擾理解為圍繞裂紋尖端周圍材料的運動[2]，則其邊界條件為：

同時滿足兩個剪切方程非零：

構件不受力的情況下，裂紋擴展方程可寫為：

其中：c2=（μ/ρ）1/2是剪切速度，μ是剪切系數，ρ是物體的質量密度。

1.3 構件裂紋擴展數學模型

將所研究的構件進行有限元離散化處理，其運動擴展方程的矩陣形式[3]如下:

其中：[M]、[C]、[K]分別為結構的總質量矩陣、結構的阻尼矩陣、結構的總剛度矩陣；{u}為結構的節點位移列陣；{F}為節點等效載荷列陣；{}為節點位移二階導數列陣，即加速度列陣；{}為節點位移的一階導數列陣。

求解線彈性有限元動力方程可采用Newmark逐步積分法，得到穩定結果。在時刻t+Δt由運動方程得

速度與位移由（10）式給出

其中α、β為控制算法精度和穩定的兩個參數。

如果t=0時的初始位移和初始速度為{u0}與{0}，則由式(8)求得初始加速度

再根據式(9)— (11)式，求出下一時刻△t的{uΔt}，{Δt}，{Δt}，由此即可得到所有時間離散點上的位移、速度與加速度，進而可以求得各個時刻的應力、應變和應力強度因子。

1.4 加載時間步長的確定

動態有限元程序，是通過逐步時間積分來解動力方程，需要選定合適的加載時間步長，加載時間步長△t選得過小，將增加計算時間，加載時間步長△t選得過大將影響計算精度，本文參考文獻[4]中所建議的公式確定時間步長。

△t≤ Δl/C，△t為單元最小尺寸，c為最大縱波波速。參考t△ ≤ Δl/C可得對不同裂紋長度計算取△t=1.0×10-5～3.0×10-4s時能得到較好的結果。

1.5 單元形式的選取和裂紋擴展運動的求解

本文采用APDL參數化語言編制三維裂紋模型，同時結合 ANSYS子模型技術進行分析[5]。由于裂紋頂端區域應變場具有r-1/2階的奇異性，為此在裂紋頂端采用四分之一中點奇異等參元。對于線性彈性問題，這種單元構造自動生成r-1/2階的奇異性。其他部位采用標準八節點等參元。

求解裂紋擴展運動的主要采用模態疊加法。有限元法是計算動態K1(t)因子的有效方法，但求解動態問題需逐步加載計算，時間步長很小，很費機時。本文引用疊加原理，將每一種動態載荷曲線分解為很多微小的分級加載疊加組成，同時由于用APDL參數化語言編制的三維裂紋程序中可以靈活地修改相關參數，這樣大大減少了有限元法的計算量。

本文按照模態疊加法步驟計算所有時間離散點上的位移、應力等信息，再利用軟件求解靜態應力強度因子，計算出某個時間點上的應力強度因子。從而可以得到應力強度因子和時間的動態關系。

2 算例分析

為了模型的普遍性，采用斷裂力學中常見的有限尺寸圓柱體中深埋圓裂紋模型[6，7]為研究對象。圓裂紋示意如圖1。

圖1 圓裂紋示意圖

圖2 八分之一圓柱體的中心裂紋模型

圖3 裂紋局部放大圖

考慮到對稱結構，建立八分之一圓柱體的中心裂紋模型，如圖2，裂紋局部放大如圖3。取有限立方的八分之一建立模型的尺寸Ro；裂紋尺寸位于a處，圓柱高度為h，材料及模型參數具體尺寸如表1。選擇單元MESH200和單元SOLID95,在單元屬性中將 MESH200設置為KEYOPT(1)=7。用ANSYS編制的APDL程序計算裂紋深度a分別為: 5 mm、9 mm、15 mm、23 mm、27 mm，時間分別為2 ms、2.5 ms、3 ms、3.5 ms、4 ms、5 ms時的應力強度因子。在裂紋尖端用退化三角形奇異等參元，其余地方用標準等參元。

表1 材料及模型參數

本文采用突加載荷，如圖4，其中σ0（t）=100 MPa。對不同裂紋長度的動態應力強度因子的有限元計算結果列入表2中。圖5所示為裂紋尖端應力場，典型的應力集中現象。

圖4 載荷示意圖

表2 突加載荷動態強度因子的結果

圖5 裂紋尖端應力場

圖6 載荷作用下裂紋尖端擴展圖

圖6為載荷作用下裂紋尖端位移擴展過程。最后利用各個時間點和相應的應力強度因子作圖，得出動態應力強度因子隨時間的變動圖，如圖7。

圖7 動態應力強度因子圖

本文模型的正確性和求解靜態應力強度因子的準確性可以通過查看《應力強度因子手冊》[8]對比得到。根據計算靜態應力強度因子的步驟，這樣可以計算出理論值為K1（t）的誤差。本文算出的應力強度因子與《應力強度因子手冊》中結果最大誤差為 4.7%，證明模型及網格劃分良好，計算結果精確。

3 結束語

(1) 利用 ANSYS建立三維模型的方法是多種多樣的，本文主要考慮的圓柱體中心圓裂紋的結構。還可以利用相同的單元和建模方法建立其他典型裂紋結構，如中心穿透裂紋，橢圓裂紋等。

(2) 通過與算例的比較，可以認為利用ANSYS的 APDL參數化語言編程建立三維裂紋模型是可行的，有限元法是計算動態K1(t)因子的有效方法，本文提供的合理計算方案，也可為其它結構件動態K1(t)因子的計算作參考。

[1] 陳傳堯. 疲勞與斷裂[M]. 武漢: 華中科技大學出版社，2001.

[2] 范天佑. 斷裂動力學原理與應用[M]. 北京:北京理工大學出版社，2006

[3] 王助成，邵敏. 有限單元法基本原理和數值方法[M].第二版. 北京: 清華大學出版社, 1997.

[4] V MURTI S VALLIAPPAN. Engineering Fracture Mechanics. 1986, 23(3): 585-614.

[5] ANSYS，Inc. ANSYS Coupled-Field Analysis Guide Release7.0.

[6] Chen Aijun. Study on Dynamic Stress Intensity Factors of Disk with a Radial Edge Crack Subjected to External Impulsive Pressure. Acta Mechanica Solida Sinica,Volume 20, Issue 1, March 2007, Pages 41-49.

[7] Yong-Dong Li, Kang Yong Lee. Dynamic Stress Intensity Factors of Two Collinear Mode-III Cracks Perpendicular to and on the Two Sides of a Bi-FGM Weak-discontinuous Interface. European Journal of Mechanics - A/Solids, Volume 27, Issue 5,September-October 2008, Pages 808-823

[8] 中國航空研究院. 應力強度因子手冊. 北京:科學出版社，1993.