MCTS試件的三維有限元計算斷裂分析

李慶芬，齊桂營，朱莉，何水清

(哈爾濱工程大學機電工程學院，黑龍江哈爾濱150001)

工程構件和機器零部件中往往含有三維復雜裂紋和其他缺陷，這些裂紋和缺陷造成受力構件的應力集中、強度下降、并會導致災難性的斷裂事故.斷裂力學作為研究含裂紋結構強度的一門學科，可以正確評價帶裂紋構件的力學特性，其中應力強度因子作為表征外力作用下彈性體裂紋尖端附近應力場強度的一個重要參量，是評價斷裂破壞和裂紋擴展規律的重要指標.因此，求解裂紋體的應力強度因子是線彈性斷裂力學中很重要的一項工作，對評估含裂紋構件剩余強度和疲勞裂紋擴展壽命等都具有重要的意義.

在斷裂力學發展的初期，為了數學處理上的簡便，往往將實際裂紋簡化為平面二維模型.目前，二維斷裂問題的研究已十分成熟，已建立了一系列應力強度因子的理論計算公式［1］，并基于不同理論確定了不同的斷裂判據［2-6］.然而工程實際構件中的裂紋，由于構件幾何形狀及承載方式的多樣化，通常呈現的是三維復合型斷裂模式，理想的二維裂紋幾乎是不存在的.為了更好地解決工程實際中的斷裂問題，研究三維復合型裂紋的斷裂特性就變得尤為重要.而三維斷裂問題雖已受到廣大研究者的關注，但有關三維復合型斷裂問題的研究還很不成熟.由于結構和裂紋的復雜性，目前只有少數典型裂紋體的應力強度因子存在理論解析解［7］，大多數三維復合型斷裂問題在求解裂紋前沿應力強度因子時存在較大困難，尤其是在裂紋平面與三維構件自由表面的交匯處，迄今尚無精確的解析解.現有的三維斷裂判據，如 S 判據［3］、POOK 判據［8］和最大主應力判據［9］等，也只給出了平面純拉伸(Ⅰ型)、平面純剪切(Ⅱ型)和反平面剪切(Ⅲ型)加載條件下的裂紋起始角和斷裂極限值，而對于三維復合型裂紋，尤其是Ⅰ+Ⅱ+Ⅲ型復合裂紋的斷裂行為，則迄今還少有研究［10-11］.

本文將采用由Richard設計的特殊加載設備［12］和帶傾斜裂紋面(傾角γ為45°)的修正緊湊拉剪試件，通過改變裂紋的幾何形狀方位構建三維復合型裂紋，應用修正的虛擬裂紋閉合積分方法，對三維裂紋體復合型裂紋的應力強度因子進行分析計算，探討其裂紋尖端的斷裂行為.該研究將為深入探討三維復合型斷裂問題做出有益的貢獻，并為解決實際工程中存在類似裂紋構件的安全可靠性評估提供參考依據.

1 數學模型與計算方法簡介

1.1 裂紋尖端應力場及應力場強度因子

裂紋尖端應力場是建立斷裂判據的理論依據.直角坐標系下裂紋尖端應力場如下:

式中:σx、σy、σz、τxy、τxz、τyz為直角坐標系下各方向的應力分量，r、θ為裂紋尖端附近點的極坐標.其中，應力強度因子Ki(i=Ⅰ、Ⅱ、Ⅲ)唯一地描述了裂紋尖端應力場的強度，是評價斷裂行為的重要參數.因此，分析計算裂紋前沿的應力強度因子的分布情況是研究裂紋體斷裂特性的基礎.

1.2 修正的虛擬裂紋閉合積分方法

Rybicki和Kanninen在虛擬裂紋閉合積分方法(virtual crack closure integral method，VCCIM)［13］的基礎上提出了修正的虛擬裂紋閉合積分方法(modified virtualcrack closure integralmethod，MVCCI)［14］.該方法用裂紋長度為a的原始裂紋尖端節點位移ΔUy，i-1(a)代替虛擬裂紋閉合積分方法所要求擴展裂紋的節點位移ΔUy，j-1(a+Δa)，避免了當裂紋擴展為a+Δa時，二次有限元分析的附加影響(見圖1).具體表達式如下:

