余海峰,吳興文,梁樹林,池茂儒

(1. 西南交通大學(xué) 牽引動力國家重點實驗室,成都 610031; 2. 西南交通大學(xué) 機械工程學(xué)院,成都 610031)

車軸作為動車組走行部中及其重要的部件,承受著各種復(fù)雜載荷,當(dāng)車軸表面出現(xiàn)裂紋后,有必要采用基于斷裂力學(xué)的損傷容限方法評估其剩余壽命。在進行斷裂安全分析時,應(yīng)力強度因子是一個重要的指標,它是判斷裂紋體在載荷作用下裂紋擴展規(guī)律和斷裂破壞的主要參數(shù)。因此,對于車軸剩余壽命的安全評估,應(yīng)力強度因子的求解極為重要。

目前針對三維應(yīng)力強度因子的求解,主要方法有實驗法、權(quán)函數(shù)法、解析法和有限元法等。趙大華等[1]和李相麟等[2]采用三維光彈性應(yīng)力凍結(jié)法對帶環(huán)形裂紋的圓軸進行了實驗分析,利用雙參數(shù)法測定I型應(yīng)力強度因子,切片逐次削去法測定III型應(yīng)力強度因子,取得了較好的結(jié)果。Wang等[3]和Jin等[4]針對有限厚度板中的半橢圓形表面裂紋,采用權(quán)函數(shù)法計算任意一維以及二維應(yīng)力場作用下的應(yīng)力強度因子,并將計算結(jié)果與有限元法對比,驗證了方法的有效性。周素霞等[5-6]通過對彎曲載荷下實心軸表面裂紋應(yīng)力強度因子解析式進行修正,得出空心軸表面裂紋應(yīng)力強度因子,并將解析式計算結(jié)果與有限元法進行比較,發(fā)現(xiàn)結(jié)果非常相近。

以上幾種方法中,實驗法受到實驗設(shè)備和條件的制約,權(quán)函數(shù)法對于車軸這種復(fù)雜結(jié)構(gòu)難于實施,解析法只能針對簡單垂向彎曲載荷下的計算,而有限元法以其強大的模擬和數(shù)值計算功能已被廣泛用于應(yīng)力強度因子的計算。一種常見的方法是通過有限元法來計算裂紋體裂紋尖端附近的位移場和應(yīng)力場,然后通過位移外推法或應(yīng)力外推法來計算應(yīng)力強度因子[7-9]。另外一種方法,更常見于有限元軟件中,是通過J積分或I積分等,利用能量與應(yīng)力強度因子之間的關(guān)系來計算應(yīng)力強度因子[10-12]。其中,第一種方法是通過將節(jié)點斷開來模擬裂紋,需要在裂紋尖端設(shè)置足夠密集的網(wǎng)格才能夠精確地模擬裂紋尖端附近的位移場或應(yīng)力場;第二種方法是通過特殊的建模手段,在裂紋尖端設(shè)置奇異單元來模擬裂紋,雖然其對網(wǎng)格的精度沒有太高的要求(由于J積分與路徑無關(guān)),但是在模擬多個不同裂紋形狀的裂紋時,需要大量的網(wǎng)格劃分工作。

綜上,雖然有限元法能夠很好地求解車軸的應(yīng)力強度因子,但由于裂紋形狀的演變,需要建立多組含裂紋車軸有限元模型來進行不同載荷下的單個應(yīng)力強度因子求解。其中,Madia等[13-14]以表格的形式給出了指定裂紋深度下的形狀因子以求解車軸應(yīng)力強度因子;Pokorny等[15]針對車軸的不同裂紋深度,基于應(yīng)力強度因子解析式的一般形式,采用多項式對形狀因子函數(shù)進行擬合?？傮w來說,目前對于車軸整個裂紋演變周期內(nèi)的任意裂紋深度的應(yīng)力強度因子的研究較少。為此,本文基于車軸應(yīng)力強度因子的一般解析式,根據(jù)有限元法中位移外推法的計算結(jié)果,使用五次多項式對形狀因子函數(shù)進行擬合,在給定裂紋形狀演變特性和載荷的情況下,推導(dǎo)出同一載荷模式下任意載荷大小和裂紋深度的應(yīng)力強度因子的解析式,并驗證了形狀因子函數(shù)和解析式的適用性。

