脈沖風(fēng)洞測力系統(tǒng)建模與載荷辨識方法研究

脈沖風(fēng)洞測力系統(tǒng)建模與載荷辨識方法研究

王鋒1,2，賀偉1，毛鵬飛1，張小慶1

(1.中國空氣動力研究與發(fā)展中心吸氣式高超聲速技術(shù)研究中心，綿陽6210002.中國空氣動力研究與發(fā)展中心高超聲速沖壓發(fā)動機技術(shù)重點實驗室，綿陽621000)

摘要：將載荷辨識技術(shù)應(yīng)用于脈沖燃燒風(fēng)洞模型測力。用子結(jié)構(gòu)綜合法建立了測力試驗系統(tǒng)的動力學(xué)模型，在時域內(nèi)將動力學(xué)方程進行離散，建立起天平測量信號與模型氣動載荷歷程之間的線性關(guān)系，作為載荷辨識的模型。采用Tikhonov正則化和子空間投影法相結(jié)合的混合正則化方法，將高維的、不適定的載荷辨識問題轉(zhuǎn)化為低維的適定問題，以利于快速求解。提出了一種新方法來確定合適的投影子空間維數(shù)，然后應(yīng)用L曲線準(zhǔn)則來尋找低維正則化問題的最優(yōu)正則化參數(shù)。最后通過算例驗證了系統(tǒng)建模方法的精度和載荷辨識算法的有效性與穩(wěn)定性。

關(guān)鍵詞：脈沖風(fēng)洞;動力學(xué)建模;載荷辨識;正則化;反問題

中圖分類號:V231;O32文獻標(biāo)志碼：A

Dynamicmodelingoftestingsysteminimpulsefacilitiesandloadidentificationmethod

WANG Feng1,2, HE Wei1,MAOPeng-fei1,ZHANGXiao-qing1(1.Air-breathingHypersonicTechnologyResearchCenter,ChinaAerodynamicsResearchandDevelopmentCenter,Mianyang621000,China;2.ScienceandTechnologyonScramjetLaboratory,ChinaAerodynamicsResearchandDevelopmentCenter,Mianyang621000,China)

Abstract:Load identification method was applied to force measurement in impulse combustion facilities. The dynamic model of the testing system was built using substructure synthesis method. Then, the differential equations were discretized to establish the linear relationship between the outputs of balance and the dynamic forces history acting on the aircraft model. A hybrid regularization method which combines the Tikhonov regularization and subspace projection method was proposed to solve the load identification problem. The hybrid method converts the ill-conditioned large scale problem into a well-posed small scale problem which can be solved easily. A new method was proposed to determine the size of the projection subspace. The L-curve criteria was employed to search the optimal regularization parameter of the dimension-reduced problem. The accuracy of the structural dynamic modeling method and the effective and stability of the load identification method were validated by a numerical example.

Keywords:impulsefacilities;structuraldynamicmodeling;loadidentification;regularization;inverseproblem

吸氣式高超聲速飛行器是當(dāng)前國際航空航天技術(shù)的熱點研究領(lǐng)域，地面風(fēng)洞試驗是最重要的研究手段，脈沖燃燒風(fēng)洞是國內(nèi)目前開展大尺度超燃發(fā)動機和機體-推進一體化高超聲速飛行器試驗研究的主力設(shè)備之一，具有口徑大、建造和運行成本較低的優(yōu)點。由于其有效的試驗時間較短，只有約200ms～300ms，給模型測力試驗帶來一定的困難。目前采用常規(guī)的測力方法，為了在短時間內(nèi)測得模型的氣動力，需盡量降低模型的質(zhì)量，同時盡量提高天平的剛度，以便提高試驗系統(tǒng)的響應(yīng)速度。吸氣式飛行器采用機身和推進系統(tǒng)一體化設(shè)計，內(nèi)外流緊密耦合，為了反映真實的氣動-推進特性，希望采用全尺寸或接近全尺寸的試驗?zāi)Ｐ?。所以，目前在風(fēng)洞開展試驗的飛行器模型尺寸和質(zhì)量越來越大，導(dǎo)致試驗系統(tǒng)的響應(yīng)速度降低，彈性振動也難以衰減，基于穩(wěn)態(tài)測量思想的常規(guī)測力方法面臨挑戰(zhàn)。因此，有必要研究新的測力方法，提高脈沖風(fēng)洞測力的準(zhǔn)確性，更好地發(fā)揮其優(yōu)勢。

本文研究采用載荷辨識技術(shù)進行間接測力的方法。載荷辨識是根據(jù)測量得到的結(jié)構(gòu)響應(yīng)(應(yīng)變、位移、速度或加速度等)，結(jié)合結(jié)構(gòu)的動力學(xué)方程(或傳遞函數(shù))反求結(jié)構(gòu)所受的載荷歷程。作為結(jié)構(gòu)動力學(xué)的第二類反問題，載荷辨識技術(shù)獲得了廣泛的研究和應(yīng)用[1-5]，其難點在于問題本身的不適定性，因此在求解時普遍使用正則化方法，文獻[6]對此進行了較全面的總結(jié)。

