王光軍 熊峰
(中國船舶重工集團(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)。
對于線性彈性均勻介質(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)度因子。
動(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ì)量密度。
將所研究的構(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)度因子。
動(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é)果。
本文采用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)系。
為了模型的普遍性,采用斷裂力學(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é)果精確。
(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.