脈沖燃燒風洞測力系統動態標定方法

武 龍, 王 鋒, 樂嘉陵

(中國空氣動力研究與發展中心 高超聲速沖壓發動機技術重點實驗室, 四川 綿陽 621000)

脈沖燃燒風洞測力系統動態標定方法

武 龍, 王 鋒*, 樂嘉陵

(中國空氣動力研究與發展中心 高超聲速沖壓發動機技術重點實驗室, 四川 綿陽 621000)

針對脈沖燃燒風洞中的測力系統，提出了一種動態標定方法。利用力錘在模型表面上不同位置，沿不同方向施加一系列集中載荷，由輸入載荷和天平輸出辨識出該表面對應的單位脈沖響應函數(UIRF)，再將各表面對應的UIRF加權得到系統的UIRF，加權系數由試驗狀態下各表面的壓力分布確定。辨識某表面對應的UIRF時，通過將其參數化使反卷積問題轉化為參數優化問題以回避問題的病態特性。求解參數優化問題時，先用遺傳算法搜索到參數全局最優解的近似值，再以此作為單純形方法的初值繼續優化得到參數最優值。在ANSYS中模擬了動態標定過程，考慮了實際試驗中輸出應變含有較大噪聲的情況，驗證了這種動態標定方法的準確性。

脈沖燃燒風洞；動態標定；參數辨識；遺傳算法；單純形方法

0 引 言

脈沖燃燒風洞的建造和運行成本較低，能夠以不同尺度、方式進行基本流動和工程問題的研究，是目前國內開展大尺度超燃沖壓發動機和機體推進一體化高超聲速飛行器試驗研究的主要設備之一[1]。然而，目前脈沖燃燒風洞的工作時間僅能達到300ms左右[2]，模型振動在試驗時間內來不及衰減，測力天平得到的信號含有較大的振動成分，傳統方法難以獲得良好的結果。如果能先獲得試驗系統的動態輸入輸出特性，就可以用天平的輸出信號反算出模型所受的載荷歷程而不需要模型達到穩定，這種載荷辨識方法已經在激波風洞內得到應用[3-4]。該方法能夠得到載荷的時變過程，這對于飛行器帶動力試驗[5-6]尤為重要。

載荷辨識的前提是已知系統的動態輸入輸出特性，動態標定試驗是獲得這一特性的可靠途徑。昆士蘭大學的Abdel-Jawad等對超高速膨脹管中的多分量應力波天平進行了動態標定，通過加載一系列標定載荷，得到了模型對分布載荷的脈沖響應矩陣[7]。牛津大學的L.J.Doherty等在Stalker管風洞中對一帶動力飛行器進行測力試驗時，用力錘敲擊方法對一臺三分量應力波天平進行了動態標定。中國空氣動力研究與發展中心的王鋒等采用瞬時卸載方式對系統進行激勵，分別用階躍函數和斜坡階躍函數描述卸載過程，對某單分量天平進行了動態標定。本文針對一個六分量測力試驗系統，給出了一種動態標定方法并對這種方法的準確性進行了驗證，為后續開展載荷辨識工作提供必要條件[8]。

1 測力試驗系統介紹

試驗系統包括蒙皮框架結構的試驗模型、天平以及固定底座3部分，如圖1所示。模型包含5個面，受載荷沖擊時通過框架將振動傳遞至天平浮動框，天平固定框與底座相連，底座固定于地面。

六分量天平由沿軸向布置的前后2個環形結構組合而成，如圖2所示，上部為浮動框，下部為固定框。浮動框與固定框之間是測量元件，其上貼有應變片，可測量應變并通過后續電路將應變信號輸出。

2 測力試驗系統模態分析

模態分析是結構動力學分析的基礎，系統在載荷作用下發生振動，可視為各階模態振動的疊加[9]。通過模態分析可近似了解系統的基本動力學特征，為開展動態標定提供重要的參考信息。利用有限元軟件計算得到系統的前6階模態，如表1所示。

表1 試驗系統的前6階模態Table 1 The first six order modes of the testing system

3 測力試驗系統動態標定方法

3.1 動態標定目的

