基于響應面優化的結構有限元模型修正

鮑 諾 王春潔

(北京航空航天大學 機械工程及自動化學院，北京100191)

在航空航天等領域存在大量的復雜結構，由于分析過程中的眾多不確定因素以及引入多種假設的建模等效簡化，如材料、邊界條件、結構連接等誤差和局部或整體的非線性等，造成有限元計算值與試驗值之間存在差異[1－2].因此，有必要利用試驗數據對模型進行修正，以便提高有限元模型的計算精度.

模型修正作為動力學的逆問題，可轉化為對設計參數的優化.張保強等以彈性模量、剛度等變量為修正參數對復雜邊界彈性梁有限元模型進行迭代修正[3].王永華等為匹配真實模型應用熵判別粒子群法對發動機性能參數進行有限元模型迭代修正[4].Coppotelli等基于靈敏度分析的模型修正方法對織女星運載火箭的上部結構系統(UCMEC，Upper Composite Mechanical Configuration)進行迭代修正，修正后模態頻率與試驗值誤差明顯減小[5].Davoodi等利用遺傳算法對典型球關節系統的參數進行有限元模型迭代修正，結果顯示模型具有良好的預測精度[6].Mthembu等以梁單元彈性模量為變量結合粒子群法研究了H梁有限元模型迭代修正中最優模型的選擇[7].可見模型參數修正總是需要面對參數篩選和有限元模型反復迭代的情況.

本文針對模型修正過程中模型迭代困難及模型參數選擇且可能存在多工況的問題，提出響應面優化的方法進行模型修正.采用統計分析中F檢驗篩選修正參數來構造響應面，用響應面來替代有限元模型結合遺傳算法進行優化，避免迭代過程每次調用有限元程序.以三自由度系統數值模型來驗證待修正模型中參數多工況時，響應面優化的模型修正能力.以GARTEUR模型來檢驗參數篩選后響應面優化的修正結果，通過三級預示水平來驗證修正后模型具有一定的普遍適用性，進而證實參數篩選后響應面優化的模型修正方法的有效性.

1 響應面優化的模型修正原理

1.1 4階多項式響應面模型

響應面法將試驗設計和數理統計相結合，通過樣本點擬合出響應面函數來模擬輸入(參數)－輸出(響應)的隱式關系.常用的試驗設計有中心復合設計、正交設計、全因子設計、拉丁方抽樣等.由于樣本點在修正參數區間的分布對響應面的近似精度影響較大，因此本文在拉丁方試驗的基礎上增加一個準則[8]，求得此準則下最優的拉丁方設計.其試驗設計矩陣每列中各個水平出現的次序和各個樣本點的因子水平分布更均勻，抽樣效果好.

采用非線性擬合精確度更高的不完全4階多項式作為響應面近似函數，將試驗設計的多組采樣點及響應值運用回歸、擬合等方法構建響應面模型，其響應面近似函數形式為

式中，xi為修正設計參數;K為修正設計參數的個數;β0，βi，βij，βii，βiii，βiiii為不完全 4 階多項式待定系數.

由試驗設計確定的樣本點組成修正設計參數矩陣X，同時可得有限元模型計算值向量Y;假設有限元計算值與響應面擬合值誤差組成的向量為ε，則各矩陣向量之間的關系如式(2)所示.根據最小二乘法，由式(4)和式(5)可得式(6)，求解出不完全4階響應面函數的系數向量β[9].

由待定系數β可得到響應面函數的表達式.

1.2 響應面有效性檢驗

判斷擬合的響應面模型是否可信，通常用均方根誤差(RMSE)相對值和決定系數(R2)兩個標準對樣本數據的擬合精度檢驗，分別如式(7)和(8)所示.RMSE→0表示響應面誤差小;R2→1表示響應面與原模型相似程度高.

1.3 F值檢驗分析

由結構的特點初步選取待修正參數，然后分析參數對各階頻率影響的顯著程度，篩選出對頻率影響顯著的參數，忽略影響小的參數.本文采用F值檢驗法對修正參數篩選.利用試驗設計在參數的設計空間內確定樣本點并進行有限元計算獲得樣本響應頻率，然后對樣本進行F值檢驗，分析多個因素(如剛度、密度、彈性模量等)對模態頻率的影響是否顯著.與傳統的靈敏度參數篩選方法相比，F值檢驗法是從全局角度出發，在整個設計空間中篩選對特征量有顯著影響的參數，避免了靈敏度分析只計算參數在某設計點處局部梯度的問題.

