范新亮，王彤，夏遵平

南京航空航天大學 機械結(jié)構(gòu)力學及控制國家重點實驗室，南京 210016

有限元模型修正首先在航空航天領域提出，發(fā)展至今日已形成了一個龐大的理論體系，并且廣泛應用于運載火箭、衛(wèi)星、航天飛機、飛機、直升機等結(jié)構(gòu)的響應與載荷預示、顫振分析、振動控制[1]。模型修正可基于不同的特征進行，主流的有基于頻域模態(tài)參數(shù)、頻響函數(shù)(Frequency Response Function,FRF)的方法以及基于時域動響應的方法。其中時域方法雖抗噪性相比前兩者稍弱，但獲取實驗數(shù)據(jù)簡便，且對于包含非線性的結(jié)構(gòu)[2]、運行狀態(tài)的結(jié)構(gòu)[3]等難以利用模態(tài)或頻響函數(shù)數(shù)據(jù)進行模型修正的情形仍適用。Modak等[4]通過數(shù)值仿真詳細對比了基于模態(tài)與基于頻響函數(shù)的方法[5-9]的差異及收斂性?；谀B(tài)參數(shù)的方法簡單有效，在多數(shù)情形下能給出較滿意的修正結(jié)果，但也存在一些缺陷：如模態(tài)參數(shù)提取過程引入了誤差及不確定性[10]，待修正參數(shù)數(shù)目受到測試模態(tài)數(shù)量的限制?；陬l響函數(shù)的方法則避免了提取模態(tài)參數(shù)帶來的誤差，且具有大量數(shù)據(jù)可供修正以解決待修正參數(shù)眾多的問題。

普遍的模型修正算法均是對整體結(jié)構(gòu)有限元模型進行修正，但有限元模型不確定性常存在于某些關(guān)鍵的局部區(qū)域，例如飛機各個部件的連接區(qū)域、機械結(jié)構(gòu)的關(guān)節(jié)點、轉(zhuǎn)子系統(tǒng)的支承等，且這類局部區(qū)域所在的子結(jié)構(gòu)(局部結(jié)構(gòu))無法從整體結(jié)構(gòu)中拆卸，或其組裝至整體結(jié)構(gòu)上的動態(tài)特性與獨立狀態(tài)有所差異，此時往往需要利用整體結(jié)構(gòu)的測試信息對其有限元模型進行迭代修正，計算效率較低。對此，Weissenburger[11]提出了一種預測局部結(jié)構(gòu)的參數(shù)變化對整體結(jié)構(gòu)振動特性影響的方法。Zhu等[12]使用主模態(tài)縮聚殘余結(jié)構(gòu)，并以局部區(qū)域?qū)淖咏Y(jié)構(gòu)有限元模型中的物理參數(shù)和殘余結(jié)構(gòu)的模態(tài)參數(shù)作為優(yōu)化變量，極小化結(jié)構(gòu)動力學特征方程的殘差來修正整體結(jié)構(gòu)的有限元模型。翁順等[13]提出了一種基于子結(jié)構(gòu)的有限元模型修正方法，當結(jié)構(gòu)僅局部參數(shù)變化時只需分析局部區(qū)域所對應子結(jié)構(gòu)即可求解整體結(jié)構(gòu)的模態(tài)特征量靈敏度矩陣。上述研究都有效提高了大型結(jié)構(gòu)有限元模型修正的效率，但均為基于模態(tài)參數(shù)的方法。在基于頻響函數(shù)的方法中，Guo等[14]推導了一種考慮殘余結(jié)構(gòu)約束的局部結(jié)構(gòu)應變頻響函數(shù)的近似計算方法，并根據(jù)其靈敏度建立了針對局部區(qū)域進行模型修正的方程，該方法與Zhu等[12]所提方法均對殘余結(jié)構(gòu)進行了縮聚而保留了完備的局部結(jié)構(gòu)，然而并未對縮聚誤差對修正結(jié)果的影響進行討論。綜上，若能得到整體結(jié)構(gòu)動響應與獨立狀態(tài)下的局部結(jié)構(gòu)的關(guān)系，則可直接利用整體結(jié)構(gòu)測試數(shù)據(jù)對獨立狀態(tài)下的局部結(jié)構(gòu)模型進行修正，從而顯著提高效率。

