高強化柴油機曲軸過渡圓角應力計算的研究

李 文，廖日東，左正興，劉麗濤

(1.北京控制工程研究所，北京 100190; 2.北京理工大學機械與車輛學院，北京 100081)

前言

隨著柴油機強化程度的不斷提高，曲軸的工作條件愈加苛刻。眾多實驗和數值計算表明［1-3］，曲軸最大應力的位置處于與發火氣缸相對應的連桿軸頸或主軸頸的兩端過渡圓角區域。為了準確進行曲軸的強度可靠性評估，有效制定曲軸過渡圓角區域的強化工藝措施，須對工作狀態下過渡圓角的應力分布情況進行詳細分析［4-6］。曲軸過渡圓角區域的應力分布情況很難通過實驗的方法準確測量［7］，有限元法是目前曲軸應力計算的常用方法。

使用有限元方法進行計算的關鍵是盡可能減小計算誤差。文獻中針對提高曲軸強度計算準確性的研究大多集中在模型的結構簡化和邊界條件施加的正確性方面。文獻［8］中使用連續梁法計算所得彎矩和支反力作為邊界條件對忽略油孔的曲軸單拐進行了靜力學計算，文獻［9］中通過使用曲軸整軸有限元模型提高了靜計算結果的精度，文獻［10］中通過使用彈簧單元模擬軸承的支撐作用，使邊界條件符合實際條件。許多文獻還進一步從曲軸結構的動力學、摩擦學角度對曲軸結構強度進行研究［11-12］。然而，這些工作并沒有深入研究有限元方法計算誤差對結果的影響。

實際上，有限元模型的離散模式也是影響計算結果收斂性的重要因素。因此，本文中以某型號高強化柴油機曲軸為研究對象，在采用合理的模型簡化和邊界條件施加方式的同時，考慮了有限元計算結果在曲軸過渡圓角應力集中區域的收斂性，以期獲得曲軸過渡圓角網格劃分形式與分析精度的關系，降低計算結果的誤差，并在此準則的基礎上建立了具有良好計算收斂性和穩定性的曲軸模型。

1 曲軸單拐及其簡化模型

1.1 曲軸單拐的有限元模型

曲軸單拐結構復雜，在工作過程中承受彎矩、轉矩和軸向力的聯合作用，但主要是爆發壓力造成的沖擊力。曲軸失效的基本形式是彎曲疲勞失效和扭轉疲勞失效，其中彎曲疲勞失效占曲軸失效的80%，因此可認為通過活塞和連桿傳遞給連桿軸頸的氣缸爆發壓力和主軸頸軸承的支反力是曲軸失效的主要原因。現階段對曲軸強度計算所采用的計算模型主要有以下兩種:(1)曲軸的整體模型，該模型能夠充分考慮相鄰曲拐之間相互影響和給定工況下各曲拐的應力情況，但計算規模的限制使軸頸過渡圓角處無法布置足夠數量的網格，以保證應力集中部位計算結果的收斂性;(2)曲軸的單拐模型，該模型通常用于分析受載最嚴重的曲拐，通常利用曲拐的結構和載荷的對稱性，取1/2單拐作為分析對象。

由于本文中研究對象是曲軸過渡圓角區域內的應力分布，因此選取某型號柴油機1/2曲軸單拐模型進行分析，主要考慮彎曲載荷的作用。

邊界條件如圖1(a)所示，在對稱面上施加對稱約束，p為作用在連桿軸頸上的爆發壓力，F1和F2是主軸頸上承受的平衡支反力。載荷沿連桿軸頸和主軸頸軸線方向按二次拋物線規律分布;沿軸頸圓周方向120°范圍內按余弦規律分布［8］。

1.2 曲軸單拐過渡圓角的簡化模型

