龔京風(fēng)，徐宗著，宣領(lǐng)寬

(1.武漢科技大學(xué)汽車與交通工程學(xué)院，湖北 武漢，430065；2.中國艦船研究設(shè)計(jì)中心船舶振動(dòng)噪聲重點(diǎn)實(shí)驗(yàn)室，湖北 武漢，430064)

活塞是柴油機(jī)的核心部件，其長期處于高溫、高壓的工作環(huán)境且承受應(yīng)力水平較高的周期性載荷，活塞的抗疲勞性能直接影響到柴油機(jī)的工作可靠性?；钊诠ぷ鬟^程中的疲勞損傷包括高周疲勞損傷和低周疲勞損傷，前者主要是指柴油機(jī)因高頻變化的溫度、燃?xì)獗l(fā)壓力及往復(fù)慣性力而導(dǎo)致的損傷，后者主要是指柴油機(jī)在啟動(dòng)、停機(jī)等工況變化中因活塞溫度大范圍變化引起的塑性變形和蠕變等損傷[1]。由于柴油機(jī)長時(shí)間處于工作循環(huán)中，而啟停工況下的運(yùn)行時(shí)間較短，因此本文主要對活塞高周疲勞進(jìn)行分析。

在計(jì)算力學(xué)領(lǐng)域，除了FEM之外，有限體積法(FVM)也得到了廣泛關(guān)注，尤其是格點(diǎn)型有限體積法(CV-FVM)。Xia等[8]將CV-FVM應(yīng)用于二維和三維懸臂梁的結(jié)構(gòu)動(dòng)力學(xué)計(jì)算，所得到的解與解析解一致。陳浩[9]基于確定性結(jié)構(gòu)動(dòng)響應(yīng)的控制方程推導(dǎo)出隨機(jī)結(jié)構(gòu)動(dòng)響應(yīng)的控制方程，并使用CV-FVM對控制方程進(jìn)行離散，給出了隨機(jī)結(jié)構(gòu)動(dòng)響應(yīng)控制方程的求解過程。井麗龍[10]通過直梁靜力彎曲分析、固有頻率分析以及集中載荷動(dòng)力響應(yīng)分析，驗(yàn)證了CV-FVM在梁、板、殼這類結(jié)構(gòu)的靜力學(xué)和動(dòng)力學(xué)分析中的準(zhǔn)確性。筆者等[11-13]使用CV-FVM對熱彈性問題進(jìn)行了研究：在文獻(xiàn)[11]中利用雙線性四邊形單元和三角形單元的混合網(wǎng)格分析功能梯度材料的熱彈性問題，發(fā)現(xiàn)相比于梯度單元和高階單元，CV-FVM可以減少計(jì)算量和內(nèi)存消耗；在文獻(xiàn)[12]中研究了各向異性材料的熱彈性問題，拓展了CV-FVM的應(yīng)用范圍；在文獻(xiàn)[13]中研究了考慮溫度熱效應(yīng)的功能梯度材料的熱彈性問題，進(jìn)一步驗(yàn)證了CV-FVM的計(jì)算精度和適用性。宣領(lǐng)寬[14]在研究結(jié)構(gòu)聲耦合問題時(shí)采用三角形單元、四邊形單元、四面體單元、三棱柱單元和六面體單元?jiǎng)澐钟?jì)算域，基于CV-FVM離散各向異性材料的彈性力學(xué)控制方程，分別使用顯式和隱式算法求解離散方程，并將該方法得到的結(jié)果與商用軟件計(jì)算結(jié)果進(jìn)行了對比。

1 高周疲勞數(shù)值模型的建立

1.1 CV-FVM控制體簡介

以圖1中節(jié)點(diǎn)N1為例對CV-FVM控制體進(jìn)行說明。C1、C2分別為N1的兩個(gè)相鄰單元1和2的中心；F1、F2、F3、F4、F5分別為N1相鄰網(wǎng)格面單元的中心，L1、L2、L3、L4分別為N1相鄰網(wǎng)格單元邊線的中心。依次連接這些中心點(diǎn)，即為N1的控制體，其中，四邊形C2F4L3F5、C2F4L4F3、C2F5L2F3為N1在單元2內(nèi)的控制體的邊界面，四邊形C1F1L1F2、C1F5L3F1、C1F2L2F5為N1在單元1內(nèi)的控制體的邊界面。

圖1 CV-FVM控制體示意圖

1.2 彈性方程的離散

設(shè)控制體的體積為V、邊界面積為S，根據(jù)力平衡條件建立彈性方程：

(1)

式中：u、v、w分別為坐標(biāo)軸x、y、z方向的位移分量；nx、ny、nz分別為控制體邊界面的法向矢量在x、y、z方向的分量；fx、fy、fz分別為體積力在x、y、z方向的分量；ρ為密度；σ表示應(yīng)力。

