基于三維有限元解的緊湊拉伸試樣應力強度因子計算公式

賈 旭，胡緒騰，宋迎東

(南京航空航天大學能源與動力學院，江蘇省航空動力系統重點實驗室，

機械結構力學及控制國家重點實驗室， 南京 210016)

?

基于三維有限元解的緊湊拉伸試樣應力強度因子計算公式

賈 旭，胡緒騰，宋迎東

(南京航空航天大學能源與動力學院，江蘇省航空動力系統重點實驗室，

機械結構力學及控制國家重點實驗室， 南京 210016)

摘要：對標準緊湊拉伸(CT)試樣的二維和三維應力強度因子有限元解進行了對比分析，基于三維有限元虛擬裂紋閉合法擬合了一種新的標準CT試樣的應力強度因子計算公式，并采用ANSYS軟件內嵌位移外推法進行了驗證。結果表明：基于二維分析所得的應力強度因子計算公式計算結果與實際三維CT試樣的具有較大的差異；在加載孔等效分布力一定的條件下，CT試樣裂紋前沿大部分區(qū)域的應力強度因子與中心點的相近，且與厚度無關；擬合得到的三維CT試樣裂紋前沿中心點的應力強度因子計算公式具有很高的精度，其計算結果與ANSYS軟件內嵌位移外推法的相對誤差在0.5%以內。

關鍵詞：三維分析；CT試樣；應力強度因子；有限元法

0引言

為提高航空發(fā)動機的安全性，20世紀80年代以來損傷容限設計已逐漸成為國內外航空發(fā)動機輪盤等關鍵件的重要設計準則之一[1-2]。損傷容限設計的關鍵是對以斷裂為主要破壞形式的關鍵件進行裂紋擴展分析，而裂紋擴展分析的基礎是應力強度因子的計算。

工程中普遍采用標準緊湊拉伸(CT)試樣進行試驗來建立材料裂紋擴展速率模型、確定材料斷裂韌性。1974年Newman[3]基于二維平面假設(平面應變假設)下的改進邊界配置法擬合得到了標準CT試樣在不同幾何配置下的應力強度因子計算公式。1976年Srawley[4]對此公式進行了改進，拓寬了標準CT試樣應力強度因子計算公式的適用范圍。目前ASTM E399-12e3《平面應變斷裂韌性測試標準》和ASTM E647-13ae1《疲勞裂紋擴展速率測試標準》仍在繼續(xù)使用Srawley改進的標準CT試樣應力強度因子計算公式。非標準試樣的裂紋擴展速率研究[5-6]也仍在使用應力強度因子手冊中的二維應力強度因子表達式。然而，二維平面假設畢竟存在一定的局限性，在實際的三維裂紋體中，即使是標準CT試樣中形狀規(guī)則、初始裂紋前緣呈直線、承受單一銷釘載荷的張開型穿透裂紋，其裂紋前沿應力場也呈現出不同的應力狀態(tài)。特別是裂紋與自由表面相交的區(qū)域，既不是平面應變狀態(tài)也不是平面應力狀態(tài)[7-8]。二維平面應力強度因子計算公式將整個裂紋前沿生硬地假設為平面應變或平面應力兩種極端情況，勢必對應力強度因子的計算造成不可避免的誤差。目前少有學者對三維標準CT試樣的應力強度因子計算公式進行研究，Towers[9]計算對比了裂紋前緣中心點、1/4厚度點、表面點的應力強度因子的不同；Neman[10]在平面應力狀態(tài)假設下研究了CT試樣的權函數。為了得到標準CT試樣三維條件下的應力強度因子計算公式，作者建立了二維平面應變狀態(tài)的CT試樣的有限元模型，通過不同方法計算了應力強度因子并與目前普遍采用的應力強度因子計算公式進行了對比；然后建立了CT試樣的三維有限元模型，得到了CT試樣三維應力強度因子分布；最后基于CT試樣的三維有限元解擬合得到了一種新的CT試樣應力強度因子表達式，并與ANSYS軟件內嵌的位移外推法進行了對比分析。

1CT試樣應力強度因子的二維有限元分析

目前在ASTM標準中廣泛采用的CT試樣應力強度因子計算公式為：

(1)

式中：K為應力強度因子；a為裂紋長度；W為CT試樣的寬度；B為CT試樣的厚度(平面狀態(tài)假設下B=1)；P為施加在加載銷釘上的力。