為獲得模型網格數量和形狀函數階數對曲軸過渡圓角區域計算結果收斂性的關系，選取圖1(a)中所示的灰色區域為研究對象，將該區域簡化為圖1(b)所示的階梯軸，其中R1=R3，R2=R4。根據曲軸主軸頸和連桿軸頸的過渡圓角的實際形狀，簡化的階梯軸模型的過渡圓角分為沉割式(見圖1(d))和非沉割式(見圖1(e))。通過計算結果的對比，兩種形式的過渡圓角具有相同的應力集中收斂性質。文中僅列出了非沉割式圓角收斂性的分析過程。在階梯軸細端端面施加沿豎直方向(z方向)的面力T，相對于對圓角位置形成彎矩，不斷調整該面力的大小，通過多次試算使應力峰值的大小與相同局部網格數量的曲軸單拐模型在最大爆發壓力作用下的圓角應力峰值相近，便于與曲軸單拐模型的計算結果進行比較。當最大爆發壓力為17MPa時，面力T的大小為42.94MPa。

2 簡化模型的收斂性分析

有限元方法的計算程度與試探空間(形狀函數)能否很好地局部逼近精確解息息相關［13］。在應力集中區域內的應力峰值和應力梯度比其它區域大很多，因此需要更多的網格數量或更高的多項式次數，才能獲得更加逼近真值的解。為獲得網格數量和形狀函數階次與應力集中區域計算結果收斂性的關系，選取圖1(b)中灰色區域為網格加密區域，劃分方式如圖1(c)所示，使 L1=L2=L3=L4，其中L1=πr/2為過渡圓角的弧長，r為過渡圓角半徑。在該區域內布置六面體網格，可以發現對于固定過渡圓角半徑和軸頸半徑的階梯軸，影響網格尺寸的變量一共有3個，分別是沿過渡圓角周向α的網格數m，過渡圓角深度方向ρ1的網格數n，軸頸周向φ的網格數l。本文中將通過一系列計算與分析對m、n、l的取值與結果收斂性進行討論。本文中所有有限元計算均使用ABAQUS standard求解器進行求解。計算中彈性模量E=210GPa，泊松比ν=0.3。

2.1 m，n取值對結果收斂性的影響

分別變化m和n值得到應力集中處的Mises應力值，見表1。由表可見:m一定時，隨著n的不斷增大，應力集中點處的應力峰值不斷增大;n一定時，改變周向網格數m，當m＜6時，應力峰值隨著m值增加而迅速增大;當m≥6時，隨著周向網格數的增加，最大應力峰值在一定范圍內波動。因此，過渡圓角周向網格數m≥6時，網格數的增加對結果收斂性的影響較小。應力值的小幅變化是由于過渡圓角周向網格數目的增加引起的節點位置變化造成的。可見過渡圓角深度方向網格數是決定應力峰值收斂性的主要因素。

表1 過渡圓角處m，n取值與最大Mises應力σAmax的關系MPa

在模型建立過程中，單元質量也是必須考慮的重要因素。低質量的網格單元會造成較大的計算誤差或結果的不收斂，因此n與m不宜相差過大。通常情況下，當n=m時能獲得較好的網格單元質量。

由表1可以得到網格數與收斂誤差的關系。將m與n均為30時的最大應力值379.30MPa作為基準，可以發現對于六面體一次單元，當m與n均為6時即可滿足工程上誤差為10%的要求;當m與n均劃分為10，即可滿足工程上誤差為5%以內的高精度要求。

2.2 l取值對結果收斂性的影響

固定m和n值(m=n=12)，改變 l值，可算得過渡圓角區域中的應力峰值σAmax，如表2所示。當該結構僅承受彎矩作用時，過渡圓角處的最大應力值隨l值的變化并不顯著。當l值由72變為6時，應力只提高了0.49%。故可認為，在保證單元質量的前提下，軸頸周向網格數對結果不產生顯著影響。

表2 m與n為12時環向網格密度變化時的最大Mises應力σAmax MPa

