邢譽(yù)峰，梁昆

(北京航空航天大學(xué) 固體力學(xué)研究所，北京100191)

飛行器結(jié)構(gòu)物理參數(shù)的主要特點(diǎn)之一是其具有時(shí)變質(zhì)量.運(yùn)載火箭除了具有時(shí)變質(zhì)量參數(shù)之外，其幾何特征是其具有類似細(xì)長梁的構(gòu)形，其柔性特征致使其彈性變形非常重要，無論對(duì)載荷分析還是飛行控制皆如此.

邢譽(yù)峰等在文獻(xiàn)［1-2］中，對(duì)過載火箭結(jié)構(gòu)系統(tǒng)的動(dòng)力學(xué)特性和變質(zhì)量火箭系統(tǒng)的分析方法進(jìn)行了研究.值得指出的是，運(yùn)載火箭尤其是捆綁運(yùn)載火箭的縱橫耦合特性是不可忽略的，這主要是由于火箭具有柔性結(jié)構(gòu)，并且同時(shí)存在縱橫載荷(如推力和風(fēng)載等).當(dāng)彈性變形不再是小變形時(shí)，系統(tǒng)將出現(xiàn)強(qiáng)耦合現(xiàn)象，否則也存在弱耦合現(xiàn)象，如推力致使橫向剛度降低.基于此，本文開展了相關(guān)的研究工作.

1995年，文獻(xiàn)［3］使用諧波增量平衡法求解了文獻(xiàn)［4-5］給出的縱橫耦合梁的主共振和倍頻主共振響應(yīng).1999年，文獻(xiàn)［6-7］推導(dǎo)了等效柔性梁結(jié)構(gòu)縱橫耦合非線性動(dòng)力學(xué)方程，并分析了自由振動(dòng)和受迫振動(dòng)響應(yīng)特性.2010年，文獻(xiàn)［8］推導(dǎo)了梁縱橫向耦合振動(dòng)的非線性單元?jiǎng)偠染仃?，用有限元方法求解了?dòng)態(tài)響應(yīng)，并分析了耦合響應(yīng)中的倍頻現(xiàn)象等.文獻(xiàn)［9］利用Galerkin法和諧波增量平衡法研究縱向運(yùn)動(dòng)梁在縱橫耦合情況下的自由振動(dòng)響應(yīng)，尤其是在橫向第1階和第2階固有頻率之比接近1∶3內(nèi)共振條件下的系統(tǒng)響應(yīng).文獻(xiàn)［10］研究了水中塔架等效梁模型的耦合振動(dòng).文獻(xiàn)［11］用模擬和試驗(yàn)方法研究了縱橫耦合振動(dòng).有些學(xué)者則研究縱向激勵(lì)對(duì)橫向振動(dòng)的影響，如文獻(xiàn)［12］.

已有文獻(xiàn)鮮有研究非線性系統(tǒng)的頻率特性和響應(yīng)頻率特性，如頻率的時(shí)變特性、系統(tǒng)頻率和響應(yīng)頻率的關(guān)系等.在這種情況下，本文開展了有關(guān)工作.除了研究頻率特性之外，還考慮了共振現(xiàn)象和橫向振動(dòng)對(duì)縱橫振動(dòng)影響的特性等.

1 基本方程

對(duì)于均勻Rayleigh梁，如圖1所示，若考慮縱橫耦合，其幾何關(guān)系為［5］

式中：ξ為非線性或耦合作用因子，ξ=1對(duì)應(yīng)幾何關(guān)系的原來形式，ξ=0對(duì)應(yīng)線性系統(tǒng)，其他情況對(duì)應(yīng)作者考慮的其他情況如弱耦合情況等;εij(i，j=x，y，z)為6 個(gè)應(yīng)變分量;位移(u1，u2，u3)為與坐標(biāo)系(x，y，z)對(duì)應(yīng)的位移場函數(shù)，其表達(dá)式為

式中：w為橫向位移;y為厚度坐標(biāo);t為時(shí)間.將式(2)代入式(1)得

系統(tǒng)應(yīng)變能函數(shù)Uε和動(dòng)能函數(shù)T分別為

式中：L為梁長;E為彈性模量;A為橫截面積;ρ為體密度;EI為彎曲剛度;u為梁的縱向位移.

