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

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

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

(1.中國空氣動(dòng)力研究與發(fā)展中心吸氣式高超聲速技術(shù)研究中心，綿陽6210002.中國空氣動(dòng)力研究與發(fā)展中心高超聲速?zèng)_壓發(fā)動(dòng)機(jī)技術(shù)重點(diǎn)實(shí)驗(yàn)室，綿陽621000)

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

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

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

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

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

1試驗(yàn)系統(tǒng)的動(dòng)力學(xué)建模

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

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

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

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

1.2建模方法

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

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

1.2.1試驗(yàn)飛行器的模型

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

(1)

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

(2)

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

(3)

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

(4)

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

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

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

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

點(diǎn)Ⅱ、點(diǎn)Ⅲ的位移向量分別記作

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

(5)

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

因?yàn)閯偠染仃囀菍?duì)稱的，所以有

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

(13)

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

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

(14)

因?yàn)辄c(diǎn)Ⅲ不受力，將ΔRⅡ、ΔRⅢ代入式(5)的第二組方程，有

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

(15)

再利用式(14)可得

(16)

KⅢ,Ⅲ=TB1KⅡ,ⅡTTB1

(17)

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

v=NB,ⅡfⅡ

(18)

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

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

(19)

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

v=SⅡηⅡ+SⅢηⅢ

(20)

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

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

(21)

所以有

(22)

代入式(20)，有

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

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

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

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

(24)

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

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

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

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

圖3　點(diǎn)Ⅰ和點(diǎn)Ⅱ的相對(duì)位置 Fig.3 Position relationship of point Ⅰ and point Ⅱ

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

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

ηⅡ=

(25)

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

(26)

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

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

(27)

其中

K=

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

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

(28)

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

(29)

2載荷辨識(shí)模型

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

(30)

y=Rx

(31)

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

3載荷辨識(shí)方法

3.1子空間投影法

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

(32)

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

(33)

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

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

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

令

(34)

(35)

其等價(jià)問題為

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

(36)

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

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

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

w0=0，β1z1=y

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

αiwi=RTzi-βiwi-1

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

式中，實(shí)數(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個(gè)向量。于是，式(38)進(jìn)一步簡化為

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

(39)

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

(40)

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

(41)

(42)

式中zα是如下問題的解

(43)

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

(44)

(45)

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

zα≈Wkζ

(46)

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

(47)

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

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

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

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

(48)

3.4求解算法

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

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

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

(3)針對(duì)每個(gè)αi求解方程(39)，由(44)、(45)式得到‖xα‖2、‖rα‖2，計(jì)算指標(biāo)函數(shù)值序列J(αi,k)，找出最小值對(duì)應(yīng)的序號(hào)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，擴(kuò)充子空間，轉(zhuǎn)到第3步。

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

4數(shù)值算例

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

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

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

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

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

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

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

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

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

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

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

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

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

4.3載荷辨識(shí)結(jié)果

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

(49)

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

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

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

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

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

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

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

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

5結(jié)論

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

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

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

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

參考文獻(xiàn)

[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]王林軍. 正則化方法及其在動(dòng)態(tài)載荷辨識(shí)中的應(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ù)精細(xì)積分法的載荷識(shí)別研究[J]. 振動(dòng)與沖擊, 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]韓旭, 劉杰, 李偉杰, 等. 時(shí)域內(nèi)多源動(dòng)態(tài)載荷的一種計(jì)算反求技術(shù)[J]. 力學(xué)學(xué)報(bào), 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.