由文獻(xiàn)[15]中的本構(gòu)方程和幾何方程得

(2)

將式(2)代入式(1)中，考慮給定面力邊界條件的影響，得

(3)

(4)

(5)

式(3～5)中：pxB、pyB、pzB分別為x、y、z方向的面力；SN為邊界面的面積。當(dāng)給定位移邊界uiB(i=x,y,z)時(shí)，邊界上的節(jié)點(diǎn)位移滿足ui=uiB，可采用文獻(xiàn)[16]中的置大數(shù)法處理給定位移邊界。

以式(3)為例，對彈性問題的控制方程進(jìn)行離散，將該式右端第一項(xiàng)和第二項(xiàng)離散為

(6)

式中：nc為當(dāng)前節(jié)點(diǎn)的相鄰單元總數(shù)；ncni為當(dāng)前節(jié)點(diǎn)第i個(gè)相鄰單元上的節(jié)點(diǎn)總數(shù)；Si為當(dāng)前節(jié)點(diǎn)在相鄰單元i內(nèi)的控制體邊界面積；Nij為當(dāng)前節(jié)點(diǎn)第i個(gè)相鄰單元上第j個(gè)節(jié)點(diǎn)處的形函數(shù)。關(guān)于不同單元類型的形函數(shù)導(dǎo)數(shù)的積分計(jì)算方法見文獻(xiàn)[14]。

將式(3)右端第三項(xiàng)離散為

(7)

式中：Vi為當(dāng)前節(jié)點(diǎn)第i個(gè)相鄰單元的體積；fxi為當(dāng)前節(jié)點(diǎn)第i個(gè)相鄰單元在x方向的體積力。

將式(3)右端第四項(xiàng)離散為

(8)

式中：nN為當(dāng)前節(jié)點(diǎn)相鄰邊界網(wǎng)格面的個(gè)數(shù)；ABi為當(dāng)前節(jié)點(diǎn)在第i個(gè)相鄰邊界網(wǎng)格面內(nèi)的控制體邊界面積。

對式(3)左端項(xiàng)采用中心差分法離散為

(9)

本文主要解決靜態(tài)彈性問題且不考慮體積力，無需考慮時(shí)間項(xiàng)和體積源項(xiàng)，因此式(3)的離散方程可簡化為

(10)

式(4)和式(5)的離散方程可用相同方式得到：

(11)

(12)

將計(jì)算域劃分為具有n個(gè)節(jié)點(diǎn)的網(wǎng)格模型，對每個(gè)節(jié)點(diǎn)聯(lián)立方程式(10)～(12)，得到一個(gè)求解彈性問題的線性方程組，表示為

AX=B

(13)

其中

A=

式中：列向量B由方程中已知邊界和體積力求得，其中下標(biāo)un、vn、wn表示各位移方向上不同節(jié)點(diǎn)對應(yīng)的已知系數(shù)；系數(shù)矩陣A中各節(jié)點(diǎn)上未知變量的系數(shù)由形函數(shù)導(dǎo)數(shù)積分以及部分給定位移邊界求得，其中下標(biāo)un1、vn1、wn1表示各位移方向不同節(jié)點(diǎn)上的第一個(gè)方程對應(yīng)的未知量系數(shù)，un2、vn2、wn2、un3、vn3、wn3定義類似。

由幾何方程以及形函數(shù)可求得各節(jié)點(diǎn)上的應(yīng)變分量：

(14)

再由式(2)計(jì)算得到計(jì)算域內(nèi)各節(jié)點(diǎn)的應(yīng)力分量，同時(shí)可由式(15)求得von Mises等效應(yīng)力

(15)

1.3 高周疲勞壽命計(jì)算方法

在疲勞壽命分析中，對于恒幅應(yīng)力循環(huán)載荷，描述載荷譜的參數(shù)主要有最大應(yīng)力σmax和最小應(yīng)力σmin，同時(shí)可以根據(jù)σmax和σmin求得應(yīng)力幅σa、平均應(yīng)力σm和應(yīng)力比R，當(dāng)R=-1時(shí)，稱為對稱循環(huán)載荷。

高周疲勞壽命分析時(shí)，通常采用S-N曲線來表征應(yīng)力幅與疲勞壽命之間的關(guān)系，恒幅應(yīng)力循環(huán)載荷作用下，采用應(yīng)力比和應(yīng)力幅描述循環(huán)應(yīng)力水平，當(dāng)已知應(yīng)力比時(shí)，應(yīng)力幅即為主要評價(jià)參量。描述材料疲勞性能的基本S-N曲線大多是根據(jù)試驗(yàn)得到的S-N數(shù)據(jù)通過擬合得到曲線表達(dá)式。本文使用分段線性插值，即每兩個(gè)數(shù)據(jù)點(diǎn)使用線性函數(shù)擬合，得到S-N曲線表達(dá)式(16)：