動(dòng)能函數(shù)中忽略了耦合項(xiàng).根據(jù)Hamilton原理得到Rayleigh梁縱橫耦合模型振動(dòng)控制方程

和邊界條件，見表1.式(6)中Fu和Fw為縱向和橫向載荷.

圖1 算例簡支梁示意圖Fig.1 Schematic diagram of simply supported beam example

表1 縱橫耦合梁邊界條件Table 1 Boundary conditions of longitudinal and transverse coupled beam

從方程(6)，初步得到如下結(jié)論：

1)若考慮非線性縱橫耦合，從方程(6a)右端項(xiàng)可以看出，縱向振動(dòng)中除了包括系統(tǒng)縱向頻率和激勵(lì)頻率信息之外，還包括橫向振動(dòng)頻率和激勵(lì)頻率的各種組合和倍頻.

2)從方程(6b)右端項(xiàng)可以看出，橫向振動(dòng)中除了包括系統(tǒng)橫向頻率的各種組合和倍頻信息之外，還包括系統(tǒng)縱橫各個(gè)頻率和激勵(lì)頻率的各種組合.

3)若沒有橫向激勵(lì)或橫向初始擾動(dòng)，縱向振動(dòng)將無法引起橫向振動(dòng)，反之不然.

在下面算例中，還將詳細(xì)討論上述結(jié)論.

2 有限元求解

2.1 位移函數(shù)

下面分別用線性拉格朗日函數(shù)和三次埃爾米特函數(shù)表達(dá)梁的縱向位移u和橫向位移場w，其有限元離散形式為［13］

式中：

式中：l為單元的長度.

2.2 結(jié)構(gòu)單元矩陣

把位移函數(shù)式(7)代入能量泛函，可以得到用結(jié)點(diǎn)位移向量表示的應(yīng)變能函數(shù)Uε和動(dòng)能函數(shù)T的標(biāo)準(zhǔn)形式，即

式中單元矩陣為

其中：

由結(jié)構(gòu)單元矩陣可以看出，單元質(zhì)量矩陣與不耦合情況下單元質(zhì)量矩陣相同，這是因?yàn)閯?dòng)能函數(shù)中沒有考慮耦合項(xiàng)，見式(5);k1和k4為通常的等應(yīng)變桿單元和三次梁單元的剛度矩陣;k2和k3為梁單元縱橫耦合剛度矩陣，其元素與振動(dòng)狀態(tài)直接相關(guān)，即剛度矩陣是時(shí)變的.在形成這兩個(gè)耦合矩陣時(shí)，本文利用前一時(shí)刻的位移向量來計(jì)算當(dāng)前時(shí)刻的剛度矩陣.在本文算例中，針對(duì)每個(gè)時(shí)間步，沒有進(jìn)行迭代.當(dāng)時(shí)間步長較小時(shí)，所得結(jié)果的精度是可以保證的.

把結(jié)構(gòu)單元組裝后，根據(jù)變分原理可得系統(tǒng)的離散形式動(dòng)力學(xué)方程：

式中：M、K和F分別為結(jié)構(gòu)的總體質(zhì)量矩陣、剛度矩陣和載荷列向量.

3 數(shù)值模擬和分析

這里所用模型的主要參數(shù)來自文獻(xiàn)［7］，表2給出了幾何模型和材料參數(shù)、有限元模型參數(shù)和邊界條件，表3給出了激勵(lì)載荷.本文作者用MATLAB［14］自編了動(dòng)力學(xué)響應(yīng)分析程序，其中用的是Newmark直接積分法.時(shí)域范圍為0～1 s，時(shí)間步長h=0.1 ms.單元類型是三次梁單元，其縱向變形是線性的，3個(gè)結(jié)點(diǎn)參數(shù)包括縱向位移、橫向位移及其導(dǎo)數(shù)(斜率).激勵(lì)作用點(diǎn)和動(dòng)態(tài)響應(yīng)觀測點(diǎn)為圖1的C點(diǎn)，為距左端距離為L/10的結(jié)點(diǎn)，見圖1.如果沒有特殊說明，非線性影響系數(shù)ξ=1.

表2 模型參數(shù)Table 2 Model parameters

表3 激勵(lì)載荷Table 3 External forces

3.1 模型的正確性

