準(zhǔn)周期響應(yīng)對稱破缺分岔點(diǎn)的一種快速計(jì)算方法1)

龔冰清 鄭澤昌 陳衍茂 劉濟(jì)科

(中山大學(xué)理論與應(yīng)用力學(xué)系,廣州 510275)

引言

周期或準(zhǔn)周期振動(dòng),是非線性動(dòng)力學(xué)系統(tǒng)出現(xiàn)的典型穩(wěn)態(tài)響應(yīng)[1-3].不同于周期振動(dòng),準(zhǔn)周期響應(yīng)不具有固定的振動(dòng)周期,狀態(tài)變量不再按某周期有規(guī)律地重復(fù)出現(xiàn).更具體地說,準(zhǔn)周期振動(dòng)的特征是其振動(dòng)頻率由多個(gè)不可公約的頻率線性組合而成[4-6].因而,這些不可約的頻率也被稱為基頻,由不可約頻率線性組合而成的準(zhǔn)周期頻率稱為組合頻率.

在不同的非線性系統(tǒng)中,準(zhǔn)周期響應(yīng)出現(xiàn)的機(jī)制機(jī)理也不盡相同.出現(xiàn)準(zhǔn)周期響應(yīng)的系統(tǒng)類型也是多種多樣,例如,受到多不可約頻率外激勵(lì)作用的系統(tǒng)[7]、非線性耦合導(dǎo)致的模態(tài)共振[8]以及參數(shù)激振與自激振動(dòng)之間的耦合[9]等情況都可能產(chǎn)生準(zhǔn)周期響應(yīng).從振動(dòng)的頻率成分看,準(zhǔn)周期振動(dòng)可以看作是不可公約頻率之間的耦合效應(yīng),而引起的復(fù)合響應(yīng).在實(shí)際工程中,準(zhǔn)周期振動(dòng)在各種科學(xué)和工程問題中均被大量報(bào)道,例如機(jī)械工程[10]、電子工程[11]、氣動(dòng)彈性[12]等領(lǐng)域.

在求解方面,數(shù)值方法由于簡單明了且易于實(shí)現(xiàn),因此被廣泛運(yùn)用于準(zhǔn)周期響應(yīng)的求解[13-15].這類基于時(shí)間離散的時(shí)程響應(yīng)跟蹤數(shù)值算法,一般無法直接獲得系統(tǒng)的不穩(wěn)定穩(wěn)態(tài)解,即不穩(wěn)定的周期或準(zhǔn)周期響應(yīng).因?yàn)?不穩(wěn)定解意味著微小的擾動(dòng)或誤差會使得解出現(xiàn)越來越大的偏差,從而偏離不穩(wěn)定解本身.眾所周知,系統(tǒng)的分岔通常會伴隨著某些解的穩(wěn)定性改變.為了進(jìn)一步研究系統(tǒng)響應(yīng)的分岔現(xiàn)象,研究人員提出了基于Fourier 級數(shù)展開的半解析半數(shù)值計(jì)算方法,用于求解周期響應(yīng),例如諧波平衡法[1]和增量諧波平衡IHB 法[16-19]諧波展開類方法、L-P 法和多尺度法等攝動(dòng)類方法[1]、平均法和KBM 法等漸進(jìn)方法[1]以及同倫分析法[20]等.這些方法有的也被改進(jìn)或推廣、用來求解準(zhǔn)周期響應(yīng),其中多時(shí)間尺度IHB 法可以計(jì)算強(qiáng)非線性、多自由度系統(tǒng),獲得了廣泛的應(yīng)用.

關(guān)于系統(tǒng)的平衡點(diǎn)和周期解的分岔問題,目前已有大量的研究和成熟的理論,例如Hopf 分岔定理[21]和Floquet 理論[22]等.相較于平衡點(diǎn)分岔,周期響應(yīng)的分岔計(jì)算是一個(gè)更加棘手的問題.20 世紀(jì)80 年代,陳予恕[1]提出了周期響應(yīng)分岔計(jì)算的C-L 方法,并在諸多科學(xué)和工程中的非線性系統(tǒng)中得以廣泛應(yīng)用.在IHB 法中,若引入諧波系數(shù)限制增量過程,可直接求解周期響應(yīng)倍周期分岔或?qū)ΨQ性破缺的臨界值[23].