動態標定的目的是獲得單位脈沖載荷作用下系統的響應函數(UIRF)。這個單位脈沖載荷是指由模型表面等效到天平中心點(或模型質心)后的單位脈沖載荷。試驗模型本身盡管整體剛度較大，但也是一個彈性體，在模型不同位置處加載對應的UIRF不同。這里近似認為模型同一個面上各位置對應的UIRF相同，因此可以在1個面上施加一系列集中載荷，通過輸入載荷和輸出應變求得該面對應的UIRF。在5個面上分別進行上述操作，然后將得到的5個UIRF加權即可得到整個試驗系統的UIRF。加權系數與試驗狀態下對應面上的平均壓力大小成正比[10]。各面上的平均壓力變化，加權系數要相應調整。仿真計算中，平均壓力由事先的CFD計算近似得到；風洞試驗時，在各面上布置測壓傳感器，每個時刻都可以測出1組壓力值進而算得1組加權系數，實際上同一車試驗中各面上平均壓力的比例關系隨時間變化不大，近似計算也可將不同時刻的加權系數取平均用于后續計算。

定義模型的5個表面為e1～e5，在模型的面e(e=1,2,3,4,5)上施加一時變載荷，其相對于天平中心有6個分量，即x、y、z方向的力和這3個方向的力矩，分別記為Fx(t)，Fy(t)，Fz(t)，Mx(t)，My(t)，Mz(t)。天平輸出Y有6個通道，將第k(k=1,2,…,6)個通道的輸出記為yk(t)。測力系統的輸入輸出關系為

*

其中，“*”表示卷積。第k行表示第k通道的輸入輸出關系：在面e上施加載荷，載荷相對天平中心分別為x、y、z方向單位脈沖力或力矩時，第k通道輸出分別為Ixek、Iyek、Izek、Imxek、Imyek、Imzek。令

則(1)式可寫為

Y=Ie*

Ie即為面e對應的單位脈沖響應函數，是一個6×6的矩陣，整個系統的單位脈沖響應函數I為，

其中，αe為面e對應的加權系數。對于某個時間點，I是一個6×6的矩陣，動態標定的目的即是獲得I。

3.2 動態標定過程

本文僅針對模型的1個表面，以天平的1個通道為例介紹動態標定方法，即只求出Ie的第k行：Iek。其他各面、各通道的動態標定方法是相同的。因此只需要考慮(1)式中的第k行

yk=Iek*

本文中僅針對e5面(e=5)，第1通道(k=1)求解對應的I51，為了簡潔，下文中略去下標e和k，將(5)式、(8)式分別簡寫為:

y=I*

下文中I均指代I51，對于某個時間點為一1×6向量，下面介紹I的求法。

動態標定采用力錘敲擊法，在e5面上安裝加載部件(圖1中有放大顯示)，為力錘敲擊提供5個不同方向的加載面。加載部件的加載面面積遠大于力錘錘頭的面積，這樣便于操作并有利于保證敲擊方向的準確性。定義加載方向為d1～d5，如圖3所示。

在點p(p=1,2,3,…)處沿方向d(d=1,2,3,4,5)敲擊，敲擊力要足夠大以保證天平有明顯輸出但又不能超出量程，力錘可記錄其輸入載荷大小的時間歷程f(t)，結合加載位置和方向計算出其相對于天平中心的載荷向量F

其中

cdp為單位載荷相對于天平中心的載荷矢量，其6個分量不獨立，由加載位置和方向決定。給定d、p后，可由簡單的幾何關系求得cdp，是已知的。

設輸出信號為ydp(t)，由(10)式

ydp=I*F=I*

定義

則

ydp=Idp*

這里，Idp即為模型表面加載點對應的脈沖響應函數，將在3.3節中給出其確定算法，此處暫認為已知。于是Idp、cdp已知，只有I待求。

由(9)式，I有6個未知分量，需要6個系數不相關的方程組成非齊次線性方程組才可解。因此，需要改變d、p進行6次不相關的敲擊。記敲擊次數為h(h=1,2,3,4,5,6)，并令:

由(14)，(16)和(17)式得到

因此，只需C可逆，即cdp_1,…,cdp_6線性無關，便可解得I

3.3Idp的求法

3.3.1 參數化Idp