為了驗(yàn)證本文工作的正確性，下面把作者自編MATLAB程序的計(jì)算結(jié)果與NASTRAN結(jié)果進(jìn)行對(duì)比.NASTRAN結(jié)果中用的是Euler梁單元(不考慮剪切變形和截面轉(zhuǎn)動(dòng)慣量).為了對(duì)比，此處MATLAB質(zhì)量矩陣中不考慮截面轉(zhuǎn)動(dòng)量的影響.

表4給出了橫向和縱向前9階線性系統(tǒng)固有振動(dòng)頻率的比較，二者與理論結(jié)果吻合.

表4 前9階線性系統(tǒng)固有振動(dòng)頻率比較Table 4 Comparisons of the first 9 frequencies Hz

用表3中Case1給出的激勵(lì)，不考慮初始條件，圖2給出了兩種方法得到的縱橫動(dòng)態(tài)響應(yīng)比較.NASTRAN用非線性瞬態(tài)響應(yīng)分析求解器SOL129;MATLAB用Newmark方法進(jìn)行計(jì)算，二者用相同的時(shí)間步長.從圖2可以看出，兩種方法響應(yīng)結(jié)果也是吻合的，尤其是橫向振動(dòng)響應(yīng).

圖2 MATLAB與NASTRAN計(jì)算的動(dòng)態(tài)響應(yīng)對(duì)比Fig.2 Comparison of dynamic responses calculated by MATLAB and NASTRAN

通過圖2和表4的比較可以得到結(jié)論，本文推導(dǎo)的公式和程序?qū)崿F(xiàn)是正確的，為下文對(duì)梁的縱橫耦合非線性特性的分析奠定了基礎(chǔ).

3.2 耦合系統(tǒng)的頻率特性

無論對(duì)于梁的線性縱向還是線性橫向振動(dòng)，參見ξ=0時(shí)的方程(6)，其響應(yīng)都是由簡諧響應(yīng)疊加而成，其頻率成分包括系統(tǒng)的固有振動(dòng)頻率和激勵(lì)頻率.對(duì)于非線性控制方程(6)，雖然其動(dòng)態(tài)響應(yīng)特性遠(yuǎn)比線性系統(tǒng)復(fù)雜，但也是由各階諧波疊加組成，諧波頻率成分可以參見在第2節(jié)中總結(jié)的結(jié)論：梁非線性系統(tǒng)動(dòng)態(tài)響應(yīng)的頻率包括系統(tǒng)縱向、橫向頻率和激勵(lì)頻率的各種組合和倍頻.

除此之外，從非線性系統(tǒng)的剛度矩陣子陣k2和k3可以看出，剛度矩陣是時(shí)變的.對(duì)于任意微小的時(shí)間段，可以假設(shè)其頻率和模態(tài)不隨著時(shí)間變化，通過式(12)可以求得對(duì)應(yīng)小時(shí)間段的頻率和模態(tài).

由此可以看出，對(duì)于本文研究的幾何非線性動(dòng)力學(xué)系統(tǒng)，系統(tǒng)頻率是時(shí)變的，隨著非線性程度降低，其頻率成分應(yīng)該逐漸趨近于對(duì)應(yīng)線性系統(tǒng)的頻率，參見圖3，其所用參數(shù)與圖2所用參數(shù)相同.對(duì)于非線性系統(tǒng)，其頻率是與響應(yīng)幅值相關(guān)的［15］.值得指出的是，對(duì)于變系數(shù)線性動(dòng)力學(xué)系統(tǒng)，例如變質(zhì)量系統(tǒng)，系統(tǒng)的頻率也是時(shí)變的;對(duì)于阻尼非線性動(dòng)力學(xué)系統(tǒng)，若阻尼較小，其對(duì)頻率的影響是可以忽略的，此時(shí)可以認(rèn)為系統(tǒng)頻率不是時(shí)變的.

圖3 非線性系統(tǒng)頻率的時(shí)變特性Fig.3 Time-variable system frequency characteristics of nonlinear system

為了進(jìn)一步考慮非線性系統(tǒng)動(dòng)態(tài)響應(yīng)頻率成分的特性，考慮如下兩種情況.

3.2.1 受迫振動(dòng)頻率耦合

系統(tǒng)初始條件為零，激勵(lì)形式為表3的Case1，但激勵(lì)頻率更換為fu=17 Hz，fw=11 Hz.