以TC4合金為試樣及銷釘材料(彈性模量E為109 GPa,泊松比ν為0.34)，建立了寬度為0.04 m的CT試樣的兩種二維平面應變有限元1/2模型。第一種(①)在裂紋尖端采用六節(jié)點奇異單元plane82網格劃分，單元個數2 837；第二種(②)在裂紋尖端采用八節(jié)點高階單元plane82分網，單元個數3 007，如圖1所示。圖1中包圍裂紋尖端的16個奇異單元的特征尺寸l為0.01W。為了說明所建立的有限元模型、加載方式、應力強度因子計算方法的準確性，分別采用了ANSYS軟件內嵌位移外推法[11](Displacement Extrapolation Method，簡稱DEXM)和虛擬裂紋閉合法[12](Virtual Crack Closure Technique，簡稱VCCT)計算了兩種不同加載方式下的應力強度因子，并分別與式(1)結果進行對比，見表1。耦合加載方法將加載孔上表面內所有節(jié)點與加載孔中心的質點進行剛性耦合后對質點進行集中力加載；接觸加載方法是通過建立銷釘模型，進行接觸求解。相對于耦合加載，接觸加載方式為非線性求解，雖然計算效率較低但最切合實際。表1中δ1，δ2，δ3分別為奇異單元劃分/耦合加載時DEXM、奇異單元劃分/接觸加載時DEXM、非奇異單元劃分/接觸加載時VCCT計算結果與式(1)的相對誤差。

圖1　CT試樣二維有限元模型Fig.1　2D finite element model about the CT specimen：(a) general graph and (b) local graph at crack front

從表1可以看出，奇異單元劃分/耦合加載時，DEXM計算結果與式(1)的相對誤差在±2%以內。奇異單元劃分/接觸加載下DEXM與非奇異元劃分/接觸加載VCCT計算結果相對式(1)的誤差較一致，均在±1%以內。這一結果不僅說明了接觸加載的方式更為合理并驗證了有限元模型和方法的正確性，也說明了以平面假設為基礎的式(1)在計算CT試樣二維平面K時具有足夠高的精度。

表1　CT試樣二維平面應變假設下的應力強度因子計算結果及相對誤差Tab.1　 Calculated stress intensity factors of CT specimens and relative errors under two-dimensional plane strain assumption

2CT試樣應力強度因子的三維有限元分析

三維模型CT試樣與銷釘的材料及參數與二維分析時相同。CT試樣的寬度W為0.04 m，寬厚比W/B分別為4，6，8，10，12，14，16，20。分三個區(qū)域對CT試樣的三維模型進行了劃分：靠近裂紋尖端采用五面體奇異單元或20節(jié)點等參元劃分，遠離裂紋尖端采用20節(jié)點等參元劃分，過渡區(qū)采用四面體等參元劃分，單元為solid 95，如圖2所示。第一種(①)網格類型的單元數目為12 210，第二種(②)網格類型的單元數目為16 194。為了說明CT試樣二維平面應力強度因子與三維模型的不同，首先以a/W為0.3繪制了不同寬厚比下裂紋前緣假設為平面應變與平面應力狀態(tài)下的應力強度因子分布，如圖3所示。其中平面應變狀態(tài)和平面應力狀態(tài)并非指三維裂紋尖端的應力狀態(tài)，而是指應力強度因子數值計算方法的前提條件。然后以a/W為0.3，W/B為10為例對比了三維應力強度因子與式(1)的相對誤差，見表2。其中坐標軸Z平行于裂紋線，Z/(B/2)=1表示CT試樣的中心位置,δ4為平面應變狀態(tài)下的結果與式(1)的相對誤差，δ5為平面應力狀態(tài)下的結果與式(1)的相對誤差。從圖3可以看出，以二維平面假設為基礎的經驗公式顯然與三維平面應變或平面應力狀態(tài)下的應力強度因子相差較大；同一裂紋尺寸、不同寬厚比下CT試樣的應力強度因子在裂紋體表面時各有差異；在靠近裂紋體中心時，逐漸匯聚到一個數值。

圖2　CT試樣的三維有限元模型Fig.2　The 3D finite element mode of CT specimen:(a)general graph and local graph and (b) mesh at crack front

表2　CT試樣三維有限元解與式(1)結果的對比Tab.2　Comparison of 3D finite element solutions and formula (1) results

圖3　a/W為0.3時不同厚度下的K分布Fig.3　Stress intensity factor distribution at differentthicknesses with a/W equal to 0.3

