張偉杰，陸秋海，緱百勇，王 波

(1.清華大學(xué) 航天航空學(xué)院 應(yīng)用力學(xué)教育部重點(diǎn)實(shí)驗(yàn)室，北京 100084；2.長(zhǎng)安福特汽車(chē)有限公司 技術(shù)開(kāi)發(fā)中心，重慶 401122)

隨著有限元計(jì)算技術(shù)的進(jìn)步，有限元建模與分析已日益成為解決工程實(shí)際問(wèn)題重要的手段，采用有限元建模分析代替費(fèi)時(shí)、費(fèi)力和昂貴的結(jié)構(gòu)試驗(yàn)分析方法，也日益成為研究者的共識(shí)。但對(duì)于實(shí)際工程結(jié)構(gòu)，尤其是復(fù)雜的多部件結(jié)構(gòu)，有限元分析模型的正確性往往難以得到保證，分析精度常不能滿(mǎn)足需求。因此，如何獲得結(jié)構(gòu)的準(zhǔn)確模型，提高模型分析和預(yù)示的精度，是結(jié)構(gòu)仿真分析成敗的關(guān)鍵，通過(guò)實(shí)驗(yàn)數(shù)據(jù)進(jìn)行有限元模型的確認(rèn)及其修正，尤其重要。

有限元模型修正目前有兩個(gè)主流方法：矩陣修正法和設(shè)計(jì)參數(shù)型法[1]。在各種有限元模型修正方法中，響應(yīng)面方法[2]的修正效率較高。現(xiàn)有文獻(xiàn)中響應(yīng)面方法較多應(yīng)用于可靠度分析問(wèn)題[3―6]，而近年來(lái)也有不少文獻(xiàn)研究基于響應(yīng)面方法的有限元模型修正[7―10]。

基于響應(yīng)面方法的有限元模型修正的基本思想是：通過(guò)擬合響應(yīng)面函數(shù)近似得到設(shè)計(jì)參數(shù)與特征參數(shù)的顯式表達(dá)式，從而快速求解設(shè)計(jì)參數(shù)的修正量。其主要步驟有：確定響應(yīng)面函數(shù)；確定試驗(yàn)設(shè)計(jì)方案；調(diào)用有限元程序；擬合響應(yīng)面；求解設(shè)計(jì)參數(shù)修正量；檢驗(yàn)修正誤差；迭代生成新響應(yīng)面。

當(dāng)響應(yīng)面回歸系數(shù)確定后，若需要求解不同特征參數(shù)對(duì)應(yīng)的設(shè)計(jì)參數(shù)修正量，則必須每次計(jì)算都通過(guò)迭代算法來(lái)求解，在實(shí)際工程應(yīng)用仍不太方便。本文受到逆系統(tǒng)方法[11]的啟發(fā)，提出了逆響應(yīng)面方法，直接建立設(shè)計(jì)參數(shù)關(guān)于特征參數(shù)的顯式關(guān)系。獲得逆響應(yīng)面模型后，可以很方便地直接計(jì)算不同特征參數(shù)對(duì)應(yīng)的設(shè)計(jì)參數(shù)值，無(wú)需迭代計(jì)算。

1 逆響應(yīng)面方法

1.1 基本思想

逆響應(yīng)面法的基本思想是對(duì)調(diào)響應(yīng)面函數(shù)中的輸入變量與輸出變量關(guān)系，即，以若干輸出變量為自變量，以若干輸入變量為因變量，通過(guò)回歸分析方法來(lái)擬合輸入變量關(guān)于輸出變量的顯式表達(dá)式。

利用逆響應(yīng)面法進(jìn)行有限元模型修正可直接得到設(shè)計(jì)參數(shù)（如密度、彈性模量等）關(guān)于特征參數(shù)（如固有頻率、振型等）的顯式表達(dá)式，那么當(dāng)逆響應(yīng)面函數(shù)的回歸系數(shù)確定后，可以根據(jù)特征參數(shù)的目標(biāo)值直接得到設(shè)計(jì)參數(shù)的修正值，而不需要進(jìn)行迭代計(jì)算，有效地提高計(jì)算速度和精度。

1.2 逆響應(yīng)面函數(shù)

逆響應(yīng)面函數(shù)可選用二次（或更高階次的）多項(xiàng)式函數(shù)，即

