基于響應面法和MOGA 的重載曲軸形貌優(yōu)化研究

馬文生，龔鈔，李方忠，劉春川

（1.重慶理工大學機械工程學院，重慶 400054；2.重慶水泵廠有限責任公司國家企業(yè)技術中心，重慶 400033；3.哈爾濱工程大學航天與建筑工程學院，黑龍江哈爾濱 150006）

0 前言

大型隔膜泵主要用于石油、化工、采礦等領域，在高磨蝕、高濃度固液兩相介質(zhì)長距離管道化輸送領域占據(jù)絕對優(yōu)勢地位［1-2］。其中重載曲軸作為大型隔膜泵動力端的關鍵零部件，承受著由活塞端帶來的連桿力、支撐軸頸處的摩擦力以及自身旋轉慣性力等引起的周期性載荷。因此，研究曲軸的強度、剛度等很有必要。

近年來，針對曲軸的優(yōu)化，國內(nèi)外學者進行了大量的相關研究。游孟平等［3］通過建立曲軸的剛柔耦合動力學模型進行動力學仿真，研究曲軸表面粗糙度對疲勞安全系數(shù)及疲勞壽命的影響，結果表明：曲軸表面粗糙度越大，曲軸的疲勞壽命越短。劉潔、盧秋霞［4］通過靜強度分析及模態(tài)分析對六缸發(fā)動機曲軸結構進行優(yōu)化，最終降低曲軸質(zhì)量12.7%。ALIAKBARI 等［5-6］、LI 等［7］對曲軸疲勞強度進行分析，發(fā)現(xiàn)曲軸的疲勞強度不足的原因是圓角處表面容易產(chǎn)生疲勞裂紋。因此，對曲軸圓角表面進行硬化處理能大大減少疲勞裂紋的產(chǎn)生。楊秀琴［8］為解決某型號曲軸在輥鍛制坯階段模具磨損嚴重導致成形時出現(xiàn)折疊的現(xiàn)象，選取過渡圓角、入模圓角、過渡斜度3個關鍵參數(shù)，以降低磨損量為目標，基于正交試驗進行模具結構優(yōu)化研究，有效減少了曲軸成形時出現(xiàn)折疊的現(xiàn)象。叢建臣等［9］通過對內(nèi)燃機曲軸進行扭轉疲勞試驗研究，發(fā)現(xiàn)連桿頸油孔是扭轉疲勞失效最常見部位，通過軸頸表面感應淬火處理降低曲軸的扭轉疲勞強度約30%，通過拋磨油孔內(nèi)壁提高曲軸的扭轉疲勞強度25%以上。

對于重載曲軸而言，其整體應力大小通常遠小于材料的強度極限，其失效往往是因為曲柄間過渡圓角處的疲勞破壞。因此，降低過渡圓角處的最大應力值可以有效提高重載曲軸的壽命。CHEN等［10］、WANG 等［11］發(fā)現(xiàn)圓角半徑對曲軸的強 度影響較大，隨著圓角半徑的減小，其應力值會顯著增大。因此，本文作者以重泵公司某型號隔膜泵曲軸為研究對象，以ANSYS 為研究工具，研究曲軸圓角處的結構參數(shù)對圓角上的最大應力值及曲軸整體變形量的影響。采用Kriging 插值法得到優(yōu)化目標P3（軸頸與曲柄過渡圓角的最大應力值）、P4（曲柄斜面處圓角的最大應力值）、P5（曲軸的最大變形量）與設計變量P1（曲柄臂斜面傾斜角）、P2（斜面處圓角半徑）之間的具體函數(shù)關系，在支持曲軸模型生成的最大尺寸范圍內(nèi)，應用多目標遺傳算法尋找最優(yōu)的結構參數(shù)。

1 參數(shù)化建模

參數(shù)化建模在形貌優(yōu)化設計過程中很常見，它通過改變參數(shù)變量自動更新產(chǎn)品的結構［12］。本文作者基于UG 建立曲軸圓角處的局部參數(shù)化模型，見圖1。如圖2 所示，變量P1、P2為圓角處結構特征參數(shù)，決定了曲柄處的形狀樣貌。在建模過程中，除了安裝工藝孔和油孔外，其他幾何特征均未做簡化，特別是計算所關心的各個圓角處。為方便后續(xù)表述，將曲拐從扭矩輸入端開始依次編號為拐1、拐2、拐3。

