平板表面裂紋應(yīng)力強(qiáng)度因子和應(yīng)力分布規(guī)律的確定

Bui Manh Tuan 陳云飛

(1東南大學(xué)機(jī)械工程學(xué)院，南京210081)

(2濉和工業(yè)學(xué)院機(jī)械工程學(xué)院，越南濉和56000)

在進(jìn)行斷裂安全分析時(shí)，應(yīng)力強(qiáng)度因子(stress intensity factor，SIF)K是判斷裂紋體在載荷作用下是否產(chǎn)生破壞的主要參數(shù).目前，確定應(yīng)力強(qiáng)度因子的方法較多，典型的有解析法、有限元法(FEM)、虛擬裂紋擴(kuò)展技術(shù)［1］(virtual crack extension techniques)、修改裂紋閉合積分技術(shù)［2］(modified crack closure integral techniques)、J-集成技術(shù)［3-4］(J-integration techniques)、節(jié)點(diǎn)位移技術(shù)［5］等.由于有限元法計(jì)算結(jié)果相對(duì)精確，因此對(duì)于結(jié)構(gòu)或裂紋形狀復(fù)雜和受復(fù)雜載荷作用的結(jié)構(gòu)件，一般采用有限單元法.常用的有限元軟件包括ANSYS，ABAQUS，NASTRAN，MSC，Hyperworks等.與其他方法相比，有限元法具有單元的布局靈活、節(jié)點(diǎn)的配置方式比較任意，對(duì)裂紋形狀、位置無(wú)特殊限制等優(yōu)點(diǎn).因此，對(duì)幾何形狀和載荷比較復(fù)雜的裂紋體，有限元法是確定應(yīng)力強(qiáng)度因子最為有效的方法.本文基于有限元法和節(jié)點(diǎn)位移方法，研究了在中間裂紋與側(cè)裂紋處裂紋尖端的應(yīng)力強(qiáng)度因子和應(yīng)力場(chǎng).使用8節(jié)點(diǎn)四邊形等參元和1/4節(jié)點(diǎn)退化單元，對(duì)網(wǎng)格密度和裂紋長(zhǎng)度對(duì)計(jì)算精度的影響進(jìn)行了研究.

1 裂紋尖端附近的應(yīng)力和位移

在斷裂力學(xué)中，按裂紋體受力(即裂紋面相對(duì)位移的方向)可把裂紋分為Ⅰ，Ⅱ，Ⅲ型，如圖1所示.

圖1 三種基本裂紋模式

設(shè) u，v，w 是裂紋沿 x，y，z方向的位移分量，將裂紋面看成位移的間斷面，即空間坐標(biāo)相同裂紋的2個(gè)表面上點(diǎn)的位移是不連續(xù)的.

1)Ⅰ型裂紋(張開(kāi)形).受垂直于裂紋面的拉應(yīng)力作用，裂紋表面在x-y面內(nèi)法向張開(kāi).

2)Ⅱ型裂紋(滑開(kāi)形).受平行于裂紋面且垂直于裂紋前緣的剪應(yīng)力作用，裂紋面在x-z面內(nèi)沿x軸方向相對(duì)滑動(dòng).

3)Ⅲ型裂紋(撕開(kāi)形).受平行于裂紋面和裂紋前緣的剪應(yīng)力作用，裂紋面在x-z面內(nèi)沿z軸方向相對(duì)滑動(dòng).

對(duì)任意裂紋前沿為曲線(xiàn)的表面裂紋，在裂紋尖端取坐標(biāo)原點(diǎn)，局部直角坐標(biāo)系的x-y平面為裂紋前緣任一點(diǎn)的法平面，y-z平面為其切平面;局部柱坐標(biāo)系的r→0平面為裂紋前緣任一點(diǎn)的法平面，如圖2所示.根據(jù)Westergaard復(fù)變函數(shù)方法可求得線(xiàn)彈性斷裂力學(xué)中Ⅰ，Ⅱ，Ⅲ型裂紋尖端(r→0)的應(yīng)力應(yīng)變場(chǎng)為［6-7］:

1)Ⅰ型

2)Ⅱ型

3)Ⅲ型

圖2 1/4節(jié)點(diǎn)單元

式中，r，θ為局部柱坐標(biāo)系中的2個(gè)坐標(biāo)分量;G為剪切彈性模量，G=E/2(1+μ)，E為彈性模量，μ為材料泊松比;σx，σy為法向應(yīng)力;τxy，τxz，τyz為剪切應(yīng)力;KⅠ，KⅡ，KⅢ為Ⅰ，Ⅱ，Ⅲ型應(yīng)力強(qiáng)度因子;k 為與材料泊松比μ有關(guān)的常數(shù)，k=(3－μ)/(1+μ).

實(shí)際結(jié)構(gòu)中的裂紋并不是受單一荷載作用的.因此，由任意外載荷作用所產(chǎn)生的裂紋尖端附近區(qū)域的應(yīng)力場(chǎng)和位移場(chǎng)，可分別由上述的3種獨(dú)立變形型式的應(yīng)力場(chǎng)和位移場(chǎng)線(xiàn)形疊加給出，即

在常見(jiàn)的裂紋尖端場(chǎng)表達(dá)式中，一般高階量O(r)，O(r0)都忽略，因?yàn)楫?dāng)r與面內(nèi)尺寸相比很小時(shí)，這些高階量與主導(dǎo)項(xiàng)相比，可以略去不計(jì).因此，通常采用主導(dǎo)項(xiàng)來(lái)表達(dá)線(xiàn)彈性裂紋尖端區(qū)的應(yīng)力位移場(chǎng).

由裂紋尖端區(qū)域的應(yīng)力場(chǎng)式(17)和位移場(chǎng)式(18)可知:如果給定KⅠ，KⅡ，KⅢ三個(gè)參數(shù)，則裂紋尖端區(qū)域的應(yīng)力場(chǎng)和位移場(chǎng)也就可以完全確定.KⅠ，KⅡ，KⅢ與坐標(biāo)(r，θ)無(wú)關(guān)，它們表征了靠近裂紋尖端區(qū)的應(yīng)力場(chǎng)與位移場(chǎng)奇異性的強(qiáng)度，其值由裂紋體的幾何和所受載荷決定，因此可以用來(lái)度量相應(yīng)的應(yīng)力應(yīng)變場(chǎng)強(qiáng)度，從而判斷裂紋是否會(huì)擴(kuò)展以及擴(kuò)展的方向［6，8］.

2 應(yīng)力強(qiáng)度因子計(jì)算方法

2.1 1/4節(jié)點(diǎn)退化單元

奇異單元(quarter-point element)是繞裂紋尖端附近構(gòu)造的一種特殊單元，裂紋尖端區(qū)域除采用奇異單元之外，其余采用常規(guī)單元.奇異單元的特點(diǎn)是位移插值函數(shù)中直接包含所需要的奇異項(xiàng)，因而能夠反映裂紋尖端附近的奇異性.目前，已有許多不同的奇異單元用于求解裂紋尖端應(yīng)力應(yīng)變場(chǎng)，其中應(yīng)用較為廣泛的是8節(jié)點(diǎn)等參數(shù)奇異單元，其節(jié)點(diǎn)退化如圖2所示［6-8］.

2.2 有限元位移法求解應(yīng)力強(qiáng)度因子

采用有限元法和節(jié)點(diǎn)位移技術(shù)來(lái)計(jì)算應(yīng)力強(qiáng)度因子K的方法，也被稱(chēng)為直接方法.

當(dāng) θ= ± π，r=ra－b時(shí)，應(yīng)力強(qiáng)度因子可以用下式計(jì)算:

對(duì)于平面應(yīng)變狀態(tài)，Ⅱ型裂紋的應(yīng)力強(qiáng)度因子KⅡ可用裂紋面的位移來(lái)表示，即