設有n個時間點t1…tj…tn，那么就存在Idp(t1)…Idp(tj)…Idp(tn)共n個未知數，若直接將f和ydp(t)代入(15)式解卷積，問題的病態特征可能使結果誤差很大。

考慮結構是線彈性的，可將脈沖響應函數Idp(t)參數化

參數化將Idp寫成m個不同頻率振動疊加的形式。其中，m為選取的模態數量，Ai是第i個模態的幅值，ωi、γi分別是第i個模態的頻率和阻尼比。參數化后，未知數變為3m個，通常，n?m，未知數的個數大大減少。

設實際測量得到的輸出信號為yodp(t)，令R(t)=yodp(t)-ydp(t)，定義殘差r

求Idp(t)的問題轉化為用Idp*f擬合yodp(t)的問題，即確定能使r取到最小值的m及其對應的Ai、ωi、γi這一優化問題。

3.3.2 優化算法

分析3.3.1節中的這一優化問題，待定參數是m階模態參數(ω1,A1,γ1)…(ωm,Am,γm)，優化目標是使殘差r取最小值。優化過程中，r逐漸減小，給定一個值rl，以r

阻尼比γ一般在10-3量級或更小，在振動頻率不太高、測量時長較短的情況下對優化結果影響很小，故先不考慮阻尼比，對每個模態僅考慮幅值和頻率2個參數，Idp(t)簡化為

求解前m未知，可預估m值，然后對3m個參數同時優化。但這對預估的要求很高，m取小了導致Idp(t)中必然缺少高階振動成分；m取大了，則優化參數過多，優化難度較大。

這里設計了一個算法，將整個優化過程分若干次進行，每次應用一定的優化算法只辨識出1個模態對應的成分。第1次優化時，待擬合曲線為yodp(t)。第i次(i=1,2…)優化終止時，計算Ri(t)并輸出辨識出的第i階頻率ωi和幅值Ai。第i+1次優化時，以Ri(t)作為待擬合曲線，重復上一步的優化過程，得到ωi+1和Ai+1。因為每次優化的目標都是使殘差盡可能減小，被識別出來的都是當次待擬合曲線中最主要的振動成分，所以各振動成分是按照幅值由大到小的順序被辨識出來的，即Ai+1≥Ai。當進行到第M次優化后，AM?A1，主要的模態成分已經被辨識出來，若達到了r

在第i次優化中，如圖4虛線框中部分，未知數為ωi和Ai，待擬合曲線為Ri-1(t)(i=1時為yodp(t))，目標函數為ri(ωi,Ai)。考慮r(ωi,Ai)是一個多峰函數，使用傳統的搜索方法極易陷入局部最優解。遺傳算法是一種概率搜索算法，可以從一個種群開始并行搜索，并且不依賴于目標函數的梯度信息，具有很強的全局優化性能[11-12]，適合提取當前幅值最大的振動成分。這里采用Matlab的遺傳算法工具箱(GA)進行優化[13-14]。收斂準則為本代殘差r與上一代殘差rlg的相對差值小于ε，即

此處ε取10-4。求解前給定初代種群的取值范圍，求解中如果某一代滿足了收斂準則，就輸出結果；如果不滿足收斂準則，就通過選擇性復制、交叉、變異

生成下一代種群，直到滿足收斂準則為止。

