肖龍江，黃建亮

(中山大學(xué)航空航天學(xué)院，廣東 廣州 510275)

準(zhǔn)周期運(yùn)動是指系統(tǒng)的振動響應(yīng)中含有多個不可約頻率的運(yùn)動[1]。當(dāng)外激勵的頻率之間不可公度或系統(tǒng)的多個模態(tài)之間存在內(nèi)共振時將會產(chǎn)生準(zhǔn)周期運(yùn)動。準(zhǔn)周期運(yùn)動的存在會導(dǎo)致響應(yīng)的振幅、頻率和相位角隨著時間連續(xù)變化，使得能量在不同模態(tài)之間實時遷移，導(dǎo)致結(jié)構(gòu)的破壞[2-4]。Afaneh和Ibrahim[5]利用多尺度法分析了1∶1內(nèi)共振條件下兩端固支屈曲梁的前兩階振動響應(yīng)，發(fā)現(xiàn)在該條件下屈曲梁可能存在準(zhǔn)周期響應(yīng)。Kreider和Nayfeh[6]通過實驗發(fā)現(xiàn)，在前兩階模態(tài)存在1∶1內(nèi)共振條件下，兩端固支屈曲梁會產(chǎn)生頻譜圖中存在等間距邊頻帶的準(zhǔn)周期運(yùn)動，Emam和Nayfeh[7]采用數(shù)值計算方法驗證了上述實驗結(jié)果。Lau等[8-9]把增量諧波平衡法和多時間尺度法結(jié)合起來，用于計算系統(tǒng)的準(zhǔn)周期運(yùn)動。Huang和Zhu[10-11]改進(jìn)了Lau等[8-9]的方法，用兩時間尺度的增量諧波平衡法研究了一端固定一端簡支的彎曲梁的準(zhǔn)周期運(yùn)動。在國內(nèi)，袁銘鴻等[12]研究了耦合系統(tǒng)的準(zhǔn)周期運(yùn)動，從他們得到的準(zhǔn)周期運(yùn)動的頻譜圖中可以觀察到等間距的邊頻帶。張丹偉等[13]用兩個時間尺度的IHB法研究了外激勵頻率不可公約時非線性系統(tǒng)的準(zhǔn)周期運(yùn)動。李小彭等[14]通過數(shù)值計算分析了多自由度內(nèi)共振系統(tǒng)的準(zhǔn)周期運(yùn)動。

本文研究了在基礎(chǔ)簡諧激勵作用下，前兩階對稱模態(tài)存在1∶3內(nèi)共振的兩端固支屈曲梁的非線性振動。先通過Galerkin方法對梁的振動方程進(jìn)行離散，然后將傳統(tǒng)的單一時間尺度的增量諧波平衡和Floquet理論[15-16]相結(jié)合，分析外激勵振幅變化時振動方程的周期響應(yīng)及其穩(wěn)定性的變化。當(dāng)梁的周期響應(yīng)發(fā)生Hopf分岔演變?yōu)闇?zhǔn)周期響應(yīng)時，采用兩時間尺度的IHB法研究準(zhǔn)周期響應(yīng)，并驗證結(jié)果的準(zhǔn)確性。

1 兩端固支屈曲梁的振動方程

(1)

圖1 基礎(chǔ)簡諧激勵作用下的兩端固支屈曲梁簡圖Fig.1 Schematic of a fixed-fixed buckled subject to base harmonic excitation

在此基礎(chǔ)諧波激勵作用下，由哈密頓原理可求得屈曲梁的振動方程為

(2)

(3)

對方程進(jìn)行無量綱化，令

(4)

帶入方程和邊界條件中可得

(5)

y(0,t)=y(1,t)=y′(0,t)=y′(1,t)=0

(6)

其中，點和撇分別表示對時間和空間求導(dǎo)，且

(7)

為了求出靜態(tài)屈曲函數(shù)v(x)，舍去方程(5)和(6)中的時間相關(guān)項和外激勵項可得

(8)

v(0)=v(1)=v′(0)=v′(1)=0

(9)

(10)

其中，b為無量綱化跨中撓度。把方程(10)帶入方程(8)中可解得

(11)

令方程的解為

y(x,t)=u(x,t)+v(x)

(12)

帶入方程(5)及邊界條件(6)可得

(13)

u(0,t)=u(1,t)=u′(0,t)=u′(1,t)=0

(14)

用文獻(xiàn)[17]中的方法求解兩端固支屈曲梁的固有頻率和模態(tài)函數(shù)[17-18]。舍去方程中的非線性項、阻尼項和外激勵項，可得到方程所對應(yīng)的線性化自由振動方程，令方程的解為

(15)

由此可得

(16)

(17)

方程(16)的通解為

Φj=c1sinλjx+c2cosλjx+c3sinhμjx

+c4coshμjx+c5cos 2πx

(18)

其中，

(19)

(20)

可得各階固有頻率所對應(yīng)的模態(tài)函數(shù)Φj，各階模態(tài)函數(shù)滿足條件