圖4給出了該受迫振動(dòng)情況的幅頻特性，從中可以看出：

1)在縱向幅頻特性中，17 Hz為縱向激勵(lì)頻率;14 Hz和36 Hz為橫向基頻與橫向激勵(lì)頻率之代數(shù)和(25±11)Hz;22 Hz為橫向激勵(lì)頻率的2倍(11×2);50 Hz為橫向基頻的 2倍(2×25)Hz;5 Hz和39 Hz分別為橫向激勵(lì)頻率的2倍與縱向激勵(lì)頻率之代數(shù)和(2×11±17)Hz;19 Hz、3 Hz和53 Hz分別為橫向基頻與縱橫激勵(lì)頻率之代數(shù)和(25+11-17、11+17±25)Hz.

2)在橫向幅頻特性中，11 Hz為橫向激勵(lì)頻率;25Hz為橫向第1階固有頻率;8 Hz和42 Hz為橫向基頻與縱向激勵(lì)力之代數(shù)和(25±17)Hz;6Hz和28Hz為縱橫激勵(lì)頻率代數(shù)和(17±11)Hz.

3.2.2 自由振動(dòng)頻率耦合

不考慮激勵(lì)，僅考慮表5給出的橫向初始位移條件，該初始條件是由橫向集中力P=4.5422×104N作用在梁上B點(diǎn)得到的，B點(diǎn)距離左端的距離為9L/10.

圖5給出了幅頻特性曲線.橫向幅頻特性中，主要頻率為前4階系統(tǒng)橫向頻率，其中并沒有包含系統(tǒng)縱向頻率成分，這是由于縱向頻率比較高的緣故.

圖4 受迫振動(dòng)的幅頻特性Fig.4 Amplitude-frequency characteristics of forced vibration

表5 橫向初始條件Table 5 Initial conditions of transverse vibration

縱向頻譜特性中，50 Hz、202 Hz、454 Hz 分別為系統(tǒng)橫向頻率的二倍頻;76 Hz、302 Hz為橫向固有頻率的三倍頻;126 Hz、378 Hz為兩個(gè)橫向頻率之差;252Hz、430Hz、328Hz、504Hz為兩個(gè)橫向頻率之和;600 Hz、650 Hz為3個(gè)橫向頻率之代數(shù)和(228+402 ±25)Hz.76 Hz、126 Hz還可以分別表示為前2階橫向頻率之差與之和、202 Hz還可以表示為第1和第3階橫向頻率之差、302 Hz還可以表示為第2和第4階橫向頻率之差.

3.3 縱橫相互影響特征

1)首先考慮表3中Case2的激勵(lì)，初始條件為零.圖6給出了縱橫動(dòng)態(tài)響應(yīng)，從圖中可以看出，如果沒有橫向激勵(lì)或橫向初始條件，縱向振動(dòng)不能引起橫向振動(dòng)，參見式(6)，這是因?yàn)榇藭r(shí)w'和w″始終為零.若沒有縱向激勵(lì)和縱向初始條件，橫向振動(dòng)卻可以引起縱向振動(dòng).

2)考慮表3中Case1中的橫向激勵(lì)，縱向激勵(lì)為0，同時(shí)考慮表5給出的橫向初始條件.圖7給出了動(dòng)態(tài)響應(yīng)結(jié)果，從圖中可以看出，縱向振動(dòng)對(duì)橫向振動(dòng)的影響較小，橫向振動(dòng)對(duì)縱向振動(dòng)的影響比較大.

圖5 自由振動(dòng)頻譜Fig.5 Amplitude-frequency characteristics of free vibration

圖6 橫向初始條件引起的振動(dòng)響應(yīng)Fig.6 Vibration responses with respect to transverse initial conditions

圖7 考慮橫向初始條件和激勵(lì)載荷時(shí)縱向、橫向振動(dòng)響應(yīng)Fig.7 Vibration response in both longitudinal and transverse directions considering transverse initial condition and excitation force

3)橫向激勵(lì)作用下縱向振動(dòng)響應(yīng)特性.