本文根據(jù)頻響函數(shù)解耦理論得到的整體結(jié)構(gòu)與獨立狀態(tài)的局部結(jié)構(gòu)兩者間頻響函數(shù)的關(guān)系式，將包含待修正參數(shù)的局部結(jié)構(gòu)動剛度矩陣與殘余結(jié)構(gòu)有限元頻響函數(shù)耦合得到整體結(jié)構(gòu)擬合頻響函數(shù)，通過極小化其與測量值的殘差得到待修正參數(shù)的估計值，從而推導了以局部結(jié)構(gòu)動力學矩陣表示的整體結(jié)構(gòu)頻響函數(shù)殘差關(guān)于待修正參數(shù)的靈敏度方程。因此，當結(jié)構(gòu)建模誤差僅發(fā)生在局部區(qū)域時，利用該方程可對分離出的局部結(jié)構(gòu)單獨進行有限元模型修正，將殘余結(jié)構(gòu)與更新的局部結(jié)構(gòu)重新裝配即得到修正后的整體結(jié)構(gòu)有限元模型。最后，通過仿真和實驗算例驗證了該方法的有效性和抗噪性。

1 理論背景

1.1 局部結(jié)構(gòu)頻響識別方法

為了利用整體結(jié)構(gòu)測試頻響函數(shù)對獨立狀態(tài)下的局部結(jié)構(gòu)有限元模型進行修正，首先需要建立整體結(jié)構(gòu)與局部結(jié)構(gòu)之間頻響函數(shù)的關(guān)系，即頻響函數(shù)解耦問題[15-19]。

圖1 局部結(jié)構(gòu)與殘余結(jié)構(gòu)Fig.1 Local structures and residual structures

整體結(jié)構(gòu)、殘余結(jié)構(gòu)及局部結(jié)構(gòu)的輸入輸出關(guān)系分別為

HAFA=XA,HRFR=XR,HLFL=XL

(1)

式中：HA、HR和HL分別為整體結(jié)構(gòu)、殘余結(jié)構(gòu)和局部結(jié)構(gòu)的頻響函數(shù)矩陣。

根據(jù)界面力協(xié)調(diào)條件及位移協(xié)調(diào)條件有

(2)

設在殘余結(jié)構(gòu)VR內(nèi)部節(jié)點中所選測試自由度對應的位移矢量為XI*=PI*TXI，其中PI*為XI*T各列在XIT中的位置矩陣。定義偽測試自由度為VR內(nèi)部測試自由度與完備界面自由度的并集，則V與VR在偽測試自由度上的位移矢量XA*、XR*為

(3)

(4)

相應載荷矢量為FA*=PA*TFA及FR*=PA*TFR。而VL在完備界面自由度上的位移矢量XL*為

(5)

其載荷矢量為FL*=PLTFL。在上述位移矢量XA*、XR*及XL*上，式(1)成為

(6)

式中：頻響函數(shù)子矩陣HA*、HR*與HL*分別為

并將式(6)按VR內(nèi)部測試自由度與完備界面自由度寫為分塊形式

(7)

(8)

(9)

(10)

其中

(11)

式中：HR*與HC*分別為由殘余結(jié)構(gòu)與局部結(jié)構(gòu)有限元模型計算得到的頻響函數(shù)。

式(10)即為有限元模型整體結(jié)構(gòu)V與殘余結(jié)構(gòu)VR、局部結(jié)構(gòu)VL間頻響函數(shù)的關(guān)系式，是后文推導針對局部結(jié)構(gòu)有限元模型修正公式的基礎。由式(3)注意到偽測試自由度包含完備的界面自由度XJ，因此計入所有界面自由度時解耦式(10) 是精確的，若忽略某些自由度則界面力協(xié)調(diào)條件不完全成立，將引起近似誤差。

1.2 局部結(jié)構(gòu)有限元模型修正方法

假設整體結(jié)構(gòu)有限元模型中局部區(qū)域存在建模誤差，其待修正參數(shù)為θ，而殘余結(jié)構(gòu)不存在或近似不存在誤差。由式(10)知局部結(jié)構(gòu)與殘余結(jié)構(gòu)頻響函數(shù)耦合所得整體結(jié)構(gòu)擬合頻響函數(shù)為

(12)

取由參數(shù)θ確定的殘差函數(shù)為

(13)

式中：ej為激勵自由度j在偽測試自由度中的位置向量；HA*ej為整體結(jié)構(gòu)偽測試自由度對應的測試頻響函數(shù)。

