摘 要：連續(xù)突風(fēng)氣彈響應(yīng)在工程中一般采用隨機理論進(jìn)行分析求解，通過將全機頻響函數(shù)，與激勵功率譜密度積分，即可得到各自由度均方根值。該方法基于頻域線性求解，無法直接進(jìn)行非線性響應(yīng)分析。為此本研究基于匹配濾波器理論，提出了一套可考慮結(jié)構(gòu)非線性的連續(xù)突風(fēng)響應(yīng)工程簡化計算方法。通過引入突風(fēng)濾波器，串聯(lián)至全機氣彈動力學(xué)系統(tǒng)，構(gòu)成了連續(xù)突風(fēng)的時域分析模型。針對所選取自由度的均方根值，通過匹配濾波器理論，計算得到一套各自由度“時間匹配”的響應(yīng)，在此基礎(chǔ)上，針對非線性結(jié)構(gòu)模型，通過調(diào)整脈沖激勵幅值，尋找最大的響應(yīng)作為非線性結(jié)果，最終建立了一套適合于工程應(yīng)用的、考慮結(jié)構(gòu)非線性的連續(xù)突風(fēng)響應(yīng)計算方法。算例表明，本研究所建立的方法針對連續(xù)突風(fēng)，可獲得滿足工程精度要求的各自由度“時間匹配”的響應(yīng)，同時可用于結(jié)構(gòu)非線性響應(yīng)計算。

關(guān)鍵詞：匹配濾波器理論；隨機理論；功率譜密度；結(jié)構(gòu)非線性；時域

中圖分類號：V211.47" 文獻(xiàn)標(biāo)志碼：A

DOI：10.11776/j.issn.1000-4939.2024.01.004

A calculation method of continuous gust aeroelastic response considering structural nonlinearity

Abstract：The continuous gust aeroelastic response is generally analyzed and solved using random theory in engineering.By integrating the frequency response functions of the airplane with the power spectral density of excitation，the root mean square value of each degree of freedom can be obtained.The method is based on the linear solution in the frequency domain，and the nonlinear response analysis cannot be directly carried out.Therefore，based on the matched filter theory，a simplified calculation method for continuous gust response considering structure nonlinearity is described in this paper.By introducing a gust filter，which is connected in series to the aeroelastic dynamic system of the whole plane，a time-domain analysis model of continuous gust is obtained.Based on the root mean square value of the selected degrees of freedom，the matched filter theory is used to calculate a set of the “time correlated\" response for the other degree of freedom.On this basis，for the nonlinear structure model，by adjusting the impulse strength amplitude，the maximum response is found as the nonlinear result，and finally a nonlinear continuous gust response calculation method suitable for engineering and structural nonlinearity is established.The simulated results showed that the method established in this paper can be used to obtain the “time-correlated” response of each degree of freedom，and can also be used for structural nonlinear response calculation.

Key words：matched filter theory；random theory；power spectral density；structural nonlinear；time domain

隨著航空技術(shù)的發(fā)展，越來越多的大型客機裝備有偏航阻尼器、載荷減緩系統(tǒng)等裝置，如何合理考慮結(jié)構(gòu)、飛控系統(tǒng)等非線性因素（各類限制器及控制面偏轉(zhuǎn)間隙等）對結(jié)構(gòu)動力學(xué)的影響，一直是型號設(shè)計中的難題。尤其是對于航路中常常遇到的連續(xù)突風(fēng)/大氣紊流[1-2]，一個準(zhǔn)確合理的非線性系統(tǒng)建模，不僅可有效避免控制律不合理導(dǎo)致連續(xù)突風(fēng)響應(yīng)-載荷過大[3]，也可極大的加強適航當(dāng)局對各類非線性減緩系統(tǒng)的審查信心。

傳統(tǒng)連續(xù)突風(fēng)設(shè)計基于von Karman功率譜密度進(jìn)行頻域求解，考慮非線性后，線性疊加已不適用，因此難以在計算中考慮這些非線性環(huán)節(jié)。針對該問題，EASA在AMC25.341[4]中給出了若干非線性求解方案，總原則是要確保非線性分析與線性分析具備相同程度的安全性。因此一般要求非線性、線性分析具有類似的計算方法，并可提供一套匹配-平衡載荷，同時要求隨著非線性程度的降低，非線性結(jié)果應(yīng)逐漸逼近線性結(jié)果。

