三維裂紋擴(kuò)展時(shí)動(dòng)態(tài)應(yīng)力強(qiáng)度因子的有限元研究

王光軍 熊峰

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

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

1 三維裂紋擴(kuò)展時(shí)動(dòng)態(tài)應(yīng)力強(qiáng)度因子有限元分析

1.1 分析思路

對于線性彈性均勻介質(zhì)，當(dāng)載荷隨時(shí)間變化，穩(wěn)定裂紋頂端的漸近位移場與靜態(tài)情況完全類似。對于平面I型問題

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

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

1.2 三維裂紋擴(kuò)展基本方程

動(dòng)態(tài)隨機(jī)載荷導(dǎo)致裂紋的隨機(jī)擴(kuò)展，如果將構(gòu)件的動(dòng)態(tài)干擾理解為圍繞裂紋尖端周圍材料的運(yùn)動(dòng)[2]，則其邊界條件為：

同時(shí)滿足兩個(gè)剪切方程非零：

構(gòu)件不受力的情況下，裂紋擴(kuò)展方程可寫為：

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

1.3 構(gòu)件裂紋擴(kuò)展數(shù)學(xué)模型

將所研究的構(gòu)件進(jìn)行有限元離散化處理，其運(yùn)動(dòng)擴(kuò)展方程的矩陣形式[3]如下:

其中：[M]、[C]、[K]分別為結(jié)構(gòu)的總質(zhì)量矩陣、結(jié)構(gòu)的阻尼矩陣、結(jié)構(gòu)的總剛度矩陣；{u}為結(jié)構(gòu)的節(jié)點(diǎn)位移列陣；{F}為節(jié)點(diǎn)等效載荷列陣；{}為節(jié)點(diǎn)位移二階導(dǎo)數(shù)列陣，即加速度列陣；{}為節(jié)點(diǎn)位移的一階導(dǎo)數(shù)列陣。

求解線彈性有限元?jiǎng)恿Ψ匠炭刹捎肗ewmark逐步積分法，得到穩(wěn)定結(jié)果。在時(shí)刻t+Δt由運(yùn)動(dòng)方程得

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

其中α、β為控制算法精度和穩(wěn)定的兩個(gè)參數(shù)。

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

再根據(jù)式(9)— (11)式，求出下一時(shí)刻△t的{uΔt}，{Δt}，{Δt}，由此即可得到所有時(shí)間離散點(diǎn)上的位移、速度與加速度，進(jìn)而可以求得各個(gè)時(shí)刻的應(yīng)力、應(yīng)變和應(yīng)力強(qiáng)度因子。

1.4 加載時(shí)間步長的確定

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

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

1.5 單元形式的選取和裂紋擴(kuò)展運(yùn)動(dòng)的求解

本文采用APDL參數(shù)化語言編制三維裂紋模型，同時(shí)結(jié)合 ANSYS子模型技術(shù)進(jìn)行分析[5]。由于裂紋頂端區(qū)域應(yīng)變場具有r-1/2階的奇異性，為此在裂紋頂端采用四分之一中點(diǎn)奇異等參元。對于線性彈性問題，這種單元構(gòu)造自動(dòng)生成r-1/2階的奇異性。其他部位采用標(biāo)準(zhǔn)八節(jié)點(diǎn)等參元。

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

本文按照模態(tài)疊加法步驟計(jì)算所有時(shí)間離散點(diǎn)上的位移、應(yīng)力等信息，再利用軟件求解靜態(tài)應(yīng)力強(qiáng)度因子，計(jì)算出某個(gè)時(shí)間點(diǎn)上的應(yīng)力強(qiáng)度因子。從而可以得到應(yīng)力強(qiáng)度因子和時(shí)間的動(dòng)態(tài)關(guān)系。

2 算例分析

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

圖1 圓裂紋示意圖

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

圖3 裂紋局部放大圖

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

表1 材料及模型參數(shù)

本文采用突加載荷，如圖4，其中σ0（t）=100 MPa。對不同裂紋長度的動(dòng)態(tài)應(yīng)力強(qiáng)度因子的有限元計(jì)算結(jié)果列入表2中。圖5所示為裂紋尖端應(yīng)力場，典型的應(yīng)力集中現(xiàn)象。

圖4 載荷示意圖

表2 突加載荷動(dòng)態(tài)強(qiáng)度因子的結(jié)果

圖5 裂紋尖端應(yīng)力場

圖6 載荷作用下裂紋尖端擴(kuò)展圖

圖6為載荷作用下裂紋尖端位移擴(kuò)展過程。最后利用各個(gè)時(shí)間點(diǎn)和相應(yīng)的應(yīng)力強(qiáng)度因子作圖，得出動(dòng)態(tài)應(yīng)力強(qiáng)度因子隨時(shí)間的變動(dòng)圖，如圖7。

圖7 動(dòng)態(tài)應(yīng)力強(qiáng)度因子圖

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

3 結(jié)束語

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

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

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

[2] 范天佑. 斷裂動(dòng)力學(xué)原理與應(yīng)用[M]. 北京:北京理工大學(xué)出版社，2006

[3] 王助成，邵敏. 有限單元法基本原理和數(shù)值方法[M].第二版. 北京: 清華大學(xué)出版社, 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] 中國航空研究院. 應(yīng)力強(qiáng)度因子手冊. 北京:科學(xué)出版社，1993.