張 宇，黃小平，閆小順

(1．上海交通大學(xué)海洋工程國家重點實驗室，上海200240;2．中國艦船研究設(shè)計中心，湖北武漢430064)

0 引言

軸類構(gòu)件在船舶結(jié)構(gòu)中有許多應(yīng)用，如用作螺旋槳軸等。這些構(gòu)件需承受周期性的應(yīng)力，并因此導(dǎo)致疲勞等結(jié)構(gòu)損傷而提前失效。在傳統(tǒng)的設(shè)計思想中，人們通常把構(gòu)件材料視作無缺陷的連續(xù)均勻體，只根據(jù)S－N曲線來估計疲勞壽命，并不考慮裂紋萌生和裂紋擴展的時間，然而在實際結(jié)構(gòu)中，初始裂紋或者缺陷無法避免。已有研究表明，由表面缺陷引起的結(jié)構(gòu)斷裂是軸類構(gòu)件失效的主要原因之一。隨著斷裂力學(xué)的發(fā)展，現(xiàn)在已經(jīng)可以通過研究疲勞裂紋的擴展行為來更準(zhǔn)確地預(yù)測構(gòu)件的疲勞壽命。應(yīng)力強度因子作為判斷結(jié)構(gòu)斷裂和裂紋擴展的重要參數(shù)，也日益引起了人們的重視。

當(dāng)前，對于軸類或者管類構(gòu)件表面裂紋的研究，多集中在拉彎載荷下的Ⅰ型應(yīng)力強度因子上［1－6］。對于許多軸類構(gòu)件來說，除了承受外力和自身重量所引起的彎曲應(yīng)力之外，還要承受扭轉(zhuǎn)載荷，因此計算彎扭組合載荷下表面裂紋的擴展壽命更有意義。由于彎扭組合載荷下的裂紋不再是單一的張開型裂紋，而是復(fù)合型裂紋，相關(guān)的研究要少得多，且對表面點應(yīng)力強度因子計算存在較大誤差［7－9］。本文在前人研究成果的基礎(chǔ)上，以含半橢圓表面裂紋的實心圓軸為研究對象，對裂紋體網(wǎng)格劃分進行了改進，運用插值方法得到了裂紋表面點處的應(yīng)力強度因子，研究了扭轉(zhuǎn)載荷下裂紋形狀對應(yīng)力強度因子的影響，以及裂縫前緣各點Ⅱ、Ⅲ型應(yīng)力強度因子沿裂縫前緣的分布規(guī)律，并探討了彎、扭組合載荷下同型應(yīng)力強度因子疊加原理的有效性。從而為評估該類裂紋的擴展壽命提供參考。

1 復(fù)雜載荷下的應(yīng)力強度因子計算

三維裂紋形狀如圖1所示，x方向是裂紋正前方，y方向是裂紋面法線方向，z方向和裂紋面平行?？紤]一個離裂端很近，位置在極坐標(biāo) (r，θ)的單元，根據(jù)彈性力學(xué)，可以將裂紋尖端的位移場表示為:

式中:u，v，w 分別為直角坐標(biāo)系中 x，y，z方向的位移分量;r和θ為極坐標(biāo)系下的2個坐標(biāo)分量，且r的值遠小于裂紋長度。G為剪切模量，G=E/2(1+μ)，κ為與材料泊松比μ有關(guān)的常數(shù)，對于平面應(yīng)變問題，κ=3－4μ;對于平面應(yīng)力問題，κ=(3－μ)/(1－μ)。KI，KII，KIII分別是 I，II，III型應(yīng)力強度因子。由于裂紋前緣絕大部分范圍為平面應(yīng)變狀態(tài)，平面應(yīng)力狀態(tài)只在緊靠自由表面的一個很小的區(qū)域內(nèi)存在。因此，本文取κ=3－4μ。

圖1 三維裂紋模型Fig.1 Three dimensional crack model

圖2 定義求解路徑Fig.2 Defined path for solving the SIFs

根據(jù)式(1)，如果裂紋表面上(θ=180°)某一點的位移已知，那么可以得到:

式中:Δu，Δv，Δw分別為2個裂紋面之間沿x，y，z方向的相對位移。在用Ansys自帶命令kcalc計算應(yīng)力強度因子時，需要定義各個節(jié)點處的局部坐標(biāo)系和計算路徑，并且要求局部坐標(biāo)系的x軸垂直于裂紋前緣，y軸垂直于裂紋所在平面，求解路徑如圖2所示。在得到2～5點的相對位移之后，代入式(2)，即可得到這些節(jié)點處的應(yīng)力強度因子。然后再利用最小二乘法，即可外推得到裂紋尖端處的各型應(yīng)力強度因子的數(shù)值解。

為了使所研究的應(yīng)力強度因子變化規(guī)律具有普適性，這里對應(yīng)力強度因子進行無量綱化處理，無量綱應(yīng)力強度因子表示如下:

式中:σb為彎曲應(yīng)力;τxy為圓軸表面最大剪切應(yīng)力;a為裂紋深度。對于實心圓軸來說，彎曲應(yīng)力σb=，剪切應(yīng)力;M為施加的彎矩;T為施加的扭矩。

2 有限元分析

2.1 模型基本參數(shù)

本文中，計算所使用的實心圓軸直徑d=50 mm，長度l=300 mm，裂紋在圓軸中部，裂紋所在橫截面如圖3所示。用無量綱參數(shù)a/d和a/c描述裂紋形狀，其中，a/d為裂紋相對深度，a/c為裂紋形狀比。a，c分別為裂紋深度和裂紋長度，β為裂紋面與圓軸表面的交角，稱為裂紋嵌入角。用x/h來表示裂紋前緣上任意一點D的相對位置，h為裂紋寬度，x為裂紋前緣上D點的橫坐標(biāo)。a/c變化范圍為0.2～1.0，a/d的變化范圍取為0.05～0.5，x/h的變化范圍為0～1.0?？紤]到扭轉(zhuǎn)載荷的不對稱性，建立的模型為全模型，建模時采用商用有限元軟件Ansys14.0。為了適應(yīng)裂紋尖端節(jié)點處的應(yīng)力奇異性，裂紋前緣采用1/4 20節(jié)點等參退化奇異單元。對于圓軸的其他部分，則采用等參三維六面體單元solid95進行自由網(wǎng)格劃分，含裂紋的實心圓軸有限元模型如圖4所示。材料彈性模量為E=210 000 MPa，泊松比μ=0.3。

在施加彎矩或者扭矩時，本文采用cerig命令的方式，在構(gòu)件中心位置建立一個節(jié)點，定義為mass21單元，將該節(jié)點和其他受力節(jié)點耦合，形成剛性區(qū)域，然后在該節(jié)點上直接施加力矩。圖5是彎矩加載示意圖。圓軸另一端施加全約束，邊界條件如圖6所示。

圖3 裂紋所在橫截面Fig.3 Geometry of semi-elliptical surface crack

圖4 實心圓軸有限元模型Fig.4 Finite element model of shaft

圖5 彎矩加載Fig.5 Applied bending moment by using mass21 element

圖6 載荷及邊界條件Fig.6 Loading and boundary conditions

2.2 裂紋體網(wǎng)格劃分

在第1節(jié)中已經(jīng)提到，用Ansys自帶命令計算裂紋前緣各節(jié)點的應(yīng)力強度因子時，需要首先定義合適的局部坐標(biāo)系和求解路徑。在目前常用的裂紋體網(wǎng)格劃分中，通常是由面單元拖拉形成整個裂紋體單元［10］，如圖7所示。然而這種方法無法保證裂紋前緣各節(jié)點處局部坐標(biāo)系中的x軸垂直于裂紋前緣。由式(2)可知，在求解Ⅰ型應(yīng)力強度因子時，僅需要y軸方向的相對位移Δv。只要局部坐標(biāo)系的y軸方向垂直于裂紋所在平面，即可得到比較準(zhǔn)確的KI值。求解Ⅱ，Ⅲ型應(yīng)力強度因子時，則需要局部坐標(biāo)系下的x，z方向的相對位移Δu，Δw，如果局部坐標(biāo)系設(shè)定不合理，KII，KIII的求解結(jié)果會產(chǎn)生很大的誤差。