假設對有限元模型的修正參數A進行F檢驗，則 A 的 F 值為[10]

式中SA和Se分別為因素A和誤差e的偏差平方和;fA和fe分別為因素A和誤差e的自由度.給定顯著水平 θ=0.05，使得

對于取定的樣本值，由式(9)計算得到FA，則

1) 若 FA≥F1－θ(fA，fe)，即 P≤0.05，則變量對響應的影響顯著;

2) 若 FA＜ F1－θ(fA，fe)，即 P ＞0.05，則變量對響應的影響不顯著.

2 響應面優化的模型修正流程

2.1 目標函數

結構有限元模型修正通過將試驗值和響應面值相結合，構造出優化目標函數，將模型修正轉換為優化問題.

式中，F(x)為實測頻率與響應面計算頻率相對誤差的平方和;為第j階實測模態頻率;fj為第j階響應面計算頻率;n為參與修正頻率的階數;xi為第i個修正參數，分別為修正參數上下限;K為修正參數的個數.

2.2 模型修正流程

由于模型的多個樣本點計算，將初始有限元模型選定的修正變量參數化，建立參數化模型.實現響應面優化的模型修正流程模塊化，見圖1.主要修正流程為:首先初始模型進行參數化，其次篩選顯著性修正參數構建響應面模型，再次對響應面參數進行優化，最后將參數優化結果代入初始有限元模型中，對修正結果進行驗證.

圖1 基于響應面優化的模型修正流程

3 三自由度系統數值算例

三自由度彈簧－質量系統[11]如圖 2所示.mi=1.0 kg(i=1，2，3)，ki=1.0 N/m(i=3，4)，其中mi為i號質量塊質量，ki為相應的彈簧剛度.模型修正中所有的誤差參數已知，假設3個待修正參數變量的初始值為:k1=k2=k5=2.0 N/m，然而3個變量真實試驗值為:k1=k2=k5=1.0N/m，與真實值相比，誤差為100%.

圖2 三自由度彈簧－質量系統

模型建模過程中，通常會有未考慮到的參數存在問題，例如真實系統中某個參數不存在或者存在誤差，在此工況下進行模型修正，驗證響應面優化的模型修正能力.假定工況1為標準工況，連接值k6=3.0 N/m(準確值)，除了初始設定的k1，k2，k53個待修正參數，不存在其他誤差參數;工況2為模型中連接值k6不存在，即k6=0 N/m，其他與工況1相同;工況3為模型中存在連接誤差，即k6=3.5 N/m，其他與工況1相同.

根據模型修正流程構建三自由度數值算例3種工況的響應面模型.3種工況的響應面模型檢驗結果，如圖3a和圖3b所示，其均方根誤差相對值和決定系數分別趨向0和1，可以準確描述三自由度數值算例的結構參數與響應的關系.

圖3 三自由度模型3種工況響應面有效性評價

三自由度數值模型的3階試驗頻率和響應面計算頻率構造的目標函數，如式(11)所示.采用遺傳算法優化，其迭代收斂曲線如圖4所示.

圖4 三自由度模型3種工況目標函數收斂曲線

3種工況修正前后結果如表1所示，工況1輸入參數k1，k2，k5(N/m)的誤差由100%降低到2%之內，代入有限元模型得到的輸出頻率f1，f2，f3(Hz)的誤差都在1%之內.

工況2和工況3由于待修正模型中k6=0 N/m(不存在)和k6=3.5 N/m(有誤差)的原因，修正后輸入參數k1，k2，k5必然出現較大偏差，但是模型修正后的輸出3階頻率f1，f2，f3誤差全部降低到2%之內.因此待修正模型中即便有參數被剔除或者存在誤差，響應面優化法的模型修正也能獲得準確的模態頻率修正值，具有良好的模型修正能力.

表1 三自由度模型修正前后參數值

4 GARTEUR模型修正與預測

4.1 初始模型構建