σa=f(N)

(16)

式中：N為材料疲勞壽命。

材料應(yīng)力和疲勞壽命數(shù)據(jù)都是利用光滑小尺寸試件在對稱循環(huán)載荷作用下得到的，對于大尺寸構(gòu)件，需要引入一個(gè)小于1的尺寸修正因子；類似地，表面粗糙度的影響也可以引入一個(gè)小于1的修正因子。將這兩種影響因素綜合考慮，定義為疲勞強(qiáng)度因子kf，代入到式(16)得

σa=kff(N)

(17)

文獻(xiàn)[17]指出，在相同的應(yīng)力幅下，隨著平均應(yīng)力的上升，構(gòu)件的疲勞壽命將縮短，因此在進(jìn)行疲勞壽命分析時(shí)需要考慮平均應(yīng)力的影響?？刹捎靡韵聝煞N方法對應(yīng)力幅進(jìn)行修正。

(1)使用Goodman理論修正：

(18)

(2)使用Gerber理論修正：

(19)

相比較來看，平均應(yīng)力增大到一定程度后，Goodman理論修正的應(yīng)力幅值要高于Gerber理論修正的應(yīng)力幅值，即Goodman理論比Gerber理論在疲勞壽命預(yù)測上更加保守。在工程實(shí)際中往往采用Goodman理論對應(yīng)力幅值進(jìn)行修正。

對于靜態(tài)彈性問題，給定最大外載荷，將計(jì)算得到的σvon_Mises作為載荷譜最大應(yīng)力，再由給定的應(yīng)力比計(jì)算得到應(yīng)力幅，即可對構(gòu)件的疲勞壽命進(jìn)行預(yù)測。

2 數(shù)值模型的驗(yàn)證

在板左側(cè)分別施加90、100、110 MPa的法向壓力載荷。板右側(cè)和上側(cè)設(shè)置為對稱邊界，下側(cè)為自由邊界，并指定應(yīng)力比R=0.1，kf=1，使用Goodman修正理論進(jìn)行應(yīng)力幅修正，驗(yàn)證模型的S-N曲線見圖3。

圖2 驗(yàn)證模型

圖3 驗(yàn)證模型的S -N曲線

分別使用FEM和CV-FVM方法計(jì)算該模型的應(yīng)力，并將最大等效應(yīng)力計(jì)算結(jié)果與文獻(xiàn)[18]中的研究結(jié)果進(jìn)行對比，見表1。

表1 最大等效應(yīng)力計(jì)算結(jié)果

從表1可以看出，在3種壓力載荷作用下，CV-FVM法與其他兩種方法的計(jì)算結(jié)果誤差在1%以內(nèi)，驗(yàn)證了CV-FVM方法的正確性。分別使用FEM和CV-FVM方法計(jì)算在110 MPa法向壓力載荷作用下矩形板的疲勞壽命，結(jié)果如圖4所示。

(a) FEM計(jì)算結(jié)果

(b) CV-FVM計(jì)算結(jié)果

使用FEM和CV-FVM方法計(jì)算得到的最小疲勞壽命分別為221 000和226 758循環(huán)周次，兩者數(shù)值極為接近，疲勞壽命云圖分布也呈現(xiàn)出相同規(guī)律，這充分驗(yàn)證了本文提出的高周疲勞數(shù)值模型的正確性。

3 活塞高周疲勞壽命分析

3.1 活塞計(jì)算模型及邊界條件

活塞材料為ZL109鋁合金，密度為2770 kg/m3，楊氏模量為69 GPa，泊松比為0.33。忽略倒角等幾何特征，采用四面體單元?jiǎng)澐钟?jì)算域?；钊?jì)算模型見圖5，該模型為完整活塞的1/4模型。

圖5 活塞計(jì)算模型

研究活塞疲勞問題時(shí)，活塞受到的燃?xì)鈮毫d荷可以近似為正弦脈沖循環(huán)載荷[19]。本文在給定最大燃?xì)鈮毫g和應(yīng)力比的基礎(chǔ)上，基于靜態(tài)應(yīng)力結(jié)果模擬計(jì)算正弦脈沖循環(huán)載荷作用下的活塞疲勞壽命。在活塞頂部、火力岸及第一道環(huán)上表面施加157.3 MPa的燃?xì)鈮毫?，第一環(huán)槽內(nèi)側(cè)及下表面按燃?xì)鈮毫Φ?5%施加，第一環(huán)岸、第二環(huán)槽上表面及下表面按燃?xì)鈮毫Φ?5%施加，第二環(huán)槽內(nèi)側(cè)按燃?xì)鈮毫Φ?0%施加[20]。同時(shí)還考慮此時(shí)的支反力FN作用在活塞銷孔上表面的表面壓力PN，支反力與表面壓力的關(guān)系滿足下式：