其中αk、αki、αkji為逆響應(yīng)面回歸系數(shù)，xi（i=1，2，…，Nx）為輸入變量（特征參數(shù)），yk（k=1，2,…，Ny）為輸出變量（待修正的設(shè)計(jì)參數(shù)），Nx是輸入變量的數(shù)量，Ny是輸出變量的數(shù)量。這個(gè)模型和通常的響應(yīng)面模型形式上完全相同，只是逆響應(yīng)面模型是以特征參數(shù)作為輸入?yún)?shù)，設(shè)計(jì)參數(shù)作為輸出參數(shù)，正好與響應(yīng)面法相反，逆響應(yīng)面法也由此得名。

1.3 試驗(yàn)設(shè)計(jì)

由于輸入變量（特征參數(shù)）不方便直接給定，所以用于逆響應(yīng)面法的試驗(yàn)設(shè)計(jì)和響應(yīng)面法相同，即由輸出變量（待修正的設(shè)計(jì)參數(shù)）的試驗(yàn)設(shè)計(jì)確定。另外，因?yàn)檩斎搿⑤敵鰯?shù)量未必相等，則試驗(yàn)設(shè)計(jì)的試驗(yàn)點(diǎn)數(shù)可能會(huì)少于或遠(yuǎn)多于逆響應(yīng)面函數(shù)所需試驗(yàn)點(diǎn)數(shù)（即逆響應(yīng)面回歸系數(shù)的個(gè)數(shù)），那么就會(huì)導(dǎo)致無(wú)法求解或明顯增加計(jì)算成本的問(wèn)題，所以計(jì)算時(shí)需要根據(jù)輸入、輸出數(shù)量的關(guān)系適當(dāng)調(diào)整試驗(yàn)點(diǎn)數(shù)。

本文選用中心復(fù)合設(shè)計(jì)（CCD）[12]，此試驗(yàn)設(shè)計(jì)由三部分組成：（1）中心點(diǎn)，即試驗(yàn)原點(diǎn)，共1個(gè)試驗(yàn)點(diǎn)；（2）軸點(diǎn)，即yi= ±a（i=1，2，…，Ny）且yj=0（j≠i），共2Ny個(gè)試驗(yàn)點(diǎn)；（3）析因設(shè)計(jì)點(diǎn)，即由Ny因素二水平（±b）正交表構(gòu)造而成的點(diǎn)，共有2Ny個(gè)試驗(yàn)點(diǎn)，即總試驗(yàn)點(diǎn)數(shù)為N=1+2Ny+2Ny。而對(duì)于二次逆響應(yīng)面函數(shù)，回歸系數(shù)的個(gè)數(shù)為M=1+2Nx+Nx(Nx-1)/2。

若N＜M，可考慮增加水平為±a/2的軸點(diǎn)或水平為±b/2的析因設(shè)計(jì)點(diǎn)；若N遠(yuǎn)大于M（如N＞2M），可考慮去掉軸點(diǎn)或析因設(shè)計(jì)點(diǎn)以便降低有限元計(jì)算成本。

1.4 逆響應(yīng)面精度

為了檢驗(yàn)?zāi)骓憫?yīng)面的擬合精度，確定結(jié)果誤差是否在允許范圍內(nèi)，利用平均相對(duì)誤差值來(lái)表示，即

其中xi*為輸入變量（特征參數(shù)）的目標(biāo)值，x^i為由有限元計(jì)算得到的特征參數(shù)的實(shí)際值。

“集體情感只有固著于某種物質(zhì)對(duì)象,它本身才能被意識(shí)到”。在灌溉分水的過(guò)程中其實(shí)是村民集體意識(shí)的加強(qiáng)與社會(huì)秩序鞏固的一個(gè)過(guò)程，通過(guò)灌溉過(guò)程中村民自覺(jué)分配用水時(shí)間和用水量，慢慢成了一種約定俗成的規(guī)矩，也是一種民間社會(huì)秩序，達(dá)成了某種共識(shí)和契約，形成了村民們之間的凝聚力。在此過(guò)程中，體現(xiàn)了村民們的集體意識(shí)和協(xié)作觀念。孫澄在水文化的固守與變遷中寫(xiě)道，通過(guò)用水分配可以“調(diào)節(jié)人與人之間的關(guān)系,調(diào)整人與社會(huì)的關(guān)系,協(xié)調(diào)漢族與周邊彝族、傣族的民族關(guān)系,維護(hù)和構(gòu)建了區(qū)域內(nèi)“和諧”取用水的公共關(guān)系和社會(huì)秩序”。芒沙村的水資源利用與管理模式就是這種文化的真實(shí)寫(xiě)照。

