基于代理模型的加工中心滑鞍有限元模型修正

陸原超，周俊榮，王瑞超，李會軍

（五邑大學 智能制造學部，廣東 江門 529020）

0 引言

伴隨著機械制造加工業(yè)的飛速發(fā)展，對加工中心的加工精度和加工效率等要求不斷上升，對加工中心的設計需求提出了新的挑戰(zhàn)，傳統(tǒng)的經(jīng)驗設計方法已經(jīng)不能滿足新的設計需求，有限元分析越來越多地應用到加工中心的設計過程之中[1]。有限元分析的第一步，就是建立一個準確的有限元模型，但由于加工中心結(jié)構(gòu)的復雜性[2]，模型修正技術應運而生[3]。有限元計算結(jié)構(gòu)輸出與輸入為隱函數(shù)關系[4]，修正過程中頻繁調(diào)用有限元模型，需要大量的計算量，通過代理模型將有限元模型用顯函數(shù)近似可以大幅減少計算量，提高模型修正效率[5]。常用的代理模型有多項式響應面[6]、支持向量回歸[7]和Kriging模型[8]等。

1 理論基礎

1.1 模型修正理論

有限元模型修正起源于20世紀60年代，可以分為矩陣法和參數(shù)法。矩陣法通過修改結(jié)構(gòu)剛度矩陣與質(zhì)量矩陣等減小與真實結(jié)構(gòu)間的誤差，但由于矩陣法缺少物理意義且不具備工程意義，已經(jīng)逐漸被淘汰。參數(shù)法通過修改密度、彈性模量、連接剛度等物理參數(shù)，保證了修正過程的物理意義，逐漸成為研究重點[2]。

通過修正有限元模型的設計參數(shù)，參數(shù)型模型修正方法逐步減小響應計算值與實測值之間的誤差，將模型修正這一參數(shù)辨識問題轉(zhuǎn)變?yōu)閮?yōu)化問題，可將優(yōu)化目標函數(shù)F(x)寫為

其中：x是待修正參數(shù)的向量；ri是第i個響應的殘差函數(shù)；

ωi是第i個殘差函數(shù)的權重系數(shù)。

1.2 Kriging模型

傳統(tǒng)參數(shù)型模型修正方法需要不斷地調(diào)用有限元模型參與計算，由于有限元模型計算的黑箱特性，通過建立一個可以模擬輸入與輸出之間函數(shù)關系的顯式代理模型可以大幅減少計算量，Kriging模型就是一種常用的代理模型[4，6]。Kriging模型可以分為回歸部分與隨機響應部分，寫為：

其中：y(x)是Kriging模型的響應預估值；f(x)是Kriging模型的回歸部分，是關于x的多項式函數(shù)，用來擬合整體趨勢；Z(x)是Kriging模型的隨機響應部分，用來擬合局部變化。Kriging模型不僅能得到預估值，還能得到方差，并利用方差在設計空間中重新取樣，對自身進行更新從而獲得較高的模型預測精度。重新取樣依據(jù)的規(guī)則被稱為加點準則，最大期望改善（Expected Improvement，EI）準則是使用較為廣泛的加點準則，其將EI函數(shù)定為評價指標來表征改善程度I(x)，并將EI函數(shù)最大值的設計點作為新設計點添加到設計點集中更新Kriging模型，如此循環(huán)直至滿足收斂準則。改善程度I(x)可以寫為

其中：Φ是標準正態(tài)累計分布函數(shù)；ψ是標準正態(tài)分布概率密度函數(shù)。通過不斷迭代添加新設計點更新Kriging模型，可以使得模型自適應進化，從而滿足工程需要。

2 滑鞍模態(tài)試驗

滑鞍是加工中心主要移動部件，直接與電主軸相連，提供加工中心Z軸移動。為了獲得滑鞍的模態(tài)參數(shù)，對其進行了力錘敲擊模態(tài)試驗。實驗過程中選用了東方所INV9832-50三向加速度傳感器與奇石樂9722A2000力錘，在滑鞍表面均布測點，共布置了20個測點。測試過程采用移動傳感器法，保持力錘在同一敲擊點不變，移動傳感器遍歷所有測點，數(shù)據(jù)采集依靠基鈦克Impaq頻譜分析儀，采樣頻率為5000 Hz，最終結(jié)果導入到計算機中，運用專業(yè)軟件ME’SCOPE進行模態(tài)分析，得到滑板模態(tài)參數(shù)，如表1所示。

表1 滑鞍模態(tài)試驗實測值

3 構(gòu)建Kriging代理模型

3.1 有限元模型建立

滑鞍整體尺寸約為320 mm×220 mm×110 mm，根據(jù)滑鞍幾何尺寸在SoildWorks中建立起三維模型，將三維模型導入到Hypermesh 中劃分結(jié)構(gòu)化網(wǎng)格，網(wǎng)格模型如圖1所示。

圖1 滑鞍有限元網(wǎng)格模型

