一種針對局部結構的有限元模型修正方法

范新亮，王彤，夏遵平

南京航空航天大學 機械結構力學及控制國家重點實驗室，南京 210016

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

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

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

1 理論背景

1.1 局部結構頻響識別方法

為了利用整體結構測試頻響函數對獨立狀態下的局部結構有限元模型進行修正，首先需要建立整體結構與局部結構之間頻響函數的關系，即頻響函數解耦問題[15-19]。

圖1 局部結構與殘余結構Fig.1 Local structures and residual structures

整體結構、殘余結構及局部結構的輸入輸出關系分別為

HAFA=XA,HRFR=XR,HLFL=XL

(1)

式中：HA、HR和HL分別為整體結構、殘余結構和局部結構的頻響函數矩陣。

根據界面力協調條件及位移協調條件有

(2)

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

(3)

(4)

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

(5)

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

(6)

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

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

(7)

(8)

(9)

(10)

其中

(11)

式中：HR*與HC*分別為由殘余結構與局部結構有限元模型計算得到的頻響函數。

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

1.2 局部結構有限元模型修正方法

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

(12)

取由參數θ確定的殘差函數為

(13)

式中：ej為激勵自由度j在偽測試自由度中的位置向量；HA*ej為整體結構偽測試自由度對應的測試頻響函數。

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

(14)

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

(15)

根據極大似然估計原理得到等價極大似然函數為

(16)

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

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

(17)

(18)

根據Newton-Gauss方法知參數θ的迭代估計式為

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

(19)

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

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

(20)

(21)

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

(22)

其中等效頻響函數定義為

(23)

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

(24)

AL(ω)=

(25)

(26)

式(24)代入式(22)得

(27)

其中

(28)

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

(29)

(30)

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

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

(31)

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

1.3 算法實現相關問題

1.3.1 頻率點的篩選準則

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

(32)

1.3.2 局部結構模型修正流程

局部結構模型修正流程如圖2所示：

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

圖2 局部結構模型修正流程Fig.2 Process of model updating of local structures

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

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

(33)

(34)

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

2 算例分析

2.1 仿真算例

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

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

圖3 噴氣式飛機局部結構與殘余結構Fig.3 Local and residual structures of a jet

圖4 修正前整體結構頻響函數對比Fig.4 Comparison of global structure FRF before updating

圖5 修正后整體結構頻響函數對比Fig.5 Comparison of global structure FRF after updating

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

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

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

2.2 實 驗

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

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

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

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

表1 有限元模型參數Table 1 Design parameters of FEM

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

圖10 修正前三角機翼飛機整體結構頻響函數對比Fig.10 Comparison of global structure FRF before updating for delta wing aircraft

圖11 修正后三角機翼飛機整體結構頻響函數對比Fig.11 Comparison of global structure FRF after updating for delta wing aircraft

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

圖14為整體結構頻響函數測試值與修正后的有限元模型計算值在全頻域的形狀[14]相關性，其相當于MAC在頻域上的推廣，定義為

(35)

圖12 三角機翼飛機修正前后頻響函數幅值相關性對比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 修正后振型相關性Fig.13 Modal assurance criterion after updating

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

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

3 結 論

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

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

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

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