載荷辨識應(yīng)用于模型測力，不需要結(jié)構(gòu)響應(yīng)達到穩(wěn)態(tài)，可以放松對天平剛度的要求和對模型質(zhì)量的限制，同時，辨識得到的是載荷隨時間的變化過程，能反映更豐富的試驗動態(tài)信息，如發(fā)動機點火時的推力建立過程，這對于精細地研究高超聲速飛行器的氣動-推力特性具有重要價值。Simmons等[7]最早將載荷辨識技術(shù)用于激波風(fēng)洞模型阻力測量，使用的應(yīng)力波天平，后來又將該技術(shù)發(fā)展用于多分量同時測量[8-10]。考慮到脈沖燃燒風(fēng)洞的試驗時間遠大于激波風(fēng)洞(5ms左右)，本文的方法在試驗系統(tǒng)建模、天平標(biāo)定方面均和基于應(yīng)力波天平的測力技術(shù)有所不同，并且提出了新的載荷辨識算法。

1試驗系統(tǒng)的動力學(xué)建模

載荷辨識需已知系統(tǒng)的動態(tài)特性，所以首先建立試驗系統(tǒng)的結(jié)構(gòu)動力學(xué)方程，并希望方程的自由度盡量少，以降低后續(xù)的計算量。

圖1　試驗系統(tǒng)示意圖 Fig.1 Schematic of the testing system

1.1試驗系統(tǒng)結(jié)構(gòu)組成

風(fēng)洞測力試驗系統(tǒng)主要由飛行器模型、天平、支架以及數(shù)據(jù)采集系統(tǒng)組成，對于結(jié)構(gòu)動力學(xué)建模，僅考慮前三部分。由于飛行器模型一般比較長、比較重，同時，為了避開其腹部的推進流道，通常采用背部支撐的方式固定在風(fēng)洞中，系統(tǒng)的剖面如圖1所示。天平的浮動框與試驗?zāi)Ｐ凸踢B，天平的固定框與支架上端固連，支架下端固定于地面。在來流作用下，試驗?zāi)Ｐ蛶犹炱礁涌蜻\動，使天平的彈性元件產(chǎn)生應(yīng)變，利用應(yīng)變片產(chǎn)生輸出信號。

1.2建模方法

可以用有限元方法直接對整個試驗系統(tǒng)進行建模，然后通過模態(tài)降階得到低階模型，但是，為了在動力學(xué)方程中保留模型的剛體自由度和直接使用天平的標(biāo)定數(shù)據(jù)，在此采用子結(jié)構(gòu)方法來建模。將系統(tǒng)分為試驗飛行器、天平和支架三個子結(jié)構(gòu)。由于試驗時天平的浮動框與試驗飛行器固連在一起，固定框與支架固連在一起，因此，在建模時把浮動框歸入試驗飛行器，把固定框歸入支架，而天平子系統(tǒng)本身只考慮連接浮動框與固定框的彈性測量元件。

子結(jié)構(gòu)采用Craig-Bampton方法[11](后面簡稱C-B方法)建模。C-B方法的子結(jié)構(gòu)模態(tài)包括界面約束模態(tài)(后面簡稱約束模態(tài))和固定界面正規(guī)模態(tài)(后面簡稱正規(guī)模態(tài))，因此，首先要確定子結(jié)構(gòu)間的連接界面，并根據(jù)需要作必要的簡化。

前面已將天平浮動框歸入試驗飛行器、將天平固定框歸入支架，所以試驗飛行器與天平的連接面就是天平測量元件與浮動框的對接面，天平與支架的連接面就是測量元件與固定框的連接面。實際中，天平的浮動框和固定框的剛度遠遠大于測量元件的剛度，相對于彈性測量元件，兩個框可近似看作剛體，于是它們與測量元件的對接界面均可簡化為一個6自由度的點。將試驗飛行器(含浮動框)的界面點記作點Ⅰ、天平(測量元件)在浮動框一側(cè)的界面點記作點Ⅱ、天平(測量元件)在固定框一側(cè)的界面點記作點Ⅲ、支架的界面點記作點Ⅳ。

1.2.1試驗飛行器的模型

試驗飛行器用有限元方法建模。為了將其與天平的對接面簡化為一個點，在適當(dāng)?shù)奈恢茫缳|(zhì)心處，設(shè)置一節(jié)點作為主節(jié)點(即點Ⅰ)，將浮動框上與天平彈性元件對接面上的節(jié)點作為從節(jié)點，用剛性單元相連。在此基礎(chǔ)上，用成熟的商業(yè)有限元軟件可以直接進行子結(jié)構(gòu)分析，輸出子結(jié)構(gòu)的質(zhì)量和剛度矩陣，從而建立其動力學(xué)方程