式中，ua和ub為y方向在點(diǎn)a和點(diǎn)b上的位移;rab為裂紋尖端到點(diǎn) b的距離(即點(diǎn) a到點(diǎn) b的距離).

有限元法計(jì)算應(yīng)力強(qiáng)度因子的準(zhǔn)確度依賴(lài)于網(wǎng)格的細(xì)度和單元類(lèi)型，如果采用1/4節(jié)點(diǎn)單元，該方法準(zhǔn)確度將得到改善.位移沿裂紋邊對(duì)1/4節(jié)點(diǎn)退化單元的值為

式中，vupper，vlower分別為1/4節(jié)點(diǎn)沿abc處和ade處的裂紋位移，如圖2所示.

根據(jù)有限元法計(jì)算結(jié)果，裂紋張開(kāi)位移(crack open displacements，COD)為

應(yīng)力強(qiáng)度因子 K 為［5，7］

KⅠ，KⅡ，KⅢ的計(jì)算公式相同，只要代入不同方向的位移量即可計(jì)算Ⅰ，Ⅱ，Ⅲ型的應(yīng)力強(qiáng)度因子.

有限元法計(jì)算的流程如圖3所示.

圖3 有限元分析流程圖

3 算例

3.1 計(jì)算模型1(中間裂紋)［9-10］

計(jì)算模型1如圖4所示，彈性模量 E=20 MN/cm2;泊松比μ=0.3;法向應(yīng)力P=30 N/cm.由于模型對(duì)稱(chēng)，只建一半模型.模型網(wǎng)格劃分采用8節(jié)點(diǎn)四邊形等參元和1/4節(jié)點(diǎn)退化單元.在裂紋尖端附近網(wǎng)格的細(xì)度劃分如圖5所示.

圖4 計(jì)算模型1

圖5 網(wǎng)格劃分與在裂紋尖端附近網(wǎng)格的細(xì)度

對(duì)于中間裂紋，應(yīng)力強(qiáng)度因子理論計(jì)算公式［9］為

式中，aL為裂紋長(zhǎng)度;bL為板跨度;Ktheory為裂紋理論計(jì)算強(qiáng)度應(yīng)子;KFEM為有限元計(jì)算強(qiáng)度應(yīng)子;Kerror為有限元計(jì)算與理論計(jì)算強(qiáng)度應(yīng)子的誤差.

應(yīng)力強(qiáng)度因子計(jì)算精度與網(wǎng)格劃分及1/4節(jié)點(diǎn)單元退化的尺寸有關(guān).采用有限元法和理論計(jì)算法計(jì)算的Ⅰ型應(yīng)力強(qiáng)度因子如圖6所示.

假設(shè)裂紋長(zhǎng)度 aL=5，10，…，50 cm，則應(yīng)力強(qiáng)度因子KⅠ與裂紋長(zhǎng)度aL之間的關(guān)系如圖7所示.

假設(shè)荷載 q=10，20，…，50 N/cm，應(yīng)力強(qiáng)度因子KⅠ與不同荷載之間的關(guān)系如圖8所示.

圖6 不同網(wǎng)格數(shù)2種計(jì)算方法得到的KⅠ比較

圖7 應(yīng)力強(qiáng)度因子KⅠ與裂紋長(zhǎng)度之間的關(guān)系

圖8 不同載荷下的應(yīng)力強(qiáng)度因子KⅠ

由圖6～圖8可以看出，精確度(收斂程度)與網(wǎng)格密度、1/4節(jié)點(diǎn)單元退化的尺寸和裂紋長(zhǎng)度有關(guān).在網(wǎng)格數(shù)為54(網(wǎng)格密度低)時(shí)，有限元法計(jì)算結(jié)果與理論值的誤差為6.861 9%，提高劃分網(wǎng)格密度值，當(dāng)網(wǎng)格數(shù)大于868時(shí)，誤差小于1.777 9%.當(dāng)裂紋長(zhǎng)度為15 cm時(shí)，有限元位移法計(jì)算結(jié)果與理論結(jié)果相比，誤差為最大(1.508 5%);當(dāng)裂紋長(zhǎng)度為35 cm時(shí)，計(jì)算結(jié)果與理論結(jié)果相比，誤差最小(0.004 4%).不同載荷下，有限元位移法計(jì)算結(jié)果與理論結(jié)果計(jì)算一致.