從上世紀(jì)80年代后期開始，國外陸續(xù)有相關(guān)研究機構(gòu)開展了大量的研究[5]，典型的有美國的FAA、NASA，歐洲的NLR等，這些方法大部分是基于隨機或確定性的時域分析來解決突風(fēng)非線性響應(yīng)求解的難題。主要有隨機模擬（stochastic simulation）[6-8]、譜突風(fēng)（spectral gust）[9]、統(tǒng)計方法（statistical method）[10-12]、匹配濾波（matched filter）[6，13]等。其中，匹配濾波作為一種典型的確定性時域分析方法，得到了廣泛的關(guān)注。匹配濾波器概念最初誕生于雷達(dá)信號處理[13]，其對已知的輸入信號（即匹配波形）可得到最大信噪比的輸出信號。將之應(yīng)用于突風(fēng)響應(yīng)計算，則問題可等效為尋找最嚴(yán)重的突風(fēng)波形，使得響應(yīng)最大?；谠摾碚摚琒COTT等[14]采用優(yōu)化的方法尋找合適的激勵波形，實現(xiàn)了非線性分析，但其未系統(tǒng)評估由于突風(fēng)濾波器的引入帶來的分析誤差。

鑒于此，本研究基于匹配濾波器理論，建立了一套可考慮結(jié)構(gòu)非線性的連續(xù)突風(fēng)氣彈響應(yīng)計算方法，系統(tǒng)評估了突風(fēng)濾波器對計算精度的影響，并開創(chuàng)性地應(yīng)用于考慮偏轉(zhuǎn)限制的、副翼連接故障場景下的連續(xù)突風(fēng)非線性模擬。首先針對條款規(guī)定的以功率譜密度表征的von Karman突風(fēng)譜，引入突風(fēng)濾波器進(jìn)行近似，其輸出作為全機氣動彈性動力學(xué)系統(tǒng)的輸入，串聯(lián)構(gòu)成了連續(xù)突風(fēng)的時域等效分析模型，為了得到目標(biāo)自由度某一確定的響應(yīng)值，通過匹配濾波器理論，反算推導(dǎo)對應(yīng)的匹配-激勵波形，將之帶入時域等效模型中計算得到一套各自由度“時間匹配”的響應(yīng)。在此基礎(chǔ)上，針對非線性模型，通過調(diào)整脈沖激勵幅值，尋找最大的響應(yīng)量作為非線性分析結(jié)果，最終建立了一套適合于工程應(yīng)用的考慮結(jié)構(gòu)非線性的連續(xù)突風(fēng)響應(yīng)計算方法。由于引入的突風(fēng)濾波器是誤差的主要來源，分別評估了2種突風(fēng)濾波器模型對計算精度的影響，結(jié)果顯示本研究建立的方法可有效地獲得滿足工程精度要求的響應(yīng)結(jié)果。最后，將副翼偏轉(zhuǎn)角度限制作為非線性環(huán)節(jié)，將所建立方法應(yīng)用于考慮副翼偏轉(zhuǎn)結(jié)構(gòu)限制的非線性連續(xù)突風(fēng)響應(yīng)模擬，并進(jìn)行了2種限制閾值的對比研究，進(jìn)一步驗證了方法的有效性。

1 計算方法

1.1 經(jīng)典連續(xù)突風(fēng)動響應(yīng)計算

在模態(tài)空間內(nèi)進(jìn)行突風(fēng)氣彈耦合計算[15-16]，控制方程為

［－Mω2+iCω+（1+ig）K－qQ］y=P（1）

式中：M、C、K分別為廣義質(zhì)量矩陣、阻尼矩陣和剛度矩陣；g為人工結(jié)構(gòu)阻尼；q為動壓；Q為廣義氣動力系數(shù)矩陣；P則為突風(fēng)引起的外部激勵，對于連續(xù)突風(fēng)來說，可表示為

P=qwgQ′w（2）

式中：Q′為突風(fēng)下洗相關(guān)廣義氣動力系數(shù)矩陣；w為突風(fēng)下洗矩陣；wg為突風(fēng)尺度因子。

最終得到氣彈系統(tǒng)各模態(tài)的頻響函數(shù)為

通過對頻響函數(shù)和激勵功率譜密度函數(shù)Φ（ω）進(jìn)行積分，即可得到各模態(tài)響應(yīng)y的響應(yīng)系數(shù)[17]。