(1)

這里，下標(biāo)B(Boundary)表示界面約束模態(tài)，下標(biāo)N(Normal)表示正規(guī)模態(tài)。ηⅠ,B是試驗飛行器約束模態(tài)的位移向量，為6個剛體位移

(2)

ηⅠ,N為正規(guī)模態(tài)位移向量，假設(shè)取前nm階正規(guī)模態(tài)，則

(3)

因此，飛行器子系統(tǒng)共有6+nm個自由度。因為約束模態(tài)均為剛體模態(tài)，所以子陣KⅠ,BB為零矩陣，而子陣MⅠ,BB為飛行器的剛體質(zhì)量矩陣(參考點為點Ⅰ)，而fⅠ,B正是飛行器所受氣動力和力矩向量(關(guān)于點Ⅰ的)。fⅠ,N是正規(guī)模態(tài)對應(yīng)的廣義力，其第Ⅰ個分量的計算式如下

(4)

1.2.2天平子系統(tǒng)的模型

天平的浮動框和固定框分別歸屬飛行器和支架后，天平子系統(tǒng)僅剩下彈性測量元件，由于測量元件尺寸相對較小，質(zhì)量也較輕，為簡化模型，忽略其慣性，從而將天平子系統(tǒng)建模為無質(zhì)量的廣義彈簧，其與試驗飛行器和支架分別有一個對接點(點Ⅱ和點Ⅲ)，所以天平子系統(tǒng)共12個自由度，其模型為12×12的剛度矩陣。簡化模型示意圖見圖2。

雖然可以用有限元計算天平的剛度矩陣，但天平是試驗系統(tǒng)的關(guān)鍵環(huán)節(jié)，實際中要經(jīng)過校準(zhǔn)，顯然利用校準(zhǔn)數(shù)據(jù)來得到其剛度矩陣將更為準(zhǔn)確。

圖2　天平子系統(tǒng)的力學(xué)模型 Fig.2 Mechanical model of the balance

點Ⅱ、點Ⅲ的位移向量分別記作

天平子系統(tǒng)的力學(xué)模型為

(5)

式中，fⅡ、fⅢ分別是作用于點Ⅱ、點Ⅲ的載荷向量，每個載荷向量包括3個力和3個力矩。在校準(zhǔn)時，是將天平固定框完全約束，通過剛性足夠的夾具與浮動框連接，然后在夾具上施加載荷，載荷施加點稱為校心，一般位于天平的幾何中心線上。在此，假設(shè)點Ⅱ即為校心。點Ⅲ選在同一中心線上位于固定框的一側(cè)，與固定框剛性連接，與點Ⅱ之間的距離為h。那么，校準(zhǔn)時的邊界條件即為：點Ⅲ完全固定，在點Ⅱ施加校準(zhǔn)載荷。于是在方程(5)中，ηⅢ=0，從而有

KⅡ,ⅡηⅡ=fⅡ

(6)

KⅢ,ⅡηⅡ=fⅢ

(7)

代入式(6)有

KⅡ,ⅡRⅡ=PⅡ

(8)

于是得到

(10)

利用式(7)和式(8)可得

KⅢ,ⅡRⅡ=-TB1PⅡ=-TB1KⅡ,ⅡRⅡ

(11)

所以有

KⅢ,Ⅱ=-TB1KⅡ,Ⅱ

(12)

因為剛度矩陣是對稱的，所以有

KⅡ,Ⅲ=KTⅢ,Ⅱ=-KTⅡ,ⅡTTB1=-KⅡ,ⅡTTB1

(13)

現(xiàn)在尚有KⅢ,Ⅲ待求。假設(shè)解除點Ⅲ的所有約束，允許天平產(chǎn)生剛體運動，在點Ⅱ沿6個自由度依次施加微位移Δx、Δy、Δz、Δθx、Δθy、Δθz，寫成位移向量矩陣

由于現(xiàn)在點Ⅲ自由，它要產(chǎn)生相協(xié)調(diào)的一組位移向量，記作ΔRⅢ，忽略二階小量，易得

(14)

因為點Ⅲ不受力，將ΔRⅡ、ΔRⅢ代入式(5)的第二組方程，有

KⅢ,ⅡΔRⅡ+KⅢ,ⅢΔRⅢ=0

(15)

再利用式(14)可得

(16)

KⅢ,Ⅲ=TB1KⅡ,ⅡTTB1

(17)

至此，利用天平校準(zhǔn)時的加載及位移測量數(shù)據(jù)可以計算出其剛度矩陣。天平校準(zhǔn)時還要獲取天平方程，即天平傳感器(一般是應(yīng)變計)輸出與校心載荷的關(guān)系，這是常規(guī)天平校準(zhǔn)的首要目的，具體過程在此不詳述，假設(shè)已獲得如下天平方程