1 有限元法求解應(yīng)力強度因子的原理

對于車軸來說,其主要受載為垂向載荷、壓裝配合載荷和橫向載荷,在運動時受到旋轉(zhuǎn)彎曲的作用,加速、制動產(chǎn)生的扭轉(zhuǎn)載荷占少數(shù)[16]。并且有實驗證明,對動車組動車空心車軸進行動應(yīng)力測試,在兩種編組運行工況下,車軸的彎曲應(yīng)力遠遠大于車軸的扭轉(zhuǎn)應(yīng)力[17-18]。在這種受載情況下,裂紋面受到拉應(yīng)力作用會張開,因此這里主要研究車軸的I型裂紋。同時許多研究表明,沿裂紋前緣絕大部分范圍為平面應(yīng)變狀態(tài),平面應(yīng)力狀態(tài)只在緊靠自由表面的一個很小的區(qū)域內(nèi)存在[17-18],因此,假設(shè)車軸裂紋前沿平面內(nèi)的應(yīng)力分布為平面應(yīng)變狀態(tài),便可以將三維裂紋簡化為二維模型來求解。

如圖1所示,考慮嵌入線彈性材料中的裂紋,對于沿三維裂紋前沿的給定位置,引入垂直于裂紋前沿的二維平面。

圖1 二維平面內(nèi)裂紋尖端坐標系的定義Fig. 1 Definition of coordinate systems at crack tip in a two-dimensional plane

在二維平面內(nèi),從裂紋尖端點引入兩個坐標系,一個為xyz直角坐標系,一個為rθ極坐標系。選取一個距離裂紋尖端很近,極坐標為(r,θ)的節(jié)點,根據(jù)線彈性斷裂力學(xué)理論,可以推得裂紋尖端的應(yīng)力表達式為

(1)

進一步,可以得到裂紋尖端的位移表達式為

(2)

式中:σx、σy和τxy分別為正應(yīng)力和切應(yīng)力分量;ux、uy分別為x和y方向上的位移分量;KI為I型應(yīng)力強度因子;E為彈性模量;υ為泊松比。

一般來說,σy和uy會使裂紋張開,因此是裂紋前沿最重要的變量。在裂紋面上,取θ=0,由式(1)可得

(3)

取θ=π,由式(2)可得

(4)

因此,裂紋尖端點的應(yīng)力強度因子KI,tip為

(5)

或

(6)

由于數(shù)值計算不可能得到r→0的應(yīng)力和位移,一般需要在裂紋尖端前取幾個距尖端很近的點,利用式(3)、式(4)計算這些點的應(yīng)力強度因子,然后通過外推獲得尖端點的應(yīng)力強度因子。

值得注意的是,由式(2)可知當(dāng)r→0時,uy=0,即假設(shè)裂紋尖端點的位移值為0。而在實際加載中,裂紋尖端點的位移不為0,即裂紋場產(chǎn)生了剛體位移,那么實際上裂紋前沿的位移場uy,real應(yīng)該為有限元計算出來的位移值uy,FEM減去裂紋尖端點的位移值uy,tip,即

uy,real=uy,FEM-uy,tip

(7)

所以,采用位移外推法時,裂紋尖端點的應(yīng)力強度因子KI,tip為

(8)

2 車軸應(yīng)力強度因子計算模型

在給定車軸材料屬性(彈性模量和泊松比)的情況下,車軸的應(yīng)力強度因子與以下3個變量有關(guān):初始裂紋的位置、初始裂紋形狀及其演變特性和車軸所受載荷。本文首先通過無裂紋車軸有限元計算結(jié)果確定了初始裂紋位置,接著給定了初始裂紋形狀及其演變特性,然后確定采用疊加法計算垂向載荷和壓裝配合載荷下的應(yīng)力強度因子,最后建立含裂紋車軸有限元模型,進行應(yīng)力強度因子的計算。

2.1 初始裂紋位置的確定