縱向激勵(lì)力為零，不考慮初始條件，分別考慮橫向激勵(lì)力幅值分別為Fw1=20 kN和Fw2=2Fw1作用時(shí)縱向動(dòng)態(tài)響應(yīng)幅值的變化規(guī)律，激勵(lì)頻率為11 Hz.圖8和表6給出了分析結(jié)果，從中可以看到兩種情況下，F(xiàn)w2作用下的縱向振動(dòng)幅值基本為Fw1作用下幅值的4倍.從方程(6a)右端第2項(xiàng)可以看出此結(jié)論的合理性.表6中14Hz近似為一階固有頻率與激勵(lì)頻率之差，22 Hz激勵(lì)頻率2倍，36Hz近似為一階固有頻率與激勵(lì)頻率之和.

圖8 橫向激勵(lì)載荷幅值加倍后縱向幅頻特性Fig.8 Longitudinal amplitude-frequency characteristics when transverse excitation forces double

表6 橫向激勵(lì)載荷加倍時(shí)縱向振動(dòng)幅頻特性Table 6 Longitudinal amplitude-frequencycharacteristics when transverse excitation forces double

3.4 共振現(xiàn)象

考慮受迫振動(dòng)情況，縱向激勵(lì)頻率fu=15 Hz，橫向激勵(lì)頻率fw=10.3125Hz，二者之和與線性情況橫向第1階固有振動(dòng)頻率相等.

由圖9可知，對(duì)于這種情況，非線性橫向和縱向振動(dòng)響應(yīng)都明顯增大.但由于非線性系統(tǒng)頻率是時(shí)變的，外部激勵(lì)頻率難與非線性系統(tǒng)頻率相等，因此也就不會(huì)出現(xiàn)因?yàn)榧?lì)頻率與系統(tǒng)固有頻率相等且無阻尼時(shí)系統(tǒng)幅值為無窮大的情況，這點(diǎn)與線性系統(tǒng)截然不同.

圖9 激勵(lì)頻率之和等于橫向基頻時(shí)的振動(dòng)響應(yīng)Fig.9 Vibration responses when summation of frequency of excitation forces equals to transverse fundamental frequency

4 結(jié)論

縱向振動(dòng)與橫向振動(dòng)耦合現(xiàn)象在空間框架和火箭結(jié)構(gòu)中普遍存在，尤其是捆綁運(yùn)載火箭更是如此.在對(duì)此類結(jié)構(gòu)，尤其是對(duì)大推力火箭系統(tǒng)，進(jìn)行力學(xué)特性分析時(shí)有必要考慮這種耦合現(xiàn)象.通過對(duì)控制方程特性和模擬結(jié)果的分析比較，得到如下結(jié)論：

1)梁的縱橫耦合系統(tǒng)頻率是時(shí)變的，若非線性程度較小時(shí)，其頻率趨近對(duì)應(yīng)線性系統(tǒng)的固有振動(dòng)頻率.

2)梁縱橫耦合非線性系統(tǒng)的響應(yīng)是由各次諧波組成的，各次諧波頻率是系統(tǒng)的縱橫頻率和激勵(lì)頻率的各種代數(shù)組合，包括倍頻等.

3)若沒有橫向激勵(lì)和橫向初始擾動(dòng)，縱向振動(dòng)不會(huì)引起橫向振動(dòng).

4)當(dāng)激勵(lì)頻率或組合和系統(tǒng)的某階頻率接近時(shí)，響應(yīng)的振幅可能增大，但通常不會(huì)像線性系統(tǒng)那樣變?yōu)闊o窮大，這是因?yàn)榧?lì)頻率是不變的，而非線性系統(tǒng)頻率是時(shí)變的，二者難以相等.當(dāng)縱向激勵(lì)力明顯大于橫向激勵(lì)力時(shí)，縱向?qū)M向振動(dòng)的耦合影響效果會(huì)顯著增強(qiáng).

5)橫向激勵(lì)的幅值增加n倍時(shí)，縱向振動(dòng)幅頻特性大約增加n2倍.

References)

［1］邢譽(yù)峰，謝珂，潘忠文.縱向過載環(huán)境下變質(zhì)量歐拉梁動(dòng)特性分析方法［J］.北京航空航天大學(xué)學(xué)報(bào)，2013，39(8)：999-1003.

Xing Y F，Xie K，Pan Z W.Dynamic analysis methods of variable-mass Euler beam under longitudinal overload［J］.Journal of Beijing UniversityofAeronauticsand Astronautics，2013，39(8)：999-1003(in Chinese).