v=NB,ⅡfⅡ

(18)

式中：v是天平輸出(一般是電壓)，其與載荷向量是等長的，系數(shù)矩陣NB,Ⅱ是方陣，且非奇異?，F(xiàn)在需要推導(dǎo)天平輸出與點Ⅱ、Ⅲ位移的關(guān)系。因為校準(zhǔn)時固定框(點Ⅲ)被完全約束，所以ηⅢ=0，式(6)成立，將其代入式(18)得

v=NB,ⅡKⅡ,ⅡηⅡ?SⅡηⅡ

(19)

這里，將SⅡ=NB,ⅡKⅡ,Ⅱ稱為點Ⅱ?qū)?yīng)的傳感矩陣。尚需求點Ⅲ對應(yīng)的傳感矩陣，將其記作SⅢ，則兩點都產(chǎn)生位移時，天平的輸出應(yīng)該是

v=SⅡηⅡ+SⅢηⅢ

(20)

如果給定點Ⅱ的位移ηⅡ，而讓點Ⅲ自由運動，則根據(jù)式(14)有ηⅢ=TB2ηⅡ，此時天平僅有剛體位移，輸出應(yīng)該為零，根據(jù)式(20)有

SⅡηⅡ+SⅢTB2ηⅡ=0

(21)

所以有

(22)

代入式(20)，有

此即天平子系統(tǒng)的傳感方程。

1.2.3支架子系統(tǒng)模型

支架子系統(tǒng)亦用有限元方法建模，其過程與試驗飛行器子系統(tǒng)一樣，需設(shè)置一節(jié)點作為主節(jié)點(即點Ⅳ)，選取固定框上與彈性元件對接處的所有節(jié)點為從節(jié)點，用剛性單元相連。在此基礎(chǔ)上，用商業(yè)有限元軟件執(zhí)行子結(jié)構(gòu)分析可以獲取子結(jié)構(gòu)的質(zhì)量和剛度矩

陣，得到其動力學(xué)方程為

(24)

支架約束模態(tài)的位移向量為

假設(shè)正規(guī)模態(tài)取前ns階，記為

則方程(24)共6+ns個自由度。

1.2.4系統(tǒng)模型綜合

圖3　點Ⅰ和點Ⅱ的相對位置 Fig.3 Position relationship of point Ⅰ and point Ⅱ

盡管各子結(jié)構(gòu)之間的界面已簡化為一個點，但是各子結(jié)構(gòu)的界面點是為建模方便而選取的，相鄰子結(jié)構(gòu)的“界面”不一定重合，需要用剛性桿將其連接，以構(gòu)成完整系統(tǒng)。

各子結(jié)構(gòu)的界面點盡管不一致，但應(yīng)該都是選在其對稱面(x-y平面，見圖1)內(nèi)。假設(shè)點Ⅰ和點Ⅱ的相對位置如圖3所示，則用剛性桿連接后，其位移(微量)向量之間存在如下關(guān)系

ηⅡ=

(25)

類似地，可建立點Ⅲ和點Ⅳ位移向量的關(guān)系

(26)

利用式(25)、(26)，可將三個子結(jié)構(gòu)的方程(1)、(5)和(24)綜合在一起，形成整個系統(tǒng)的動力學(xué)方程，其中獨立的自由度向量只有ηⅠ,B、ηⅠ,N、ηⅣ,B、ηⅣ,N，列在一起記作

系統(tǒng)的動力學(xué)方程為

(27)

其中

K=

實際中，天平置于飛行器內(nèi)部，是不受外載荷作用的，因此fⅡ和fⅢ均為零。支架外部可以加整流外罩進行遮擋，也不受氣動載荷，因此fⅣ,B和fⅣ,N為零。非零載荷只有試驗飛行器的氣動載荷，其被分解為約束模態(tài)載荷和正規(guī)模態(tài)載荷，根據(jù)它們的定義，實際中前者是主要的，也是人們通常所指的氣動載荷，后者則相對較小，而且不直接引起天平輸出，因此可將fⅠ,N忽略，僅保留fⅠ,B。于是，試驗系統(tǒng)的力載荷向量簡化為

實際結(jié)構(gòu)總是有阻尼的，考慮阻尼后系統(tǒng)動力學(xué)方程為

(28)

傳感方程(23)也需要用η來表示，易得

(29)

2載荷辨識模型

式(28)、(29)構(gòu)成完整的試驗系統(tǒng)動力學(xué)輸入輸出方程。通過擴維降階將其變?yōu)闋顟B(tài)空間形式

(30)

y=Rx

(31)

式中：輸出y已知，輸入x待求。

3載荷辨識方法

3.1子空間投影法