圖1 曲軸的三維模型Fig.1 3D model of crankshaft

圖2 曲柄參數(shù)化結構Fig.2 Crank parametric structure

2 靜強度計算

本文作者采用靜強度應力分析方法來確定曲軸的危險工況及應力分布情況。為獲取更加精確的曲軸交變應力水平，綜合考慮計算精度和計算時間因素，將曲軸一個旋轉周期內(nèi)的受力情況劃為72 個載荷工況。當曲柄轉角θ＝0°時，曲軸受力情況為第一載荷工況。當曲柄轉角θ＝5°時，曲軸受力情況為第二載荷工況，以此類推。

2.1 曲軸參數(shù)

文中研究對象為三拐四支撐結構曲軸，各曲拐之間相隔120°，曲軸三維結構模型見圖1。曲軸材料參數(shù)及曲軸-連桿機構部分參數(shù)如表1 所示。

表1 曲軸參數(shù)（部分）Tab.1 Crankshaft parameters（partial）

2.2 載荷分析

隔膜泵工作時，作用在曲軸上的力有：經(jīng)連桿傳遞到曲軸的連桿力、曲軸與支撐軸承間的摩擦力、曲軸自身的旋轉慣性力。通過設置摩擦接觸的方式添加曲軸與支撐軸承間的摩擦力，通過添加初始旋轉速度來考慮旋轉慣性力對曲軸的影響，而連桿力則通過推導計算得到具體數(shù)值。

曲軸-十字頭系統(tǒng)的機構運動簡圖及受力見圖3。θ為曲軸曲拐與x軸夾角，當曲拐轉角θ＝0～180°時，對應工作行程，十字頭所受活塞力F＝1 530 kN；當曲拐轉角θ＝180°～360°時，對應非工作行程，根據(jù)文獻［13］可以將十字頭所受活塞力近似看為F＝0 kN。

圖3 曲軸-十字頭系統(tǒng)的機構運動簡圖Fig.3 Schematic of mechanism movement of crankshaftcrosshead system

在工作行程時，根據(jù)正弦定理可以得到

則作用于曲軸上的連桿力Fi可表示為

式中：θ為曲拐轉角，（°）；F為活塞力，F(xiàn)＝1 530 kN；AB為曲拐到軸旋轉中心距離，AB＝260 mm；AC為連桿長度，AC＝1 700 mm。根據(jù)上述公式，可得到不同曲拐轉角θ下的連桿力Fi的具體數(shù)值。需要注意的是，在靜強度計算過程中，曲軸某一工況位置上的載荷情況應為該工況下3 個曲拐受力情況的疊加狀態(tài)。

2.3 邊界條件

在支撐軸頸處增加軸套來模擬軸承與曲軸之間的摩擦影響，在軸套與曲軸之間添加摩擦因數(shù)為0.12的摩擦接觸，并設置接觸間隙為0.25 mm。此外，如圖4 所示，約束曲軸鍵槽側面的z向旋轉自由度；約束兩端支撐軸套的徑向跳動（即x、y方向上的位移）；約束另外兩個支撐軸套表面的徑向和軸向跳動（即x、y、z方向上的位移）；約束曲軸底端面的軸向跳動（即z方向上的位移）。

圖4 曲軸靜強度計算的邊界條件Fig.4 Boundary conditions for static strength calculation of crankshaft

2.4 網(wǎng)格劃分

為了兼顧強度計算的準確性和計算時間，應在批量計算之前，對計算模型的網(wǎng)格進行無關性驗證。由于曲軸形狀的不規(guī)則，仿真過程中采用四面體網(wǎng)格劃分法對曲軸整體進行體網(wǎng)格劃分，細化各過渡圓角表面上的面網(wǎng)格。

如表2 所示，經(jīng)過多次試驗，當網(wǎng)格尺寸由體網(wǎng)格為20 mm、面網(wǎng)格為5 mm 轉變?yōu)轶w網(wǎng)格為10 mm、面網(wǎng)格為5 mm 時，P3和P4值的變化幾乎忽略不計。因此，確定以曲軸體網(wǎng)格大小20 mm、圓角處面網(wǎng)格5 mm 作為后續(xù)計算時的網(wǎng)格尺寸。此時，模型包括1 925 088 個節(jié)點、1 313 433 個單元，如圖5 所示。

