鮑 諾 王春潔

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

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

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

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

1 響應(yīng)面優(yōu)化的模型修正原理

1.1 4階多項式響應(yīng)面模型

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

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

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

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

由待定系數(shù)β可得到響應(yīng)面函數(shù)的表達式.

1.2 響應(yīng)面有效性檢驗

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

1.3 F值檢驗分析

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

假設(shè)對有限元模型的修正參數(shù)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，則變量對響應(yīng)的影響顯著;

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

2 響應(yīng)面優(yōu)化的模型修正流程

2.1 目標(biāo)函數(shù)

結(jié)構(gòu)有限元模型修正通過將試驗值和響應(yīng)面值相結(jié)合，構(gòu)造出優(yōu)化目標(biāo)函數(shù)，將模型修正轉(zhuǎn)換為優(yōu)化問題.

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

2.2 模型修正流程

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

圖1 基于響應(yīng)面優(yōu)化的模型修正流程

3 三自由度系統(tǒng)數(shù)值算例

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

圖2 三自由度彈簧－質(zhì)量系統(tǒng)

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

根據(jù)模型修正流程構(gòu)建三自由度數(shù)值算例3種工況的響應(yīng)面模型.3種工況的響應(yīng)面模型檢驗結(jié)果，如圖3a和圖3b所示，其均方根誤差相對值和決定系數(shù)分別趨向0和1，可以準(zhǔn)確描述三自由度數(shù)值算例的結(jié)構(gòu)參數(shù)與響應(yīng)的關(guān)系.

圖3 三自由度模型3種工況響應(yīng)面有效性評價

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

圖4 三自由度模型3種工況目標(biāo)函數(shù)收斂曲線

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

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

表1 三自由度模型修正前后參數(shù)值

4 GARTEUR模型修正與預(yù)測

4.1 初始模型構(gòu)建

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

圖5 GARTEUR結(jié)構(gòu)模型

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

圖6 GARTEUR有限元模型

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

4.2 參數(shù)篩選及響應(yīng)面有效性評價

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

圖7 12個初選參數(shù)的顯著性P值

設(shè)定P=0.05為顯著性水平臨界值，從初始模型的x1～x12的12個參數(shù)中篩選出7個對各階模態(tài)頻率顯著性的參數(shù)，如表2所示.將篩選出的7個顯著參數(shù)構(gòu)建響應(yīng)面模型，響應(yīng)面檢驗結(jié)果如圖8所示，其均方根誤差相對值和決定系數(shù)分別趨向0和1，準(zhǔn)確地描述了GARTEUR結(jié)構(gòu)參數(shù)與響應(yīng)的關(guān)系.

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

圖8 GARTEUR模型的響應(yīng)面有效性評價

4.3 響應(yīng)面優(yōu)化的模型修正結(jié)果

驗證修正的模型是否具有一定的普適性，需檢驗修正模型復(fù)現(xiàn)及預(yù)示能力，高質(zhì)量的修正模型應(yīng)達到三級預(yù)示水平[14]:

1)能夠精確復(fù)現(xiàn)參與修正頻段內(nèi)的試驗數(shù)據(jù);

2)能夠準(zhǔn)確預(yù)測參與修正頻段外的試驗數(shù)據(jù);

3)對結(jié)構(gòu)局部修改，模型不重新修正只需相應(yīng)的改動就能預(yù)測結(jié)構(gòu)修改后的動力學(xué)特性.

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

圖9 GARTEUR模型目標(biāo)函數(shù)收斂曲線

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

表3 GARTEUR模型的模態(tài)頻率修正結(jié)果及誤差對比

修正頻段內(nèi)的前6階模態(tài)頻率的平均誤差由修正前10.13%降低到0.652%;預(yù)測頻段內(nèi)的后4階模態(tài)頻率的平均誤差由4.48%降低到1.12%;總平均誤差由7.87%降低到0.839%.可見，修正后有限元模型不但能復(fù)現(xiàn)修正頻段的頻率而且還能準(zhǔn)確預(yù)測修正頻段外的頻率，滿足一級和二級預(yù)測標(biāo)準(zhǔn).

4.4 模型結(jié)構(gòu)修改的預(yù)測性檢驗

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

表4 GARTEUR模型局部修改后的頻率預(yù)測結(jié)果及誤差對比

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

5 結(jié)論

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

2)三自由度系統(tǒng)的3種不同工況模型修正結(jié)果表明，待修正模型中有參數(shù)不存在或者有誤差的情況下，響應(yīng)面優(yōu)化的模型修正也能得到準(zhǔn)確的模態(tài)頻率修正值，具有良好的模型修正能力.

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

References)

[1]丁繼鋒，韓增堯，馬興瑞.航天器動力學(xué)模型試驗驗證技術(shù)研究進展[J].力學(xué)進展，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]張保強，陳國平，郭勤濤.基于模態(tài)頻率和有效質(zhì)量的有限元模型修正[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]王永華，楊欣毅，蘇珉，等.熵判別粒子群優(yōu)化算法在發(fā)動機模型修正中的應(yīng)用[J].航空動力學(xué)報，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，F(xiàn)riswell 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，F(xiàn)riswell 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，F(xiàn)regolent A.Dynamic model updating using virtual antiresonances[J].Shock and Vibration，2004，11(3/4):351－363

[14]Link M，F(xiàn)riswell 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