2.3 收斂準則的適用性分析

為討論上述收斂性準則對不同圓角幾何尺寸的適用性，保持邊界條件和網格劃分方式不變，改變過渡圓角半徑r，令m=n且將其數值從2逐漸增加至30，可以獲得r取值不同時網格數量與應力峰值的關系，如表3所示。由表可見:圓角半徑r取值的不同改變了應力峰值的大小，但通過對表3的數據進行歸一化處理發現，隨著r取值不同網格劃分數量與結果收斂性的關系并沒有多大變化，說明所獲得的收斂性準則適用于不同尺寸的過渡圓角。

表3 m與n相等時，圓角半徑與最大Mises應力σAmax 的關系 MPa

3 階梯軸與曲軸單拐的應力對比

應力集中造成的數值計算誤差主要是由于高應力梯度造成的。為了證明由階梯軸模型得到的過渡圓角區域承受彎矩作用時收斂準則能夠適用于曲軸單拐，將階梯軸與曲軸單拐的過渡圓角處應力分布與應力梯度進行對比。其中階梯軸模型的過渡圓角分為沉割式和普通式兩種。

曲軸單拐過渡圓角采取與階梯軸相同的區域選擇和劃分方式，且階梯軸的軸頸半徑與過渡圓角半徑之比和曲軸單拐主軸頸與過渡圓角之比相同。為了獲得較好的收斂性，曲軸單拐和階梯軸的過渡圓角周向網格數和深度方向網格數均取12，并且嚴格控制過渡圓角區域處的節點位置，在過渡圓角區域分別建立局部坐標系O3和O4，如圖2所示。

選取應力峰值點處的沿過渡圓角周向、深度方向和軸頸周向的應力分布進行對比，可以發現以下規律。

(1)無論是沉割式還是普通過渡圓角，在過渡圓角表面，階梯軸和曲軸單拐沿過渡圓角周向的應力分布曲線形狀相同，但應力峰值的位置不同，如圖3(a)所示。

(2)在由應力峰值點為起點的沿過渡圓角深度方向ρ1的路徑上，沉割式過渡圓角與普通過渡圓角的應力分布曲線基本重合，如圖3(b)所示。應力值由過渡圓角表面沿ρ1方向迅速降低，曲軸單拐過渡圓角應力峰值的應力分布曲線與階梯軸的應力分布曲線整體趨勢相同，應力的下降速度稍快于階梯軸。

(3)在以應力峰值點為起點的沿軸頸周向φ的路徑上，沉割式過渡圓角與普通過渡圓角的應力分布曲線完全重合，與曲軸單拐過渡圓角的應力分布差別較大，如圖3(c)所示。

(4)將階梯軸和曲軸單拐過渡圓角區域的應力梯度進行對比，如圖3(d)所示，選取各模型過渡圓角應力峰值點處沿軸頸周向、過渡圓角周向和深度方向的應力曲線，求其沿各自方向的梯度。應力梯度的計算公式為

式中:q為路徑方向，qi為路徑q上的單元節點，σqi為節點qi處的Mises應力值，li為節點i至路徑原點處的路徑長度。由于各個方向的路徑長度相差較大，為便于將各方向的應力梯度進行對比，選取lqi/lq作為橫坐標，Kqi值為縱坐標，其中lq為q方向的路徑總長度。通過對比可以發現，沿軸頸周向，各模型的應力梯度遠遠小于過渡圓角周向和深度方向的應力梯度;沿過渡圓角周向，各模型的應力梯度值為幅值在251.4～-231.18MPa/mm之間連續變化的曲線;沿過渡圓角深度方向，各模型的應力梯度曲線均為一條單調遞減曲線，且其應力峰值基本相同。

由上述規律可以得到以下結論。

(1)曲軸單拐模型和其過渡圓角區域的簡化模型在彎矩作用下，除了在過渡圓角周向上應力集中點位置不同之外，沿各方向的應力分布曲線較為接近，且沿各方向應力梯度水平分別相同。

