基于振動試驗數據的有限元模型修正技術研究

王澤宇，劉 闖，馮咬齊

（北京衛星環境工程研究所，北京 100094）

0 引言

利用試驗數據對航天器結構有限元模型進行修正，可以為航天器的動態仿真試驗奠定基礎。這對于提高模型的動力學預示精度，制定準確的試驗條件有著重要的意義[1]。目前，模型修正的方法主要有兩類：模態法和頻響函數法[2]。模態法利用模態試驗數據提取出來的固有頻率和振型對數學模型中的參數進行修正，這種方法出現得較早，也比較成熟[3-6]。頻響函數是結構的固有屬性，也是系統頻域內的重要特征量，利用頻響數據對模型進行修正的方法可以回避模態辨識帶來的困難和誤差，其在模態密集和高阻尼的結構中仍然可以獲得。從20世紀90年代開始，基于頻響函數的有限元模型修正方法取得了較大的發展[7-9]。

然而無論是利用模態法還是頻響函數法對航天器模型進行修正，都需要進行額外的力學試驗。如果能利用航天器的振動試驗數據完成模型修正工作，必然能節約成本，縮短研制周期。鑒于此，本文提出了利用基于基礎激勵傳遞特性的模型修正方法[10]來修正航天器結構的有限元模型。

1 基于基礎激勵傳遞特性的模型修正方法

航天器振動試驗是將試驗件由夾具緊固在振動臺上，控制振動臺給試驗件提供一個加速度邊界條件，測量在此條件下一些關鍵部位的響應。其主要目的是為了暴露航天器元器件、結構在工藝和材料方面的缺陷，考核產品耐受力學環境的能力。

航天器結構在振動臺上的運動方程為

假設結構做無阻尼振動，上述方程可簡化為

這樣方程(2)可寫成

方程(3)又可以寫成

根據節點是否與外界連接，剛度矩陣和質量矩陣也可以相應地寫成

式中下標1代表對應于與基礎連接的節點的自由度，下標2代表與外界沒有連接的節點的自由度。將式(5)代入式(4)中，可以得到如下方程：

假定臺面各處的加速度相同，方程(6)兩邊同時除以0u，那么

其中： ｛e｝ 為單位矩陣， ｛e｝=｛1， 1， 1，…｝T。這里定義一個新的傳遞函數｛H(ω)｝ = ｛u｝/u0，這個傳遞函數反映了結構傳導基礎激勵響應的特性，是結構本身的固有屬性，可以用此特性作為數學模型與實際物理模型是否一致的標準并用來修正模型，它具體如下式所示：

對于理論模型和試驗模型，由式(7)可以分別得到：

式中下標A代表理論分析模型，X代表試驗模型。理論模型與試驗模型之間存在著偏差，設：

同時

將式(11)與式(12)代入式(10)，并減去式(9)，得到

假設

將式(14)代入式(13)，就可以得方程：

將式(16)代入式(15)，并定義S(ω)、｛u｝和｛Δα(ω)｝分別表示如下關系：

最終得到如下的修正方程：

可以根據不同的ω建立多個方程，將式(20)變成一個超靜定的方程組，這個方程的解就是設計參數的變化量。

在實際測量中，不可能把有限元模型所有節點的響應都測量出來，測量數據是不完備的。因此在實際應用的過程中，必須首先對理論模型進行縮聚[11]。結構的動剛度矩陣與頻響函數矩陣之間的關系： × =ZHI，即：

根據實際測量的自由度t與未測量的自由度n對模型進行分塊。此外，與基礎連接的那些節點的自由度也要保留。通過方程(21)可以提取出如下兩個方程：

并可以得到

縮聚后的動剛度矩陣與頻響函數矩陣應該還是互逆的關系，所以有

由此，修正方程相應的變為

此時的｛e｝base只包括模型中與基礎運動狀態相同的自由度，｛HX｝test與｛ΔH｝test也只對應于進行了測量的節點自由度。

2 模型修正方法的程序實現

模型修正的算法雖然出現得比較早，而且比較容易編寫計算機程序，然而實際應用中的問題是如何將修正的程序與現有的有限元分析軟件連接起來，使得一些商用軟件在作為有限元求解器的同時，還能夠輸出模型修正所需數據并將修正的結果賦予到有限元模型中[12]。只有解決了這一問題，模型修正技術才能真正應用在構型復雜、規模大、材料豐富的結構中。基于這樣的考慮，本文在商用有限元分析軟件NASTRAN的基礎上進行了二次開發，使其可以在進行有限元求解的同時，輸出模型修正所需要的信息，進行模型修正工作。

2.1 模型修正程序的設計思路

模型修正程序結構見圖 1。整個修正程序分兩部分：一部分是利用NASTRAN接口語言編程，使之輸出進行修正工作所需要的信息；另一部分為編寫 Matlab修正程序。接口程序 1 可以從NASTRAN的輸出文件中讀取出動剛度矩陣和材料參數等信息，并將此信息輸送給修正程序。修正程序將處理這些信息并得到動剛度矩陣對材料參數的靈敏度矩陣，完成修正參數和頻率點的選取等工作，組建修正方程并求解，然后將修正過的參數通過接口程序2傳遞給NASTRAN。如果計算結果不滿足控制條件，再調用NASTRAN進行計算，進而組建新的修正方程，得到新的參數。整個修正過程有兩個控制條件：迭代次數和參數的變化范圍。當滿足其中一個條件時，修正程序運行結束。

圖1 模型修正程序框圖Fig.1 The program block diagram of model modification