3.2 計(jì)算模型 2(側(cè)裂紋)［9，11-12］

計(jì)算模型2如圖9所示，彈性模量 E=20 MN/cm2;泊松比 μ=0.3，切向應(yīng)力Q=30 N/cm.

圖9 計(jì)算模型2

對(duì)于側(cè)裂紋理論計(jì)算公式為［7］

應(yīng)力強(qiáng)度因子計(jì)算精度與網(wǎng)格劃分及1/4節(jié)點(diǎn)退化單元的尺寸有關(guān).采用有限元位移法和理論計(jì)算法計(jì)算的Ⅱ型應(yīng)力強(qiáng)度因子如圖10所示.

假設(shè)裂紋長(zhǎng)度aL為10，20，…，80 cm，則應(yīng)力強(qiáng)度因子KⅡ與裂紋長(zhǎng)度aL之間的關(guān)系如圖11所示.

假設(shè)荷載 Q=10，20，…，50 N/cm，則應(yīng)力強(qiáng)度因子KⅡ與不同載荷分布力之間的關(guān)系如圖12所示.

由圖8～圖12可以看出，KⅡ有限元法計(jì)算結(jié)果與理論值的誤差大部分小于0.773 2%.在網(wǎng)格數(shù)為454時(shí)，誤差小于0.266 0%.不同裂紋長(zhǎng)度應(yīng)力強(qiáng)度因子KⅡ的有限元計(jì)算結(jié)果與理論結(jié)果相比，誤差很小.當(dāng)裂紋長(zhǎng)度為80 cm時(shí)，有限元位移法計(jì)算結(jié)果與理論結(jié)果相比誤差為最大(1.049 1%).當(dāng)裂紋長(zhǎng)度為25 cm時(shí)，有限元位移法計(jì)算結(jié)果與理論結(jié)果相比誤差為最小(0.061 8%).不同載荷下，有限元位移法計(jì)算結(jié)果與理論結(jié)果計(jì)算一致.

圖10 不同網(wǎng)格數(shù)2種計(jì)算方法應(yīng)力強(qiáng)度因子KⅡ比較

圖11 應(yīng)力強(qiáng)度因子KⅡ與裂紋長(zhǎng)度之間的關(guān)系

圖12 不同載荷下的應(yīng)力強(qiáng)度因子KⅡ

4 裂紋尖端附近的位移場(chǎng)與應(yīng)力場(chǎng)計(jì)算

除了在裂紋尖端附近網(wǎng)格細(xì)度加密外，網(wǎng)格的劃分采用8個(gè)節(jié)點(diǎn)單元類(lèi)型和1/4節(jié)點(diǎn)單元退化.

位移場(chǎng)與應(yīng)力場(chǎng)的中間裂紋有限元模擬結(jié)果如圖13所示.位移場(chǎng)與應(yīng)力場(chǎng)的側(cè)裂紋有限元模擬結(jié)果如圖14所示.

圖13 位移場(chǎng)和應(yīng)力場(chǎng)的中間裂紋

由位移場(chǎng)和應(yīng)力場(chǎng)計(jì)算結(jié)果可見(jiàn)，在裂紋尖端附近存在嚴(yán)重退化，且應(yīng)力最大.

5 結(jié)語(yǔ)