為此，本文對裂紋體網(wǎng)格劃分方法作了如下改進:將裂紋體切割成兩部分，再分別對這2個部分進行網(wǎng)格劃分，如圖8所示。從而使除了表面點附近區(qū)域外的裂紋前緣各節(jié)點處局部坐標(biāo)系中的x軸盡量垂直于裂紋前緣。網(wǎng)格劃分方法改進前后裂紋前緣各節(jié)點的局部坐標(biāo)系對比如圖9所示。當(dāng)a/c=0.4，a/d=0.1時，分別對裂紋體網(wǎng)格劃分改進前后的模型施加彎曲和扭轉(zhuǎn)載荷，計算得到KI，KII，KIII沿著裂紋前緣的分布規(guī)律如圖10所示?？梢钥吹?，網(wǎng)格劃分改進與否對于KI的影響很小，而對于裂紋前緣上的各點來說，愈靠近表面點處，網(wǎng)格劃分對局部坐標(biāo)系中x軸的方向影響愈大，從而對KII，KIII的計算結(jié)果的影響也愈大。需要指出的是，改進后的網(wǎng)格劃分方法仍然無法保證每一個局部坐標(biāo)系的x軸都完全垂直于裂紋前緣，但和改進前的網(wǎng)格劃分方法相比，相對而言要更加合理，下文中均采用改進后的網(wǎng)格劃分方法。

圖7 常規(guī)裂紋體網(wǎng)格劃分Fig.7 Regular crack meshing

圖8 改進網(wǎng)格劃分方法后的裂紋體模型Fig.8 Crack model after improving meshing method

圖9 網(wǎng)格劃分改進前后裂紋前緣各節(jié)點局部坐標(biāo)系對比Fig.9 Local coordinate system comparison between regular crack meshing and changed crack meshing

圖10 K I、K II、K III沿裂縫前緣分布對比(a/c=0.4，a/d=0.1)Fig.10 The comparison of K I，K II，K III along crack surface

2.3 模型驗證

為了驗證上述模型以及網(wǎng)格劃分的正確性，本文計算了在彎曲載荷下，裂紋形狀比a/c=1.0時最深點A處和近表面點B'處的Ⅰ型應(yīng)力強度因子，并 和 Murakami［3］， Shih ＆ Chen［4］， Couroneau ＆Royer［11］等的計算結(jié)果進行對比，如圖11所示。對比結(jié)果顯示，最深點的無量綱應(yīng)力強度因子KI與Murakami，Shih＆Chen以及Couroneau＆Royer等人的計算結(jié)果基本相符，誤差極小，只有裂紋深度很小時誤差才大于3.0%;對于近表面點B'處，本文的結(jié)果和已有文獻的平均值相比，誤差基本不超過5%，且變化趨勢和上述文獻相同，從而證明了本文建模以及應(yīng)力強度因子求解方法的正確性。

圖11 彎曲載荷下應(yīng)力強度因子對比結(jié)果Fig.11 Comparison of SIFs under bending loading

2.4 表面點應(yīng)力強度因子處理

有必要指出，線彈性斷裂力學(xué)的主要假設(shè)之一是裂紋前緣的連續(xù)性，對于裂紋前緣和幾何體自由表面的交點，即裂紋的表面端點來說，這一假設(shè)不再適用［7］。已經(jīng)有人從斷裂能量的角度出發(fā)，在理論上證明了在裂紋表面端點附近，應(yīng)力奇異性的次數(shù)取決于材料泊松比 μ 以及裂紋嵌入角 β［12－13］。Pook指出，在用三維有限元方法分析含有裂紋的幾何體時，所得到的裂紋表面端點附近的斷裂力學(xué)參數(shù)不可信。對于裂紋表面端點來說，運用退化1/4節(jié)點單元來模擬應(yīng)力和應(yīng)變的奇異性并不能得到準(zhǔn)確的結(jié)果。當(dāng)裂紋嵌入角β很小時，表面點處的應(yīng)力強度因子會突然趨向于0，而在裂紋嵌入角β很大時，表面點處的應(yīng)力強度因子會突然趨向于無窮大。但是這種突變只發(fā)生于表面端點附近很小的區(qū)域之中［12］。在已有的文獻里，人們在研究Ⅱ、Ⅲ型應(yīng)力強度因子沿裂紋前緣的分布時，通常避開這一區(qū)域，將x/h的最大值取在0.8～0.9之間［9］。然而對于構(gòu)件疲勞壽命預(yù)報來說，裂紋表面點的應(yīng)力強度因子非常重要。考慮到應(yīng)力強度因子沿裂縫前緣分布的連續(xù)性，本文采用插值方法外插得到表面點處的各型應(yīng)力強度因子，從而將x/h的取值范圍拓展到0～1.0。

