王澤宇，劉 闖，馮咬齊

（北京衛(wèi)星環(huán)境工程研究所，北京 100094）

0 引言

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

然而無論是利用模態(tài)法還是頻響函數(shù)法對(duì)航天器模型進(jìn)行修正，都需要進(jìn)行額外的力學(xué)試驗(yàn)。如果能利用航天器的振動(dòng)試驗(yàn)數(shù)據(jù)完成模型修正工作，必然能節(jié)約成本，縮短研制周期。鑒于此，本文提出了利用基于基礎(chǔ)激勵(lì)傳遞特性的模型修正方法[10]來修正航天器結(jié)構(gòu)的有限元模型。

1 基于基礎(chǔ)激勵(lì)傳遞特性的模型修正方法

航天器振動(dòng)試驗(yàn)是將試驗(yàn)件由夾具緊固在振動(dòng)臺(tái)上，控制振動(dòng)臺(tái)給試驗(yàn)件提供一個(gè)加速度邊界條件，測(cè)量在此條件下一些關(guān)鍵部位的響應(yīng)。其主要目的是為了暴露航天器元器件、結(jié)構(gòu)在工藝和材料方面的缺陷，考核產(chǎn)品耐受力學(xué)環(huán)境的能力。

航天器結(jié)構(gòu)在振動(dòng)臺(tái)上的運(yùn)動(dòng)方程為

假設(shè)結(jié)構(gòu)做無阻尼振動(dòng)，上述方程可簡(jiǎn)化為

這樣方程(2)可寫成

方程(3)又可以寫成

根據(jù)節(jié)點(diǎn)是否與外界連接，剛度矩陣和質(zhì)量矩陣也可以相應(yīng)地寫成

式中下標(biāo)1代表對(duì)應(yīng)于與基礎(chǔ)連接的節(jié)點(diǎn)的自由度，下標(biāo)2代表與外界沒有連接的節(jié)點(diǎn)的自由度。將式(5)代入式(4)中，可以得到如下方程：

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

其中： ｛e｝ 為單位矩陣， ｛e｝=｛1， 1， 1，…｝T。這里定義一個(gè)新的傳遞函數(shù)｛H(ω)｝ = ｛u｝/u0，這個(gè)傳遞函數(shù)反映了結(jié)構(gòu)傳導(dǎo)基礎(chǔ)激勵(lì)響應(yīng)的特性，是結(jié)構(gòu)本身的固有屬性，可以用此特性作為數(shù)學(xué)模型與實(shí)際物理模型是否一致的標(biāo)準(zhǔn)并用來修正模型，它具體如下式所示：

對(duì)于理論模型和試驗(yàn)?zāi)Ｐ?，由?7)可以分別得到：

式中下標(biāo)A代表理論分析模型，X代表試驗(yàn)?zāi)Ｐ?。理論模型與試驗(yàn)?zāi)Ｐ椭g存在著偏差，設(shè)：

同時(shí)

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

假設(shè)

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

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

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

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

在實(shí)際測(cè)量中，不可能把有限元模型所有節(jié)點(diǎn)的響應(yīng)都測(cè)量出來，測(cè)量數(shù)據(jù)是不完備的。因此在實(shí)際應(yīng)用的過程中，必須首先對(duì)理論模型進(jìn)行縮聚[11]。結(jié)構(gòu)的動(dòng)剛度矩陣與頻響函數(shù)矩陣之間的關(guān)系： × =ZHI，即：

根據(jù)實(shí)際測(cè)量的自由度t與未測(cè)量的自由度n對(duì)模型進(jìn)行分塊。此外，與基礎(chǔ)連接的那些節(jié)點(diǎn)的自由度也要保留。通過方程(21)可以提取出如下兩個(gè)方程：

并可以得到

縮聚后的動(dòng)剛度矩陣與頻響函數(shù)矩陣應(yīng)該還是互逆的關(guān)系，所以有

由此，修正方程相應(yīng)的變?yōu)?/p>

此時(shí)的｛e｝base只包括模型中與基礎(chǔ)運(yùn)動(dòng)狀態(tài)相同的自由度，｛HX｝test與｛ΔH｝test也只對(duì)應(yīng)于進(jìn)行了測(cè)量的節(jié)點(diǎn)自由度。

