馮美君，蘇國(guó)營(yíng)

(1.中國(guó)航發(fā)哈爾濱東安發(fā)動(dòng)機(jī)有限公司，黑龍江 哈爾濱 150000；2.北京航空航天大學(xué) 機(jī)械工程及自動(dòng)化學(xué)院，北京 100191)

在研究實(shí)際工程問(wèn)題時(shí)，有很多問(wèn)題的運(yùn)動(dòng)方程和數(shù)學(xué)模型可以用含有參數(shù)激勵(lì)的非線性動(dòng)力系統(tǒng)來(lái)描述，例如，在橫浪或縱向波沖擊下船舶的操縱穩(wěn)定性和傾覆機(jī)理問(wèn)題[1-3]，繩系衛(wèi)星的振動(dòng)與控制問(wèn)題[4]，內(nèi)燃機(jī)氣門(mén)機(jī)構(gòu)和離心擺式減振器的振動(dòng)問(wèn)題[5]，以及機(jī)械柔性結(jié)構(gòu)系統(tǒng)的屈曲問(wèn)題[6]等。而螺旋錐齒輪傳動(dòng)系統(tǒng)實(shí)質(zhì)上是一個(gè)高維強(qiáng)參數(shù)周期系統(tǒng)，研究動(dòng)力學(xué)特性本質(zhì)上是對(duì)參數(shù)激勵(lì)非線性系統(tǒng)的響應(yīng)和穩(wěn)定性分析，因而研究包含參數(shù)激勵(lì)的非線性動(dòng)力系統(tǒng)對(duì)解決實(shí)際工程問(wèn)題具有十分重要的意義。

Bogoliubov和Mitroplsky等[7]用平均值的方法研究分析了Mathieu方程的漸近解，推動(dòng)了含有參數(shù)激勵(lì)非線性振動(dòng)理論的發(fā)展。Sethna和Bajaj等[8-9]把分岔思想引入到非線性振動(dòng)理論研究中，利用分岔思想研究含有參數(shù)激勵(lì)的非線性動(dòng)力系統(tǒng)的分岔，結(jié)果反饋良好。Malhotra和SriNamachchivaya等[10]研究分析了非線性可逆系統(tǒng)在非半單共振條件下的全局動(dòng)力學(xué)問(wèn)題。國(guó)內(nèi)陳予恕等[11-13]利用L-S(Liapunov-Schmidt)方法和奇異性理論研究分析了含有參數(shù)激勵(lì)的非線性動(dòng)力系統(tǒng)的分岔。

上述各位學(xué)者的研究方法雖然各有所長(zhǎng)，但在對(duì)類似螺旋錐齒輪的高維強(qiáng)參數(shù)周期系統(tǒng)的處理上仍然十分困難，現(xiàn)有的非線性振動(dòng)理論對(duì)實(shí)際工程中具有周期運(yùn)動(dòng)特性傳動(dòng)系統(tǒng)的動(dòng)態(tài)設(shè)計(jì)和分析有很多不足。

本文從狀態(tài)轉(zhuǎn)移矩陣、瞬態(tài)和穩(wěn)態(tài)受迫響應(yīng)等方面推導(dǎo)并構(gòu)建了偽解析方法，對(duì)包含無(wú)限多頻率成分的時(shí)變參數(shù)在規(guī)定頻譜內(nèi)進(jìn)行連續(xù)階躍調(diào)制，并進(jìn)行了螺旋錐齒輪嚙合振動(dòng)的Hill-Meissner算例驗(yàn)證。

1 參數(shù)激勵(lì)系統(tǒng)的受迫響應(yīng)

對(duì)于N自由度廣義參數(shù)激勵(lì)系統(tǒng)的動(dòng)力學(xué)問(wèn)題，可由一個(gè)2N階矩陣方程來(lái)描述：

(1)

(2)

構(gòu)建狀態(tài)轉(zhuǎn)移矩陣Φ(t,0)=W(t)W-1(0)，則式2可以轉(zhuǎn)化為：

(3)

根據(jù)Floquet理論可知，時(shí)變周期為T(mén)的狀態(tài)轉(zhuǎn)移矩陣Φ(T,0)可轉(zhuǎn)換為關(guān)于矩陣Γ的矩陣指數(shù)函數(shù)：

Φ(T,0)=eΓT

(4)

式中，Γ是一個(gè)固定矩陣。然后定義1個(gè)周期為T(mén)的矩陣為：

P(t)=Φ(t,0)e-ΓT

(5)