3 結(jié)果分析與討論

3.1 扭轉(zhuǎn)載荷

對模型施加扭轉(zhuǎn)載荷后，本文計算了裂紋前緣各點的Ⅱ、Ⅲ型無量綱應(yīng)力強度因子KII，KIII。當(dāng)裂紋形狀比a/c=1.0，a/d=0.05～0.5時，KII，

圖12 各型應(yīng)力強度因子沿裂縫前緣分布Fig.12 K I，K II variation along the crack surface

KIII沿著裂縫前緣的分布如圖12所示。從圖中可以看到，除最深點外，裂紋前緣上其他各點均同時存在KII和KIII，說明扭轉(zhuǎn)載荷下裂紋并不是單一型裂紋，而是II－III型復(fù)合裂紋。KII在裂縫最深點為0，沿裂縫前緣隨著x/h的增大而逐漸增大，且隨著裂紋相對深度的增加，其增長速度逐漸變快;KIII在裂紋最深處達到最大值，沿著裂縫前緣隨x/h的增加而逐漸減小。

對于裂紋最深點，在扭轉(zhuǎn)載荷下，只有KIII存在，隨著相對深度a/d的增加，KIII逐漸減小;隨著裂紋形狀比a/c的增加，KIII逐漸減小。其變化規(guī)律如圖13和圖14所示。

圖13 最深點K III隨a/d變化Fig.13 K III variation with a/d at deepest point

圖14 最深點K III隨a/c變化Fig.14 K III variation with a/c at deepest point

3.2 彎扭組合載荷

在裂紋形狀比a/c=1.0，相對深度比a/d=0.1時，本文分別計算了彎曲、扭轉(zhuǎn)和彎扭組合載荷這3種載荷情況下的各型應(yīng)力強度因子。3種載荷下的KI，KII，KIII沿裂紋前緣分布對比分別如表1～表3所示?？梢钥吹?，在計算表面裂紋的應(yīng)力強度因子時，彎扭組合載荷可以視作彎曲載荷和扭轉(zhuǎn)載荷的簡單疊加，兩者之間相互影響可以忽略。如果已知彎曲載荷和扭轉(zhuǎn)載荷下的應(yīng)力強度因子經(jīng)驗公式，可以直接用于求解彎扭組合載荷下的各型應(yīng)力強度因子，從而避免復(fù)雜耗時的有限元計算。

表1 K I沿裂紋前緣分布對比(a/c=1.0，a/d=0.1)Tab.1 K I variation along the crack surface(a/c=1.0，a/d=0.1)

表2 K II沿裂紋前緣分布對比(a/c=1.0，a/d=0.1)Tab.2 K II variation along the crack surface(a/c=1.0，a/d=0.1)

表3 K III沿裂紋前緣分布對比(a/c=1.0，a/d=0.1)Tab.3 K III variation along the crack surface(a/c=1.0，a/d=0.1)

4 結(jié)語

本文采用三維有限元方法，對彎曲、扭轉(zhuǎn)和彎扭組合載荷下實心圓軸表面裂紋應(yīng)力強度因子進行了研究，主要研究內(nèi)容和結(jié)論如下:

1)對裂紋體網(wǎng)格劃分進行了改進，使得裂紋前緣各點的Ⅱ，Ⅲ型應(yīng)力強度因子計算結(jié)果更加準(zhǔn)確。

2)基于應(yīng)力強度因子沿著裂縫前緣分布的連續(xù)性，運用插值方法外插得到表面點處的應(yīng)力強度因子。