本文基于有限元法和節(jié)點(diǎn)位移方法采用8節(jié)點(diǎn)四邊形等參元和1/4節(jié)點(diǎn)退化單元，對(duì)中間裂紋和側(cè)裂紋尖端的強(qiáng)度應(yīng)子進(jìn)行了計(jì)算.采用整體較稀疏的網(wǎng)格劃分，但在裂紋尖端附近區(qū)域增加網(wǎng)格和節(jié)點(diǎn)的密度，能夠大大減少網(wǎng)格數(shù)量，便于方程組的求解，節(jié)省計(jì)算時(shí)間，同時(shí)確保了問(wèn)題的準(zhǔn)確性.

通過(guò)有限元模型計(jì)算結(jié)果與理論計(jì)算結(jié)果相比表明:有限元模型計(jì)算精度與網(wǎng)格密度、1/4節(jié)點(diǎn)單元退化的尺寸和裂紋尖端的裂紋長(zhǎng)度有關(guān)，其中網(wǎng)格劃分是最主要的影響因數(shù)，提高劃分密度可以有效提高計(jì)算精度.位移場(chǎng)和應(yīng)力場(chǎng)計(jì)算結(jié)果表明，中間裂紋和側(cè)裂紋的裂紋尖端附近存在嚴(yán)重退化，且應(yīng)力最大.

圖14 位移場(chǎng)和應(yīng)力場(chǎng)的側(cè)裂紋

References)

［1］Liu P F，Hou S J，Chu J K，et al.Finite element analysis of postbuckling and delamination of composite laminates using virtual crack closure technique［J］.Composite Structures，2011，93(6):1549-1560.

［2］Muthu N，F(xiàn)alzon B G，Maiti S K，et al.Modified crack closure integral technique for extraction of SIFs in meshfree methods［J］.Finite Elements in Analysis and Design，2014，78(5):25-39.

［3］Okada Hiroshi，Ohata Shogo.Three-dimensional J-integral evaluation for cracks with arbitrary curvatures and kinks based on domain integral method for quadratic tetrahedral finite element［J］.Engineering Fracture Mechanics，2013，109(9):58-77.

［4］劉明堯，柯夢(mèng)龍，周祖德，等.裂紋尖端應(yīng)力強(qiáng)度因子的有限元計(jì)算方法分析［J］.武漢理工大學(xué)學(xué)報(bào)，2011，33(6):116-121.Liu Mingyao，Ke Menglong，Zhou Zude，et al.Analysis of finite element calculation methods for crack-tip stress intensity factor［J］.Journal of Wuhan University of Technology，2011，33(6):116-121.(in Chinese)

［5］Fehl Barry D，Truman Kevin Z.An evaluation of fracture mechanics quarter-point displacement techniques used for computing stress intensity factors［J］.Engineering Structures，1999，21(5):406-415.

［6］張俊清.高速列車(chē)空心車(chē)軸表面裂紋應(yīng)力強(qiáng)度因子研究［D］.北京:北京交通大學(xué)交通運(yùn)輸學(xué)院，2011.

［7］Tracy Dennis M.Finite elements for determination of crack tip elastic stress intensity factors［J］.Engineering Fracture Mechanics，1971，3(3):255-266.

［8］Wu Zhixue.The shape of a surface crack in a plate based on a given stress intensity factor distribution［J］.International Journal of Pressure Vessels and Piping，2006，83(3):168-180.

［9］Tada H，Paris P，Irwin G.The stress analysis of cracks handbook［M］.Washington DC，USA:Washington U-niversity，1957:41-52.

［10］Brennan F P，Teh L S.Determination of crack tip stress intensity factors in complex geometries by composition of weight function solutions［J］.Fatigue ＆Fracture of Engineering Materials＆Structures，2004，27(1):1-7.

［11］Muthu N，Maiti S K，F(xiàn)alzon B G，et al.A comparison of stress intensity factors obtained through crack closure integral and other approaches using extended element-free Galerkin method［J］.Computational Mechanics，2013，52(3):587-605.

［12］Salari R H，Rahimi D M，Rahimi P S，et al.Meshless EFG simulation of linear elastic fracture propagation under various loadings［J］.Arabian Journal for Science and Engineering，2011，36(7):1381-1392.