1.2 匹配濾波器理論

基于氣彈系統(tǒng)各模態(tài)的頻響函數(shù)，針對輸入波形為x（t）的信號，可以通過傅里葉變換及逆變換，得到時域的響應(yīng)y（t），即

另外通過單位脈沖激勵δ（t）的響應(yīng)hy（t），時域響應(yīng)也可通過對其卷積計算得到，有

假設(shè)輸入信號x（t）滿足如下規(guī)律

x（t）=hy（－t+t0）/K（8）

式中，K為任意常數(shù)，則其頻域形式有

此y（t）可表示為

考慮到X*（ω）X（ω）是功率譜密度函數(shù)，用Sxx（ω）表示，由于積分函數(shù)是一個對于頻率的偶函數(shù)，所以可得

當(dāng)達(dá)到最大響應(yīng)，即t=t0時，

當(dāng)激勵波形x（t）的均方根值為1時，有

從而K值即等于時域響應(yīng)的最大值，且在t=t0時刻，有

ymax=y（t0）=K=σhy（14）

而對于均方根值等于1的任意形狀信號x′（t），其最大響應(yīng)值y′（t）不會超過σhy，證明如下。

根據(jù)施瓦爾茲不等式可得

將式（13）代入，則可得

y′（t）2≤σ2hy（17）

1.3 基于MFT的線性連續(xù)突風(fēng)計算

經(jīng)典連續(xù)突風(fēng)響應(yīng)計算基于隨機振動進(jìn)行線性求解，難以直接處理各類非線性如舵面偏轉(zhuǎn)角度及速率限制、載荷減緩等問題，一種工程處理方法是將隨機振動轉(zhuǎn)為確定性時域分析求解，為后續(xù)非線性處理建立基礎(chǔ)。

由于連續(xù)突風(fēng)一般以頻域功率譜密度的形式表征，典型的如von Karman譜，為了開展時域分析，可將之以傳遞函數(shù)的形式進(jìn)行近似，兩者關(guān)系為

Φ（ω）=G（ω）2（18）

式中：G（ω）為連續(xù)突風(fēng)的傳遞函數(shù)表征形式。

式中，Hy（ω）是氣動彈性系統(tǒng)對連續(xù)突風(fēng)的頻響函數(shù)。已知突風(fēng)濾波器G（ω），突風(fēng)形狀頻域成分可以表示為

Wg（ω）=G（ω）X（ω）（20）

應(yīng)用匹配濾波器理論，將X（ω）代入，有

從而得到使得設(shè)計量y數(shù)值最大的臨界突風(fēng)形狀，即

考慮到Φ（ω）=G（ω）G*（ω），臨界突風(fēng)形狀也可直接簡化采用Φ（ω）進(jìn)行計算，即

式（22）～（23）給出的突風(fēng)形狀可直接用于后續(xù)線性及非線性的相關(guān)計算中。

1.4 考慮結(jié)構(gòu)非線性的連續(xù)突風(fēng)計算

對于存在結(jié)構(gòu)非線性環(huán)節(jié)的連續(xù)突風(fēng)動態(tài)響應(yīng)計算，傳統(tǒng)線性疊加原理已無法適用，因此需在上述線性模型基礎(chǔ)上進(jìn)行拓展，一種可行的方法是如實的進(jìn)行時域響應(yīng)分析，通過調(diào)整脈沖激勵δ（t）的幅值k，反復(fù)進(jìn)行基于MFT的時域響應(yīng)計算，使得響應(yīng)量y的幅值最大，此時該幅值即為y的設(shè)計值，非線性計算流程如圖1所示。

2 計算模型

2.1 線性模型介紹

以GTA模型[18]進(jìn)行對比分析，采用梁單元、偶極子網(wǎng)格法建立突風(fēng)動響應(yīng)氣彈耦合分析模型，如圖2、圖3所示，副翼、方向舵和升降舵分別以單獨的梁單元和氣動面進(jìn)行模擬。各舵面與主翼面的連接剛度通過一維彈簧元模擬。

2.2 非線性環(huán)節(jié)