載荷辨識是典型的反問題，方程(31)通常是病態(tài)的，輸出向量y中的噪聲會使得解x嚴重偏離真解，失去利用價值，而實際中測量噪聲(誤差)是不可避免的。為了得到有意義的解，必須對其進行正則化，Tikhonov正則化方法[15-16]是最常用的，它用如下最小化問題的解作為問題(31)的近似解[17]

(32)

式中α(>0)是正則化參數(shù)。上式與如下問題等價

(33)

實際中，由于結(jié)構(gòu)頻率較高，測量的采樣率也很高，若待辨識的時段稍長則矩陣R的規(guī)模會很大，而在優(yōu)化選擇正則化參數(shù)α的過程中，要反復(fù)求解式(33)，因此很難直接基于式(33)來求解。

為了使實際問題的求解更加可行，這里采用子空間投影法，將問題(32)的解限制在某個維數(shù)較低的解域上，從而降低求解規(guī)模。假設(shè)選定了一組基向量w1、w2、…、wk，將近似解x限制在由其張成的k維空間Wk上，即

x∈Wk=span{w1,w2,…,wk}

令

(34)

(35)

其等價問題為

(WTkRTRWk+αWTkWk)cα=(RWk)Ty

(36)

方程(36)的系數(shù)矩陣大小為k×k，通常k遠遠小于R的行數(shù)和列數(shù)，求解規(guī)模大大降低。

顯然，投影子空間的選擇，也就是投影基向量w1、w2、…、wk的選擇是影響近似解好壞的關(guān)鍵，為了逼近原問題的解，該子空間需盡量包含系數(shù)矩陣R前k個右奇異向量的信息。Krylov子空間具備這樣的特性[18]。k維Krylov子空間的定義如下[16, 18]

Krylov子空間的基向量可以按定義直接計算，但采用Golub-Kahan雙對角迭代法[19]可以得到性質(zhì)更好的基向量，并為問題的求解帶來方便?；舅惴ㄈ缦拢?/p>

w0=0，β1z1=y

對i=1,2,3,…,k，執(zhí)行以下迭代：

αiwi=RTzi-βiwi-1

βi+1zi+1=Rwi-αkzi

式中，實數(shù)αi、βi的作用是分別使向量wi和zi單位化。迭代的結(jié)果得到如下關(guān)系式

RWk=Zk+1Bk

(37)

式中，

并且，Wk、Zk+1矩陣各自的列向量是相互正交的單位向量。將式(37)代入式(36)并利用Wk、Zk+1的正交性質(zhì)，式(36)可簡化為

(BTkBk+αI)cα=BTkZTk+1y

(38)

根據(jù)前面的迭代算法，y=β1z1，所以BTkZTk+1y=β1BTkZTk+1z1=β1BTke1，其中e1標(biāo)準(zhǔn)歐氏空間的第1個向量。于是，式(38)進一步簡化為

(BTkBk+αI)cα=β1BTke1

(39)

注意BTkBk+αI為三對角矩陣，因此上式可以用追趕法快速求解。得到cα后，最終解由下式給出

(40)

3.2確定正則化參數(shù)的L曲線準(zhǔn)則

(41)

(42)

式中zα是如下問題的解

(43)

當(dāng)采用子空間投影后,有

(44)

(45)

為減小計算量，zα也用子空間投影法來近似計算，即令

zα≈Wkζ

(46)

從而將問題(43)變?yōu)槿缦氯菀浊蠼獾膯栴}

(47)

于是，曲率的κ(α) 計算量大大降低，使得L曲線準(zhǔn)則容易實施。

3.3子空間維數(shù)的確定方法

現(xiàn)在討論如何確定適當(dāng)?shù)淖涌臻g維數(shù)k。顯然，對同一個正則化參數(shù)α，取不同的子空間維數(shù)，得到的正則化解是不同的，因此，不同子空間維數(shù)下作出的L曲線是不重合的。但是，由于k維Krylov子空間的基向量具有逼近矩陣R前k個右奇異向量的性質(zhì)，隨著k的增加，后面的基向量對正則解的貢獻逐漸減弱，所以L曲線本身隨k的增加將趨于收斂。在此，取如下指標(biāo)作為判斷L曲線收斂的依據(jù)

J=‖xα‖2·‖rα‖2

(48)

3.4求解算法

根據(jù)前面的討論，對問題(31)完整的求解算法包括兩步，首先用3.3節(jié)的方法確定合適的子空間維數(shù)，然后在該子空間下用3.2節(jié)的L曲線準(zhǔn)則確定最優(yōu)的正則化參數(shù)，并得到對應(yīng)的正則化解。第一步的詳細算法流程為：

(1)給定初始子空間維數(shù)k0和維數(shù)的增量Δk；給定i0和nα，計算出αi序列；給定參數(shù)ε；