GARTEUR飛機模型由法國國家航空航天研究院設計制造，該模型具有低剛度、高柔度、模態頻率低且密集的特點，可以用來評估動力學試驗與模型修正技術[12－13].飛機模型由矩形截面的鋁制梁構成，其中機身長1.5 m，翼展2.0 m，機翼上表面附著一層黏彈性阻尼材料.主要包括機身、機翼、水平尾翼、垂直尾翼、兩邊小翼等部件，如圖5所示.

圖5 GARTEUR結構模型

由于GARTEUR模型各部件細長的特點，主體結構采用梁單元模擬，其有限元模型如圖6所示.機翼連接處為兩個對稱剛性單元，垂直尾翼與機身連接為剛性單元，機身與機翼連接為BUSH單元.在各部件的螺栓等連接處相應地增加集中質量單元來模擬連接件質量.初始模型中各部件的幾何和材料參數相對準確，但結構連接處模型簡化較多、誤差較大，因此需要對模型進行修正.

圖6 GARTEUR有限元模型

根據建模過程中連接結構的特點，初步選出x1～x12共12個待修正參數，分別為:機身/機翼連接處扭轉剛度、機身彎曲剛度z向和y向、垂直尾翼彎曲剛度x向、垂直尾翼扭轉剛度、機翼彎曲剛度x向和z向、機翼水平剛性單元長度、兩邊小翼彎曲剛度、垂直尾翼剛性單元長度、機翼對機身垂直偏置距離、機翼扭轉剛度.

4.2 參數篩選及響應面有效性評價

對GARTEUR模型的初選參數進行F值檢驗，獲得12個參數對10階模態頻率的顯著性P值，如圖7所示.

圖7 12個初選參數的顯著性P值

設定P=0.05為顯著性水平臨界值，從初始模型的x1～x12的12個參數中篩選出7個對各階模態頻率顯著性的參數，如表2所示.將篩選出的7個顯著參數構建響應面模型，響應面檢驗結果如圖8所示，其均方根誤差相對值和決定系數分別趨向0和1，準確地描述了GARTEUR結構參數與響應的關系.

表2 篩選出的顯著性參數及修正前后變化

圖8 GARTEUR模型的響應面有效性評價

4.3 響應面優化的模型修正結果

驗證修正的模型是否具有一定的普適性，需檢驗修正模型復現及預示能力，高質量的修正模型應達到三級預示水平[14]:

1)能夠精確復現參與修正頻段內的試驗數據;

2)能夠準確預測參與修正頻段外的試驗數據;

3)對結構局部修改，模型不重新修正只需相應的改動就能預測結構修改后的動力學特性.

以文獻[14]試驗測試數據為參考，選取前6階模態頻率作為修正頻段并構造目標函數(一級檢驗)，將模型修正轉化為優化算法尋優的問題，如式(11)所示.后4階模態頻率作為預測頻段，用于修正模型的二級預測檢驗.對構建的響應面模型采用遺傳算法進行優化，獲得參數最優解，代入有限元模型得到修正頻率.GARTEUR模型的目標函數迭代收斂曲線，如圖9所示.

圖9 GARTEUR模型目標函數收斂曲線

模型修正前后的有限元分析結果和實測10階模態頻率值，如表3所示.

表3 GARTEUR模型的模態頻率修正結果及誤差對比

修正頻段內的前6階模態頻率的平均誤差由修正前10.13%降低到0.652%;預測頻段內的后4階模態頻率的平均誤差由4.48%降低到1.12%;總平均誤差由7.87%降低到0.839%.可見，修正后有限元模型不但能復現修正頻段的頻率而且還能準確預測修正頻段外的頻率，滿足一級和二級預測標準.

4.4 模型結構修改的預測性檢驗

將GARTEUR模型局部修改，如圖6所示.機翼尖上的配重由0.15 kg改為0.72 kg，考察修正模型對局部結構修改后模態特性預測能力(三級預測)[14]，分析預測結果與實測模態頻率對比，如表4所示.

表4 GARTEUR模型局部修改后的頻率預測結果及誤差對比