采用極大似然估計來識別參數(shù)θ，使殘差vj取極小值，即參數(shù)θ所決定的整體結(jié)構(gòu)頻響函數(shù)擬合值與測量值誤差最小。極大似然估計需要測量數(shù)據(jù)的“先驗”噪聲信息，當滿足互不相關(guān)的白噪聲假設時，“先驗”信息可由測量數(shù)據(jù)的標準方差來代替[20]。對于具有No個輸出、Ni個輸入的頻響函數(shù)，共包含NoNi個隨機變量。假設所有頻率點處的頻響函數(shù)測量值相互獨立，則此NoNi個隨機變量的聯(lián)合概率密度函數(shù)為

(14)

式中：H為共軛轉(zhuǎn)置符號；ωk為第k個頻率點(為行文簡潔，推導中非必要處略去ωk)；v由Ni個殘差向量vj組合而得；C為頻響函數(shù)噪聲的協(xié)方差矩陣，當各列噪聲相互獨立時其為Cj組成的對角矩陣，且

(15)

根據(jù)極大似然估計原理得到等價極大似然函數(shù)為

(16)

式中：Nf為頻率點數(shù)目。將式(16)簡寫為

L(θ)=Q(θ)HQ(θ)

(17)

(18)

根據(jù)Newton-Gauss方法知參數(shù)θ的迭代估計式為

J(θ(r))JH(θ(r))d(r+1)=-J(θ(r))Q(θ(r))

(19)

式中：d(r+1)為第r+1個迭代步的參數(shù)增量，即

d(r+1)=θ(r+1)-θ(r)

(20)

(21)

式(12)、式(18)代入式(21)左端雅可比矩陣可化簡得

(22)

其中等效頻響函數(shù)定義為

(23)

式中：I為單位矩陣；DL(θ)為局部結(jié)構(gòu)的動剛度矩陣，其參數(shù)化的表達式[12]為

(24)

AL(ω)=

(25)

(26)

式(24)代入式(22)得

(27)

其中

(28)

將式(27)代入式(21)得與ωk及ej對應的迭代方程

(29)

(30)

式(30)即為求解待修正參數(shù)的迭代方程，簡寫為

S(θ(r))d(r+1)=f(θ(r))

(31)

式中:f(θ)為待修正參數(shù)θ所確定的整體結(jié)構(gòu)頻響函數(shù)擬合值與測量值的加權(quán)殘差；S(θ)為以局部結(jié)構(gòu)的動力學矩陣MLe等表示的殘差f(θ)關(guān)于參數(shù)θ的靈敏度。綜上，對整體結(jié)構(gòu)測試得到HA*ej后，再由殘余結(jié)構(gòu)計算得到HR*，即可由式(31)對獨立狀態(tài)下的局部結(jié)構(gòu)有限元模型進行修正。

1.3 算法實現(xiàn)相關(guān)問題

1.3.1 頻率點的篩選準則

實際測試中，反共振區(qū)頻響函數(shù)易受到噪聲污染，而共振區(qū)不僅噪聲小，對參數(shù)變化的靈敏度高，并且其幅值能有效反映阻尼參數(shù)。因此在進行模型修正的頻率范圍內(nèi)選擇共振區(qū)的頻率點，更容易得到穩(wěn)健的參數(shù)識別結(jié)果[21]。所提方法以頻響函數(shù)測試值與擬合值的殘差為最小化目標函數(shù)，因此選取殘差較小的頻率點參與參數(shù)估計可有效減小測試噪聲的擾動影響，故此處以頻響函數(shù)幅值相關(guān)性系數(shù)[14]作為篩選頻率點的準則

(32)

1.3.2 局部結(jié)構(gòu)模型修正流程

局部結(jié)構(gòu)模型修正流程如圖2所示：

1) 將含建模誤差的局部結(jié)構(gòu)VL從整體結(jié)構(gòu)V中分離出來。確定完備界面自由度及其對應的位置矩陣PR,J與PL、整體結(jié)構(gòu)偽測試自由度及其對應位置矩陣PA*以及激勵自由度及其對應位置向量ej。需要注意的是偽測試自由度并非真實進行測試的自由度，因為完備的界面自由度信息幾乎不可能由試驗得到。并且偽測試自由度為真實測試自由度與完備界面自由度的并集。

圖2 局部結(jié)構(gòu)模型修正流程Fig.2 Process of model updating of local structures

6) 重復步驟3)～5)，直至滿足停止準則。迭代結(jié)束后將局部結(jié)構(gòu)的修正參數(shù)代回整體結(jié)構(gòu)，即得到修正后的整體有限元模型。