計算應(yīng)力強度因子的第一步是要確定初始裂紋的位置。通常在確定車軸剩余壽命時,所考慮的初始裂紋位置位于可能導(dǎo)致車軸最早疲勞失效的地方,由于此處研究的是I型裂紋,軸向主應(yīng)力對其裂紋擴展起決定性作用,因此本文以車軸最大軸向主應(yīng)力為判據(jù)來確定初始裂紋的位置。為此首先建立了無裂紋車軸有限元模型,來計算車軸在垂向載荷和壓裝配合載荷下的應(yīng)力分布。

輪對模型如圖2所示,包括空心車軸、車輪和齒輪箱。

圖2 無裂紋車軸有限元模型Fig. 2 Finite element model of a crack-free axle

由于幾何和載荷的對稱性,使用ABAQUS建模時僅建立輪對一半的結(jié)構(gòu),網(wǎng)格類型使用C3D8三維線性六面體單元。設(shè)置車軸為線彈性體,其材料的彈性模量為206 GPa,泊松比為0.3。加載設(shè)置如下:

1) 第一個載荷步,壓裝配合工況,一側(cè)軸端全約束,選取齒輪箱和車輪內(nèi)表面為主面,空心車軸外表面為從面,設(shè)置摩擦因數(shù)為0.6,過盈量為0.28 mm,以逐漸遞增的方式分多個載荷子步施加,保證接觸的良好收斂。

2) 第二個載荷步,垂向載荷工況,車輪與鋼軌接觸部分全約束,在車軸兩側(cè)軸肩軸承安裝處施加垂向載荷(F=34 kN模擬車輛靜載荷,該載荷為車軸受到的主要載荷68 kN的一半)。

在施加過盈接觸載荷和垂向載荷后,計算的有限元結(jié)果如圖3所示,車軸的軸向主應(yīng)力在輪座和齒輪箱座之間的過渡槽(S型過渡槽)處達到最大。分析原因: 1) 由于S型過渡槽處存在缺口效應(yīng),導(dǎo)致此處容易產(chǎn)生應(yīng)力集中; 2) 由于此處垂向載荷產(chǎn)生的彎曲應(yīng)力較大; 3) 由于車輪、齒輪箱與車軸的壓裝配合,改變了此處的應(yīng)力比,使得車軸輪座和齒輪箱座表面處受到殘余壓應(yīng)力,而S型過渡槽表面處受到殘余拉應(yīng)力?？紤]軸向拉應(yīng)力對I型裂紋的擴展起主要作用,因此將初始裂紋位置確定為S型過渡槽處。

圖3 有限元計算結(jié)果Fig. 3 Calculation results of finite element model

2.2 初始裂紋形狀及其演變特性

影響車軸應(yīng)力強度因子的另一個重要因素是裂紋前沿的形狀,即初始裂紋形狀及其演變特性。文獻中通常使用半橢圓形表面裂紋來描述車軸的裂紋形貌[15-16,19],橢圓軸之間的比率在裂紋擴展過程中不斷變化,符合車軸真實裂紋前沿的演變過程。因此,本文采用半橢圓形裂紋,用于車軸的裂紋擴展仿真。

根據(jù)經(jīng)驗和文獻[20],疲勞裂紋垂直于最大主應(yīng)力擴展,而最大主應(yīng)力的方向幾乎與軸向主應(yīng)力相同,因此,可以假設(shè)疲勞裂紋垂直于車軸擴展。由于文獻[21]模擬的車軸裂紋初始位置以及考慮的載荷種類均與本文一致,因此本文引用文獻中的裂紋形狀演變曲線來進行應(yīng)力強度因子的計算。如圖4所示,考慮無損探傷的精度,設(shè)置初始裂紋深度為1 mm,裂紋深度a從1、1.5、2、2.5、3、4、5、6.5、9、12、15、20、25、32、38、45、55 mm逐漸遞增,裂紋形狀比a/c的比率隨著裂紋深度a的增長從0.814逐漸減小到0.281。

圖4 裂紋形狀演變曲線Fig. 4 Evolution curve of crack shape

2.3 疊加法求解不同載荷下的應(yīng)力強度因子