［2］邢譽(yù)峰，謝珂，潘忠文.變質(zhì)量系統(tǒng)振動(dòng)分析的兩種方法［J］.北京航空航天大學(xué)學(xué)報(bào)，2013，39(7)：858-862.

Xing Y F，Xie K，Pan Z W.Two methods for vibration analysis of variable-mass systems［J］.Journal of Beijing University of Aeronautics and Astronautics，2013，39(7)：858-862(in Chinese).

［3］夏品奇，Pun D.縱橫向耦合梁的諧波響應(yīng)分析［J］.振動(dòng)工程學(xué)報(bào)，1995，8(1)：67-72.

Xia P Q，Pun D.Harmonic responses of beams with longitudinal and transversal coupling［J］.Journal of Vibration Engineering，1995，8(1)：67-72(in Chinese).

［4］奈弗A H，穆克D T.非線性振動(dòng)(上)［M］.北京：高等教育出版社，1990：215-242.

Nayfeh A H，Mook D T.Nonlinear oscillations volumel 1［M］.Beijing：Higher Education Press，1990：215-242(in Chinese).

［5］奈弗A H，穆克D T.非線性振動(dòng)(下)［M］.北京：高等教育出版社，1990：93-140.

Nayfeh A H，Mook D T.Nonlinear oscillations volume 2［M］.Beijing：Higher Education Press，1990：93-140(in Chinese).

［6］ Han S M，Benaroya H.Non-linear coupled transverse and axial vibration of a compliant structure，Part 1：Formulation and free vibration［J］.Journal of Sound and Vibration，2000，237(5)：837-873.

［7］ Han S M，Benaroya H.Non-linear coupled transverse and axial vibration of a compliant structure，Part 2：Forced vibration［J］.Journal of Sound and Vibration，2000，237(5)：875-900.

［8］胡義，楊建國.梁縱橫耦合振動(dòng)研究［J］.武漢理工大學(xué)學(xué)報(bào)：交通科學(xué)與工程版，2010，34(3)：537-541.

Hu Y，Yang J G.Studies on the longitudinal and lateral coupled vibration of beam［J］.Journal of Wuhan University of Technology：Transportation Science ＆ Engineering，2010，34(3)：537-541(in Chinese).

［9］黃建亮，陳樹輝.縱橫向耦合軸向運(yùn)動(dòng)梁的自由振動(dòng)響應(yīng)研究［J］.動(dòng)力學(xué)與控制學(xué)報(bào)，2010，8(4)：316-321.

Huang J L，Chen S H.Study on vibration of an axially moving beam with coupled transverse and longitudinal motions［J］.Journal of Dynamics and Control，2010，8(4)：316-321(in Chinese).

［10］ Kuchnicki S N，Benaroya H.Coupled transverse and axial motion of a compliant structure in response to vortex-shedding loads［J］.Journal of Sound and Vibration，2002，257(5)：903-929.

［11］ Xu X X，Jiao S J，Cheng J L，et al.Study and test of coupled axial and transverse vibrations on drillpipe for rotary drilling rig［J］.Applied Mechanics and Materials，2014，472：73-78.

［12］張雷，王永.尾部有推力自由梁的振動(dòng)分析［J］.彈道學(xué)報(bào)，2010，22(1)：83-86.

Zhang L，Wang Y.Vibration analysis of free-free beam with end rocket thrust［J］.Journal of Ballistics，2014，22(1)：83-86(in Chinese).

［13］邢譽(yù)峰，李敏.計(jì)算固體力學(xué)原理與方法［M］.北京：北京航空航天大學(xué)出版社，2011：51-59.

Xing Y F，Li M.The principle and method of computational solid mechanics［M］.Beijing：Beihang University Press，2011：51-59(in Chinese).

［14］徐榮橋.結(jié)構(gòu)分析有限元法與MATLAB程序設(shè)計(jì)［M］.北京：人民交通出版社，2011：32-40.

Xu R Q.Structural analysis and MATLAB program design of the finite element method［M］.Beijing：China Communications Press，2006：32-40(in Chinese).

［15］邢譽(yù)峰，李敏.工程振動(dòng)基礎(chǔ)［M］.2版.北京：北京航空航天大學(xué)出版社，2011：199-210.

Xing Y F，Li M.Foundation of engineering vibration［M］.2nd ed.Beijing：Beihang University Press，2011：199-210(in Chinese).