(2)沉割式階梯軸模型和普通階梯軸模型在彎矩作用下，除了在過渡圓角周向上應力集中點位置不同之外，沿各方向的應力分布和應力梯度分布基本相同。因此，在研究過渡圓角應力集中區域收斂性的過程中，將沉割式階梯軸簡化為普通階梯軸是合理的。

(3)各模型沿軸頸周向的應力梯度遠遠小于其它兩個方向，因此在該方向的網格分布數量對結果的收斂性基本不產生影響，與2.2節中的數值試驗結果相符。

(4)各模型沿過渡圓角深度方向的應力梯度峰值遠遠大于其它兩個方向，且應力梯度本身的變化率也較大，因此在該方向的網格分布數量是影響結果的收斂性的主要因素，與2.1節中的數值試驗結果相符。

因此，由階梯軸模型得到的彎曲載荷作用下的過渡圓角區域網格收斂性準則適用于曲軸單拐主軸頸過渡圓角區域和連桿軸頸過渡圓角區域。

4 曲軸單拐應力分布特征

由上節結論可知，過渡圓角區域周向和深度方向網格數均為12時的1/2曲軸單拐有限元模型具有較高的計算精度，因此，采用該網格劃分方式對上文所研究曲軸的主軸頸和連桿軸頸過渡圓角區域的應力分布情況進行分析。

4.1 主軸頸和連桿軸頸過渡圓角應力比較

圖4為該曲軸單拐在載荷作用下的Mises應力分布情況。由圖可見，過渡圓角區域的應力水平較高，其中主軸頸過渡圓角位置承受壓應力，連桿軸頸過渡圓角位置承受拉應力。

圖5(a)中曲線是在主軸頸過渡圓角圓心處的局部坐標系O3下，0＜β＜160°范圍內，距過渡圓角表面0～3mm范圍內(3mm＜ρ2＜6mm)的應力分布情況;圖5(b)中曲線是在連桿軸頸過渡圓角圓心處的局部坐標系O4下，0＜γ＜90°范圍內，距過渡圓角表面0～3mm范圍內(3mm＜ρ3＜6mm)的應力分布情況，其應力梯度的分布趨勢與圖5(a)中的應力分布趨勢類似，但是其應力梯度與應力峰值均小于圖5(a)中的應力分布曲線。由此可以得出，在過渡圓角周向方向，主軸頸過渡圓角的應力集中較為嚴重。

為了進一步研究應力沿曲軸過渡圓角表面深度方向的變化情況，分別對主軸頸過渡圓角表面和連桿軸頸過渡圓角表面沿深度方向選取13條路徑，得到過渡圓角不同位置處的沿深度方向的應力分布曲線，見圖5(c)和圖5(d)。經過對比可以發現:主軸頸過渡圓角區域在表面處的應力值最大，隨著深度增加應力迅速減小，在距離表面3～4mm處降至較低水平，隨后應力值變化趨于平緩。連桿軸頸過渡圓角區域與主軸頸過渡圓角區域的應力分布情況類似，但是其應力峰值較低。

圖5(e)中的曲線為過渡圓角表面Mises應力沿主軸頸周向的應力分布情況。應力值在φ=0～90°范圍內較大，在φ=90°～180°范圍內較小。圖5(f)中的曲線為過渡圓角表面Mises應力沿連桿軸頸周向的應力分布情況。應力集中發生在γ=37.5°，θ=180°處，在 θ=0～110°范圍內，由于外載荷的作用，存在一定程度的應力變化，在θ=110°～180°范圍內應力較大，并且在θ=180°處達到應力最大值。

4.2 主軸頸和連桿軸頸過渡圓角處的應力梯度