表2 網(wǎng)格無關性驗證試驗數(shù)據(jù)Tab.2 Grid independence verification test data

圖5 曲軸整體網(wǎng)格Fig.5 The grid of crankshaft

2.5 計算結果

在solver 中將載荷步增至72 個，每個載荷步代表一個曲軸工況，將前面計算得到的72 個工況下的連桿力分解到x、y軸方向上，并以bear loading［14］的方式依次施加到對應曲拐上，實現(xiàn)對72 個工況進行批量計算。經(jīng)驗證，此方法與單獨計算各個工況所得結果一致。

圖6 展示了72 個工況下的各圓角處最大應力值。其中，P3的最大值出現(xiàn)在第25 個工況下，為80.79 MPa。此時拐1 曲柄轉角為0°，所受連桿力為1 530 kN；拐2 曲柄轉角為120°，所受連桿力為1 543.6 kN；拐3 曲柄轉角為240°，所受連桿力為0 kN。而P4的最大值出現(xiàn)在第13 個工況下，為66.50 MPa。此時拐1 曲柄轉角為300°，所受連桿力為0 kN；拐2曲柄轉角為60°，所受連桿力為1 543.6 kN；拐3 曲柄轉角為180°，所受連桿力為1 530 kN。

圖6 曲軸在72 個工況下的圓角處最大應力值Fig.6 Maximum stress of crankshaft at fillet under 72 working conditions

很顯然，曲軸的整體受力大小在第13 工況和第25 工況下一樣。但拐1 和拐3 受力情況發(fā)生交換，而拐2 沒有變化。考慮到P3在所有工況下均大于P4，且P4在第13 工況下和第25 工況下的最大值相差不大，僅為0.03 MPa。因此，確定第25 工況作為后續(xù)曲柄形貌優(yōu)化的計算工況。第25 工況下的曲軸圓角處等效應力云圖如圖7 所示。

3 曲柄形貌優(yōu)化

串聯(lián)UG 和Workbench 的數(shù)據(jù)接口，結合Design Exploration 模塊搭建曲柄形貌優(yōu)化的計算環(huán)境。

圖8 所示為整個優(yōu)化過程的流程。首先，建立參數(shù)集，其中包含設計變量和優(yōu)化目標。確定設計變量的約束范圍，在約束范圍內(nèi)均勻地設計試驗方案，并計算得到不同方案下的計算結果。然后，基于計算結果構建各目標關于變量的響應面，并通過多目標遺傳算法在約束條件下尋找目標的最優(yōu)值。最后，將最優(yōu)結果對應的設計變量返回到幾何模型中進行計算，得到仿真結果。如果仿真結果與優(yōu)化結果之間的誤差較小，說明優(yōu)化成功，并選取一個最優(yōu)候選點，對其取整作為優(yōu)化最終結果。否則，重新構建各響應面，并再次尋優(yōu)、驗證。

圖8 曲柄形貌優(yōu)化流程Fig.8 Flow of crank shape optimization

3.1 構建響應面模型

響應面模型是利用合理的試驗設計方法并通過得到對應的試驗數(shù)據(jù)，采用插值或者擬合的方法來構建變量與目標之間的數(shù)學函數(shù)關系。它可以清晰地表現(xiàn)出各變量對目標的影響關系，便于尋找最優(yōu)目標。響應面的構建包括試驗設計和模型構建這兩部分［15］。

3.1.1 試驗設計

通過設計一系列的試驗方案，運用盡量少的資源和計算時間，獲得可反映設計變量與優(yōu)化目標的高度相關信息數(shù)據(jù)。試驗設計中，好的試驗方案，既能降低計算成本又能提高響應面精度［16］。根據(jù)第3.2 節(jié)中的不等式（9）（10），P2的取值范圍會隨著P1的變化而變化，為了更加均勻地設計試驗樣本點，將P1設置為6 水平因素，P2為多水平因素，一共31 組試驗樣本點。

部分試驗方案及結果見表3，圖9 給出了P1處于不同水平下，曲軸最危險工況下各優(yōu)化目標隨P2變化的規(guī)律曲線。

表3 部分試驗方案及結果Tab.3 Partial test scheme and results

圖9 各優(yōu)化目標隨P1、P2變化的規(guī)律曲線Fig.9 Regularity curves of optimization objectives changing with P1 and P2：（a） P3；（b） P4；（c） P5