式(2)適用于承受平面純拉伸載荷即二維張開型(Ⅰ型)裂紋模式，圖1(a)為二維模型裂紋尖端的有限元網格.式中:t為試件厚度，F為節點力，ΔUy為節點的相對位移，Wy為節點力作用產生節點位移所做的功，G1c(a)為二維裂紋尖端的能量釋放率.

圖1 MVCCI計算方法的低階單元有限元網格Fig.1 FE low order standard elements for MVCCI-method

MVCCI這種高精度的解析方法還可以直接結合低階的實體單元來計算復雜的三維斷裂問題.圖1(b)為三維裂紋前沿的有限元網格.具體的三維MVCCI解析計算公式如下:

MVCCI方法是一種近似的方法，其值收斂于精確解.該方法已被用于計算典型的二維和三維斷裂問題［15］，并取得了較好的計算精度.

2 帶傾斜裂紋面的MCTS試件模型及有限元分析

2.1 MCTS試件及加載設備

本文選用帶傾斜裂紋面(傾角γ為45°)的修正緊湊拉剪試件(modified compact tension shear，MCTS).試件及其特定的加載設備如圖2所示.

裂紋平面的傾斜導致裂紋方位關于加載方向的不對稱分布，形成復合型裂紋，裂紋在不同加載條件下將會產生不同的復合型斷裂模式.隨加載角度α的改變，試件將從承受純拉伸(α=0°)加載，拉伸/剪切(0°＜α＜90°)復合型加載，變化到純剪切(α=90°)加載.

圖2 帶傾斜裂紋面的MCTS試件及其特定加載設備Fig.2 MCTS-specimen with inclined crack plane and the special loading device

試件的主要幾何參數如下:寬度W=90 mm，高度h=145 mm，厚度t=10 mm，裂紋長度a=50 mm，裂紋面傾角γ為45°.所選材料參數為:彈性模量E=80 kN/mm2，泊松比 ν =0.3，材料為鋁合金，試件所承受的靜態載荷F=40 kN.

2.2 有限元模型的建立

采用HyperMesh軟件提供的模型變形技術對已有的緊湊拉剪試件(compact tension shear specimen，CTS)［16］的計算模型進行改進，通過傾斜其裂紋平面構建具有復合型裂紋的MCTS計算分析模型，揭示裂紋面方位對斷裂行為的影響.為了形成傾斜裂紋平面并保持各個網格之間的關聯關系，建模過程中將對CTS試樣上下1/2有限元模型分別建立修改域，修改域中的操縱點沿總體坐標系y軸上下移動，使裂紋平面所在的整個橫截面發生傾斜，然后通過節點耦合模擬傾斜裂紋平面的直裂紋前沿，并合并非裂紋面區域的雙重節點，最終形成傾斜裂紋平面MCTS試樣的有限元模型.由于上述數值計算方法無需用奇異單元來模擬裂紋尖端的奇異性，所以整個三維有限元模型均采用八節點的solid45單元，并對裂紋前沿附近進行了網格細化，以保證計算精度，最精細(Δa=0.25 mm)的三維有限元模型共4 104個單元，4 886個節點.

圖3給出了帶傾斜裂紋面MCTS試件的整體有限元模型及裂紋前沿區域網格.

圖3 帶傾斜裂紋面MCTS試件的有限元模型及網格Fig.3 FE-model and mesh of the CTS-specimen with inclined crack plane

2.3 有限元計算方法

本文采用商用ANSYS軟件分別對Ⅰ型、Ⅱ型、及Ⅰ+Ⅲ復合型(α=45°)加載條件下的MCTS試樣進行了有限元計算.

首先根據MVCCI方法，計算得出三維裂紋體沿裂紋前沿各個節點處的能量釋放率(strain energy release rate，SERR)Gi(a，z/t)，i=Ⅰ、Ⅱ、Ⅲ.然后應用Irwin關系式［13］將能量釋放率轉換為應力強度因子如下:

將結果標準化如下:

式中:σ0=F/wt，σ0為裂紋面正應力.

最后，繪制出具有傾斜裂紋平面的MCTS試樣直裂紋前沿的應力強度因子分布曲線.進而探討傾斜裂紋平面MCTS試樣在不同加載條件下的斷裂行為.

3 計算結果與分析

3.1 Ⅰ型載荷作用下的計算結果與分析

Ⅰ型載荷作為脆性斷裂最危險的加載方式一直是研究重點.圖4給出了I型載荷(α=0°)作用下帶傾斜裂紋面(γ=45°)和直裂紋面(γ=0°)試樣裂紋前沿的應力強度因子的分布曲線.由圖可知:直裂紋面CTS的裂紋前沿上只有Ⅰ型應力強度因子的存在;而MCTS試樣的裂紋前沿上卻同時存在著Ⅰ型、Ⅱ型和Ⅲ型應力強度因子.說明在純張開型載荷作用下，直裂紋面試樣只能發生張開型斷裂，而傾斜裂紋面試樣不僅發生了張開型斷裂，同時還產生了滑開型和撕開型斷裂，即產生了Ⅰ+Ⅱ+Ⅲ復合型斷裂模式.說明I型載荷下，傾斜裂紋面試樣的斷裂模式不再滿足單一載荷與斷裂模式唯一對應的關系.

圖4 Ⅰ型載荷作用下MCTS試樣裂紋前沿的標準化SIFFig.4 Normalized SIF along crack front of MCTS-specimen under mode-Ⅰloading

MCTS試樣在I型加載下的變形情況如圖5(a)所示.由圖可看出，試件整體變形以張開型斷裂為主，但可以清楚地看到傾斜裂紋平面同時發生了滑開型和撕開型位移.CTS試樣的外形尺寸及加載條件與MCTS試樣完全相同，而MCTS試樣裂紋平面的傾斜，導致裂紋方位關于加載方向的不對稱分布，形成了復合型裂紋，即傾斜裂紋平面所形成的復雜裂紋幾何形狀是產生三維復合型斷裂模式的主要原因.

由計算結果還可看出，Ⅱ型應力強度因子在z/t=0.5處具有最大值，由此可知試件將在兩側自由表面附近首先發生失穩擴展.

圖5 MCTS試樣在不同加載下的變形情況Fig.5 Deformed FE-model of the MCTS-specimen under different loading conditions

3.2 Ⅱ型載荷作用下的計算結果與分析

圖6給出了Ⅱ型載荷(α =90°)作用下帶傾斜裂紋面(γ=45°)MCTS試樣裂紋前沿的應力強度因子的分布曲線.

由圖6可知:帶傾斜裂紋面MCTS試樣的裂紋前沿上雖然仍以Ⅱ型應力強度因子為主，但在試樣兩側自由表面附近區域可以看到少量Ⅰ型和Ⅲ型應力強度因子的存在.說明在純滑開型載荷作用下，傾斜裂紋面試樣不僅發生了滑開型斷裂，同時還產生了張開型和撕開型斷裂，即產生了Ⅰ+Ⅱ+Ⅲ復合型斷裂模式.這與直裂紋面試樣在純Ⅱ型載荷下往往只發生滑開型斷裂模式是不同的.由計算結果還可看出，應力強度因子在兩側自由表面附近具有最大值，由此可知試件將在兩側自由表面附近首先發生失穩擴展.

圖6 Ⅱ型載荷作用下MCTS試樣裂紋前沿的標準化SIFFig.6 Normalized SIF along crack front of MCTS-specimen under mode-Ⅱloading

MCTS試樣在 II型加載下的變形情況，如圖5(b)所示.由圖可看出，試件整體變形以滑開型斷裂為主，但也可以清楚的看到傾斜裂紋平面同時發生了張開型和撕開型位移.

3.3 復合型載荷作用下作用下的計算結果與分析