裂紋萌生位置是靠近輪軸壓裝配合的S型過渡槽處,而該處的載荷主要來源于壓裝配合載荷和垂向載荷。對于車軸的給定位置,壓裝配合載荷在車軸裝配時已經(jīng)確定,僅與過盈量和摩擦因數(shù)有關(guān),在整個服役周期內(nèi)基本保持不變。而垂向載荷會由于列車的不同運行狀態(tài)(通過曲線、道岔等)而發(fā)生改變。針對垂向載荷,可通過試驗或仿真獲得車軸在服役期間的載荷,再通過雨流計數(shù)法將其轉(zhuǎn)變成不同等級的載荷譜。

由于應(yīng)力強度因子的概念是建立在線彈性斷裂力學(xué)基礎(chǔ)上的,由式(3)可得應(yīng)力強度因子與裂紋尖端附近的應(yīng)力呈線性關(guān)系,那么在相同幾何形狀(包括裂紋幾何形狀)的前提下,可以采用疊加原理,將復(fù)雜載荷下的應(yīng)力強度因子求解轉(zhuǎn)變?yōu)閱蝹€載荷下的應(yīng)力強度因子的求解,然后疊加求和。

因此,有

KTol=KBending+KPF

(9)

式中:KTol為垂向載荷和壓裝配合載荷共同作用下的應(yīng)力強度因子;KBending為垂向載荷作用下的應(yīng)力強度因子;KPF為壓裝配合載荷作用下的應(yīng)力強度因子。

為了驗證疊加法求解應(yīng)力強度因子的可行性,本文比較了2種情況下無裂紋車軸S型過渡槽處沿厚度方向上的應(yīng)力結(jié)果:1) 同時施加壓裝配合載荷和垂向載荷;2) 分別單獨施加壓裝配合載荷(車軸與車輪及齒輪箱設(shè)置過盈接觸,過盈量為0.28 mm,摩擦因數(shù)為0.6)和垂向載荷(F=34 kN,此時不考慮過盈配合,車軸與車輪及齒輪箱共節(jié)點),然后疊加。

應(yīng)力計算結(jié)果如圖5所示,從圖5中可以看出兩種方法得到的應(yīng)力結(jié)果沒有明顯差異,而采用疊加法計算應(yīng)力強度因子不僅能夠?qū)⒑愣ǖ膲貉b配合載荷和變化的垂向載荷(載荷譜中不同的載荷等級)分開,方便計算應(yīng)力強度因子,還能夠顯著降低計算時間,提高計算效率,因此本文采用疊加法計算應(yīng)力強度因子。

圖5 兩種情況下沿厚度方向上的軸向主應(yīng)力比較Fig. 5 Comparison of axial principal stress along the thickness direction in two cases

2.4 含裂紋車軸有限元模型建模

在應(yīng)力強度因子的計算過程中,若裂紋深度發(fā)生改變,只需要將裂紋截面附近的網(wǎng)格重構(gòu),以實現(xiàn)不同裂紋形狀的模擬,上述部分稱為裂紋模塊。此外,車軸其它部分的網(wǎng)格保持不變,使用Tie綁定功能將裂紋模塊與車軸剩余部分網(wǎng)格進行連接,以保持應(yīng)力應(yīng)變傳遞的準確性。為了保證應(yīng)力強度因子的計算精度,裂紋前沿的網(wǎng)格尺寸應(yīng)當(dāng)足夠小以保證應(yīng)力外推或位移外推的準確性,文獻[19]采用應(yīng)力外推法,將裂紋前沿網(wǎng)格尺寸設(shè)置為裂紋深度a的10%。本文經(jīng)過多次仿真,設(shè)置不同的裂紋前沿網(wǎng)格尺寸來計算對應(yīng)的應(yīng)力強度因子,發(fā)現(xiàn)當(dāng)網(wǎng)格尺寸取為裂紋深度的4%時,應(yīng)力強度因子的計算趨于收斂,因此將裂紋前沿網(wǎng)格尺寸設(shè)置為裂紋深度a的4%,具體模型如圖6所示。

圖6 裂紋模塊示意圖Fig. 6 Schematic diagram of crack module