選取主軸頸和連桿軸頸過渡圓角應力峰值點處沿軸頸周向、過渡圓角周向和過渡圓角深度方向的應力曲線，使用式(3)求其沿各自方向的梯度，如圖6所示。

Kφ和Kθ為主軸頸和連桿軸頸周向的應力梯度。其中|Kφ|的最大值為15.034MPa/mm，|Kθ|最大值為15.788MPa/mm。Kβ和Kγ分別為主軸頸和連桿軸頸過渡圓角周向的應力梯度，其中Kβ為一條幅值在200.73～-231.18MPa/mm之間連續變化的曲線;Kγ是一條單調遞減的曲線，其最大值為165.09MPa/mm，最小值為 -155.13MPa/mm，兩條曲線的應力梯度峰值相差不大，但是由于受力狀態不同，其應力梯度分布形式不同。Kρ2和Kρ3為兩條單調遞增曲線，其中Kρ2最小值為-526.75MPa/mm，Kρ3最小值為-349.26MPa/mm，兩條曲線快速增加至-100MPa/mm之后緩慢上升并趨向于0。

由于連桿軸頸沿各方向的應力梯度與主軸頸的應力梯度水平相差不大，因此得到的各方向網格數與收斂性的關系也適用于連桿軸頸。

進一步對應力梯度的研究可以發現，沿路徑方向，Kρ2和 Kρ3逐漸趨近于零，而 Kβ和 Kγ在路徑兩端的幅值較大，其中Kγ在路徑的兩端分別取最大值和最小值，說明圖2中所示的網格區域并未完全包括該結構的大應力梯度區域，因此須擴大網格細化區域。由于正六面體網格的收斂性最好，因此選取圖1(c)中與圓角區域相接的正方形區域進行網格細化，保證各方向網格數相同，使網格為正方形。經過分析發現，該區域在曲軸內部的邊界上應力梯度很小，因此可以認為，本文中選取的網格細化區域是合理有效的。

4.3 曲軸過渡圓角各方向應力分布表達式

由圖5(e)可知，沿主軸頸表面周向，φ=0截面上的應力水平最高，且此處沿軸頸周向的應力梯度最小，因此對該截面的應力分布進行研究。在極坐標系 O3下，10°≤β≤160°，3mm≤ρ2≤6mm 范圍內，假設法向力和切向力的分布可表示為如下的傅里葉級數:

使用MATLAB對以上各式進行數據擬合，通過多次試算，確定此處n=1時，以上各式能夠較好地表現過渡圓角區域內的法向力和切向力的分布。因此該模型過渡圓角處的法向力和切向力分布為

式中:A0，A1，B1，C0，C1，D1，E0，E1，F1為 r/ρ2和 σm的函數，ωβ，ωρ2，ωβρ2為 r/ρ2的函數;σm為假設連桿=軸頸承受均布載荷時的等效外載荷，即σm通過計算可得該型曲軸的σm=70.4887MPa/mm2;r=3mm為主軸頸過渡圓角半徑。

當 ρ2分別取3，3.5，4，4.5，5，5.5，6mm 時，取過渡圓角周向均勻分布的12個節點上的應力值，對式(7)～式(9)進行數值擬合，得到ρ2不同時各式的參數取值。數據點與擬合公式的相關程度使用相關系數平方δ2表示，如表4所示。由表4可以看出，擬合的應力函數與數據點具有很好的相關性。因此，可以認為采用傅里葉級數展開的形式表示曲軸過渡圓角法向力和切向力的分布是合理的。對式(7)～式(9)中的各參數以r/ρ2

表4 主軸頸應力分布的相關性系數平方

為變量采用多項式形式進行擬合，得到各參數的多項式形式的表達式為

擬合所得表達式與數據點之間的關系通過相關系數平方δ2表示，如表5所示。由表5可知，擬合的參數表達式與數據點具有很好的相關性，表達式與數據點間的偏差較小。同樣的方法可以表示連桿軸頸過渡圓角的應力分布。