對于準(zhǔn)周期響應(yīng),其分岔計(jì)算方法和分析理論就更加棘手,也更加缺乏.目前來看,對于非線性動(dòng)力系統(tǒng)的準(zhǔn)周期響應(yīng),在大多數(shù)情況下,其分岔點(diǎn)的確定是通過反復(fù)的數(shù)值模擬來實(shí)現(xiàn)的.這種方法的原理簡單、過程直接,但一方面計(jì)算量往往很大,另一方面如前面所述、數(shù)值方法也無法提供不穩(wěn)定的解,而且存在不利于分析參數(shù)的影響規(guī)律等問題.為此,本文將周期響應(yīng)對稱破缺分岔計(jì)算方法[23]進(jìn)一步推廣,使之對于準(zhǔn)周期同樣有效,該方法是基于增量諧波平衡過程的.利用多時(shí)間尺度的IHB 法,設(shè)計(jì)迭代格式,通過追蹤特定的諧波系數(shù),從而快速求解準(zhǔn)周期響應(yīng)對稱破缺分岔點(diǎn).以多頻激勵(lì)Duffing系統(tǒng)以及Duffing-van der Pol 耦合系統(tǒng)為例,闡述方法的主要步驟,驗(yàn)證方法的正確性和有效性.

1 增量諧波平衡法

增量諧波平衡法是求解非線性振動(dòng)問題的一種半數(shù)值、半解析方法,是由Lau 等[24]在20 世紀(jì)80 年代提出并發(fā)展起來的.IHB 法的本質(zhì)是Newton-Raphson 增量過程和諧波平衡過程的結(jié)合[25].經(jīng)過了幾十年的發(fā)展,IHB 不僅可用于周期或準(zhǔn)周期響應(yīng)的求解,已應(yīng)用于動(dòng)力學(xué)響應(yīng)分析的諸多領(lǐng)域,如混沌控制[26]、反問題識別[27]等.

考慮含不可約頻率雙外激勵(lì)的非線性振動(dòng)系統(tǒng)

其中M,C,K分別是質(zhì)量、阻尼和剛度矩陣,G(x,x′)是非線性函數(shù).f1和f2是外激勵(lì)向量幅值,ω1和 ω2是外激勵(lì)的不可約頻率.不失一般性,將參數(shù) ω1作為分岔控制參數(shù).

IHB 法的第一步是增量過程,即Newton-Raphson 過程.讓xi表示一種振動(dòng)狀態(tài),對應(yīng)的控制參數(shù)為 ω1=ω1i,其臨近狀態(tài)可以通過增量的方式表示為

將式(2)代入式(1),展開并忽略高階小量,得到線性化的增量方程

其中

稱為不平衡力,當(dāng)xi和ω1i為方程(1) 的精確解時(shí),R(xi,ω1i)=0.

IHB 法的第二步是諧波平衡過程,也即Ritz-Garlerkin 平均過程,令

把式(4)和式(5)代入式(3),并應(yīng)用諧波平衡過程,可以得到線性方程組

其 中 ?u=[?c0,?c10,?s10,···,?cMN,?sMN]T,方 陣Ku是 (M+N)(M+N+1)+1 維的雅可比矩陣,Kω是控制參數(shù)的梯度向量,KR是剩余向量.這些矩陣和向量可通過初始解xi和 ωi確定.在IHB 法中,通常應(yīng)選擇一個(gè)主動(dòng)變化的參數(shù)(稱為主動(dòng)增量)來控制解的連續(xù)延拓.如果 ? ω 被給出,方程(6)描述的每一步增量為 ?u的一組方程,可以迭代求解.

2 對稱破缺分岔點(diǎn)的快速算法

下面引入補(bǔ)充方程,并將 ? ω 看成被動(dòng)增量參與迭代,進(jìn)而獲得系統(tǒng)發(fā)生對稱性破缺的分岔值.作為一般性的討論,可將準(zhǔn)周期解嚴(yán)格表示為廣義的Fourier 級數(shù)

其中y(t)也可作為x(t)的一個(gè)分量.在實(shí)際計(jì)算中,對于單自由度系統(tǒng),y(t)即是x(t);對于多自由度系統(tǒng),y(t)可取為x(t)的第一個(gè)分量.對稱的準(zhǔn)周期響應(yīng)的相圖指的是,對于任意一點(diǎn) (y(t),y′(t)),總能找到關(guān)于原點(diǎn)與之對稱的點(diǎn) (?y(t),?y′(t)).也即是表明,如果y(t) 是方程(1)的解,那么 ?y(t) 也是一個(gè)解.

引入雙時(shí)間尺度 τ1=ω1t和τ2=ω2t,和周期解對稱性質(zhì)類似[23],對于關(guān)于原點(diǎn)對稱的準(zhǔn)周期解,其對稱性由以下等式描述