本文算例中，由于已知設(shè)計(jì)參數(shù)的目標(biāo)值，故可通過(guò)計(jì)算修正值的平均誤差來(lái)檢驗(yàn)?zāi)骓憫?yīng)面的擬合精度。

1.5 求解步驟

（1）選擇合適階次的多項(xiàng)式作為逆響應(yīng)面函數(shù)；

（2）根據(jù)N、M的值確定中心復(fù)合設(shè)計(jì)的試驗(yàn)點(diǎn)；

（3）調(diào)用有限元程序計(jì)算試驗(yàn)點(diǎn)的輸出值；（4）利用最小二乘法求解逆響應(yīng)面回歸系數(shù)；（5）把輸入變量的目標(biāo)值代入逆響應(yīng)面函數(shù)計(jì)算修正值；

（6）調(diào)用有限元程序計(jì)算修正值的誤差；

（7）若誤差在允許范圍內(nèi)，則求解結(jié)束；若誤差較大，可適當(dāng)提高逆響應(yīng)面函數(shù)的階次或增加試驗(yàn)點(diǎn)，并重復(fù)步驟（3）—（6），直到誤差在允許范圍內(nèi)。

2 算例驗(yàn)證

2.1 算例信息

本文利用一個(gè)簡(jiǎn)支梁（如圖1）作為算例驗(yàn)證逆響應(yīng)面方法，梁長(zhǎng)為6 m，橫截面形狀為矩形，尺寸為0.2 m×0.2 m，密度為2 500 kg/m3，彈性模量為30 GPa，泊松比為0.3。建立相應(yīng)的由15個(gè)梁?jiǎn)卧M成的有限元模型，并以其中的單元3、單元8和單元9的彈性模量作為待修正設(shè)計(jì)參數(shù)，以前4階固有頻率作為特征參數(shù)。

2.2 試驗(yàn)數(shù)據(jù)

給定設(shè)計(jì)參數(shù)的可變化量為50%，即試驗(yàn)點(diǎn)中彈性模量最小值為15 GPa，最大值為45 GPa，利用中心復(fù)合設(shè)計(jì)建立試驗(yàn)點(diǎn)并進(jìn)行有限元分析，得到如表1的試驗(yàn)數(shù)據(jù)，其中，E3、E8和E9分別表示單元3、單元8和單元9的彈性模量，f1、f2、f3和f4表示前四階的固有頻率。

為了對(duì)比不同特征參數(shù)數(shù)量的結(jié)果，逆響應(yīng)面的擬合將分成四輸入三輸出、三輸入三輸出和二輸入三輸出三種情況，分別以前4、3和2階固有頻率作為特征參數(shù)（輸入）。對(duì)于四輸入三輸出情況，N=M=15，為了保證計(jì)算精度，增加6個(gè)軸點(diǎn)作為試驗(yàn)點(diǎn)，即表1中序號(hào)16—21點(diǎn)，即，所用到的試驗(yàn)數(shù)據(jù)為表1中1—21；對(duì)于三輸入三輸出情況，試驗(yàn)點(diǎn)數(shù)合適，則所用到的試驗(yàn)數(shù)據(jù)為點(diǎn)1—15；對(duì)于二輸入三輸出情況，N=15，M=6，可考慮僅用部分試驗(yàn)點(diǎn)計(jì)算以縮短有限元計(jì)算時(shí)間，如用點(diǎn)1及點(diǎn)8—15。

圖1 簡(jiǎn)支梁及其橫截面示意圖Fig.1 Simply supported beam and its cross section

表1 試驗(yàn)數(shù)據(jù)Tab.1 Experimental data

2.3 四輸入三輸出結(jié)果

由表1中1—21數(shù)據(jù)，采用最小二乘擬合得到三個(gè)修正參數(shù)的逆響應(yīng)面函數(shù)分別為式（3）（4）（5），從而得到E3、E8和E9的修正值分別為 21.719 GPa、42.527 GPa和26.162 GPa，其與響應(yīng)面法計(jì)算（所用的試驗(yàn)數(shù)據(jù)一致）結(jié)果對(duì)比如表2。

2.4 三輸入三輸出結(jié)果

由表1中點(diǎn)1—15數(shù)據(jù)，通過(guò)最小二乘擬合得到三個(gè)修正參數(shù)的逆響應(yīng)面函數(shù)（限于篇幅，表達(dá)式從略），得到E3、E8和E9的修正值分別為19.669 GPa、42.295 GPa和26.072 GPa，其與響應(yīng)面法計(jì)算（所用的試驗(yàn)數(shù)據(jù)一致）結(jié)果對(duì)比如表3。