則系統(tǒng)任一時(shí)刻的狀態(tài)轉(zhuǎn)移矩陣及其可逆陣可以表述為：

Φ(t,0)=P(t)eΓ t
Φ-1(τ,0)=e-Γ τP-1(τ)

(6)

另外，系統(tǒng)的受迫響應(yīng)解也可以由若干個(gè)余弦激勵(lì)響應(yīng)解構(gòu)成，此時(shí)強(qiáng)迫激勵(lì)項(xiàng)可表述為：

(7)

式中，r是余弦激勵(lì)的幅值向量；ωf是強(qiáng)迫激勵(lì)頻率。將式7代入式3中可得：

(8)

由于Γ為對(duì)角化且其矩陣冪級(jí)數(shù)收斂，則式8中的矩陣指數(shù)函數(shù)應(yīng)用Lagrange-Sylvester定理將變化成Jordan規(guī)范形的變換陣和Γ的特征值的函數(shù)，即：

e±Γ τ=Vdiag(eμit)V-1(i=1,2,…,2N)

(9)

式中，μi是矩陣Γ的特征值也即特征指數(shù)；V是特征向量矩陣。因?yàn)橹芷诰仃嘝(t)和P-1(t)都是連續(xù)且非奇異矩陣，因而求響應(yīng)的偽閉合解時(shí)，可以使用傅里葉級(jí)數(shù)展開(kāi)方法，因此可以表示為：

(10)

式中，ωp是參數(shù)激勵(lì)頻率，且ωp=2π/T。

至此，為消去式8中積分運(yùn)算所需的變換都已經(jīng)引入，可以得出最終的求解受迫響應(yīng)的方程式。首先，將式9和式10代入式8中替換積分內(nèi)部項(xiàng)；其次，將求和部分移出再展開(kāi)積分；再將式9和式10代入上步，求解式中替換求和運(yùn)算外部項(xiàng)，消去可乘的互逆矩陣，進(jìn)行合并求和運(yùn)算；最后，將每項(xiàng)中2個(gè)對(duì)角矩陣的相應(yīng)元素相乘，可得：

(11)

式11就是參數(shù)激勵(lì)系統(tǒng)在強(qiáng)迫激勵(lì)下受迫響應(yīng)的廣義級(jí)數(shù)形式，右邊第1項(xiàng)表示穩(wěn)定響應(yīng)部分，第2項(xiàng)表示瞬態(tài)響應(yīng)部分。

如果系統(tǒng)是穩(wěn)定的，則矩陣?；騀loquet指數(shù)μ的特征值具有全部負(fù)實(shí)數(shù)，隨著時(shí)間的增加，受迫響應(yīng)的瞬態(tài)部分將接近零；相反，如果任何Floquet指數(shù)具有正實(shí)數(shù)，則解的瞬態(tài)部分將發(fā)散到無(wú)窮大，這表明不穩(wěn)定性。

本文中，計(jì)算狀態(tài)轉(zhuǎn)換矩陣的方法是將時(shí)間段劃分為更小的子區(qū)間，假設(shè)每個(gè)子區(qū)間上的時(shí)間是恒定不變的，例如，將周期T分為n個(gè)剛度恒定的子區(qū)間，n個(gè)子區(qū)間終點(diǎn)對(duì)應(yīng)的狀態(tài)轉(zhuǎn)移矩陣分別為Φ(t1,0)、Φ(t2,0)直至Φ(tn-1,0)、Φ(T,0)，則周期狀態(tài)轉(zhuǎn)移矩陣Φ(T,0)可以用子區(qū)間的狀態(tài)轉(zhuǎn)移矩陣合成，如下：

Φ(T,0)=Φ(t1,0)Φ(t2,t1)…Φ(T,tn-1)

(12)

另外，矩陣Γ和Φ(T,0)共享特征向量矩陣，且Φ(T,0)的特征值λi和Γ的特征值μi間存在關(guān)系μi=ln(λi)/T；因此，可以求得式11中特征指數(shù)μi和V矩陣。這樣，式11實(shí)現(xiàn)了受迫響應(yīng)的級(jí)數(shù)形式解析計(jì)算，是一種偽閉合解。

2 時(shí)變參數(shù)的連續(xù)階躍調(diào)制

對(duì)于受迫響應(yīng)的偽閉合解，其精度完全取決于周期矩陣的傅里葉級(jí)數(shù)展開(kāi)和周期狀態(tài)轉(zhuǎn)移矩陣。