(2)執(zhí)行k0步Golub-Kahan雙對角迭代，得到式(37)的各矩陣；

(3)針對每個αi求解方程(39)，由(44)、(45)式得到‖xα‖2、‖rα‖2，計算指標(biāo)函數(shù)值序列J(αi,k)，找出最小值對應(yīng)的序號ik和最小值Jk；

(4)將當(dāng)前子空間下的結(jié)果ik、Jk與上一步的結(jié)果ik-Δk、Jk-Δk作比較，若滿足條件

終止，得到子空間維數(shù)k；否則往下執(zhí)行；

(5)k←k+Δk，繼續(xù)執(zhí)行Δk步Golub-Kahan，擴充子空間，轉(zhuǎn)到第3步。

確定了子空間維數(shù)后，對應(yīng)低維正則化問題的L曲線就確定了，現(xiàn)在需要尋找最大曲率點的正則化參數(shù)。這時，L曲線上各點的曲率κ是參數(shù)α的函數(shù)，即式(41)，可以采用一維搜索算法[23]求函數(shù)κ(α)的近似極大值解，起點可以取前面確定子空間維數(shù)時得到的αik，具體算法在此不再詳細列出。得到問題(39)的正則解cα后代入式(40)，最終得到問題(31)的正則化解。

4數(shù)值算例

考慮一個虛擬的試驗飛行器及其測力系統(tǒng)，其模型、天平及支架如圖4所示。飛行器長2.1m，材料為超硬鋁，重255.24kg，天平和支架為鋼材。

圖4　試驗系統(tǒng)三維圖 Fig.4 View of the testing system

用有限元法直接對全系統(tǒng)建模，稱其為全階模型，作為簡化模型的對比基準(zhǔn)。天平也用有限元建模，通過靜力計算，模擬實際校準(zhǔn)過程，用得到的數(shù)據(jù)建立簡化模型的剛度矩陣。每個測量通道取四個點的應(yīng)變組成惠氏電橋，根據(jù)經(jīng)驗設(shè)置應(yīng)變片的靈敏度，將應(yīng)變轉(zhuǎn)化為電壓，據(jù)此計算天平傳感矩陣。

4.1模態(tài)計算結(jié)果

用有限元軟件分別對試驗飛行器和支架進行子結(jié)構(gòu)分析并輸出質(zhì)量、剛度矩陣，得到子結(jié)構(gòu)動力學(xué)模型。其中試驗飛行器取7個正規(guī)模態(tài)，支架取6個正規(guī)模態(tài)，加上各自的6個約束模態(tài)，系統(tǒng)簡化模型共25個自由度，而全階模型的節(jié)點數(shù)超過30萬，自由度超過100萬。對全階模型和簡化模型進行模態(tài)分析，得到的前10階固有頻率列于表1，可見簡化模型縱向和橫側(cè)向的前兩階頻率精度都是比較高的。

表1　固有頻率計算結(jié)果對比

4.2瞬態(tài)響應(yīng)結(jié)果

某脈沖燃燒風(fēng)洞典型的來流動壓變化曲線如圖5所示(對原始數(shù)據(jù)作低通濾波的結(jié)果)?，F(xiàn)在計算試驗系統(tǒng)在其作用下的天平輸出。對于全階模型，用工程方法[24]計算試驗飛行器表面氣流壓力分布，模擬真實的受力狀態(tài)。對于簡化模型，先將飛行器的壓力分布系數(shù)積分，得到合力和力矩(關(guān)于質(zhì)心)系數(shù)，再作為集中力加載。計算時，兩種模型的阻尼均以模態(tài)阻尼方式施加，各階阻尼比均取0.005。

假設(shè)側(cè)滑角為零，則氣動載荷只有阻力、升力和俯仰力矩。在圖5所示來流作用下，全階模型和簡化模型得到的天平三個通道的輸出對比如圖6所示，低頻信號吻合良好。高頻響應(yīng)的誤差來自兩個方面，一是簡化模型本身的高頻特性誤差較大，二是兩種模型的加載方式不同，全階模型加載的是分布力，而簡化模型加載的是等效集中力，忽略了式(1)中的fⅠ,N。

由模態(tài)分析和瞬態(tài)響應(yīng)結(jié)果的對比，說明低階簡化模型反映了試驗系統(tǒng)主要的動力學(xué)特性。

圖5　典型的風(fēng)洞來流動壓曲線 Fig.5 The history of dynamic pressure of the flow

圖6　瞬態(tài)響應(yīng)對比 Fig.6 Comparison of transient responses

圖7　載荷辨識結(jié)果對比 Fig.7 Comparison of the identified loads

4.3載荷辨識結(jié)果