本研究在上述GTA模型的基礎(chǔ)上，通過降低左側(cè)副翼的連接剛度來近似模擬故障情況，同時對副翼偏轉(zhuǎn)角度的限制，作為本研究基于MFT方法進(jìn)行結(jié)構(gòu)非線性分析的算例。為了方便的開展非線性突風(fēng)響應(yīng)分析，采用了簡化分析方法，在傳統(tǒng)線性環(huán)節(jié)之后，串接上非線性環(huán)節(jié)，將結(jié)構(gòu)非線性分析作為線性分析的后處理。當(dāng)左側(cè)副翼超過一定閾值時，認(rèn)為已滿偏，對其大于閾值上的幅值進(jìn)行一定的折減，從而得到減緩系數(shù)，對于其他的變量值，則直接乘以該減緩系數(shù)。計算框圖如圖4所示。

2.3 突風(fēng)參數(shù)

針對適航條款[4]規(guī)定的von Karman譜來定義連續(xù)突風(fēng)形狀，其功率譜密度函數(shù)如式（24）給出，即

式中：σw為突風(fēng)速度均方根值；L為突風(fēng)尺度;Lv為突風(fēng)尺度除飛行速度。

為了對比分析突風(fēng)濾波器對MFT方法計算結(jié)果的影響，給出了2種對von Karman譜近似擬合的濾波器表達(dá)式。

其系數(shù)見表1。

2.4 模態(tài)參數(shù)

全機氣彈系統(tǒng)的頻響函數(shù)頻率計算間隔為0.1Hz，頻率截止到40Hz。采用單位模態(tài)質(zhì)量法對模態(tài)振型歸一化，并截取

前20階模態(tài)構(gòu)成模態(tài)空間進(jìn)行計算，依次為6階組合剛體模態(tài)、左側(cè)副翼近似剛體模態(tài)以及13階彈性模態(tài)，本研究以第1階彈性模態(tài)，即機翼對稱1階彎曲模態(tài)為基準(zhǔn)，采用MFT方法進(jìn)行計算得到與該模態(tài)均方根值匹配的激勵波形、時域突風(fēng)形狀，及對應(yīng)各模態(tài)的時域歷程，在此基礎(chǔ)上，開展非線性響應(yīng)分析。

3 計算結(jié)果

3.1 突風(fēng)濾波器評估

首先基于目前計算參數(shù)，對比評估2種濾波器對von Karman突風(fēng)模型的近似結(jié)果，均方根值積分范圍為0～40Hz，間隔0.1Hz，共401個頻率點，前20階模態(tài)響應(yīng)的均方根值見表2。從表中可看出，采用第二種近似系數(shù)誤差較大，對于本模型最大達(dá)到29.2%，而第一種近似方法，最大誤差僅為5.9%。

3.2 線性MFT方法評估

以第8階模態(tài)（機翼對稱1階彎曲）為基準(zhǔn)，開展線性MFT方法研究。其中t0時刻設(shè)置為3s，計算時長為8s，圖5給出基于2種濾波器計算得到的歸一化后的匹配濾波波形。從圖中可以看出，由于該波形是氣彈系統(tǒng)脈沖響應(yīng)時域歷程上的時間求逆，因此其初始時刻趨近于0，而隨著時間的推進(jìn)，逐漸振蕩增大直到峰值，且2種濾波器的峰值略有差異。

對于上述求得的匹配濾波波形，輸入至突風(fēng)濾波器即可得到滿足von Karman譜的突風(fēng)時域曲線，2者對比如圖6所示，具體峰值對比見表3，可以看出，基于第1種濾波器傳遞函數(shù)計算得到的突風(fēng)速度正向最大值相比第2種濾波器偏大。

基于圖5的濾波波形或基于圖6的突風(fēng)速度輸入，都可進(jìn)行時域仿真且可得到一致的結(jié)果，典型模態(tài)的響應(yīng)結(jié)果見圖7，其中第1階為剛體模態(tài)，為沉浮模態(tài)（向上為正）；第7階為近似剛體模態(tài)，對應(yīng)為左側(cè)副翼偏轉(zhuǎn)模態(tài)；第8階為機翼對稱1階彎曲。考慮到第1階模態(tài)是剛體模態(tài)，可以近似地看出飛機姿態(tài)的變化，在飛機進(jìn)入負(fù)向的突風(fēng)場后，飛機產(chǎn)生向下的運動，而在3s時刻前后，由于突風(fēng)速度瞬間從負(fù)向最小值穿越直至正向最大值，飛機逐漸開始向上運動，而對于第8階機翼彎曲模態(tài)來說，其為匹配波形輸入的基準(zhǔn)模態(tài)，嚴(yán)格滿足在t0=3s時刻達(dá)到最大值。

