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

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

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

0 引言

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

1 理論基礎

1.1 模型修正理論

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

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

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

ωi是第i個殘差函數的權重系數。

1.2 Kriging模型

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

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

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

2 滑鞍模態試驗

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

表1 滑鞍模態試驗實測值

3 構建Kriging代理模型

3.1 有限元模型建立

滑鞍整體尺寸約為320 mm×220 mm×110 mm，根據滑鞍幾何尺寸在SoildWorks中建立起三維模型，將三維模型導入到Hypermesh 中劃分結構化網格，網格模型如圖1所示。

圖1 滑鞍有限元網格模型

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

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

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

3.2 參數選擇與靈敏度分析

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

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

圖2 設計參數靈敏度分析

3.3 自適應Kriging模型構建

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

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

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

4 模型修正

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

表3 修正參數在模型修正前后的對比

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

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

5 結語

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