對于給定裂紋形狀的車軸模型,由于采用疊加法求解應(yīng)力強度因子,根據(jù)載荷的不同設(shè)置2組模型:1) 考慮裂紋(使用Tie綁定),不考慮過盈配合載荷(車軸與車輪及齒輪箱共節(jié)點),施加垂向力載荷F=34 kN; 2) 考慮裂紋(使用Tie綁定),施加壓裝配合載荷(車軸與車輪及齒輪箱設(shè)置過盈接觸,過盈量為0.28 mm,摩擦因數(shù)為0.6)。

3 計算結(jié)果分析

依據(jù)2.4節(jié)含裂紋車軸建模方法,建立了不同裂紋形狀的車軸有限元模型,進而求解給定載荷下的應(yīng)力強度因子。本節(jié)首先設(shè)置裂紋形狀比a/c=0.5,裂紋相對深度a/T=0.102,即a=6.426 mm,c=12.852 mm的裂紋,采用位移外推法和應(yīng)力外推法計算垂向載荷和壓裝配合載荷下的應(yīng)力強度因子KBending和KPF;然后將計算結(jié)果與文獻[13-14]對比,驗證了模型和方法的有效性;接著根據(jù)2.2節(jié)設(shè)置的裂紋形狀演變曲線,采用位移外推法計算出隨著裂紋深度變化的應(yīng)力強度因子,基于應(yīng)力強度因子解析式的一般形式,使用五次多項式來擬合給定載荷下的形狀因子函數(shù);最后設(shè)置不同裂紋深度及載荷,將該解析式計算結(jié)果與位移外推法比較,驗證了形狀因子函數(shù)與載荷的無關(guān)性。

3.1 位移外推法及應(yīng)力外推法計算結(jié)果

位移外推法擬合曲線如圖7所示,應(yīng)力外推法擬合曲線如圖8所示。

圖7 位移外推法擬合曲線Fig. 7 Fitted curve by using displacement extrapolation method

圖8 應(yīng)力外推法擬合曲線Fig. 8 Fitted curve by using stress extrapolation method

由于靠近裂紋尖端的計算結(jié)果誤差較大,故舍棄靠近裂紋尖端較近的幾個點,位移外推法從第5個點開始外推,應(yīng)力外推法從第6個點開始外推,采用最小二乘法進行外推,得到直線方程,其截距即為計算得到的應(yīng)力強度因子,如表1所示。

表1 應(yīng)力強度因子計算結(jié)果Tab. 1 Calculation results of stress intensity factors

3.2 計算結(jié)果驗證與對比

文獻[13-14]給出了裂紋位于S型過渡槽處的A點應(yīng)力強度因子表達式,將其變形如下

(10)

式中:YA,Bending和YA,PF為垂向載荷和壓裝配合載荷作用下的形狀因子,與裂紋深度a以及裂紋形狀比a/c有關(guān),可通過查表得到;σBending和σPF為無裂紋車軸垂向載荷和壓裝配合載荷作用下外表面的軸向主應(yīng)力,計算得σBending為22.27 MPa,σPF為66.33 MPa;θ為車軸旋轉(zhuǎn)角度。

根據(jù)文獻[14]查表可得a/c為0.5,a/T為0.102時的形狀因子YA,Bending為0.98,YA,PF為0.629。取θ=0,帶入式(10)可得KA,Bending為3.100 86 MPa·m0.5,KA,PF為5.928 25 MPa·m0.5,KA,Tol為9.029 11 MPa·m0.5,將位移外推法與應(yīng)力外推法的計算結(jié)果和文獻[13-14]所提方法進行對比,結(jié)果如圖9所示。

圖9 應(yīng)力強度因子計算結(jié)果對比Fig. 9 Comparison of stress intensity factor calculation results

從圖9中可以看出,基于位移外推法和應(yīng)力外推法計算得到的KA,Tol,和通過文獻形狀因子公式推導(dǎo)出的結(jié)果相比差別不大,其相對誤差分別為2.92%、-6.26%,驗證了方法和模型的有效性。從中還可以看出基于位移外推法的應(yīng)力強度因子具有更高的精度,這是由于在有限元中,一般是先通過剛度矩陣求解得到位移值,再通過偏導(dǎo)求解得到應(yīng)力值,同時應(yīng)力值對網(wǎng)格大小有很強的敏感性,所以應(yīng)力值通常不如位移值精確[7]。

3.3 形狀因子函數(shù)及應(yīng)力強度因子擬合