對于連續(xù)突風(fēng)響應(yīng)來說，由于其基于隨機理論分析開展，直接得到的是均方根值，用各模態(tài)的均方根值進(jìn)行評估過于保守，因此工程中常常采用相關(guān)系數(shù)[16]進(jìn)行合理的考慮其他自由度的響應(yīng)，而本研究采用的MFT方法屬于確定性的時域分析，各響應(yīng)量在物理時間上是嚴(yán)格對應(yīng)一致的，因此找出t0時刻的響應(yīng)即可。表4給出了3種模型的相關(guān)系數(shù)計算誤差分析，其中第二列為通過頻響函數(shù)直接求得的相關(guān)系數(shù)[16]，并以此為基準(zhǔn)，第三、四列為通過MFT方法計算得到的誤差，第五列PSD特指通過在式（23）的突風(fēng)激勵下計算得到的誤差，可以看出，2種濾波器是對von Karman突風(fēng)譜的近似，因此二者誤差整體都比PSD方法大，同時對于相關(guān)系數(shù)絕對值較小的模態(tài)，3種方法得到的誤差也相對較大。需要注意的是，對于PSD法，由于其缺乏對應(yīng)的匹配濾波波形，無法進(jìn)行后續(xù)非線性計算。

圖8給出了對應(yīng)誤差隨相關(guān)系數(shù)變化的關(guān)系圖，可直觀的看出，對于第1種突風(fēng)濾波器，其誤差隨相關(guān)系數(shù)增大而明顯減小，考慮到工程中通常只關(guān)注相關(guān)系數(shù)較大的自由度，因此該誤差趨勢可接受。

3.3 非線性計算評估

為驗證本研究非線性計算模型的合理性，對左副翼偏轉(zhuǎn)幅值進(jìn)行限制，即對第7階模態(tài)設(shè)置閾值，若其數(shù)值大于該值，則認(rèn)為其已飽和，繼續(xù)增加需疊加折減系數(shù)，為評估方法魯棒性，考慮到其最大值為0.005466，因此閾值分別設(shè)置為0.005和0.003，定義為輕微、嚴(yán)重非線性，折減系數(shù)統(tǒng)一設(shè)置為0.01，確保超過閾值后的超調(diào)量較小，突風(fēng)濾波器采用第1種形式。

圖9給出了2種閾值下的非線性模型與線性模型的對比圖，其中第7階模態(tài)位移為控制模態(tài)，第8階為被控制模態(tài)。從圖中可以看出，由于本研究采用的非線性模型僅在線性模型后串聯(lián)了非線性環(huán)節(jié)，因此2種模型在線性段基本一致，而在非線性段差別較大，而且隨著設(shè)置閾值增大，2種模型在非線性段的差異增大，符合預(yù)期。

圖10給出了2種閾值設(shè)置下，基于優(yōu)化算法得到的脈沖激勵幅值k與非線性響應(yīng)之間的迭代關(guān)系，從圖中可以看出，本研究建立的算法迭代次數(shù)較少，可高效地開展非線性連續(xù)突風(fēng)載荷分析。

4 結(jié) 論

本研究基于匹配濾波器理論，推導(dǎo)了適合連續(xù)突風(fēng)的線性時域分析模型，在此基礎(chǔ)上發(fā)展了一套適合結(jié)構(gòu)非線性分析的連續(xù)突風(fēng)響應(yīng)的計算方法，通過算例計算，得出以下結(jié)論。

1）本研究建立的連續(xù)突風(fēng)計算模型涉及到von Karman譜的濾波器近似，在結(jié)構(gòu)動力學(xué)響應(yīng)計算中存在一定的誤差，誤差可接受，且該誤差隨相關(guān)系數(shù)增大有減小趨勢。

2）本研究建立的計算模型可有效地開展考慮結(jié)構(gòu)非線性的連續(xù)突風(fēng)響應(yīng)分析，并提供一套可供強度分析的“時間匹配”響應(yīng)。

3）本研究建立的方法可直接用于帶非線性環(huán)節(jié)的突風(fēng)載荷減緩分析中。

參考文獻(xiàn)：

[1] HOBLIT F M.Gust loads on aircraft：Concepts and applications[M].Washington：AIAA，1988.