2 模型修正方法的程序?qū)崿F(xiàn)

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

2.1 模型修正程序的設(shè)計(jì)思路

模型修正程序結(jié)構(gòu)見圖 1。整個(gè)修正程序分兩部分：一部分是利用NASTRAN接口語言編程，使之輸出進(jìn)行修正工作所需要的信息；另一部分為編寫 Matlab修正程序。接口程序 1 可以從NASTRAN的輸出文件中讀取出動(dòng)剛度矩陣和材料參數(shù)等信息，并將此信息輸送給修正程序。修正程序?qū)⑻幚磉@些信息并得到動(dòng)剛度矩陣對(duì)材料參數(shù)的靈敏度矩陣，完成修正參數(shù)和頻率點(diǎn)的選取等工作，組建修正方程并求解，然后將修正過的參數(shù)通過接口程序2傳遞給NASTRAN。如果計(jì)算結(jié)果不滿足控制條件，再調(diào)用NASTRAN進(jìn)行計(jì)算，進(jìn)而組建新的修正方程，得到新的參數(shù)。整個(gè)修正過程有兩個(gè)控制條件：迭代次數(shù)和參數(shù)的變化范圍。當(dāng)滿足其中一個(gè)條件時(shí)，修正程序運(yùn)行結(jié)束。

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

2.2 NASTRAN接口語言的編程實(shí)現(xiàn)方法

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

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

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

3 模型修正方法的仿真驗(yàn)證

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

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

表1 模型對(duì)應(yīng)位置的參數(shù)誤差Table 1 Parameter errors of model related locations

1）利用全部節(jié)點(diǎn)響應(yīng)數(shù)據(jù)進(jìn)行修正

首先利用全部節(jié)點(diǎn)的x、y、z三個(gè)方向的響應(yīng)信息對(duì)模型進(jìn)行修正，參數(shù)修正結(jié)果見表 2?？梢娎盟械墓?jié)點(diǎn)信息對(duì)模型進(jìn)行修正，能得到非常好的效果，迭代次數(shù)也較少。

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

圖5 參數(shù)迭代歷程Fig.5 Parameter iteration history

表2 參數(shù)修正結(jié)果Table 2 The resulting parameters of model modification

2）利用含噪聲的部分節(jié)點(diǎn)響應(yīng)數(shù)據(jù)進(jìn)行修正

在實(shí)際的振動(dòng)試驗(yàn)中，往往只是測(cè)量了少部分節(jié)點(diǎn)的響應(yīng)，另外測(cè)量信號(hào)還會(huì)不同程度地遭受噪聲的污染。為了模擬真實(shí)的試驗(yàn)數(shù)據(jù)，此次修正只利用了部分節(jié)點(diǎn)的響應(yīng)信息，同時(shí)在測(cè)量信號(hào)中加入 5%的高斯白噪聲來模擬真實(shí)數(shù)據(jù)。測(cè)量點(diǎn)分布情況如圖3所示。以節(jié)點(diǎn)43的基礎(chǔ)激勵(lì)傳遞函數(shù)為例，觀察修正過程，見圖6和圖7。參數(shù)修正結(jié)果見表3。

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

圖7 參數(shù)迭代歷程Fig.7 Parameter iteration history

表3 參數(shù)修正結(jié)果Table 3 The resulting parameters of model modification

仿真結(jié)果表明，此種方法在測(cè)量數(shù)據(jù)不足、試驗(yàn)數(shù)據(jù)受到噪聲污染的情況下，仍然得到了非常好的效果。

4 結(jié)論

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

（References）

[1]馬興瑞, 于登云, 韓增堯, 等.星箭力學(xué)環(huán)境預(yù)示與試驗(yàn)技術(shù)研究進(jìn)展[J].宇航學(xué)報(bào), 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]秦仙蓉.基于靈敏度分析的結(jié)構(gòu)計(jì)算模型修正技術(shù)及相關(guān)問題研究[D].南京航空航天大學(xué)學(xué)位論文, 2001

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

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