由式(7),可知

根據(jù)方程(7)～方程(9)可知,當(dāng)i+j為奇數(shù)時(shí),對稱條件自然滿足;當(dāng)i+j為偶數(shù)時(shí),若要滿足對稱條件,則意味著

注意到i+j的奇偶性和 |i|+|j| 相同.那么,當(dāng)偶次諧波系數(shù)取到非零的小量時(shí),準(zhǔn)周期解便失去對稱性,也可以看成是系統(tǒng)發(fā)生對稱破缺分岔的臨界狀態(tài).為此,在迭代格式(6)中,補(bǔ)充方程

其中 0 <λ0?1,那么通過迭代過程,構(gòu)造以下迭代格式

其中Kc=(1 0 ··· 0) 是補(bǔ)充方程(11) 關(guān)于諧波系數(shù)的Jacobi 矩陣.只要 λ0足夠小時(shí),便可以近似確定關(guān)于參數(shù) ω1的對稱破缺分岔點(diǎn).理論上,λ0越小則 ω1越可能接近真實(shí)的分岔值.在實(shí)際計(jì)算中,可以取 λ0為非常小的正數(shù),例如λ0=1.0×10?8等.當(dāng)然,最終的近似程度也還受到解的截?cái)嘀C波系數(shù)的影響.這點(diǎn)在下面的算例部分將予以討論.

3 結(jié)果驗(yàn)證與討論

首先,考慮含不可約頻率雙外激勵(lì)的Duffing 振子系統(tǒng)[20,28-29]

圖1 是所考慮Duffing 系統(tǒng)的準(zhǔn)周期響應(yīng)相圖,其中圖1(a1)～圖1(c1)為IHB 法結(jié)果,圖1(a2)～圖1(c2)為4 階Runge-Kutta (RK)法結(jié)果.在求解過程中,IHB 解截取了最高階M+N=10 階諧波,數(shù)值解采用了IHB 解作為初始條件,通過RK 法積分獲得.準(zhǔn)周期解的穩(wěn)定性判別目前還是一個(gè)非常棘手的問題.一般的數(shù)值積分方法如RK 法和Newmark 法等無法追蹤不穩(wěn)定解,通過這一直觀特點(diǎn)簡單判斷解的穩(wěn)定性.分別對比圖1(a1)、圖1(a2)(ω1=0.8)及圖1(c1)、圖1(c2)(ω1=1),IHB 結(jié)果與RK 結(jié)果是完全吻合,說明這些準(zhǔn)周期解都是穩(wěn)定的.進(jìn)一步觀察到,在圖1(a1)中,準(zhǔn)周期響應(yīng)相圖是關(guān)于原點(diǎn)中心對稱的,因此是對稱的.

圖1 Duffing 系統(tǒng)(13)準(zhǔn)周期響應(yīng)相圖,(a1) ω1=0.8,穩(wěn)定IHB 法解;(b1) ω1=1.0,不穩(wěn)定的IHB 法解;(c1) ω1=1.0,穩(wěn)定的IHB 法解;(a2),(b2),(c2)為相應(yīng)的RK 數(shù)值解.其中,IHB 法的截取諧波數(shù)為 M+N=10Fig.1 The quasi-periodic responses of Duffing system (13).(a1) Stable with ω1=0.8,(b1) unstable with ω1=1.0 and (c1) stable with ω1=1.0,all obtained by the IHB method with the truncated number of harmonics as M+N=10 ;(a2),(b2),(c2) the corresponding solutions provided by the RK method

當(dāng) ω1=1.0 時(shí),而圖1(b1),圖1(b2)中的結(jié)果并不完全相同.其中,由IHB 法解作為初始值所得的RK法結(jié)果(圖1(b2)所示),從對稱的IHB 法解(圖1(b1)所示)開始,逐漸收斂到另一個(gè)非對稱的準(zhǔn)周期解,即圖1(b2)所示非對稱部分.實(shí)際上,圖1(b2)所示的非對稱部分、在相圖上和圖1(c1)及圖1(c2)所示解是一樣的.IHB 法是半解析半數(shù)值方法,可以得到不穩(wěn)定的對稱準(zhǔn)周期解和穩(wěn)定的非對稱準(zhǔn)周期解.由此可知系統(tǒng)在參數(shù)范圍 0.8<ω1<1.0 發(fā)生了對稱破缺分岔,且對稱的準(zhǔn)周期解失去了穩(wěn)定性.