滑鞍依靠導軌滑塊和螺母絲桿與滑板相連，用彈簧阻尼單元對滑鞍邊界條件進行模擬，用彈簧的剛度值代替接合面剛度值，其中導軌滑塊可以簡化為法向彈簧與軸向彈簧，絲杠螺母可以簡化為徑向彈簧和軸向彈簧，其接觸剛度值可以通過廠家手冊查得，滑塊的法向剛度k1為8.62×108N/m，切向剛度k2為5.64×108N/m；絲桿的徑向剛度k3為6.32×108N/m，軸向剛度k4為7.32×108N/m。滑鞍材料為灰鑄鐵，材料參數(shù)取自ANSYS材料庫，彈性模量E為110 GPa，密度ρ為7200 kg/m3，泊松比μ為0.28。邊界條件設定模擬模態(tài)試驗時的外部條件，將有限元模型導入到ANSYS Workbench中進行有限元模態(tài)分析，與模態(tài)試驗的結(jié)果進行對比，如表2所示。

表2 滑鞍固有頻率實測值與計算值的對比

從對比結(jié)果可以看出，計算固有頻率和實測固有頻率最大相對誤差為19.37%，發(fā)生在第六階，前六階的平均誤差為11%左右，表明有限元模型雖然能一定程度上模擬真實模型，但仍需要進一步修正以滿足工程需求。

3.2 參數(shù)選擇與靈敏度分析

由于滑鞍材料參數(shù)選用的是ANSYS材料庫中灰鑄鐵的材料參數(shù)，與滑鞍實際材料可能有所偏差，故將滑鞍材料參數(shù)納入修正參數(shù)的選擇當中。導軌滑塊和絲桿螺母的剛度參數(shù)是從廠家手冊上查得，但在有限元模型中將其簡化為彈簧阻尼單元，因此將彈簧剛度也納入修正參數(shù)當中。綜上所述，初步選擇了滑鞍密度參數(shù)ρ、滑鞍彈性模量E、滑鞍泊松比μ、法向彈簧剛度k1、切向彈簧剛度k2、徑向彈簧剛度k3和軸向彈簧剛度k4等7個結(jié)構(gòu)參數(shù)作為設計參數(shù)。

對上述7個待修正參數(shù)進行模態(tài)頻率靈敏度計算，計算結(jié)果如圖2所示。可以看出滑鞍密度ρ、彈性模量E、法向彈簧剛度k1、徑向彈簧剛度k3和軸向彈簧剛度k4對滑鞍模態(tài)頻率的靈敏度較高，因此，有限元模型修正的修正參數(shù)選 擇 ρ、E、k1、k3和k4等5 個設計參數(shù)。結(jié)合現(xiàn)實考慮，修正參數(shù)的取值范圍在其初值的0.8～1.2倍之間。

圖2 設計參數(shù)靈敏度分析

3.3 自適應Kriging模型構(gòu)建

實驗設計是一種基于數(shù)理統(tǒng)計和概率論的科學方法，可以通過較少的試驗次數(shù)得到輸入與輸出間的關系。拉丁超立方實驗設計是蒙特卡羅方法的一種，具有較高的采樣效率，在實際應用中使用較廣。五設計參數(shù)的拉丁超立方實驗設計，取樣方式選為中心復合設計，初始設計點共有27個。將選定的參數(shù)提交到有限元模型中參與計算，得到對應的前六階固有頻率數(shù)值，根據(jù)此輸入與輸出關系可以構(gòu)建出初始的Kriging模型，然后可以根據(jù)EI加點準則增加細化點，以0.01為收斂值，最終新增了4個細化點，構(gòu)建成立自適應Kriging模型。以法向彈簧剛度k1和徑向彈簧剛度k3與第一階固有頻率f1的關系為例，標準和自適應Kriging響應面模型如圖3所示。

從圖3可以看出，自適應Kriging模型在函數(shù)的極大值點附近插入了新設計點，從單點極大值變成了兩點極大值。從擬合優(yōu)度來看，標準Kriging模型的擬合優(yōu)度為96.52%，自適應Kriging 模型的擬合優(yōu)度為99.24%，相較于標準Kriging 模型擬合精度得到了提高，是對有限元模型更高精度的近似。

圖3 標準Kriging和自適應Kriging一階固有頻率響應面

4 模型修正

基于自適應Kriging模型，采用遺傳算法搜尋設計空間內(nèi)目標函數(shù)的最小值，遺傳算法初代種群大小設定為1000，子代種群大小設定為500，最大迭代代數(shù)為20，允許的最大帕累托百分比設定為70%，穩(wěn)定收斂閾設定為2%，修正前后各參數(shù)與結(jié)果的對比列在表3。從中可以看出，變動量較大的為滑鞍的彈性模量、滑塊法向接觸剛度與絲桿軸向接觸剛度，但均未接近設計空間邊緣，有較高的置信度。

表3 修正參數(shù)在模型修正前后的對比

表4給出了修正之后的有限元模型前六階模態(tài)計算頻率與實測頻率的對比結(jié)果。從中可以看出，與未修正的有限元模型相比，修正后的模型與實際模型的相似度更高，最大相對誤差從原先的19.37%降低為4.52%，最大相對誤差在5%以內(nèi)，可以滿足工程計算需求。

表4 模型修正后計算頻率與實測頻率的對比

5 結(jié)語

基于自適應Kriging模型，本文對有限元模型進行了參數(shù)修正，較好地平衡了計算效率與計算精度，最終修正后的有限元模型與實測模型相對誤差降至5%以下，可以滿足工程計算需求。