從圖9 可以看出：在P1不變的情況下，優(yōu)化目標P3隨著P2的增大而增大，而P4和P5則隨P2的增大而減小。其中，P4的下降趨勢遠大于P5。對P4來講，P2即其對應處的圓角半徑是影響其大小的主要因素，這與文獻［10-11］的結論一致。此外，當P1＞55°之后，它對P4的影響將變得很??；而P3、P5受P1的變化影響卻十分明顯。因此，想要通過直接觀察計算結果找出最優(yōu)的設計變量是不可能的，需要通過構建響應面模型來精確表達設計變量與優(yōu)化目標之間的數(shù)學函數(shù)關系，并應用優(yōu)化算法來尋找最優(yōu)解。

3.1.2 基于Kriging 法構建響應面模型

Kriging 法是一種元建模算法，可提供改進的響應質(zhì)量并擬合輸出參數(shù)的更高階變化。它是一種精確的多維插值，可以插值試驗設計樣本點，結合了類似于標準響應曲面的多項式模型（提供設計空間的“全局” 模型）加上局部偏差，其表達式為

其中：y（x）是待擬合的響應函數(shù)；β是基函數(shù)回歸系數(shù)；f（x）是變量x的多項式響應面模型，提供設計空間的全局近似模型；z（x）是期望為0、方差為σ2的高斯隨機函數(shù)。其協(xié)方差矩陣為

式中：R為相關矩陣。選擇高斯相關函數(shù)為

式中：θk是相關函數(shù)參數(shù)，用于衡量兩樣本點之間的相關性隨兩點間的距離增加的衰減度，相關性越小，生成的響應面越光滑；xi，k和xj，k分別是樣本點xi和xj的第k個元素。

基于上述Kriging 響應面模型理論以及前面得到的試驗設計樣本點，建立Kriging 響應面，得到設計變量對優(yōu)化目標的響應面模型，如圖10 所示。其中各響應面模型的均方根誤差分別為2.044 9×10-7、2.326×10-6、2.538 2×10-11，這表示響應面擬合效果良好。

圖10 各優(yōu)化目標的響應面模型Fig.10 Response surface models of each optimization objective：（a） P3；（b） P4；（c） P5

3.2 多目標遺傳算法尋優(yōu)

MOGA（Multi-objective Genetic Algorithm）是當前應用較為廣泛的多目標優(yōu)化算法之一，相比較其他優(yōu)化算法，其算法復雜性低、解集的收斂性好、運行速度快。MOGA 賦予多個目標函數(shù)不同的權重并將它們組合成一個標量適應度函數(shù)。與單目標遺傳算法不同的是，MOGA 中多個目標函數(shù)的權重不是恒定的，而是隨機指定的。因此，MOGA 的搜索方向是不固定的，如圖11 所示。當多目標優(yōu)化問題具有凹Pareto前沿時，恒定權重的單目標遺傳算法往往無法找到整個Pareto 前沿［17］（即所有Pareto 最優(yōu)解）。因此，作者采用多目標遺傳算法尋找最優(yōu)解。

圖11 多目標遺傳算法的尋優(yōu)方向Fig.11 Optimization direction of multi-objective genetic algorithm

3.2.1 建立多目標優(yōu)化模型

進行多目標優(yōu)化需要先確定此次優(yōu)化的設計變量、目標函數(shù)和約束條件，以建立優(yōu)化模型。本文作者以P1、P2作為設計變量，以P3、P4、P5作為優(yōu)化目標，以P1、P2的尺寸約束范圍作為約束條件，以響應面模型表達優(yōu)化目標與設計變量之間的函數(shù)關系。優(yōu)化模型用下面的形式表示：

式中：f（xi）為優(yōu) 化目標；xi為設 計變量；xi、L、xi，U分別為設計變量xi的下限值和上限值；wi為各優(yōu)化目標所占權重，表示對應目標對于整體優(yōu)化的重要性。每個權重表示如下：

式中：rj為非負隨機數(shù)。MOGA 中每次選擇字符串的時候都會重新賦予新的不同權重值，以保證尋優(yōu)方向的隨機性。

如圖12 所示，將曲柄外形輪廓之間的關系單獨提取出來，此時的DE長度是P2的極限值。如果P2值超過極限值，UG 中將不能生成對應三維模型。

圖12 P1、P2 之間的關系Fig.12 Relationship between P1 and P2