2.5 二輸入三輸出結(jié)果

由表1中點(diǎn)1及點(diǎn)8—15數(shù)據(jù)，通過(guò)最小二乘擬合得到三個(gè)修正參數(shù)的逆響應(yīng)面函數(shù)（表達(dá)式從略），從而得到E3、E8和E9的修正值分別為20.618 GPa、36.362 GPa和35.109 GPa，其與響應(yīng)面法計(jì)算（對(duì)二輸入來(lái)說(shuō)，9個(gè)試驗(yàn)點(diǎn)不足以求解回歸系數(shù)，故使用點(diǎn)1—15作為試驗(yàn)數(shù)據(jù)）結(jié)果對(duì)比如表4。

2.6 計(jì)算結(jié)果分析

對(duì)于四輸入三輸出情況，響應(yīng)面法與逆響應(yīng)面法所用數(shù)據(jù)量一樣時(shí)，逆響應(yīng)面法的精度較高。這表明對(duì)于特征參數(shù)多于待修正設(shè)計(jì)參數(shù)的情況，相比響應(yīng)面法，運(yùn)用逆響應(yīng)面法可有效提高計(jì)算精度，但同時(shí)也需要增加相應(yīng)的試驗(yàn)數(shù)據(jù)，增大計(jì)算成本。

表2 計(jì)算結(jié)果對(duì)比（四輸入三輸出）Tab.2 Results comparison(4 inputs and 3 outputs)

表3 計(jì)算結(jié)果對(duì)比（三輸入三輸出）Tab.3 Results comparison(3 inputs and 3 outputs)

表4 計(jì)算結(jié)果對(duì)比（二輸入三輸出）Tab.4 Results comparison(2 inputs and 3 outputs)

對(duì)于三輸入三輸出情況，由于輸入變量與輸出變量的數(shù)量一致，計(jì)算精度差別不大。

對(duì)于二輸入三輸出情況，響應(yīng)面和逆響應(yīng)面的擬合精度都較差，但此時(shí)逆響應(yīng)面所需的試驗(yàn)數(shù)據(jù)較少，有限元計(jì)算成本較低。這說(shuō)明了，對(duì)于特征參數(shù)少于待修正設(shè)計(jì)參數(shù)的情況，相比響應(yīng)面法，運(yùn)用逆響應(yīng)面法可降低計(jì)算成本，對(duì)于規(guī)模較大的有限元模型來(lái)說(shuō)，運(yùn)用逆響應(yīng)面法可顯著縮短計(jì)算時(shí)間。為了進(jìn)一步提高此情況下的擬合精度，作者嘗試過(guò)利用郭勤濤等[5]的方法增加試驗(yàn)數(shù)據(jù)并提高多項(xiàng)式階次，或利用Ren等[9]的方法縮小設(shè)計(jì)參數(shù)的可變范圍并生成新的響應(yīng)面，但無(wú)論是響應(yīng)面法還是逆響應(yīng)面法，修正精度均無(wú)法大幅提高。因此，這可能說(shuō)明特征參數(shù)數(shù)量需大于或等于待修正設(shè)計(jì)參數(shù)數(shù)量才能對(duì)有限元模型進(jìn)行有效的修正，Deng等[10]也提到同樣的觀點(diǎn)。而在工程實(shí)際中，待修正設(shè)計(jì)參數(shù)通常是未知的，很有可能碰到特征參數(shù)少于待修正設(shè)計(jì)參數(shù)的問(wèn)題，解決這個(gè)問(wèn)題的有效方法仍需要作進(jìn)一步探討。

3 工程算例驗(yàn)證

3.1 算例信息

某實(shí)際工程結(jié)構(gòu)幾何模型（經(jīng)簡(jiǎn)化）如圖2所示。本文將零件1和零件2的彈性模量E p1和E p2作為待修改設(shè)計(jì)參數(shù)，前兩階固有頻率f p1和f p2作為特征參數(shù)。特征參數(shù)的目標(biāo)值（試驗(yàn)實(shí)測(cè)值）為12.862 Hz和210.372 Hz，零件1和零件2的彈性模量初始值分別為200 GPa和68 GPa，給定設(shè)計(jì)參數(shù)的可變化量為50%，利用中心復(fù)合設(shè)計(jì)建立試驗(yàn)點(diǎn)并進(jìn)行有限元分析，得到如表5的試驗(yàn)數(shù)據(jù)。

