彎扭組合載荷下圓軸表面裂紋應力強度因子計算

張 宇，黃小平，閆小順

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

0 引言

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

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

1 復雜載荷下的應力強度因子計算

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

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

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

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

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

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

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

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

2 有限元分析

2.1 模型基本參數

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

在施加彎矩或者扭矩時，本文采用cerig命令的方式，在構件中心位置建立一個節點，定義為mass21單元，將該節點和其他受力節點耦合，形成剛性區域，然后在該節點上直接施加力矩。圖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 裂紋體網格劃分

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

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

圖7 常規裂紋體網格劃分Fig.7 Regular crack meshing

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

圖9 網格劃分改進前后裂紋前緣各節點局部坐標系對比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 模型驗證

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

圖11 彎曲載荷下應力強度因子對比結果Fig.11 Comparison of SIFs under bending loading

2.4 表面點應力強度因子處理

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

3 結果分析與討論

3.1 扭轉載荷

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

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

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

對于裂紋最深點，在扭轉載荷下，只有KIII存在，隨著相對深度a/d的增加，KIII逐漸減小;隨著裂紋形狀比a/c的增加，KIII逐漸減小。其變化規律如圖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時，本文分別計算了彎曲、扭轉和彎扭組合載荷這3種載荷情況下的各型應力強度因子。3種載荷下的KI，KII，KIII沿裂紋前緣分布對比分別如表1～表3所示。可以看到，在計算表面裂紋的應力強度因子時，彎扭組合載荷可以視作彎曲載荷和扭轉載荷的簡單疊加，兩者之間相互影響可以忽略。如果已知彎曲載荷和扭轉載荷下的應力強度因子經驗公式，可以直接用于求解彎扭組合載荷下的各型應力強度因子，從而避免復雜耗時的有限元計算。

表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 結語

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

1)對裂紋體網格劃分進行了改進，使得裂紋前緣各點的Ⅱ，Ⅲ型應力強度因子計算結果更加準確。

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

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

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

［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］瞿偉廉，魯麗君，李明．帶三維穿透裂紋結構的有限元實體建模方法［J］．武漢理工大學學報，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］王效貴，郭乙木，許金泉．與界面相交的裂紋尖端的應力奇異性分析［J］．固體力學學報，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．