邢譽(yù)峰，梁昆

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

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

邢譽(yù)峰等在文獻(xiàn)［1-2］中，對過載火箭結(jié)構(gòu)系統(tǒng)的動力學(xué)特性和變質(zhì)量火箭系統(tǒng)的分析方法進(jìn)行了研究.值得指出的是，運(yùn)載火箭尤其是捆綁運(yùn)載火箭的縱橫耦合特性是不可忽略的，這主要是由于火箭具有柔性結(jié)構(gòu)，并且同時存在縱橫載荷(如推力和風(fēng)載等).當(dāng)彈性變形不再是小變形時，系統(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)縱橫耦合非線性動力學(xué)方程，并分析了自由振動和受迫振動響應(yīng)特性.2010年，文獻(xiàn)［8］推導(dǎo)了梁縱橫向耦合振動的非線性單元剛度矩陣，用有限元方法求解了動態(tài)響應(yīng)，并分析了耦合響應(yīng)中的倍頻現(xiàn)象等.文獻(xiàn)［9］利用Galerkin法和諧波增量平衡法研究縱向運(yùn)動梁在縱橫耦合情況下的自由振動響應(yīng)，尤其是在橫向第1階和第2階固有頻率之比接近1∶3內(nèi)共振條件下的系統(tǒng)響應(yīng).文獻(xiàn)［10］研究了水中塔架等效梁模型的耦合振動.文獻(xiàn)［11］用模擬和試驗(yàn)方法研究了縱橫耦合振動.有些學(xué)者則研究縱向激勵對橫向振動的影響，如文獻(xiàn)［12］.

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

1 基本方程

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

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

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

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

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

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

和邊界條件，見表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)可以看出，縱向振動中除了包括系統(tǒng)縱向頻率和激勵頻率信息之外，還包括橫向振動頻率和激勵頻率的各種組合和倍頻.

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

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

在下面算例中，還將詳細(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ε和動能函數(shù)T的標(biāo)準(zhǔn)形式，即

式中單元矩陣為

其中：

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

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

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

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

這里所用模型的主要參數(shù)來自文獻(xiàn)［7］，表2給出了幾何模型和材料參數(shù)、有限元模型參數(shù)和邊界條件，表3給出了激勵載荷.本文作者用MATLAB［14］自編了動力學(xué)響應(yīng)分析程序，其中用的是Newmark直接積分法.時域范圍為0～1 s，時間步長h=0.1 ms.單元類型是三次梁單元，其縱向變形是線性的，3個結(jié)點(diǎn)參數(shù)包括縱向位移、橫向位移及其導(dǎo)數(shù)(斜率).激勵作用點(diǎn)和動態(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 激勵載荷Table 3 External forces

3.1 模型的正確性

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

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

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

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

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

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

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

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

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

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

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

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

3.2.1 受迫振動頻率耦合

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

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

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

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

3.2.2 自由振動頻率耦合

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

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

圖4 受迫振動的幅頻特性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為兩個橫向頻率之差;252Hz、430Hz、328Hz、504Hz為兩個橫向頻率之和;600 Hz、650 Hz為3個橫向頻率之代數(shù)和(228+402 ±25)Hz.76 Hz、126 Hz還可以分別表示為前2階橫向頻率之差與之和、202 Hz還可以表示為第1和第3階橫向頻率之差、302 Hz還可以表示為第2和第4階橫向頻率之差.

3.3 縱橫相互影響特征

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

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

圖5 自由振動頻譜Fig.5 Amplitude-frequency characteristics of free vibration

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

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

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

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

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

表6 橫向激勵載荷加倍時縱向振動幅頻特性Table 6 Longitudinal amplitude-frequencycharacteristics when transverse excitation forces double

3.4 共振現(xiàn)象

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

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

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

4 結(jié)論

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

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

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

3)若沒有橫向激勵和橫向初始擾動，縱向振動不會引起橫向振動.

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

5)橫向激勵的幅值增加n倍時，縱向振動幅頻特性大約增加n2倍.

References)

［1］邢譽(yù)峰，謝珂，潘忠文.縱向過載環(huán)境下變質(zhì)量歐拉梁動特性分析方法［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)振動分析的兩種方法［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］.振動工程學(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.非線性振動(上)［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.非線性振動(下)［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］胡義，楊建國.梁縱橫耦合振動研究［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)動梁的自由振動響應(yīng)研究［J］.動力學(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］張雷，王永.尾部有推力自由梁的振動分析［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ù)峰，李敏.工程振動基礎(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).