3.2 利用逆響應(yīng)面法修正

圖2 某實(shí)際工程結(jié)構(gòu)幾何模型Fig.2 Engineering structure model

利用表5的試驗(yàn)數(shù)據(jù)擬合相應(yīng)的逆響應(yīng)面函數(shù)，將特征參數(shù)的目標(biāo)值（12.862 Hz，210.372 Hz）代入逆響應(yīng)面中，從而得到設(shè)計(jì)參數(shù)的修正值分別為172.5 GPa和85.9 GPa。

將此修正值代入有限元程序中進(jìn)行驗(yàn)證，得到前兩階固有頻率為12.821 Hz和210.759 Hz，與目標(biāo)值誤差小于1%。

因此，基于逆響應(yīng)面法的有限元模型修正能有效地應(yīng)用于實(shí)際工程項(xiàng)目中。

4 結(jié)語(yǔ)

本文提出了有限元模型修正的逆響應(yīng)面法，建立了2階多項(xiàng)式逆響應(yīng)面模型，探討了試驗(yàn)設(shè)計(jì)問(wèn)題，并對(duì)本文方法進(jìn)行了仿真驗(yàn)證。算例表明，逆響應(yīng)面方法可用于有限元模型修正，修正精度與響應(yīng)面法相當(dāng)或略?xún)?yōu)。工程算例證明，基于逆響應(yīng)面法的有限元模型修正能有效地應(yīng)用于實(shí)際工程項(xiàng)目中。

獲得響應(yīng)面模型后，如果改變特征參數(shù)的目標(biāo)值，則響應(yīng)面法必須通過(guò)迭代重新求解修正值；而逆響應(yīng)面法得到的是設(shè)計(jì)參數(shù)關(guān)于特征參數(shù)的顯式表達(dá)式，可以直接得到相應(yīng)的設(shè)計(jì)參數(shù)值，無(wú)需復(fù)雜的迭代計(jì)算，這在實(shí)際工程應(yīng)用上有很大的優(yōu)勢(shì)。

表5 試驗(yàn)數(shù)據(jù)Tab.5 Experimental data

為了保證逆響應(yīng)面法的穩(wěn)定性，回歸系數(shù)計(jì)算過(guò)程可能需要較響應(yīng)面法多的試驗(yàn)數(shù)據(jù)，這在一定程度上增加了有限元計(jì)算的開(kāi)銷(xiāo)。

[1]Mottershead J E,Friswell M I.Model updating in structural dynamics:A survey[J].Journal of Sound and Vibration,,1993,167(2):347-375.

[2]Myers R C.Response surface methodology[M].Boston:Allyn and Bacon,Inc.1971.

[3]Bucher C G,Bourgund U.A fast and efficient response surface approach for structural reliability problems[J].Structural Safety,,1990,7(1):57-66.

[4]Kim S H,Na S W.Response surface method using vector projected sampling points[J].Structural Safety,,1997,19:3-19.

[5]郭勤濤，張令彌，費(fèi)慶國(guó).用于確定性計(jì)算仿真的響應(yīng)面法及其試驗(yàn)設(shè)計(jì)研究[J].航空學(xué)報(bào)，2005，27(1)：55-61.

[6]Nguyen X S,Sellier A,Duprat F,Pons G.Adaptive response surface method based on a double weighted regression technique[J]. Probabilistic Engineering Mechanniiccss,,2009,24(2):135-143.

[7]鄧苗毅，任偉新.基于響應(yīng)面方法的結(jié)構(gòu)有限元模型修正研究進(jìn)展[J].鐵道科學(xué)與工程學(xué)報(bào)，2008，5(3)：42-45.

[8]Marwala T.Finite element model updating using response surface method[A].Paper AIAA-2004-2005[C].California:45 th AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics and Materials Conference.2004,7:5165-5173.

[9]Ren W X,Fang S E,Deng M Y.Response surface-based finite-element-model updating using structural static responses[J].Journal of Engineering Mechanics-ASCE 2011,137(4):248-257.

[10]Deng L,Cai C S.Bridge model updating using response surface method and genetic algorithm[J].Journaall ooff Bridge Engineering,2010,15(5):553-564.

[11]Bendat J S.Nonlinear system techniques and applications[M].New York:John Wiley ＆ Sons.Inc.1998.

[12]Montgomery D C.Design and analysis of experiments[M].New York:John Wiley ＆ Sons.Inc..1991.