于是可以得到：

當P1＞90°后，曲柄明顯向內(nèi)凹陷，會導致結果剛度顯著降低。因此P1的范圍為

以式（9）和式（10）作為優(yōu)化的約束條件，得到具體的優(yōu)化模型為

式中：P3、P4、P5代表對應的響應面模型。

3.2.2 優(yōu)化結果

設置初始樣本數(shù)為2 000，每次迭代的樣本數(shù)為400，允許Pareto 的最大百分比值為70%，收斂穩(wěn)定性設置為2%。得到優(yōu)化后的Pareto 前沿見及其對應的Pareto 解集見圖13。

圖13 優(yōu)化結果Fig.13 Optimization results：（a）Pareto frontier；（b）solution set of Pareto

根據(jù)文獻［18］可得曲軸的剛度條件為

式中：［y］為軸的允許撓度，mm；L為軸的跨度，L＝3 574 mm。

由式（13）可得曲軸允許的最大變形量為0.714 8 mm。對比Pareto 前沿可知，所有的非支配解都滿足曲軸的最大變形量要求。從圖13（b）可以發(fā)現(xiàn)，Pareto 解集中P1都接近于90°。針對這一發(fā)現(xiàn)，作者做出的合理解釋是：雖然隨著曲柄斜面傾斜角P1的增大，曲軸的最大變形量也將增大，但在滿足曲軸剛度條件的情況下，曲軸的曲柄斜面傾斜角接近90°可能是一種能夠最大程度地降低曲柄斜面處圓角最大應力值和軸頸與曲柄過渡圓角處的最大應力值的極限狀態(tài)。這將在后續(xù)研究中進一步驗證。

選擇P3＝54.94 MPa、P4＝54.55 MPa、P5＝0.535 mm 作為最終優(yōu)化結果，此時，P1＝89.20°、P2＝73.71 mm。將結果圓整為P1＝90°、P2＝75 mm返回到幾何模型中，得到仿真結果為P3＝54.89 MPa、P4＝54.31 MPa、P5＝0.535 mm，與優(yōu)化結果幾乎一致。最終結果與優(yōu)化前相比，P3下降了32.06%，P4下降了18.29%，P5上升了9.41%，但仍滿足剛度條件要求。仿真結果見圖14，曲軸曲柄優(yōu)化前、后對比見圖15。

圖14 仿真結果Fig.14 Simulation results：（a） P3；（b） P4；（c） P5

圖15 曲軸曲柄優(yōu)化前后對比Fig.15 Comparison of crankshaft crank before and after optimization：（a）before optimization；（b）after optimization

4 結論

本文作者從圓角處結構參數(shù)優(yōu)化的角度出發(fā)，對曲軸進行形貌優(yōu)化。首先，分析曲軸的受力情況，將曲軸運轉過程分為72 個計算工況。采用靜強度分析方法確定曲軸圓角處于最大應力值時的計算工況，并以此工況作為后續(xù)形貌優(yōu)化的計算工況。隨后，根據(jù)曲軸三維模型的生成條件確定合理的試驗樣本點，并得到各樣本點的計算結果。利用響應面模型定量地表示設計變量P1、P2與優(yōu)化目標P3、P4、P5之間的函數(shù)關系。最后，構建多目標優(yōu)化模型，使用多目標遺傳算法尋找P1、P2的非支配解集?？梢缘玫揭韵陆Y論：

（1）曲軸的危險工況為第25 個工況，此時曲軸圓角處的最大應力值為80.79 MPa。

（2）在P1不變的情況下，P2的增大會導致P3的增大，P4、P5的減小。其中，P4的減小趨勢顯然大于P5的減小趨勢。

（3）對于P4來講，其對應的圓角半徑P1是影響其大小的主要因素，并且當曲柄斜面傾斜角大于55°之后，其對于P4的影響將逐漸變小；而P3、P5隨P1的變化影響始終表現(xiàn)得較為明顯。

（4）得到的Pareto 解集中P1都接近于90°。這可能是在滿足曲軸剛度條件的情況下，P1接近90°能夠最大程度地降低P3和P4的極限狀態(tài)。

（5）確定最優(yōu)結果P1＝90°、P2＝75 mm，此時P3＝54.89 MPa、下降了32.06%，P4＝54.31 MPa、下降了18.29%，P5＝0.535 mm、上升了9.41%（但仍滿足要求）。