對(duì)周期矩陣的傅里葉級(jí)數(shù)展開(kāi)，一般要求最大諧波次數(shù)對(duì)應(yīng)的頻率大于時(shí)變參數(shù)激勵(lì)最高頻率成分的2倍。當(dāng)時(shí)變參數(shù)激勵(lì)包含無(wú)限多頻率成分或最大頻率成分不可知時(shí)，應(yīng)用離散譜分析可以確定其主要頻率成分和幅值，進(jìn)而選定合適的諧波次數(shù)以構(gòu)建周期矩陣的傅里葉級(jí)數(shù)。

對(duì)采用周期分段法求得的時(shí)變剛度系統(tǒng)周期狀態(tài)轉(zhuǎn)移矩陣，其精度與周期分段數(shù)n以及時(shí)變剛度在各個(gè)段上的取值有關(guān)。最簡(jiǎn)單直接的時(shí)變剛度分段采樣方法是根據(jù)其曲線的變化形式，取在每個(gè)子區(qū)間上的平均值；但這種采樣往往包含無(wú)限多頻率成分，造成的結(jié)果是在響應(yīng)計(jì)算時(shí)P(t)傅里葉級(jí)數(shù)展開(kāi)時(shí)會(huì)出現(xiàn)混疊效應(yīng)，進(jìn)而使式11出現(xiàn)失真。

針對(duì)時(shí)變剛度參數(shù)各個(gè)分段上的采樣值進(jìn)行連續(xù)階躍調(diào)制[14]。時(shí)變參數(shù)的連續(xù)階躍調(diào)制采樣如圖1所示。

圖1 時(shí)變參數(shù)的連續(xù)階躍調(diào)制采樣

取時(shí)變剛度的采樣值(如平均值)的指數(shù)展開(kāi)傅里葉系數(shù)為Km，則調(diào)制后的時(shí)變參數(shù)可表示為：

(13)

階躍調(diào)制后，時(shí)變參數(shù)在各個(gè)子區(qū)間上的采樣值為：

(14)

則有：

(15)

且其指數(shù)展開(kāi)的傅里葉系數(shù)滿足：

(16)

另外，將式16代入式14，可將調(diào)制后時(shí)變參數(shù)在各個(gè)子區(qū)間上的采樣值表示為有限項(xiàng)的和，可得連續(xù)階躍調(diào)制后時(shí)變參數(shù)采集值的一般方程形式為：

(17)

3 含不規(guī)律時(shí)變參數(shù)的二自由度Hill-Meisssner方程

在航空事業(yè)中廣泛應(yīng)用的螺旋錐齒輪，其單齒嚙合接觸軌跡不能簡(jiǎn)單的簡(jiǎn)化為直線，且每個(gè)瞬時(shí)接觸點(diǎn)的載荷呈不均勻性，其齒面幾何形狀還與機(jī)床加工調(diào)整參數(shù)有密切關(guān)系，這些因素的疊加使得螺旋錐齒輪嚙合剛度呈現(xiàn)無(wú)規(guī)律的時(shí)變特性。不規(guī)律周期時(shí)變參數(shù)如圖2所示。

圖2 不規(guī)律周期時(shí)變參數(shù)

考慮兩端帶支承剛度的二自由度螺旋錐齒輪傳動(dòng)系統(tǒng)動(dòng)力學(xué)Hill方程：

(18)

將式10中周期矩陣傅里葉級(jí)數(shù)展開(kāi)的諧波次數(shù)分別取32和64，將周期π劃分成80個(gè)子區(qū)間，采用上述方法求解得到x2的受迫響應(yīng)時(shí)間歷程如圖3所示。由圖3可以看出，在含有余弦強(qiáng)迫激勵(lì)和不規(guī)律時(shí)變參數(shù)情況下，該系統(tǒng)方程是穩(wěn)定的，可以獲得穩(wěn)定周期解。

圖3 含不規(guī)律時(shí)變參數(shù)方程穩(wěn)定解的時(shí)間歷程

若不考慮一端的支承剛度，則含時(shí)變參數(shù)剛度矩陣的右下角項(xiàng)由2+δ變?yōu)棣模匦虑蠼獾玫侥Ｗ畲蟮腇loquet乘子從(1,0)點(diǎn)穿出單位圓，出現(xiàn)鞍結(jié)分岔，x1和x2的失穩(wěn)受迫響應(yīng)相軌跡如圖4所示，其中，x1為典型的極限環(huán)振動(dòng)[15]。

圖4 含不規(guī)律時(shí)變參數(shù)方程失穩(wěn)解的相軌跡