從表2可以發(fā)現，雖然三維CT試樣中心區(qū)域接近為平面應變狀態(tài)，但中心區(qū)域內的應力強度因子與式(1)的結果也有很大的差別，最大相對誤差達到10.36%。雖然三維CT試樣表面點為平面應力狀態(tài)，但與式(1)的結果也有很大的差別，最大相對誤差達到-16.12%。

3基于三維有限元解的CT試樣應力強度因子計算公式

由式(1)可以得到：

(2)

由式(2)中可以看出,P/(0.125W×2×B)為P等效在加載孔上的分布力，當加載孔上的等效分布力一定時，CT試樣的應力強度因子只與W和a有關，與試樣的厚度無關。為考證這一特性，并定性地描述不同厚度下三維CT試樣裂紋前沿的應力強度因子分布特性，選用了計算效率更高的耦合加載方式以及ANSYS軟件內嵌的DEXM法，在同一加載孔上等效分布力下，分別采用平面應變狀態(tài)計算了三維CT試樣W/B分別為4，6，8，10，12，14，16，18，20，a/W分別為0.2，0.3，0.4，0.5，0.6，0.7，0.8下的應力強度因子分布，如圖4所示。

圖4　三維CT試樣裂紋前沿的應力強度因子分布Fig.4　Stress intensity factor distribution at the crack frontof the 3D CT specimen

由圖4發(fā)現，對于同一裂紋尺寸、不同寬厚比的CT試樣，位于試樣中心區(qū)域的應力強度因子基本分布在同一直線上，直線區(qū)域約占整個厚度的2/3。由此可見，在加載孔等效分布力一定的條件下，CT試樣裂紋前沿大部分區(qū)域的應力強度因子與中心點的應力強度因子接近，并且與厚度無關。

基于以上結論，作者在同一加載孔等效分布力的基礎上，采用了更精確的接觸加載方式，并采用了對網格密度依賴性較低的非奇異單元VCCT法計算了CT試樣中心點的應力強度因子。以a/W為0.5，W/B為10為例，采用裂紋尖端不同單元尺寸(l)進行網格劃分驗證了應力強度因子計算的網格無關性，見表3。其中δ′表示表中后一項K值相對于前一項的誤差。

由表3可知，l/a為0.01的K值相對于0.02時的誤差僅有0.158%，0.005的相對于0.01時的誤差僅有0.079%，相對誤差較小并且隨l/a的減小而減小。說明當l/a為0.01時，有限元法已經達到預期收斂效果，其計算結果與網格尺寸無關，為了方便劃分網格和提高計算效率，選用l/a為0.01。

表3　網格無關性驗證Tab.3　Mesh-independent verification

隨后，采用既定劃分網格方式計算了不同裂紋尺寸下的應力強度因子，并與式(1)的進行了對比，如圖5所示。從圖5中可以看出三維CT試樣中心部位的應力強度因子隨裂紋尺寸的變化趨勢與式(1)的基本一致。

圖5　VCCT法計算結果與式(1)結果的對比Fig.5　Comparison of VCCT calculation results andformula (1) results

鑒于目前采用的應力強度因子計算公式不適用于三維CT試樣，為了建立一種準確的適用于三維CT試樣的應力強度因子計算公式，在式(1)的基礎上對VCCT法的計算結果進行了應力強度因子的擬合，結果如下：

(3)

為了驗證式(3)的精度，對比了VCCT法(擬合數據點)的結果和式(1)與式(3)的相對誤差，并采用接觸加載下的奇異單元DEXM法對式(3)進行了驗證，結果見表4。表中δ6為式(1)的結果相對式(3)的誤差，δ7為VCCT法(擬合數據點)的結果相對式(3)的誤差，δ8為DEXM數據相對式(3)的誤差。

從表4可以看出，式(1)結果與式(3)的相對誤差在-7.79%～-11.75%，式(1)結果比式(3)的小；式(3)結果與VCCT法的誤差在±0.04%以內，說明式(3)具有很高的擬合精度；DEXM法的結果與式(3)的相對誤差在±0.5%以內，說明新的CT試樣應力強度因子表達式(3)在計算三維CT試樣應力強度因子時具有足夠高的精度。

表4　三維CT試樣裂紋中心部位的應力強度因子的計算結果與相對誤差Tab.4　Calculated stress intensity factors and relative errors at the crack center of the 3D CT specimen

4結論

(1) 目前廣泛采用的應力強度因子計算公式是通過二維假設分析所得，其計算結果與實際三維CT試樣的有很大差異，最大誤差達到10%以上。