為準(zhǔn)確找到系統(tǒng)發(fā)生對稱破缺的分岔點(diǎn),根據(jù)上述方法,將 λ0逐步減小,得到收斂結(jié)果,其中 ω1的變化如圖2 所示.隨著 λ0的取值越來越小,ω1的值也逐漸收斂.當(dāng)諧波數(shù)取M+N=8 時(shí),可以計(jì)算得到分岔值約為0.933 02,當(dāng)諧波數(shù)取M+N=10 時(shí),可以計(jì)算得到分岔值約為0.937 58.

圖2 采用不同的 λ0 所對應(yīng)的對稱破缺分岔點(diǎn)Fig.2 The symmetry breaking point predicted by the presented method with different value of λ0

圖3 是由IHB 法獲得的Duffing 系統(tǒng)準(zhǔn)周期響應(yīng)的分岔圖.在P點(diǎn)(約為 ω1=0.935)附近,對稱的準(zhǔn)周期解產(chǎn)生非對稱的分支.此時(shí),系統(tǒng)出現(xiàn)了對稱破缺分岔現(xiàn)象.和周期解分岔類似,經(jīng)歷了分岔后、準(zhǔn)周期解也一般會穩(wěn)定性改變.周期解的穩(wěn)定性一般可通過經(jīng)典的Floquet 理論準(zhǔn)確判別[1].然而,對于準(zhǔn)周期解,其穩(wěn)定性判別是一個(gè)相當(dāng)棘手的問題.為此,通過分析系統(tǒng)穩(wěn)態(tài)響應(yīng)的相圖變化,如圖1 所示,進(jìn)而直觀地判別IHB 法所得結(jié)果的穩(wěn)定性.研究發(fā)現(xiàn),對稱準(zhǔn)周期響應(yīng)由于對稱破壞而失去穩(wěn)定性,同時(shí)出現(xiàn)的非對稱解是穩(wěn)定的.

圖3 由IHB 法獲得的Duffing 系統(tǒng)(13)準(zhǔn)周期響應(yīng)分岔圖,其中P 為對稱破缺分岔點(diǎn)Fig.3 The bifurcation chart of the quasi-periodic response of Duffing system (13),where P denotes the symmetry breaking point

進(jìn)一步,用數(shù)值法驗(yàn)證準(zhǔn)周期響應(yīng)對稱破缺的發(fā)生.圖4 給出由RK 法得到的解的最大和最小振幅隨 ω1的變化圖.可以觀察到,在 ω1=0.937 5 前,準(zhǔn)周期解的最大值和最小值完全重合,即解是對稱的;而隨著 ω1增大,最大值和最小值之間逐漸出現(xiàn)差值.也即是表明,系統(tǒng)已經(jīng)失去了對稱性,并且大約在ω1=0.937 5 時(shí)發(fā)生了對稱破缺分岔.這也與IHB 法結(jié)果基本是一致的.

圖4 RK 法所得Duffing 系統(tǒng)(13)穩(wěn)定的準(zhǔn)周期解的最大和最小振幅變化圖Fig.4 The maximal and minimal of the amplitude of the stable quasiperiodic response of Duffing system (13),provided by the RK method

圖5 給出了分別由快速Fourier 變換和IHB 法獲得的準(zhǔn)周期響應(yīng)頻譜圖.首先,IHB 法結(jié)果和數(shù)值結(jié)果吻合良好.數(shù)值結(jié)果的部分峰值沒有IHB 法結(jié)果與之對應(yīng),因?yàn)镮HB 法中組合諧波進(jìn)行了截?cái)?取了M+N=7,即 |i|+|j|≥8 的組合諧波沒有包含在IHB 法結(jié)果中.當(dāng)參數(shù)取 ω1=0.937 時(shí),圖5(a)所示,響應(yīng)只包含奇數(shù)次諧波成分.對應(yīng)于解的表達(dá)式(7),出現(xiàn)了當(dāng) |i|+|j| 是奇數(shù)(即1,3,5,7 等)時(shí)的若干諧波成分.當(dāng)參數(shù)為 ω1=0.938 時(shí),除了占主要部分的奇數(shù)次諧波之外,還出現(xiàn)了幅值較小的偶數(shù)次諧波,圖5(b)所示的 2 ω1,即i+j=2.圖6 是準(zhǔn)周期響應(yīng)的頻率譜瀑布圖,由此可更直觀地觀察到,當(dāng)參數(shù) ω1逐漸增加,響應(yīng)出現(xiàn)了偶次諧波成分.根據(jù)基于關(guān)系式(8)～式(10)的討論,偶數(shù)次諧波分量的出現(xiàn),表明準(zhǔn)周期響應(yīng)逐漸失去了對稱性.由此推斷,對稱性破缺分岔點(diǎn)出現(xiàn)在0.937 和0.938之間.