[2] FULLER J R.Evolution of airplane gust loads design requirements[J].Journal of aircraft，1995，32（2）：235-246.

[3] JONES J G.Nonlinear aircraft loads in severe atmospheric turbulence[J].Journal of aircraft，2009，46（5）：1627-1633.

[4] EASA.Certification specifications and acceptable means of compliance for large aeroplanes CS-25[M].[S.l.]：European Aviation Safety Agency，2017.

[5] SILVA W A，PERRY B III.Maximized gust loads of a closed-loop，nonlinear aeroelastic system using nonlinear systems theory[C]//40th Structures，Structural Dynamics，and Materials Conference and Exhibit.Reston，VA，USA：AIAA，1999：AIAA 1999-1474.

[6] SCOTT R C，POTOTZKY A S，PERRY B III.Matched-filter and stochastic-simulation-based methods of gust loads prediction[J].Journal of aircraft，1995，32（5）：1047-1055.

[7] SCOTT R C，POTOTZKY A S，PERRY B III.Further studies using matched filter theory and stochastic simulation for gust loads prediction[C]//34th Structures，Structural Dynamics and Materials Conference.Reston，VA，USA：AIAA，1993：AIAA 1993-1365.

[8] 劉瀅瀅，劉景光.連續(xù)突風(fēng)頻域響應(yīng)過程的時域等效模擬方法[J].裝備環(huán)境工程，2020，17（9）：21-25.

LIU Yingying，LIU Jingguang.Time domain equivalence simulation method in continuous gust frequency response process[J].Equipment environmental engineering，2020，17（9）：21-25（in Chinese）.

[9] LUSEBRINK H，BRINK-SPALINK J.Treatment of non-linear systems by timeplane-transformed CT methods-the spectral gust method[C]∥

Proceeding of AGARD Structures and Materials Panel Workshop Lillehammer：Canada Communication Group，1994：29-36.

[10]BARBATI M，RICCI S，CLIMENT H，et al.Continuous turbulence loads calculations in presence of nonlinear controls[C]//International Forum on Aeroelasticity and Structural Dynamics（IFASD 2013）.[S.l.]：Royal Aeronautical Society，2013：1-13.

[11]JONES J G.Statistical-discrete-gust method for predicting aircraft loads and dynamic response[J].Journal of aircraft，1989，26（4）：382-392.

[12]JONES J G，WATSON G H，F(xiàn)OSTER G W.Non-gaussian statistics of atmospheric turbulence and related effects on aircraft loads[J].AIAA journal，2004，42（12）：2438-2447.

[13]ZEILER T A.Matched filter concept and maximum gust loads[J].Journal of aircraft，1997，34（1）：101-108.

[14]SCOTT R C，POTOTZKY A S，PERRY B III.Maximized gust loads for a nonlinear airplane using matched filter theory and constrained optimization[C]//18th Atmospheric Flight Mechanics Conference.Reston，VA，USA：AIAA，1991：AIAA 1991-2848.

[15]楊秋艷，王育林.飛機連續(xù)陣風(fēng)載荷計算方法應(yīng)用研究[J].飛行力學(xué)，2008，26（4）：65-68.

YANG Qiuyan，WANG Yulin.Research and application of the continuous gust loads analysis method for aircraft[J].Flight dynamics，2008，26（4）：65-68（in Chinese）.

[16]嚴(yán)德，楊超，肖志鵬.彈性飛機平衡的陣風(fēng)外載荷計算與分析[J].北京航空航天大學(xué)學(xué)報，2012，38（10）：1321-1325.

YAN De，YANG Chao，XIAO Zhipeng.Balanced external gust loads computation and analysis for elastic aircraft[J].Journal of Beijing University of Aeronautics and Astronautics，2012，38（10）：1321-1325（in Chinese）.

[17]吳志剛，楊超.彈性導(dǎo)彈的連續(xù)與離散陣風(fēng)響應(yīng)[J].北京航空航天大學(xué)學(xué)報，2007，33（2）：136-140.

WU Zhigang，YANG Chao.Continuous and discrete gust responses of elastic missiles[J].Journal of Beijing University of Aeronautics and Astronautics，2007，33（2）：136-140（in Chinese）.

[18]KARPEL M，MOULIN B，CHEN P C，Dynamic response of aeroservoelatic systems to guest excitation［J］.Journal of aircraft，2005，42（5）：1264-1272.