在上述過程中，若rl取得較小則會出現ωi≈ωj,Ai?Aj，(i

前面已經得到了(ω1,A1)…(ωm,Am)和Idp(t)。現加入對阻尼比γ的考慮，若繼續基于上述的遺傳算法進行優化，減小ε和rl，是可以得到預期的Idp(t)的，但是求解效率很低。進一步優化時采用Matlab優化工具箱的fminsearch函數，它是基于單純形算法的優化工具。將前面得到的(ω1,A1)…(ωm,Am)作為頻率和幅值的初值，0作為阻尼比γ的初值。由于這些初值已經很接近最優解，fminsearch函數只需在最優解附近搜索即可，這就避免了其易于陷入局部最優解的缺點，而充分發揮了其收斂速度更快的優勢[15]。經過若干次迭代，殘差迅速下降并達到穩定，此時終止優化，輸出(ω1,A1,γ1)…(ωm,Am,γm)，并代入(20)式計算出Idp(t)。

4 測力試驗系統動態標定方法驗證

第2節中介紹了動態標定試驗的操作方法和求系統單位脈沖響應函數的算法。本節將利用有限元軟件模擬試驗過程，利用上述算法求得I，并檢驗I的準確性。

4.1 模擬動態標定試驗

在Ansys Workbench中進行瞬態動力學分析，模擬上述的動態標定過程。分別在點p1(d1,d2,d4方向)、點p2(d1,d2方向)、點p3(d1方向)模擬力錘加載，各加載點位置如圖3所示。由幾何關系分別算得這6次加載的載荷矢量：c11、c21、c41、c12、c22、c13，代入(17)式得到矩陣C

C的秩r(C)=6，載荷矢量滿足線性無關的要求。

載荷歷程為f1(t)

測得對應的輸出，記為yo11(t)、yo21(t)、yo41(t)、yo12(t)、yo22(t)、yo13(t)。

4.2 求I(t)

4.2.1 求6次敲擊分別對應的Idp(t)

利用3.3中的算法求這6次敲擊對應的Idp(t)，即I11(t)、I21(t)、I41(t)、I12(t)、I22(t)、I13(t)。

以I11(t)為例，利用遺傳算法得到前10個振動成分，如表2所示。

表2 前10個振動成分Table 2 The first ten frequencies and amplitudes

可見，幅值遞減，前5個幅值明顯大于后續幅值，符合預期。前10個頻率都在5.8Hz、9.4Hz、12.3Hz、16.9Hz、23.0Hz附近，這也與第2節中模態分析所得到的系統前幾階固有頻率相符。進一步利用fminsearch函數優化時只需考慮這5個振動成分，按照(24)、(25)式算得優化初值，如表3所示。

表3 優化初值Table 3 Initial values of optimization

由于初值已經比較準確，殘差r在前1000步迭代中便快速下降，經過約4000步穩定在10-9量級，得到模態參數，如表4所示。

表4 優化結果Table 4 Results of optimization

將表4中的模態參數代入(20)式，即可得到I11(t)。同理，可求得I21(t)、I41(t)、I12(t)、I22(t)、I13(t)。再按照(16)中的定義即可得到

4.2.2 求I(t)

將(26)、(28)式代入(19)式，即可求得I(t)

I(t)=S·

由上文，此處I(t)省略了下標e和k，其中

(29)式即為模型e5面上天平1通道對應的UIRF，其他各面、各通道對應的UIRF求法相同。

得到其他4個面上天平1通道對應的UIRF，均與(29)式相近，其中表征頻率和阻尼比的S相同。表征幅值的系數矩陣略有不同，但其中較大的系數相對變化很小，說明在不同面上敲擊激發的起主導作用的模態是相同的。由CFD計算算得5個面上的平均壓力進而算出各面對應的加權系數分別為：0.3865，0.1725，0.0062，0.1150，0.3197。由(7)式可得整個模型對天平1通道的UIRF,

I(t)=S·

測力系統在天平1通道的動態標定完成。

4.3 驗證I(t)的準確性

在p4點施加d5方向的載荷，載荷歷程為f2(t)

測得輸出應變yo54(t)，并利用4.2中得到的I(t)計算出理論輸出y54(t)。將yo54(t)、y54(t)繪制在圖5中。兩條曲線幾乎完全重合，說明動態標定得到的I(t)是準確的。

Fig.5 Comparison between computed strainy54(t)andmeasuredstrainyo54(t)

4.4 驗證標定方法在測得應變含噪聲時的準確性

實際試驗中，輸出應變中含有噪聲，這里驗證上述動態標定方法在ydp(t)含較大噪聲時的準確性。

以p1點處沿d1方向的敲擊為例，在yo11(t)中加入以該信號標準差15%為幅值的隨機信號，記為yn11(t)，模擬含有較大噪聲的輸出信號，局部圖如圖6所示。

利用3.3中的算法由yn11(t)計算出In11(t)，對比In11(t)和I11(t)。先利用遺傳算法得到In11(t)的前10階振動成分，發現與表2中的結果十分接近。然后利用fminsearch函數繼續優化，殘差r穩定在約1.77，得到In11(t)的各振動成分，如表5所示。

表5 優化結果Table 5 Results of optimization

表5與表4結果非常接近，只有低頻成分的阻尼比稍有差異。這是因為本節為了驗證算法，加入的噪聲很大，低頻阻尼比的影響遠遠小于噪聲的影響。實際試驗中噪聲并不會這么大，低頻阻尼比的誤差會更小。繪制I11(t)、In11(t)，局部圖如圖7所示。兩條曲線幾乎重合，說明即使yn11(t)中含有明顯的噪聲，利用3.3中的算法，得到的In11(t)仍然準確。

可見，yn11(t)中的噪聲基本沒有被擬合，而是保留在了殘差r中。

5 結 論

本文給出了一種針對脈沖燃燒風洞測力系統的動態標定方法并利用ANSYS仿真對其進行了驗證，得到了以下結論：

(1) 求解UIRF時通過參數化將解卷積問題轉化為參數優化問題，可以有效回避問題的病態特性。

(2) 求解參數優化問題時先利用遺傳算法搜索到全局最優解的近似值，再以其作為單純形方法的初值繼續優化，可以迅速求得全局最優解。

(3) 即使實際試驗測得的應變信號含有較大的噪聲，這種動態標定方法仍具有很高的精度。

(4) 本動態標定方法在實際應用中可能存在的誤差主要由力錘敲擊方向的偏差引起，試驗人員在操作前需要進行適當的練習以盡量減小敲擊方向的偏差。

[1]樂嘉陵, 劉偉雄, 賀偉, 等. 脈沖燃燒風洞及其在火箭和超燃發動機研究中的應用[J]. 實驗流體力學, 2005, 19(1): 1-10.

Le J L, Liu W X, He W, et al. Impulse combustion wind tunnel and its application in rocket and scramjet research[J]. Journal of Experiments in Fluid Mechanics, 2005, 19(1): 1-10.

[2]劉偉雄, 譚宇, 毛雄兵, 等. 一種新運行方式脈沖燃燒風洞研制及初步應用[J]. 實驗流體力學, 2007, 21(4): 59-64.

Liu W X, Tan Y, Mao X B, et al. The development and preliminary application of a pulse combustion wind tunnel with new running way[J]. Journal of Experiments in Fluid Mechanics, 2007, 21(4): 59-64.

[3]Robinson M J, Mee D J, Tsai C Y, et al. Three-component force measurements on a large scramjet in a shock tunnel[J]. Journal of Spacecraft and Rockets, 2004, 41(3): 416-425.

[4]Robinson M J, Hannemann K, Schramm J M. Design and implementation of an internal stress wave force balance in a shock tunnel[J]. CEAS Space Journal, 2011, 1(1): 45-57.

[5]賀偉, 于時恩, 李宏斌. 高超聲速一體化飛行器推阻特性測量研究[J]. 實驗流體力學, 2010, 24(2): 65-68.

He W, Yu S E, Li H B. Experimental investigation on thrust drag performance of hypersonic integrative vehicle[J]. Journal of Experiments in Fluid Mechanics, 2010, 24(2): 65-68.

[6]賀偉, 童澤潤, 李宏斌. 單模塊超燃發動機推力測量天平研制[J]. 航空動力學報, 2010, 25(10): 2285-2289.

He W, Tong Z R, Li H B. Investigation of thrust balance for the single module scramjet[J]. Journal of Aerospace Power, 2010, 25(10): 2285-2289.

[7]Abdel-Jawad M M, Mee D J, Morgan R G. New calibration technique for multiple-component stress wave force balances[J]. Review of Scientific Instruments, 2007, 78(6): 065101-1-065101-7.

[8]王鋒, 任虎, 周正, 等. 載荷辨識方法用于脈沖風洞模型阻力測量研究[J]. 振動與沖擊, 2015, 34(24): 202-208.

Wang F, Ren H, Zhou Z, et al. Drag force measurement in impulse facilities by using load identification method[J]. Journal of Vibration and Shock, 2015, 34(24): 202-208.

[9]李東旭. 高等結構動力學[M]. 北京: 科學出版社, 2010: 273-298.

Li D X. Advanced structural dynamics[M]. Beijing: The Science Publishing Company, 2010: 273-298.

[10]Doherty L J, Smart M K, Mee D J. Measurement of three-components of force on an airframe integrated scramjet at Mach 10[R]. AIAA-2015-3523, 2015.

[11]劉國春, 費強, 趙武云, 等. 基于Matlab 遺傳算法優化工具箱的應用[J]. 機械研究與應用, 2014, 27(2): 71-73.

Liu G C, Fei Q, Zhao W Y, et al. Application of genetic algorithm optimization toolbox based on matlab[J]. Mechanical Research and Application, 2014, 27(2): 71-73.

[12]林鴻彬. 基于遺傳算法的數據擬合在 MATLAB 環境中的實現[J]. 湖南農機, 2010, 37(3): 92-97.

Lin H B. Data fitting based on genetic algorithm implementation in MATLAB environment[J]. Hunan Agricultural Machinery, 2010, 37(3): 92-97.

[13]羅述全. 傳統優化算法與遺傳算法的比較[J]. 湖北工業大學學報, 2007, 22(3): 32-35.

Luo S Q. Comparison between traditional optimized algorithm and heredity algorithm[J]. Hubei University of Technology Journal, 2007, 22(3): 32-35.

[14]郭海雙, 梁佳雯, 張劭昀. MATLAB 遺傳算法工具箱GADS優化及應用[J]. 電子設計工程, 2015, 23(10): 27-30.

Guo H S, Liang J W, Zhang S Y. Optimization and examples in Matlab GA toolbox GADS[J]. Electronic Design Engineering, 2015, 23(10): 27-30.

[15]楊改強, 霍麗娟, 楊國義, 等. 利用MATLAB擬合van Genuchten方程參數的研究[J]. 土壤, 2010, 42(2): 268-274.

Yang G Q, Huo L J, Yang G Y, et al. Research on fitting van genuchten equation parameter with MATLAB software[J]. Soils, 2010, 42(2): 268-274.

(編輯：張巧蕓)

A dynamic calibration method for a dynamometric systemin impulse combustion facilities

Wu Long, Wang Feng*, Le Jialing

(Science and Technology on Scramjet Laboratory, China Aerodynamics Research and Development Center, Mianyang Sichuan 621000, China)

A new dynamic calibration method for a dynamometric system in impulse combustion facilities is proposed. The calibration involved uses an instrumented impact hammer to apply individual calibration forces in different directions at different positions on a face of the model and calculates the unit impulse response function (UIRF) of the face from input loads and output strains. UIRFs of different faces are weighted to obtain the UIRF of the dynamometric system and the weighting coefficients are determined by the pressure on each face under the experimental condition. By parameterization, the problem is converted into a parameter optimization problem to solve the UIRF. Using a genetic algorithm to obtain the approximation of the global optimal solution of the parameters and setting it as the initial value of a simplex algorithm, the exact solution is obtained by the simplex algorithm then. ANSYS simulation of the dynamic calibration is presented. Input loads and output strains are recorded and noises are added to the output strains to simulate the actual experimental situation. The simulation validates the accuracy and feasibility of the dynamic calibration method.

impulse facilities;dynamic calibration;parameter identification;genetic algorithm;simplex algorithm

1672-9897(2017)04-0051-08

10.11729/syltlx20160158

2016-10-21;

2017-01-09

國家自然科學基金項目(11372339)；高超聲速沖壓發動機技術重點實驗室基金項目(STSKFKT2012001)

WuL,WangF,LeJL.Adynamiccalibrationmethodforadynamometricsysteminimpulsecombustionfacilities.JournalofExperimentsinFluidMechanics, 2017, 31(4): 51-58. 武 龍, 王 鋒, 樂嘉陵. 脈沖燃燒風洞測力系統動態標定方法. 實驗流體力學, 2017, 31(4): 51-58.

V211.7

A

武 龍(1991-)，男，黑龍江省大慶市人，碩士研究生。研究方向：載荷辨識。通信地址：四川省綿陽市涪城區劍門路西段278號(621000)。E-mail: 779483196@qq.com。

*通信作者 E-mail: wfscholar@163.com