相對于初始有限元模型，修正后模型的預測頻率平均誤差從8.415%降低到1.492%，預測結果的最大誤差從初始模型的14.474%降低到2.077%，最小誤差從1.716%降低到1.382%，表明修正模型的三級預測有效性.響應面優化后的修正模型不僅是具有試驗數據復現能力的等效模型，還能基本反映結構的真實動態特性，模型具有一定普遍適用性.

5 結論

1)以優化拉丁方試驗設計樣本點，F值檢驗篩選修正參數來構造不完全4階多項式響應面，最小二乘法確定多項式系數的方法建立的響應面能夠代替有限元模型進行模型修正，避免修正過程中有限元模型的多次迭代，提高了分析效率.

2)三自由度系統的3種不同工況模型修正結果表明，待修正模型中有參數不存在或者有誤差的情況下，響應面優化的模型修正也能得到準確的模態頻率修正值，具有良好的模型修正能力.

3)通過F值檢驗篩選出GARTEUR模型對模態頻率顯著的參數來構造響應面進行修正.結果表明修正模型滿足三級預示水平，即修正頻段和預測頻段具有良好的復現和預測能力，還能預測結構局部修改后的頻率.修正模型具有一定的普遍適用性，驗證了經過參數篩選后響應面優化的模型修正有效性.

References)

[1]丁繼鋒，韓增堯，馬興瑞.航天器動力學模型試驗驗證技術研究進展[J].力學進展，2012，42(4):395－405

Ding Jifeng，Han Zengyao，Ma Xingrui.Research evolution on the test verification of spacecraft dynamic model[J].Advances in Mechanics，2012，42(4):395 － 405(in Chinese)

[2]Xiong Y，Chen W，Tsui K L，et al.A better understanding of model updating strategies invalidating engineering models[J].Computer Methods in Applied Mechanics and Engineering，2009，198(15/16):1327 －1337

[3]張保強，陳國平，郭勤濤.基于模態頻率和有效質量的有限元模型修正[J].振動與沖擊，2012，31(24):69 －73

Zhang Baoqiang，Chen Guoping，Guo Qintao.Finite element model updating based on modal frequency and effective modal mass[J].Journal of Vibration and Shock，2012，31(24):69 －73(in Chinese)

[4]王永華，楊欣毅，蘇珉，等.熵判別粒子群優化算法在發動機模型修正中的應用[J].航空動力學報，2013，28(1):74 －81

Wang Yonghua，Yang Xinyi，Su Min，et al.Engine model correction based on entropy criterion PSO[J].Journal of Aerospace Power，2013，28(1):74 －81(in Chinese)

[5]Coppotelli G，Rinaldi R，Crema L B.Structural updating of the VEGA－UCMEC FE model using vibration test data[R].AIAA－2008－1850，2008

[6]Davoodi M R，Amiri J V，Gholampour S，et al.Determination of nonlinear behavior of a ball joint system by model updating[J].Journal of Constructional Steel Research，2012，71:52 －62

[7]Mthembu L，Marwala T，Friswell M I，et al.Finite element model selection using particle swarm optimization[C]//Conference Proceedings of the Society for Experimental Mechanics Series.New York:Springer，2011，4:41 －52

[8]Jin R C，Chen W，Sudjianto A.An efficient algorithm for constructing optimal design of computer experiments[J].Journal of Statistical Planning and Inference，2005，134(1):268 －287

[9]Khuri A I，Mukhopadhyay S.Response surface methodology[J].WIREs Computational Statistics，2010，2(2):128 － 149

[10]Rice J A.Mathematical statistics and data analysis[M].3rd ed.Florence:Brooks Cole，2007

[11]Mares C，Mottershead J E，Friswell M I.Stochastic model updating:part1theory and simulated example[J].Mechanical Systems and Signal Processing，2006，20(7):1674 －1695

[12]Kozak M T，?ztürk M，?zgüven H N.A method in model updating using miscorrelation index sensitivity[J].Mechanical Systems and Signal Processing，2009，23(6):1747 －1758

[13]D’Ambrogio W，Fregolent A.Dynamic model updating using virtual antiresonances[J].Shock and Vibration，2004，11(3/4):351－363

[14]Link M，Friswell M I.Working group 1:generation of validated structural dynamic models－results of a benchmark study utilizing the CARTEUR SM－AG19 test－bed[J].Mechanical Systems and Signal Processing，2003，17(1):9 －20