3)扭轉(zhuǎn)載荷之下的裂紋是Ⅱ－Ⅲ型復(fù)合裂紋，在裂紋最深點只有Ⅲ型應(yīng)力強度因子存在，且最深點Ⅲ型應(yīng)力強度因子隨著裂紋相對深度a/d的增大而減小，隨著裂紋形狀比a/c的增大而減小。裂縫前緣除最深點外均同時存在Ⅱ型和Ⅲ型應(yīng)力強度因子，隨著x/h的增大，Ⅱ型應(yīng)力強度因子占比逐漸增大。

4)彎扭組合載荷下的裂紋是Ⅰ－Ⅱ－Ⅲ型復(fù)合裂紋，裂紋前緣各點的Ⅰ、Ⅱ、Ⅲ型應(yīng)力強度因子可以視作分別施加彎曲載荷和扭轉(zhuǎn)載荷下裂紋前緣相同位置的同型應(yīng)力強度因子的疊加。本結(jié)論可為彎扭組合載荷下軸類構(gòu)件疲勞壽命預(yù)報提供參考。

［1］ RAJU I S，NEWMAN J C．Stress intensity factors for circumferential surface cracks in pipes and rods under tension and bending loads［J］．Fracture Mechanics，1986，17:789－805．

［2］ SHIRATORI M，MIYOSHIT，SAKAY Y，et al．Analysis and application of influence coefficients for round bar with a semi-elliptical surface crack［M］．In:Murakami Y，editor．Stress Intensity Factors Handbook， Vol II．Oxford:Pergamon Press，1987:659 － 65．

［3］ MURAKAMI Y，TSURU H．Stress intensity factor equations for semi-elliptical surface crack in a shaft under bending［M］．In:Murakami Y，editor．Stress Intensity Factors Handbook，Vol II．Oxford:Pergamon Press，1987:657 －8．

［4］ SHIH Y S，CHEN JJ．The stress intensity factor study of an elliptical cracked shaft［J］．Nuclear Engineering and Design，2002，214:137 －145．

［5］ CARPINTERI A，BRIGHENTI R．Fatigue propagation of surface flaws in round bars:a three-parameter theoretical model［J］．Fatigue ＆ Fracture of Engineering Materials ＆Structures，1996，19(12):1471 －1480．

［6］ TOBIRIO J．A critical review of stress intensity factor solutions for surface cracks in round bars subjected to tension loading［J］．Engineering Fatigue Analysis，2009，16(3):794－809．

［7］ MANUAL DA F，MANUEL DE F．Stress intensity factor for semi－elliptical surface cracks in round bars under bending and torsion［J］．Journal of Fatigue，1999，21(5):457－463．

［8］ SHAHANI R，HABIBI S E．Stress intensity factors in a hollow cylinder containing a circumferential semi-elliptical crack subjected to combined loading［J］．International Journal of Fatigue，2007，29:128 －140．

［9］ ISMAIL A E，ARIFFIN A K，ABDULLAH S，et al．Stress intensity factors under combined tension and torsion loadings［J］．Indian Journal of Engineering ＆ Materials Sciences，2012，19:5 －16．

［10］瞿偉廉，魯麗君，李明．帶三維穿透裂紋結(jié)構(gòu)的有限元實體建模方法［J］．武漢理工大學(xué)學(xué)報，2008，30(1):87－90．QUE Wei-lian，LU Li-jun，LI Ming．Solid modeling method for structure with 3 － D straight through crack［J］．Journal of Wuhan University of Technology，2008，30(1):87 －90．

［11］ COURONEAU N，ROYER J．Simplified model for the fatigue growth analysis of surface cracks in round bars under mode I［J］．Journal of Fatigue，1998，20(10):711－718．

［12］ POOK L P．On fatigue crack paths［J］．Journal of Fatigue，1995，17(1):5 －13．

［13］王效貴，郭乙木，許金泉．與界面相交的裂紋尖端的應(yīng)力奇異性分析［J］．固體力學(xué)學(xué)報，2002，23(4):412 －418．WANG Xiao-gui，GUO Yi-mu，XU Jin-quan．On the stress singularity of crack terminating at the interface［J］．Acta Mechanica Solida Sinica，2002，23(4):412 －418．