(20)

式中：l為銷座長度；r*為銷孔半徑。

所研究的活塞中，一側(cè)銷孔的支反力FN=177 kN，l=59.7 mm，r*=65 mm，可得PN=12.83 MPa。為約束活塞的剛體位移，在兩個(gè)剖面分別施加對稱邊界條件，銷孔上表面約束活塞軸線方向位移。

活塞疲勞計(jì)算中給定應(yīng)力比R=0.01，極限強(qiáng)度σu=239.3 MPa，S -N曲線上部分?jǐn)?shù)據(jù)見表2。

表2 活塞S -N曲線上部分?jǐn)?shù)據(jù)

活塞尺寸與材料疲勞性能試件的尺寸存在較大差異，需考慮尺寸效應(yīng)，對于圓截面構(gòu)件，可由式(21)計(jì)算疲勞強(qiáng)度因子[17]?；钊诩庸み^程進(jìn)行了拋光處理，表面粗糙度修正因子為1，即可忽略表面粗糙度對活塞疲勞壽命的影響[21]。

(21)

式中：d為活塞直徑，本文中d=0.0835 m。計(jì)算得到kf=0.77。

3.2 活塞高周疲勞計(jì)算結(jié)果

按前述邊界條件，分別使用FEM和CV-FVM計(jì)算活塞應(yīng)力場，得到活塞的總位移云圖(見圖6)。從圖6可以看出，兩種方法計(jì)算得到的云圖在活塞位移分布規(guī)律上高度吻合?；钊畲笪灰浦饕性陬^部區(qū)域，活塞銷孔位置由于y方向的位移被約束，最小位移出現(xiàn)于此處。

(a) FEM計(jì)算結(jié)果

(b) CV-FVM計(jì)算結(jié)果

圖7為兩種方法計(jì)算得到的等效應(yīng)力云圖，其中FEM計(jì)算得到最大等效應(yīng)力為86.1 MPa，CV-FVM計(jì)算得到的最大等效應(yīng)力為86.5 MPa，兩種方法的計(jì)算結(jié)果均在同一范圍內(nèi)變化，云圖分布規(guī)律也一致，最大等效應(yīng)力主要集中在活塞銷孔上表面內(nèi)側(cè)區(qū)域。

(a)FEM計(jì)算結(jié)果

(b) CV-FVM計(jì)算結(jié)果

在銷孔上側(cè)邊沿取特征點(diǎn)(見圖8)上的等效應(yīng)力值，比較兩種方法在數(shù)值上的精度，見圖9。從圖9中可以看出，CV-FVM和FEM計(jì)算得到的等效應(yīng)力數(shù)值吻合良好，驗(yàn)證了使用CV-FVM計(jì)算活塞應(yīng)力場的準(zhǔn)確性。

圖8 特征點(diǎn)示意圖

圖9 特征點(diǎn)上的等效應(yīng)力

在上述計(jì)算基礎(chǔ)上，使用Goodman和Gerber平均應(yīng)力修正理論計(jì)算活塞高周疲勞壽命，結(jié)果分別為8.88×1014循環(huán)周次和3.69×1015循環(huán)周次。為保證構(gòu)件的使用安全，一般對疲勞壽命進(jìn)行保守估計(jì)，因此這里使用Goodman修正方法的計(jì)算結(jié)果。圖10分別為FEM和CV-FVM計(jì)算得到的疲勞壽命云圖。FEM所得最小疲勞壽命為1.01×1015循環(huán)周次，CV-FVM所得最小疲勞壽命為8.88×1014循環(huán)周次，二者基本一致。在機(jī)械載荷的作用下，活塞最小疲勞壽命出現(xiàn)在銷孔上表面內(nèi)側(cè)區(qū)域，與最大等效應(yīng)力位于同一區(qū)域，即活塞等效應(yīng)力與疲勞壽命成反比，考察機(jī)械載荷作用下的活塞疲勞破壞，要重點(diǎn)關(guān)注其銷孔處的應(yīng)力情況。

(a)FEM計(jì)算結(jié)果

(b) CV-FVM計(jì)算結(jié)果

4 結(jié)語

使用本文數(shù)值模型預(yù)測活塞高周疲勞壽命，其結(jié)果與有限元法計(jì)算結(jié)果誤差很小。在機(jī)械載荷作用下，活塞等效應(yīng)力與疲勞壽命成反比，最小疲勞壽命出現(xiàn)在活塞銷孔上表面內(nèi)側(cè)區(qū)域，與最大等效應(yīng)力出現(xiàn)的區(qū)域一致。