以全階模型計算的瞬態(tài)響應(yīng)當(dāng)作試驗測量結(jié)果，對飛行器所受載荷歷程進行辨識。輸出采樣率取5 000Hz，采樣時長0.8s，有3個輸出，3個待辨識載荷，原始問題式(31)中矩陣R的大小為12 000×12 003。根據(jù)對真實試驗測量信號的統(tǒng)計分析，測量噪聲的標(biāo)準(zhǔn)差約為0.002 3mV，在此按0.005mV給計算輸出加白噪聲。計算時，初始子空間維數(shù)取550，增量取10，子空間收斂判定參數(shù)ε取0.001，正則化參數(shù)序列取10-30～10-20之間15個對數(shù)均布的點。最終確定的子空間維數(shù)為920，遠低于原始問題的維數(shù)，最優(yōu)正則化參數(shù)為α=1.453 8×10-24。3個載荷的辨識結(jié)果與真實值的對比見圖7，二者吻合良好。定義辨識誤差為

(49)

對于辨識結(jié)果中的噪聲，可以通過適當(dāng)?shù)钠交幚韥碛行У亟档?。圖8是對阻力曲線用20點平均法進行平滑處理后的結(jié)果，噪聲大大降低，但載荷本身的低頻信息仍得以保留。

圖8　平滑后的阻力曲線 Fig.8 Smoothed drag force

現(xiàn)在考察測量信號中不同噪聲強度對辨識結(jié)果的影響。分別取噪聲的標(biāo)準(zhǔn)差為0.002mV、0.004mV、0.008mV、0.016mV、0.032mV和0.064mV時得到的辨識結(jié)果如表2所示。其中，辨識誤差取的是3個載荷的平均誤差。可見，隨著噪聲強度的增加，子空間維數(shù)在減小、正則化參數(shù)在增大，即對問題的正則化作用不斷加強，相應(yīng)地，辨識誤差也在逐漸變大。均符合預(yù)期的規(guī)律。

表2　噪聲對辨識結(jié)果的影響

表3　不同計算條件下的結(jié)果

4.4算法的穩(wěn)定性分析

輸出測量信號中的噪聲是隨機的，因此算法每次執(zhí)行得到的結(jié)果亦可能有所不同，包括子空間大小、正則化參數(shù)以及辨識結(jié)果的精度。但是噪聲的統(tǒng)計特性基本不變，所以，對于一個穩(wěn)定的算法，計算結(jié)果不應(yīng)有大幅波動。另外，子空間維數(shù)增量Δk和正則化參數(shù)序列αi的設(shè)置也會對計算過程和結(jié)果產(chǎn)生影響。

現(xiàn)分別取Δk等于5、10、15，每種情況執(zhí)行5次辨識算法，每次的噪聲隨機生成，正則化參數(shù)序列也在10-30～10-20之間隨機生成15個點。計算結(jié)果如表3所示，可見盡管子空間維數(shù)有所波動，但變化不大，而辨識結(jié)果的精度也變化很小，表明算法是穩(wěn)定、魯棒的。

5結(jié)論

本文應(yīng)用子結(jié)構(gòu)綜合法，并充分利用天平校準(zhǔn)數(shù)據(jù)，建立了試驗系統(tǒng)的低維動力學(xué)模型，可以比較準(zhǔn)確地反映系統(tǒng)的低頻動態(tài)特性。通過對動力學(xué)方程進行時域離散，建立起試驗飛行器氣動力(力矩)與天平響應(yīng)之間的線性關(guān)系方程，用來進行氣動載荷的辨識。

將Tikhonov正則化方法與基于Krylov子空間的投影法結(jié)合，把高維的不適定問題變?yōu)榈途S的適定問題，并充分利用Golub-Kahan雙對角迭代算法的性質(zhì)，大大提高了求解效率。

提出了確定投影子空間大小的方法，該方法僅需計算低維問題的正則解及剩余向量的范數(shù)，計算簡單。子空間維數(shù)確定后，問題的規(guī)模大大減小，再應(yīng)用L曲線準(zhǔn)則確定最優(yōu)的正則化參數(shù)，計算量大為降低，用一維搜索算法可快速得到最優(yōu)參數(shù)。數(shù)值試驗表明，算法具有良好的穩(wěn)定性。

在對試驗系統(tǒng)進行建模時，結(jié)構(gòu)的阻尼參數(shù)是假定已知的，而實際中通常難以準(zhǔn)確知道結(jié)構(gòu)的阻尼參數(shù)，因此，利用測量響應(yīng)，對載荷和阻尼參數(shù)同時進行辨識是值得進一步研究的課題。

參考文獻

[1]HollandsworthPE,BusbyHR.Impactforceidentificationusingthegeneralinversetechnique[J].InternationalJournalofImpactEngineering, 1989, 8(4): 315-322.

[2]ChoiK,ChangFK.Identificationofimpactforceandlocationusingdistributedsensors[J].AIAAJournal, 1996, 34(1): 136-142.