4 結(jié)語(yǔ)

本文通過(guò)采用分段法解析計(jì)算了狀態(tài)轉(zhuǎn)移矩陣，對(duì)周期矩陣序列的傅里葉級(jí)數(shù)的開(kāi)展和矩陣指數(shù)函數(shù)的Lagrange-Sylvester定理變換，獲得了只依賴于狀態(tài)轉(zhuǎn)移矩陣的瞬態(tài)響應(yīng)和穩(wěn)態(tài)響應(yīng)的偽閉合解。對(duì)包含無(wú)限多頻率成分的時(shí)變參數(shù)在規(guī)定頻譜內(nèi)進(jìn)行連續(xù)階躍調(diào)制，避免了周期分段法計(jì)算狀態(tài)轉(zhuǎn)移矩陣引起的響應(yīng)失真解。通過(guò)對(duì)包含不規(guī)律時(shí)變嚙合剛度的螺旋錐齒輪嚙合振動(dòng)二自由度Hill方程的算例分析，驗(yàn)證了本文所建立數(shù)學(xué)分析方法的正確性和穩(wěn)定性。對(duì)于航空事業(yè)中廣泛應(yīng)用的螺旋錐齒輪，可以有效分析其不規(guī)律時(shí)變嚙合剛度引起的振動(dòng)，為螺旋錐齒輪傳動(dòng)系統(tǒng)動(dòng)力學(xué)研究建立了數(shù)學(xué)層面的基礎(chǔ)。

[1] Sanchez N E, Nayfeh A H. Nonlinear rolling motions of ships in longitudinal waves[J]. International Shipbuilding Progress, 1990, 37(4): 247-272.

[2] Sanchez N E, Nayfeh A H. Prediction of bifurcations in a parametrically excited duffing oscillator[J]. International Journal of Non-Linear Mechanics, 1990, 25(2/3): 163-176.

[3] 陳予恕, 張偉. 現(xiàn)代數(shù)學(xué)和力學(xué)(MMM)VI：船舶橫搖運(yùn)動(dòng)的全局分岔和混沌動(dòng)力學(xué)[M]. 蘇州: 蘇州大學(xué)出版社, 1995.

[4] 彭建華, 劉延柱. 繩系衛(wèi)星的混沌運(yùn)動(dòng)[J]. 上海交通大學(xué)學(xué)報(bào), 1996, 30(11): 32-36.

[5] Zhang W, Ye M. Local and global bifurcations of valve mechanism[J]. Nonlinear Dynamics, 1994, 6(3): 301-316.

[6] 張清杰, 劉士光, 鄭際嘉, 等. 結(jié)構(gòu)動(dòng)力屈曲研究進(jìn)展[J]. 力學(xué)進(jìn)展, 1993, 23(4): 530-538.

[7] Bogoliubov N N, Mitropolsky Y M. Asymptotic methods in the theory of nonlinear oscillations[M]. New York: Gordon and Breach, 1961.

[8] Sethna P R, Bajaj A K. Bifurcations in dynamical systems with internal resonances[J]. ASME Journal of Applied Mechanics, 1978, 45(4): 895-902.

[9] Bajaj A K, Shaw S W. Foreword[J]. Nonlinear Dynamics, 1993, 4(6): 527-530.

[10] Malhotra N, Sri Namachchivaya N. Global dynamics of parametrically excited nonlinear reversible systems with nonsemisimple 1∶1 resonance[J]. Physica D: Nonlinear Phenomena, 1995, 89(1/2): 43-70.

[11] 陳予恕, Langford W F. 非線性馬休方程的亞諧分岔解及歐拉動(dòng)彎曲問(wèn)題[J]. 力學(xué)學(xué)報(bào), 1988, 20(6): 522-532.

[12] 陳予恕. 非線性振動(dòng)系統(tǒng)的分岔和混沌理論[M]. 北京: 高等教育出版社, 1993.

[13] 陳予恕, Han R P S. 現(xiàn)代工程非線性動(dòng)力學(xué)的現(xiàn)狀展望[J]. 力學(xué)與實(shí)踐, 1995, 17(6): 11-19.

[14] Franks L E, Sandberg I W. An alternative approach to the realization of network transfer functions: The N-path filter[J]. Bell System Technical Journal, 1960, 39(5): 1321-1350.

[15] Schleimer J H, Stemmler M. Coding of information in limit cycle oscillators[J]. Physical Review Letters, 2009, 103(24): 248105.1-248105.4.