圖5 由快速Fourier 變換得到的Duffing 系統(tǒng)(13)準(zhǔn)周期響應(yīng)的頻譜圖,其中藍(lán)線是數(shù)值方法結(jié)果,紅點(diǎn)是IHB 法結(jié)果,組合諧波截?cái)鄶?shù)為 M+N=7Fig.5 The FFT spectrum for the quasi-periodic response of Duffing system (13),where the blue line denotes the numerical result and red the IHB result,with the number of truncated harmonics as M+N=7

圖6 Duffing 系統(tǒng)(13)準(zhǔn)周期響應(yīng)的頻率譜瀑布圖Fig.6 The frequency waterfall for the quasi-periodic response of Duffing system (13)

Duffing 系統(tǒng)和van der Pol 系統(tǒng)都是典型的非線性振動(dòng)系統(tǒng).由它們耦合組成的多自由度系統(tǒng),更是蘊(yùn)含了豐富的動(dòng)力學(xué)行為,是非線性研究廣泛采用的模型之一[30-31].為了進(jìn)一步檢驗(yàn)所提方法的有效性,考慮Duffing-van der Pol 耦合振子

其中,Duffing 振子受雙諧波外激勵(lì),參數(shù)選擇和系統(tǒng)(13)相同;耦合系數(shù)取 δ=0.2 ;ρ=1 是van der Pol振子的非線性系數(shù).

對于含強(qiáng)非線性環(huán)節(jié)的van der Pol 系統(tǒng),一些依賴于小參數(shù)的攝動(dòng)類方法有時(shí)難以獲得準(zhǔn)確的周期或準(zhǔn)周期解.由于不依賴于小參數(shù),IHB 和同倫分析等方法對含van der Pol 振子的強(qiáng)非線性系統(tǒng)有效,因而被廣泛采用[20,31-33].

圖7 給出了本文算法所得的分岔值與預(yù)設(shè)小量λ0之間的關(guān)系.當(dāng) λ0逐漸減小時(shí),所得分岔值快速收斂于固定值,約為1.000 4.圖8 所示準(zhǔn)周期響應(yīng)的最大和最小振幅,當(dāng) ω1超過1 且接近1.001 時(shí),最大和最小振幅出現(xiàn)了微小的偏差.這表明,在此參數(shù)區(qū)間,響應(yīng)逐漸失去了對稱性,也即是出現(xiàn)了對稱性破缺分岔,驗(yàn)證了計(jì)算結(jié)果的準(zhǔn)確性.

圖7 Duffing-van der Pol 系統(tǒng)(14)的對稱破缺分岔值關(guān)于 λ0 的收斂圖,其中截取諧波數(shù)為 M+N=8Fig.7 The symmetry breaking point predicted by the presented method with different value of λ0,with the number of truncated number of harmonics as M+N=8 for Duffing-van der Pol system (14)

圖8 RK 法所得Duffing-van der Pol 系統(tǒng)(14)準(zhǔn)周期解的最大和最小振幅Fig.8 The maximal and minimal of the amplitude of the stable quasiperiodic response of Duffing-van der Pol system (14),provided by the RK method

4 結(jié)論

本文基于IHB 法,提出了一種快速計(jì)算準(zhǔn)周期響應(yīng)對稱破缺分岔的方法.方法的基本依據(jù)是,當(dāng)發(fā)生對稱破缺分岔前,準(zhǔn)周期響應(yīng)的某些諧波系數(shù)為0;發(fā)生分岔的過程中,這些系數(shù)逐漸變?yōu)榉? 的小量.那么,依據(jù)這一性質(zhì),在IHB 法迭代格式中,預(yù)先令某個(gè)諧波系數(shù)為給定的小量,可以將控制參數(shù)視為被動(dòng)變量,通過收斂的IHB 法結(jié)果,直接獲得響應(yīng)的分岔值.

研究了含不可公約頻率雙外激勵(lì)Duffing 振子系統(tǒng)、以及Duffing-van der Pol 耦合系統(tǒng),驗(yàn)證了所提方法.結(jié)果表明,方法能準(zhǔn)確得到系統(tǒng)準(zhǔn)周期響應(yīng)的對稱破缺分岔點(diǎn),得到了數(shù)值法結(jié)果的驗(yàn)證.本文方法不僅可以找到不穩(wěn)定解,更重要的是可以直接提供控制參數(shù),而不需要反復(fù)的跟蹤和試錯(cuò)過程.