(2) 在加載孔等效分布力一定的條件下，CT試樣裂紋前沿大部分區(qū)域的應力強度因子與中心點的相近，并且與厚度無關。

(3) 擬合得到的三維CT試樣裂紋前沿中心點的應力強度因子計算公式具有很高的精度，相對于其他方法的誤差在0.5%以內。

參考文獻：

[1]原航空工業(yè)部第三零一研究所.航空發(fā)動機結構完整性指南:GJB/Z 101-97[S].北京:[出版者不詳],1997：36-44.

[2]呂文林，陳俊粵，田德義. 航空渦噴、渦扇發(fā)動機結構設計準則(研究報告)第二冊：輪盤[R]. 北京: 中國航空工業(yè)總公司發(fā)動機系統工程局，1997.

[3]NEWMAN J C. Stress analysis of compact specimens including the effects of pin loading[J].ASTM STP,1974，560:105-105.

[4]SRAWLEY J E. Wide range stress intensity factor expressions for ASTM method E399 standard fracture toughness specimens[J]. International Journal of Fracture, 1976，12(3): 475-476.

[5]李旭東，穆志韜，賈明明. 加載頻率對航空鋁合金腐蝕疲勞裂紋擴展速率的影響[J].機械工程材料,2014，38(7)：50-52.

[6]余圣甫, 王鐵琦, 楊其良, 等. ZG20MnSi鋼斷裂韌度和疲勞裂紋擴展速率試驗研究[J].機械工程材料,2005,29(6):17-19.

[7]GARCIA-MANRIQUE J, CAMAS D, LOPEZ-CRESPO P, et al. Stress intensity factor analysis of through thickness effects[J]. International Journal of Fatigue, 2013, 46: 58-66.

[8]BAZANT Z P, ESTENSSORO L F. Surface singularity and crack propagation[J]. International Journal of Solids and Structures, 1979, 15(5): 405-426.

[9]TOWERS O L, SMITH A P. Stress intensity factors for curved crack fronts in compact tension specimens[J]. International Journal of Fracture, 1984, 25(2): 43-48

[10]NEWMAN J C, YAMADA Y, JAMES M A. Stress-intensity-factor equations for compact specimen subjected to concentrated forces[J]. Engineering Fracture Mechanics, 2010, 77(6): 1025-1029.

[11]SLAVIK D C, MCCLAIN R D, LEWIS K. Stress intensity predictions with ANSYS for use in aircraft engine component life prediction[J]. Fatigue and Fracture Mechanics, 2000, 31: 371-390.

[12]解德，錢勤，李長安. 斷裂力學中的數值計算方法及工程應用[M]. 北京: 科學出版社, 2009.

Calculation Formula for Stress Intensity Factors of CT Specimens

based on Three Dimensional Finite Element Solutions

JIA Xu, HU Xu-teng, SONG Ying-dong

(Jiangsu Province Key Laboratory of Aerospace Power Systems,

State Key Laboratory of Mechanics and Control of Mechanical Structures,

College of Energy and Power Engineering, Nanjing University of Aeronautics and Astronautics, Nanjing 210016, China)

Abstract:Two-dimensional and three-dimensional finite element calculations of stress intensity factors for standard compact tension(CT) specimens were compared and analyzed. And a new calculation formula of stress intensity factors for standard CT specimens was established based on the three-dimensional finite element virtual crack closure method and verified using displacement extrapolation method of ANSYS. The results show that great differences existed between the stress intensity factors from calculation formula based on the two-dimensional analysis and those of actual three-dimensional CT specimens. Under the certain equivalent distributed force on the load-hole, the stress intensity factors at the most crack front regions were similar to those at the center of the CT specimen and were independent of the thickness. The fitting calculation formula of stress intensity factors at the crack front center of three-dimensional CT specimens was of high precision. And the relative errors between the fitting calculation values and the results of the displacement extrapolation method of ANSYS were less than 0.5%.

Key words:three dimension analysis; CT specimen; stress intensity factor; finite element method

通訊作者：李落星教授

作者簡介：王冠(1985-)，男，河南鄭州人，講師，博士。

基金項目：國家自然科學基金面上資助項目(51475156)；寧夏大學自然科學研究基金資助項目(ZR1403)；寧夏大學人才引進科研啟動基金資助項目(BQD2014018)

收稿日期:2015-06-02；

修訂日期:2015-10-23

DOI:10.11973/jxgccl201512009

中圖分類號：O346.1

文獻標志碼：A

文章編號：1000-3738(2015)12-0030-05