注意到當殘余結(jié)構(gòu)有限元模型規(guī)模較大時，HR*的計算量較大，因此HR*僅計算一次且在迭代中不重復該計算。由此可見，所提方法將無需進行模型修正的殘余結(jié)構(gòu)的計算工作置于迭代之外，效率相比傳統(tǒng)方法得到了提高。進一步地，可對殘余結(jié)構(gòu)有限元模型進行縮聚后計算HR*，但這將引入額外誤差。此外，若直接從整體結(jié)構(gòu)測試數(shù)據(jù)中提取局部結(jié)構(gòu)信息[23]，即

(33)

(34)

而后利用所得HL*對獨立狀態(tài)下的局部結(jié)構(gòu)有限元模型進行修正，將存在以下缺陷：需得到完整的測試頻響函數(shù)HA*，在實際中工程量巨大甚至難以實現(xiàn)；且提取HL*的過程需求解大量逆矩陣，使噪聲的擾動放大而導致抗噪性較差[24]。

2 算例分析

2.1 仿真算例

采用某噴氣式飛機模型驗證本文方法的有效性及抗噪性。模型整體結(jié)構(gòu)V及局部結(jié)構(gòu)VL、殘余結(jié)構(gòu)VR如圖3所示，采用Shell63板單元建模。其頻帶內(nèi)主要彈性模態(tài)為機翼一階、二階、三階彎曲及尾翼的彎曲。設置建模誤差存在于機翼所在局部區(qū)域，而機身與尾翼所在殘余結(jié)構(gòu)VR的模型參數(shù)無誤差。其中VL又劃分為蒙皮、襟翼、內(nèi)骨架等5個組(圖3)，每個組的彈性模量、密度及阻尼系數(shù)作為待修正參數(shù)θ。

以某組參數(shù)η*對應的整體結(jié)構(gòu)有限元模型作為真實結(jié)構(gòu)，在所選激勵自由度及偽測試自由度上計算其頻響函數(shù)后添加10%噪聲(白噪聲與有色噪聲各5%)得到整體結(jié)構(gòu)測試頻響函數(shù)HA*ej；VR參數(shù)與η*中相應分量ζ*(η*=[ζ*θ*])取相同，計算得偽測試自由度上殘余結(jié)構(gòu)頻響函數(shù)HR*；VL的參數(shù)θ相對η*中相應分量θ*設置一定的偏差作為待修正參數(shù)初始值，計算得初始局部結(jié)構(gòu)的頻響函數(shù)HL(θ(0))，并與HR*重構(gòu)得整體結(jié)構(gòu)頻響函數(shù)初始擬合值，與測量值對比如圖4所示。

圖3 噴氣式飛機局部結(jié)構(gòu)與殘余結(jié)構(gòu)Fig.3 Local and residual structures of a jet

圖4 修正前整體結(jié)構(gòu)頻響函數(shù)對比Fig.4 Comparison of global structure FRF before updating

圖5 修正后整體結(jié)構(gòu)頻響函數(shù)對比Fig.5 Comparison of global structure FRF after updating

圖6 修正前后頻響函數(shù)幅值相關(guān)性對比Fig.6 FRF magnitude correlation coefficient plot before and after updating

圖7為θ(r)相對真實參數(shù)比值的迭代歷程以及修正前后對比，E1至E5、ρ1至ρ5分別為5個分組的彈性模量及密度。由圖知修正后大多參數(shù)收斂于真實值附近，僅ρ3及ρ5識別精度較差，推測其原因為目標殘差函數(shù)對其靈敏度較小，使得該參數(shù)易受噪聲干擾。也正因此，在所分析頻段內(nèi)，該參數(shù)對整體結(jié)構(gòu)的動態(tài)特性作用較小，其誤差不會影響修正后的整體模型預測動力學響應的準確性。通過該仿真算例表明所提方法能在噪聲干擾下有效地對含建模誤差的局部區(qū)域進行修正，使得整體結(jié)構(gòu)模型的動力學特性與真實情形一致。

圖7 迭代歷程及修正前后參數(shù)值對比Fig.7 Iterative process and comparison of parameters before and after updating

2.2 實 驗

圖8 三角機翼飛機實驗系統(tǒng)Fig.8 Experimental system for delta-winged aircraft

