邢譽峰，梁昆

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

飛行器結構物理參數的主要特點之一是其具有時變質量.運載火箭除了具有時變質量參數之外，其幾何特征是其具有類似細長梁的構形，其柔性特征致使其彈性變形非常重要，無論對載荷分析還是飛行控制皆如此.

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

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

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

1 基本方程

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

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

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

系統應變能函數Uε和動能函數T分別為

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

動能函數中忽略了耦合項.根據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)，初步得到如下結論：

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

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

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

在下面算例中，還將詳細討論上述結論.

2 有限元求解

2.1 位移函數

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

式中：

式中：l為單元的長度.

2.2 結構單元矩陣

把位移函數式(7)代入能量泛函，可以得到用結點位移向量表示的應變能函數Uε和動能函數T的標準形式，即

式中單元矩陣為

其中：

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

把結構單元組裝后，根據變分原理可得系統的離散形式動力學方程：

式中：M、K和F分別為結構的總體質量矩陣、剛度矩陣和載荷列向量.

3 數值模擬和分析

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

表2 模型參數Table 2 Model parameters

表3 激勵載荷Table 3 External forces

3.1 模型的正確性

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

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

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

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

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

通過圖2和表4的比較可以得到結論，本文推導的公式和程序實現是正確的，為下文對梁的縱橫耦合非線性特性的分析奠定了基礎.

3.2 耦合系統的頻率特性

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

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

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

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

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

3.2.1 受迫振動頻率耦合

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

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

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

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

3.2.2 自由振動頻率耦合

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

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

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

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

縱向頻譜特性中，50 Hz、202 Hz、454 Hz 分別為系統橫向頻率的二倍頻;76 Hz、302 Hz為橫向固有頻率的三倍頻;126 Hz、378 Hz為兩個橫向頻率之差;252Hz、430Hz、328Hz、504Hz為兩個橫向頻率之和;600 Hz、650 Hz為3個橫向頻率之代數和(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)響應，從圖中可以看出，如果沒有橫向激勵或橫向初始條件，縱向振動不能引起橫向振動，參見式(6)，這是因為此時w'和w″始終為零.若沒有縱向激勵和縱向初始條件，橫向振動卻可以引起縱向振動.

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

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

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

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

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

縱向激勵力為零，不考慮初始條件，分別考慮橫向激勵力幅值分別為Fw1=20 kN和Fw2=2Fw1作用時縱向動態(tài)響應幅值的變化規(guī)律，激勵頻率為11 Hz.圖8和表6給出了分析結果，從中可以看到兩種情況下，Fw2作用下的縱向振動幅值基本為Fw1作用下幅值的4倍.從方程(6a)右端第2項可以看出此結論的合理性.表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 共振現象

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

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

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

4 結論

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

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

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

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

4)當激勵頻率或組合和系統的某階頻率接近時，響應的振幅可能增大，但通常不會像線性系統那樣變?yōu)闊o窮大，這是因為激勵頻率是不變的，而非線性系統頻率是時變的，二者難以相等.當縱向激勵力明顯大于橫向激勵力時，縱向對橫向振動的耦合影響效果會顯著增強.

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

References)

［1］邢譽峰，謝珂，潘忠文.縱向過載環(huán)境下變質量歐拉梁動特性分析方法［J］.北京航空航天大學學報，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］邢譽峰，謝珂，潘忠文.變質量系統振動分析的兩種方法［J］.北京航空航天大學學報，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.縱橫向耦合梁的諧波響應分析［J］.振動工程學報，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］.武漢理工大學學報：交通科學與工程版，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］黃建亮，陳樹輝.縱橫向耦合軸向運動梁的自由振動響應研究［J］.動力學與控制學報，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］.彈道學報，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］邢譽峰，李敏.計算固體力學原理與方法［M］.北京：北京航空航天大學出版社，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］徐榮橋.結構分析有限元法與MATLAB程序設計［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］邢譽峰，李敏.工程振動基礎［M］.2版.北京：北京航空航天大學出版社，2011：199-210.

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