2.2 NASTRAN接口語言的編程實現方法

NASTRAN從讀取模型文件信息到形成用于計算的質量矩陣，剛度矩陣的流程如圖 2所示，虛線箭頭代表中間省略了部分流程。整個模型的初始化都是在PHASE0.DAT這個功能模塊下完成的。在執行該模塊的過程中，數據從模型文件中讀出并存儲到數據庫中，以便后面模塊的利用[13]。

圖2 NASTRAN SOL111部分流程Fig.2 Part of NASTRAN SOL111 sequence flowchart

這樣，就得到了進行模型修正工作所需要的信息。值得一提的是，為了節省計算機空間和文件的讀取速度，NASTRAN所輸出的文件都以二進制的形式給出。

3 模型修正方法的仿真驗證

為驗證修正方法的可行性，在圖 3所示的桁架模型上開展仿真試驗。此模型節點較多，模態也比較密集，常作為模型修正方法的范例[8-10]。本文對此結構采用矩形截面梁建模，共劃分75個單元，81個節點，每個節點有6個自由度。其材料參數分別如下：E=7.5×1011N/m2，密度ρ=2 800.0 kg/m3。水平、豎直以及對角線梁的截面同為A=0.15 m2，截面慣性積Ixx=0.003 125 m4，Iyy=0.001 125 m4。給定這些參數誤差，以此作為模型修正的對象（參見表1）。修正過程僅使用所有節點x、y、z三個方向的響應。因為在實際試驗中，加速度傳感器僅能測得這三個方向的響應，這樣做比較接近實際。

圖3 桁架模型Fig.3 The truss structure model

表1 模型對應位置的參數誤差Table 1 Parameter errors of model related locations

1）利用全部節點響應數據進行修正

首先利用全部節點的x、y、z三個方向的響應信息對模型進行修正，參數修正結果見表 2。可見利用所有的節點信息對模型進行修正，能得到非常好的效果，迭代次數也較少。

圖4 修正前后模型傳遞特性對比Fig.4 Comparison of response function curves before and after model modification

圖5 參數迭代歷程Fig.5 Parameter iteration history

表2 參數修正結果Table 2 The resulting parameters of model modification

2）利用含噪聲的部分節點響應數據進行修正

在實際的振動試驗中，往往只是測量了少部分節點的響應，另外測量信號還會不同程度地遭受噪聲的污染。為了模擬真實的試驗數據，此次修正只利用了部分節點的響應信息，同時在測量信號中加入 5%的高斯白噪聲來模擬真實數據。測量點分布情況如圖3所示。以節點43的基礎激勵傳遞函數為例，觀察修正過程，見圖6和圖7。參數修正結果見表3。

圖6 修正前后模型傳遞特性對比Fig.6 Comparison of response function curves before and after model modification

圖7 參數迭代歷程Fig.7 Parameter iteration history

表3 參數修正結果Table 3 The resulting parameters of model modification

仿真結果表明，此種方法在測量數據不足、試驗數據受到噪聲污染的情況下，仍然得到了非常好的效果。

4 結論

在振動試驗過程中，往往不能測得傳統的頻響函數或者進行模態辨識工作，這給航天器有限元模型的修正帶來了困難。本文提出利用基礎激勵傳遞特性的修正方法對航天器結構有限元模型進行修正，并以復雜桁架結構為對象進行了仿真分析。分析結果初步驗證了該方法的可行性，并確認了該方法在測量數據不完備且受到噪聲污染的情況下，仍然可以修正有限元模型中的誤差。此方法在工程中的應用還會遇到誤差定位、阻尼參數的修正等問題，今后還需就這些問題展開深入研究。

（References）

[1]馬興瑞, 于登云, 韓增堯, 等.星箭力學環境預示與試驗技術研究進展[J].宇航學報, 2006, 27(3)： 323-331

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

[3]Farhat C, Hemez F M.Updating finite element dynamic model using an element-by-element sensitivity methodology[J].AIAA Journal, 1993, 31(9)： 1702-1711

[4]Alvin K F.Finite element model update via bayesian estimation and minimization of dynamic residuals[J].AIAA Journal, 1997, 35 (5)： 879-886

[5]Hua H.On a stctistical optimization method used in finite element model updating[J].Journal of Sound and Vibration, 2000, 231(4)： 1071-1078

[6]Marwala T, Sibisi S.Finite element model updating via bayesian framework and modal properties[J].Journal of Aircraft, 2005, 42 (1)： 275-278

[7]Lin R, Ewins D J.Model updating using FRF data [C]//Proc.15thInternational Seminar on Modal Analysis, 1990

[8]Imregun M, Visser W J, Ewins D J.Finite element model updating using frequency response function data-1 theory and initial investigation[J].Mechanical Systems and Signal Processing, 1995, 9： 187-202

[9]Lin R, Ewins D J.Analytical model improvement using frequency reponse functions[J].Mechanical Systems and Signal Processing, 1994, 8(4)： 437-458

[10]Lin R, Zhu J.Finite element model updating using vibration test data under base excitation[J].Journal of Sound and Vibration, 2007, 303： 596-613

[11]秦仙蓉.基于靈敏度分析的結構計算模型修正技術及相關問題研究[D].南京航空航天大學學位論文, 2001

[12]陳德成, 魏震松, 等.有限元模型修正技術的工程應用[J].中國工程科學, 2001, 3(10)： 59-63

[13]Reymond M.NASTRAN 2004 reference manual [M].Los Angeles： MSC Software Corporation, 2003