[3]MaoYM,GuoXL,ZhaoY.Astatespaceforceidentificationmethodbasedonmarkovparametersprecisecomputationandregularizationtechnique[J].JournalofSoundandVibration, 2010, 329(15): 3008-3019.

[4]HwangJS,KareemA,KimH.Windloadidentificationusingwindtunneltestdatabyinverseanalysis[J].JournalofWindEngineeringandIndustrialAero-dynamics, 2011, 99(1): 18-26.

[5]RonasiH,JohanssonH,LarssonF.Anumericalframeworkforloadidentificationandregularizationwithapplicationtorollingdiscproblem[J].Compu-ters&Structures, 2011, 89: 38-47.

[6]王林軍. 正則化方法及其在動態(tài)載荷辨識中的應(yīng)用[D]. 長沙：湖南大學(xué), 2011.

[7]SandersonSR,SimmonsJM.Dragbalanceforhypervelocityimpulsefacilities[J].AIAAJournal, 1991, 29(12): 2185-2191.

[8]SmithAL,MeeDJ,DanielWJT,etal.Design,modellingandanalysisofasixcomponentforcebalanceforhypervelocitywindtunneltesting[J].ComputersandStructures, 2001, 79: 1077-1088.

[9]RobinsonMJ,MeeDJ,TsaiCY,etal.Three-componentforcemeasurementsonalargescramjetinashocktunnel[J].JournalofSpacecraftandRockets, 2004, 41(3): 416-425.

[10]RobinsonMJ,HannemannK,SchrammJM.Designandimplementationofaninternalstresswaveforcebalanceinashocktunnel[J].CEASSpaceJournal, 2011, (1): 45-57.

[11]CraigRR,Jr,BamptonMCC.Couplingofsubstructuresfordynamicanalysis[J].AIAAJournal, 1968, 6: 1313-1319.

[12]郭杏林, 毛玉明, 趙巖, 等. 基于Markov參數(shù)精細積分法的載荷識別研究[J]. 振動與沖擊, 2009, 28(3): 27-31.

GUOXin-lin,MAOYu-ming,ZHAOYan,etal.Loadidentificationbasedonprecisetime-stepintegrationforMarkovparameters[J].JournalofVibrationandShock, 2009, 28(3): 27-31.

[13]ZhongWX,WilliamsFW.Aprecisetimeintegrationmethod[J].JournalofMechanicalEngineeringScience, 1994, 208: 427-430.

[14]韓旭, 劉杰, 李偉杰, 等. 時域內(nèi)多源動態(tài)載荷的一種計算反求技術(shù)[J]. 力學(xué)學(xué)報, 2009, 41(4): 595-602.

HANXu,LIUJie,LIWei-jie,etal.Acomputationalinversetechniqueforreconstructionofmultisourceloadsintimedomain[J].ChineseJournalofTheoreticalandAppliedMechanics, 2009, 41(4): 595-602.

[15]CalvettiD,MorigiS,ReichelL,etal.Tikhonovregulariza-tionandtheL-curveforlargediscreteill-posedproblems[J].JournalofComputationalandAppliedMathematics, 2000, 123(1-2): 423-446.

[16]LampeJ,ReichelL,VossH.Large-scaletikhonovregulari-zationviareductionbyorthogonalprojection[J].LinearAlgebraanditsApplications, 2012, 436(8): 2845-2865.

[17]BauerF,LukasMA.Comparingparameterchoicemethodsforregularizationofill-posedproblems[J].MathematicsandComputersinSimulation, 2011, 81(9): 1795-1841.

[18]HansenPC.Discreteinverseproblems:Insightandalgorithms[M].Philadelphia:SocietyforIndustrialandAppliedMathematics, 2010.

[19]GolubG,KahanW.Calculatingthesingularvaluesandpseudo-inverseofamatrix[J].JournaloftheSocietyforIndustrial&AppliedMathematics,SeriesB:NumericalAnalysis, 1965, 2(2): 205-224.

[20]HansenPC,O’learyD.Theuseofl-curveintheregulariza-tionofdiscreteill-posedproblems[J].SIAMJournalonScientificComputing, 1993, 14: 1487-1503.

[21]VogelCR.Computationalmethodsforinverseproblems[M].Philadelphia:SIAM, 2002.

[22]RegińskaT.Aregularizationparameterindiscreteill-posedproblems[J].SIAMJournalonScientificComputing, 1996, 17(3): 740-749.

[23]袁亞湘, 孫文瑜. 最優(yōu)化理論與方法[M]. 北京; 科學(xué)出版社. 1997: 56-107.

[24]BonelliF,CutroneL,VottaR,etal.Preliminarydesignofahypersonicair-breathingvehicle[C]. 17thAIAAInternationalSpacePlanesandHypersonicSystemsandTechnologiesConference,SanFrancisco,California, 11-14April, 2011.