(21)

其中，δij為狄拉克函數(shù)。用Galerkin方法求解方程，為此令

(22)

將式(22)帶入方程(13)，在方程兩邊左乘Φi(x)，并對方程兩邊從0到1對x積分，可得

=γficos(ωt),i= 1,2,…,n

(23)

其中，

取

ξi=ξ,i=1,2,…,n

(25)

ξ為無量綱化模態(tài)阻尼系數(shù)，令

(26)

可將方程(23)化為

=γficos(ωt),i=1,2,…,n

(27)

2 兩個時間尺度的IHB法

用兩個新的無量綱時間變量

τ1=ωt,τ2=ωdt

(28)

代替變量t，因此qi(t)變?yōu)樽兞喀?和τ2的函數(shù)。

qi(t)=qi(ωt,ωdt)=qi(τ1,τ2)

(29)

由此可得到

(30)

把式(28)、(29)、(30)帶入式(27)可得

=γficosτ1,i=1,2,…,n

(31)

式(31)可以寫為矩陣形式

+(K+K(2)+K(3))q=γFcosτ1

(32)

(33)

兩時間尺度IHB法的第一步是增量過程。令q0、ωd0和γ0表示一個振動狀態(tài)，其鄰近振動狀態(tài)以增量形式表示為

q=q0+Δq,

ωd=ωd0+Δωd,

γ=γ0+Δγ

(34)

其中，q0=[q10,q20,…,qn0]T，Δq=[Δq1,Δq2,…,Δqn]T。把式(34)帶入式(32)，并且略去高于一階的小量，可得

(35)

其中

(36)

R是誤差向量。當(dāng)q0、ωd0和γ0是方程(32)的準(zhǔn)確解時，R=0。

兩時間尺度的IHB法的第二步是諧波平衡過程。因為準(zhǔn)周期運(yùn)動的頻譜圖在載波頻率(頻率的整數(shù)倍)的周圍存在著等間距的邊頻帶，所以可以用Galerkin方法把qi0展開為多重Fourier級數(shù)。

=Cs·Ai

(38)

其中，nc和ns分別為載波頻率的余弦項項數(shù)和正弦項項數(shù)，nf是一個載波頻率左側(cè)或右側(cè)所攜帶的邊頻帶頻率的數(shù)目，同時

Cs=[C0C1C2…CncS0S1…Sns]

(39)

(40)

其中

C0=[1cosτ2… cosnfτ2]nf+1,

S0=[sinτ2sin2τ2… sinnfτ2]nf,

(41)

Ci1=[cos(i1τ1-nfτ2) … cosi1τ1… cos(i1τ1+nfτ2)]2nf+1,i1=1,2,…,nc

(42)

Si1=[sin(i1τ1-nfτ2) … sini1τ1… sin(i1τ1+nfτ2)]2nf+1,i1=1,2,…,ns

(43)

(44)

(45)

把增量Δqi也展開為多重Fourier級數(shù)

Δqi=Cs·ΔAi

(47)

其中，

由此可得

q0=SA, Δq=SΔA

(48)

其中，S=diag(Cs,Cs, … ,Cs)，A=[A1A2…An]T，ΔA=[ΔA1ΔA2… ΔAn]T。

對式(35)運(yùn)用Galerkin 平均過程可得

(49)

把式(48)帶入式(49)可以產(chǎn)生關(guān)于變量ΔA、 Δωd和Δγ的線性方程

(50)

其中

(51)

(52)

其中

方程(50)的未知數(shù)數(shù)量比方程數(shù)量多2，在計算時應(yīng)先指定兩個增量作為控制變量。由于 Δωd無法預(yù)先得知，因此應(yīng)該在向量ΔA中指定兩個增量作為控制變量(為了得到準(zhǔn)周期解，兩個增量不能都是載波頻率對應(yīng)的諧波系數(shù))，指定控制變量后，先令兩個增量為0，由式(50)可求解得到 Δωd、Δγ以及ΔA中除了控制變量之外的其他量。把各增量加到原來的解上得到新解，判斷新解是否能使誤差向量R的模小于給定精度(在本文中，精度為1.0×10-10)；如果不能，把新解帶入式(50)，重復(fù)上述過程，不斷循環(huán)直至找到符合精度要求的解。之后，給控制變量加上給定的增量，使循環(huán)得以更新，重復(fù)上一循環(huán)，計算出新情況下的解。

3 計算結(jié)果與分析

為了準(zhǔn)確反映屈曲梁非線性振動的特性，需要選取足夠多的模態(tài)數(shù)量對振動方程進(jìn)行Galerkin離散。Emam和Nayfeh[7, 19]的計算表明，對于兩端固支屈曲梁，選取前四階模態(tài)進(jìn)行計算是必要的。本文選取前四階模態(tài)進(jìn)行計算。

取nc=ns=5，nf=0，Δωd=0帶入式(50)可得傳統(tǒng)的單一時間尺度的IHB法，以Δγ為增量可計算出周期響應(yīng)。把周期解的形式由式(38)變?yōu)?/p>