為進一步驗證所提方法對實際結(jié)構(gòu)局部區(qū)域進行模型修正的有效性，對三角機翼飛機進行了實驗，如圖8所示。厚3 mm的機翼與尾翼通過螺栓與厚6 mm的機身連接，采用Shell63板單元建模并劃分為8個組(圖9)，且螺栓連接簡化為固接。由于該結(jié)構(gòu)在低頻范圍內(nèi)主要為機翼的振動，且機翼與機身連接處存在不確定的建模誤差，因此將機翼及其連接部分分離出來作為待修正的

圖9 三角機翼飛機有限元模型Fig.9 FE model of delta-winged aircraft

局部結(jié)構(gòu)VL，剩余部分則為近似無誤差的殘余結(jié)構(gòu)VR。實驗中測試頻率范圍設置為100 Hz，取800條譜線。整體結(jié)構(gòu)采用彈性繩懸掛方式模擬自由—自由邊界條件，傳感器位于兩側(cè)機翼上，在均勻分布的41個測試點上進行敲擊，獲得其頻響函數(shù)HA*ej。為抑制噪聲[23]，須合理選擇包含較多結(jié)構(gòu)信息的測試自由度集合。

表1 有限元模型參數(shù)Table 1 Design parameters of FEM

根據(jù)所提方法利用HA*ej和HR*對獨立狀態(tài)下的VL進行修正后，將更新的VL重新與VR裝配即得到修正后的整體結(jié)構(gòu)有限元模型，其計算的頻響函數(shù)與測量值吻合較好(圖11)，且兩者幅值相關(guān)性相比修正前有了明顯的提高(圖12)。由此證明了所提方法能有效地對局部區(qū)域存在誤差的實際結(jié)構(gòu)進行模型修正。

圖10 修正前三角機翼飛機整體結(jié)構(gòu)頻響函數(shù)對比Fig.10 Comparison of global structure FRF before updating for delta wing aircraft

圖11 修正后三角機翼飛機整體結(jié)構(gòu)頻響函數(shù)對比Fig.11 Comparison of global structure FRF after updating for delta wing aircraft

檢驗修正前后的模態(tài)特征量變化，如表2所示，修正前后平均頻率誤差由8.34%降低至0.39%， 振型相關(guān)性(MAC)平均值為0.95(圖13)。

圖14為整體結(jié)構(gòu)頻響函數(shù)測試值與修正后的有限元模型計算值在全頻域的形狀[14]相關(guān)性，其相當于MAC在頻域上的推廣，定義為

(35)

圖12 三角機翼飛機修正前后頻響函數(shù)幅值相關(guān)性對比Fig.12 FRF magnitude correlation coefficient plot before and after updating for delta wing aircraft

表2 修正前后各階頻率對比Table 2 Frequency comparison before and after updating

圖13 修正后振型相關(guān)性Fig.13 Modal assurance criterion after updating

式中：ωTk、ωFi分別為測試頻響函數(shù)與有限元模型計算頻響函數(shù)的頻率點。顯然，有限元模型與實際結(jié)構(gòu)越接近，αs的對角元素越接近1。由圖14知修正后的整體結(jié)構(gòu)有限元模型的準確度相比修正前明顯改善。

圖14 修正前后頻響函數(shù)形狀相關(guān)性對比Fig.14 FRF shape correlation coefficient plot before and after updating

3 結(jié) 論

本文提出了一種針對局部結(jié)構(gòu)的有限元模型修正方法。通過數(shù)值算例與實驗算例表明：

1) 對于一個復雜結(jié)構(gòu)，當結(jié)構(gòu)建模誤差可以確定為僅發(fā)生在某些關(guān)鍵的局部區(qū)域時，利用整體結(jié)構(gòu)測試頻響函數(shù)及近似無誤差的殘余結(jié)構(gòu)頻響函數(shù)可直接對獨立狀態(tài)下的局部結(jié)構(gòu)進行有限元模型修正，將更新的局部結(jié)構(gòu)模型與殘余結(jié)構(gòu)重新裝配即得修正后的整體結(jié)構(gòu)有限元模型。

2) 在測試噪聲較大、頻響函數(shù)初始擬合值與測試值殘差較大、待修正參數(shù)較多等復雜情形下，所提方法仍具有較好的收斂性。

3) 對三角機翼飛機的機翼所在局部區(qū)域進行修正后，所計算的整體有限元模型動態(tài)特性與實測數(shù)據(jù)吻合較好，表明所提方法能適用于實際結(jié)構(gòu)的模型修正，并提高效率。