圖7給出了復合型載荷(α=45°)作用下帶傾斜裂紋面(γ=45°)MCTS試樣裂紋前沿的應力強度因子的分布曲線.由圖可知:在Ⅰ+Ⅱ復合型載荷作用下，裂紋前沿產生了Ⅰ+Ⅱ+Ⅲ復合型斷裂.裂紋前沿強度因子的分布曲線關于中心平面不對稱，在z/t=0.5處應力強度因子具有最大值，由此可知自由表面附近最易發生裂紋的失穩擴展.

圖7 復合型載荷作用下MCTS試樣裂紋前沿的標準化SIFFig.7 Normalized SIF along crack front of MCTS-specimen under mixed-mode loading

MCTS試樣在Ⅰ+Ⅱ復合型加載下的變形情況如圖5(c)所示.從圖7中也可看出，試件上下傾斜裂紋平面在Ⅰ+Ⅱ復合型載荷作用下同時產生了張開型、滑開型和撕開型位移.

上述數值分析結果表明:帶傾斜裂紋平面的MCTS試樣在不同加載條件下，其裂紋前沿均產生了Ⅰ+Ⅱ+Ⅲ復合型變形模式，即裂紋失穩擴展的臨界條件將受到3種主要斷裂模式共同作用的影響，其不再沿原裂紋平面擴展.

裂紋平面傾斜導致的裂紋方位不對稱及斷裂模式之間的耦合作用是產生復合型斷裂模式的主要原因.根據Richard判據，該復合型裂紋在3種基本斷裂模式的共同影響下，其裂紋前沿在初始擴展階段發生偏折的同時也產生了扭轉變形，斷裂破壞的不穩定區域則隨載荷的變化而變化，MCTS試樣在I型或Ⅱ型載荷作用下將在兩側自由表面附近首先發生失穩擴展，而對于平面內任意復合型加載條件，z/t=0.5處自由表面附近區域是該復合型裂紋斷裂破壞的危險區域.

縱觀不同加載條件下裂紋前沿應力強度因子的分布曲線，I型載荷作用時裂紋前沿各節點應力強度因子值均大于其他載荷情況，由此可知，傾斜裂紋平面MCTS試樣在平面純拉伸載荷作用下最易發生失穩擴展，說明I型載荷對于三維復合型裂紋來說仍是最具破壞的加載模式.

4 實驗結果

為了對數值計算結果進行驗證，采用透明有機玻璃脆性材料(便于觀察裂紋開裂情況)加工了帶有45°傾斜平面裂紋的MCTS試樣，并分別在Ⅰ型和Ⅱ型加載模式下進行了實驗(見圖8、9).

由實驗結果可知，在I型加載模式下(圖8)，裂紋在初始擴展階段同時發生扭轉和偏折變形，形成一個新的“螺旋”狀的裂紋面，隨著裂紋的不斷擴展，扭轉的裂紋面逐漸匯合成一個新的光滑平面，新裂紋面與加載方向垂直，裂紋沿這個新的裂紋表面方向擴展直至斷裂，即傾斜裂紋平面MCTS試樣最終發生了張開型斷裂.這與傾斜平面裂紋在初始變形階段，試件整體變形以張開型斷裂為主，但同時發生了滑開型和撕開型斷裂的計算結果是一致的.

圖8 在Ⅰ型加載下的MCTS試樣Fig.8 MCTS-specimen under mode Ⅰ loading

在Ⅱ型加載模式下(圖9)，裂紋面在初始擴展階段發生偏折的同時也產生了輕微的扭轉變形，扭轉變形在試樣兩側自由表面附近變形幅度較大，然而隨著擴展的深入，扭轉現象逐漸消失，裂紋最終發生了滑開型斷裂.這與裂紋在初始變形階段，試件整體變形以滑開型斷裂為主，但也同時發生了張開型和撕開型位移的計算結果是一致的.

圖9 在Ⅱ型加載下的MCTS試樣Fig.9 MCTS-specimen under mode Ⅱ loading

根據對比分析，實驗中試樣的裂紋開裂擴展特征與本文的計算結果是一致的，一定程度上驗證了本文采用MVCCI數值計算方法所做的計算分析.

5 結論

本文對帶有傾斜裂紋面的MCTS試件的應力強度因子進行了有限元分析及實驗驗證.研究結果表明:

1)試件在不同加載條件下均產生了Ⅰ+Ⅱ+Ⅲ復合型斷裂，裂紋平面傾斜導致的裂紋方位不對稱及斷裂模式之間的耦合作用是產生三維復合型斷裂模式的主要原因.

2)實驗中試樣的裂紋開裂擴展特征與有限元分析的結果一致，一定程度上驗證了本文的計算分析.

本文的研究拓展了三維斷裂研究的空間，為深入探討三維復合型斷裂問題做出了有益的貢獻，并為解決實際工程中存在類似裂紋構件的安全可靠性評估提供了參考依據.

［1］中國航空研究院.應力強度因子手冊［M］.北京:科學出版社，1993:7-438.Aeronautical Academy of China.Academy handbook of stress intensity factors［M］.Beijing:Science Press，1993:7-438.

［2］ERDOGAN F，SIH G C.On crack extension in plates under plane loading and transverse shear［J］.Trans ASME，Journal of Basic Engineering，1963，85:519-527.

［3］SIH G C.Strain energy density factor applied to mixed mode crack problems［J］.Int J Fracture，1974，10(3):305-321.

［4］陳增濤，王鐸.復合型裂紋的混合開裂準則［J］.哈爾濱工業大學學報，1994，26(4):127-130.CHEN Zengtao，WANG Duo.A mixed fracture criterion for mixed mode crack propagation［J］.Journal of Harbin Institute of Technology，1994，26(4):127-130.

［5］崔德渝，吳紹富，張行.復合型裂紋的能量釋放率計算與斷裂判據［J］.北京航空航天大學學報，1990，3:43-48.CUI Deyu，WU Shaofu，ZHANG Xing.Energy release rate and failure criterion for mixed mode cracking［J］.Journal of Beijing University of Aeronautics and Astronautics，1990，3:43-48.

［6］蔣玉川，胡興福.復合型裂紋擴展的最大塑性區尺度準則［J］.四川大學學報，2005，37(4):11-14.JIANG Yuchuan，HU Xingfu.The fracture criterion of the maximum size of plastic zone of mixed mode crack growth［J］.Journal of Sichuan University，2005，37(4):11-14.

［7］范天佑.斷裂理論基礎［M］.科學出版社，2003:138-186.FAN Tianyou.Foundation of fracture theory［M］.Beijing:Science Press，2003:138-186.

［8］POOK L P.A finite element analysis of the angle crack specimen［M］.London:Mechanical Engineering Publications，1993:285-302.

［9］SCHOLLMANN M，RICHARD H A，KULLMER G，FULLAND M.A new criterion for the prediction of crack development in multiaxially loaded structures［J］.Int J Fract，2002，117(2):129-141.

［10］QI Guiying，BUCHHOLZ F G，YAN Shengyuan，et al.Computational fracture analyses of a compact tension shear(CTS)specimen with an inclined crack front［J］.Key Engineering Materials，2010，417-418:629-632.

［11］CITARELLA R，CRICRì G.Comparison of DBEM and FEM crack path predictions in a notched shaft under torsion ［J］.Engineering Fracture Mechanics，2010，77(11):1730-1749.

［12］RICHARD H A，BENITZ K.A loading device for the creation of mixed mode in fracture mechanics［J］.Int J Fracture，22，1983:R55-R58.

［13］IRWIN G R.Analysis of stresses and strains near the end of a crack traversing a plate［J］.J of Applied Mechanics，1957，24:361-364.

［14］RYBICKI E F，KANNINEN M F.A finite element calculation of stress intensity factors by a modified crack closure integral［J］.Engng Fract Mech，1977，9:931-938.

［15］FRIEDRICH-G B，AZZEDINE C，HANS A R.Computational fracture analyses by means of virtual crack closure integral methods［C］//CD-ROM-Proceedings of the 6th Congreso Argentino de Mecanica Computational(MECOM99).Mendoza，Argentina，1999:29.1-29.19.

［16］RICHARD H A.A new compact shear specimen［J］.Int J Fracture，1981，17:R105-R107.