基于2.2節(jié)的裂紋形狀演變曲線,采用精度更高的位移外推法計算了隨裂紋深度a變化的應(yīng)力強度因子值。應(yīng)力強度因子的一般解析式為:

(11)

(12)

根據(jù)計算結(jié)果,采用五次多項式來擬合形狀因子函數(shù),即:

(13)

(14)

式中:a為裂紋深度;YA,Bending(a)和YA,PF(a)為垂向載荷和壓裝配合載荷作用下的形狀因子函數(shù);pi和qi為五次多項式擬合常數(shù)(i=0～5)。

將擬合的形狀因子函數(shù)代入到式(11)、式(12)中,可求得任意裂紋深度下的應(yīng)力強度因子,形狀因子函數(shù)擬合常數(shù)如表2所示。

表2 形狀因子函數(shù)擬合常數(shù)Tab. 2 Constants for fitting the shape factor function

形狀因子函數(shù)擬合結(jié)果如圖10所示,應(yīng)力強度因子擬合結(jié)果如圖11所示。

圖10 形狀因子函數(shù)擬合曲線Fig. 10 Fitted curve of shape factor function

圖11 應(yīng)力強度因子擬合曲線Fig. 11 Fitted curve of stress intensity factors

3.4 形狀因子函數(shù)驗證

一般來說,裂紋深度和載荷大小都對應(yīng)力強度因子的計算有重大影響,而式(11)與式(12)的優(yōu)點在于將裂紋深度和載荷這兩項影響因素分隔開。

在給定裂紋形狀演變曲線的基礎(chǔ)上,保持載荷不變(垂向載荷F=34 kN,過盈量為0.28 mm),依次增大裂紋深度,推導(dǎo)出了隨裂紋深度變化的形狀因子函數(shù)。理論上,形狀因子函數(shù)與載荷值無關(guān),即對于不同大小的載荷,式(11)、 式(12)依然適用,僅需要變動與載荷有關(guān)的變量σBending和σPF,為了驗證形狀因子函數(shù)的適用性,保持裂紋形狀演變曲線不變,設(shè)置垂向載荷為F=41 kN,過盈量為0.30 mm,分別使用式(11)、式(12)和位移外推法計算裂紋深度a為3 mm和9 mm下的應(yīng)力強度因子,計算結(jié)果如圖12所示。

圖12 不同載荷下的應(yīng)力強度因子對比Fig. 12 Comparison of stress intensity factors under different loads

從圖12中可以看出,設(shè)置裂紋深度a為3 mm和9 mm時,與位移外推法相比,使用式(11)、式(12)計算的KA,Tol相對誤差分別為-1.91%和-1.39%,這驗證了形狀因子函數(shù)與載荷的無關(guān)性,也就是說,式(11)、式(12)適用于同一載荷模式下任意載荷大小和裂紋深度的應(yīng)力強度因子的計算。

4 結(jié)論

1) 考慮壓裝配合載荷和垂向載荷,S型過渡槽處存在嚴重的應(yīng)力集中,易于萌生裂紋,確定S型過渡槽處為車軸臨界安全位置。

2) 采用疊加法求解車軸應(yīng)力強度因子,不僅能夠?qū)⒑愣ǖ膲貉b配合載荷和變化的垂向載荷分開,還能夠顯著減少計算的時間,提升計算效率。

3) 與文獻形狀因子公式法進行對比,位移外推法和應(yīng)力外推法均能很好地計算車軸應(yīng)力強度因子,其相對誤差分別為2.92%和-6.26%,即位移外推法的精度更高。

4) 根據(jù)位移外推法的計算結(jié)果,采用五次多項式來擬合形狀因子函數(shù),進而推導(dǎo)出給定載荷下的任意裂紋深度的應(yīng)力強度因子的一般解析式,結(jié)果表明該擬合效果良好。

5) 設(shè)置不同的裂紋深度和載荷大小,驗證了形狀因子函數(shù)與載荷的無關(guān)性,表明該解析式能夠?qū)⑤d荷大小和裂紋深度對應(yīng)力強度因子的影響分開,適用于同一載荷模式下任意載荷大小和裂紋深度的應(yīng)力強度因子的計算。