(55)

其中

(56)

由于

(57)

基礎(chǔ)諧波激勵只加載在梁的對稱模態(tài)而未加載在反對稱模態(tài)上。因此在外激勵振幅較小時，系統(tǒng)的反對稱模態(tài)并未被激發(fā)。直至γ增加到0.22時，由于第一和第三模態(tài)之間存在1∶3內(nèi)共振，能量向反對稱模態(tài)轉(zhuǎn)移，反對稱模態(tài)被激發(fā)，被激發(fā)的反對稱模態(tài)和對稱模態(tài)所含的頻率成分相同。繼續(xù)增加r值，追蹤反對稱模態(tài)被激發(fā)之后的周期響應(yīng)，并用Floquet理論分析周期解的穩(wěn)定性?？梢园l(fā)現(xiàn)在r∈[0.27,0.67]時，周期解不穩(wěn)定，F(xiàn)loquet乘子中有一對復(fù)共軛特征值與單位圓相交。因此，在這一段內(nèi)周期解發(fā)生Hopf分岔，演變?yōu)闇?zhǔn)周期解，如圖2所示。

圖2 ω=37，b=0.002 62時的振幅響應(yīng)曲線Fig.2 Amplitude response curves whileω=37, b=0.002 62

(58)

其中

(59)

由此，可計算得到準(zhǔn)周期響應(yīng)曲線，如圖2所示。準(zhǔn)周期運(yùn)動頻譜圖的邊頻帶間隔ωd的值也隨之計算出來。在本例中，隨著γ的增加，邊頻帶間隔ωd逐漸增加，如圖3所示。

圖3 屈曲梁的準(zhǔn)周期響應(yīng)中ωd隨 γ的變化Fig.3 Variation of ωdwithγin the quasi-periodic responses of buckled beam

以γ=0.27時的準(zhǔn)周期解為例進(jìn)行進(jìn)一步的探究，其第一模態(tài)的頻譜圖如圖4所示。從圖中可以看出，在載波頻率(頻率ω的整數(shù)倍)周圍，存在著等間距的邊頻帶，相鄰兩個邊頻帶的間隔為ωd≈3.65。ω與ωd的比值是一個無理數(shù)，正是這兩個頻率不可公約導(dǎo)致系統(tǒng)的響應(yīng)變?yōu)闇?zhǔn)周期運(yùn)動，使它的時間歷程圖表現(xiàn)出明顯的“拍現(xiàn)象”，如圖5所示。從圖5還可以看出用IHB法得到的時間歷程圖與用四階Runge-Kutta(RK)法的數(shù)值結(jié)果完全重合。為了進(jìn)一步驗證兩時間尺度的IHB法的準(zhǔn)確性，圖6和圖7分別比較了用兩時間尺度的IHB法和四階Runge-Kutta法得到的準(zhǔn)周期運(yùn)動的相圖和龐加萊截面圖，結(jié)果也是完全一致的。從圖7可以看出，準(zhǔn)周期運(yùn)動的龐加萊截面圖是一條封閉曲線。

圖4 ω=37，b=0.002 62，γ=0.27時q1的頻譜圖Fig.4 Spectrum of q1 while ω=37, b=0.002 62, γ=0.27

圖5 ω=37，b=0.002 62，γ=0.27時準(zhǔn)周期響應(yīng)的時間歷程圖Fig.5 Time histories of the quasi-periodic response while ω=37, b=0.002 62, γ=0.27

圖6 ω=37，b=0.002 62，γ=0.27時準(zhǔn)周期響應(yīng)的相圖Fig.6 Phase plane portraits of the quasi-periodic response while ω=37, b=0.002 62, γ=0.27

圖7 ω=37，b=0.002 62，γ=0.27時準(zhǔn)周期響應(yīng)的龐加萊截面圖Fig.7 Poincaré sections of the quasi-periodic response while ω=37, b=0.002 62, γ=0.27

4 結(jié) 論

本文用多時間尺度法與IHB法相結(jié)合，得到了適用于求解系統(tǒng)的準(zhǔn)周期運(yùn)動的兩時間尺度的IHB法。對于前兩階對稱模態(tài)存在1∶3內(nèi)共振的兩端固支屈曲梁，當(dāng)作用在其上的基礎(chǔ)簡諧激勵的振幅處于一定范圍內(nèi)時，系統(tǒng)的周期響應(yīng)會發(fā)生Hopf分岔，演變?yōu)闇?zhǔn)周期響應(yīng)。文中用兩時間尺度的IHB得到了屈曲梁的準(zhǔn)周期響應(yīng)曲線圖，同時比較了兩時間尺度的IHB法和四階Runge-Kutta法所得到的結(jié)果。對比結(jié)果表明，兩時間尺度的IHB法有較高的精度，適用于解決頻譜圖中存在的等間距邊頻帶的準(zhǔn)周期運(yùn)動。