表5 主軸頸應力分布各參數擬合的相關性系數平方

由于曲柄臂厚度、軸頸重合度等幾何參數的影響，不同曲軸的連桿軸頸過渡圓角和主軸頸過渡圓角的相互影響程度各不相同。本文中雖然僅針對某一特定曲軸進行了討論，但所得到的規律具有一定的參考價值，針對不同幾何特征尺寸的曲軸均可以采用上述方法進行具體的研究，從而獲得準確的應力分布表達式。

5 結論

(1)對于曲軸過渡圓角區域應力值的有限元計算，網格模型應采用以下收斂性準則:過渡圓角及其周圍一定區域必須進行網格細化;過渡圓角周向網格數取值應該大于等于6，過渡圓角深度方向網格數應大于或等于過渡圓角周向網格數，周向與深度方向網格數均劃分為10時即可滿足工程上誤差為5%以內的高精度要求;曲軸軸頸周向網格數對過渡圓角區域的應力逼近特性不產生顯著影響，在保證單元質量的前提下，可使用較少的網格數。

(2)對于本文中使用的高強化柴油機曲軸模型，主軸頸和連桿軸頸過渡圓角區域應力峰值所在截面的各方向應力分布可使用以過渡圓角角度為變量的一次傅里葉級數的形式表示，其表達式和數據點的相關性較好。

(3)該研究方法能夠運用于其他結構應力集中位置的收斂性分析，通過擬合得到的危險截面過渡圓角區域各方向應力的表達式的方法可廣泛應用于曲軸過渡圓角區域強化工藝方面的研究。

［1］Osman A.Failure Analysis of a Crankshaft Made from Ductile Cast Iron［J］.Engineering Failure Analysis，2006，13:1260-1267.

［2］Pandey R K.Failure of Diesel Engine Crankshafts［J］.Eng Fail Anal，2003，10:165-175.

［3］Yu Z，Xu X.Failure Analysis of a Diesel Engine Crankshaft［J］.Eng Fail Anal，2005，12:487-495.

［4］Choi K S，Pan J.Simulations of Stress Distributions in Crankshaft Sections Under Fillet Rolling and Bending Fatigue Tests［J］.Int J Fatigue，2009，31:544-557.

［5］Ren W，Li K，Lee Y.Optical Measurement of Residual Stress at the Deep-rolled Crankshaft Fillet［C］.SAE Paper 2004-01-1500.

［6］Guo Y B，Barkey M E.FE-Simulation of the Effects of Machininginduced Residual Stress Profile on Rolling Contact of Hard Machined Components［J］.Int J Mech Sci，2004，46:371-388.

［7］Lee Y L，Morrissey W.Uncertainties of Experimental Crankshaft Fatigue Strength Assessment［J］.Int J Mater.Prod.Tech，2001，16:379-392.

［8］尹建民，王德海，袁銀南.X6135柴油機曲軸強度的三維有限元研究［J］.內燃機工程，1997(2):71-77.

［9］馮國勝，張幽彤，張玉申.柴油機曲軸靜動特性的三維有限元分析［J］.內燃機工程，2003，24(2):73-77.

［10］沈海濤，鄭水英，李志海.基于彈簧支承的柴油機曲軸強度有限元分析［J］.機械強度，2007，29(1):161-164.

［11］Zissimos P M.A Crankshaft System Model for Structure Dynamic Analysis of Internal Combustion Engines［J］.Computers and Structures，2001，79:2009-2027.

［12］何芝仙，桂長林，李震，等.計入軸瓦變形的曲軸動應力和疲勞強度計算［J］.機械工程學報，2009，45(11):91-97.

［13］Kelly D W，De J P，Gago S R.A Posteriori Error Analysis and Adaptive Processes in the Finite Element Method:Part I-Error A-nalysis［J］.Int J Num Meth Eng，1983，19:1621-1656.