重構(gòu)諧波平衡法及其求解復(fù)雜非線(xiàn)性問(wèn)題應(yīng)用1)

代洪華 王其偲 嚴(yán)子樸 岳曉奎

(西北工業(yè)大學(xué)航天學(xué)院,西安 710072)

(航天飛行動(dòng)力學(xué)技術(shù)國(guó)家級(jí)重點(diǎn)實(shí)驗(yàn)室,西安 710072)

引言

工程中多數(shù)動(dòng)力學(xué)問(wèn)題,其數(shù)學(xué)模型都是非線(xiàn)性的,線(xiàn)性系統(tǒng)只是在一定假設(shè)及限制條件下對(duì)非線(xiàn)性系統(tǒng)的理想化近似[1].隨著控制科學(xué)與航空宇航任務(wù)日益復(fù)雜,強(qiáng)非線(xiàn)性的影響越發(fā)不容忽視.諧波平衡法(harmonic balance method,HB)在求解非線(xiàn)性系統(tǒng)周期解中應(yīng)用廣泛,但作為一種半解析半數(shù)值方法,隨著系統(tǒng)自由度、方法階數(shù)的增加,公式推導(dǎo)工作將變得困難[2].借助計(jì)算機(jī)數(shù)學(xué)軟件可輔助推導(dǎo)工作,但對(duì)計(jì)算效率及性能的提升作用仍然有限.采用7 階HB 法推導(dǎo)立方非線(xiàn)性項(xiàng),Mathematica代碼長(zhǎng)達(dá)1321 行.Hall 等[3]提出了高維諧波平衡法(high-dimensional harmonic balance,HDHB),通過(guò)“頻域非線(xiàn)性量的時(shí)域近似計(jì)算”來(lái)簡(jiǎn)化公式推導(dǎo),但是,由于引入近似關(guān)系導(dǎo)致非物理假解問(wèn)題[4-5].Dai 等[6]發(fā)現(xiàn)了頻域和時(shí)域非線(xiàn)性項(xiàng)之間的條件等價(jià)恒等式,基于此,首創(chuàng)了重構(gòu)諧波平衡法(reconstruction harmonic balance,RHB)實(shí)現(xiàn)超高階(N＞ 100)高精計(jì)算,并給出時(shí)域配點(diǎn)計(jì)算的最優(yōu)采樣定理,從理論上消除非物理解.

HB 類(lèi)方法(HB 法及其改進(jìn)方法)不僅用于計(jì)算簡(jiǎn)單非線(xiàn)性系統(tǒng)(杜芬、范德波方程等),還發(fā)展到航空航天、深空探測(cè)等前沿領(lǐng)域.哈爾濱工業(yè)大學(xué)陳毅等[7]使用諧波平衡?交變頻域/時(shí)域法(HBalternating frequency-time,HB-AFT)用于求解航空發(fā)動(dòng)機(jī)中雙轉(zhuǎn)子?軸承?機(jī)匣系統(tǒng)動(dòng)力學(xué)方程.中山大學(xué)陳衍茂等[8]使用增量諧波平衡法對(duì)帶機(jī)外掛載的二元機(jī)翼進(jìn)行動(dòng)力學(xué)特性研究.RHB 法也被作為三體軌道的設(shè)計(jì)依據(jù),計(jì)算結(jié)果符合鵲橋中繼衛(wèi)星實(shí)際飛行數(shù)據(jù)[6].

但是該方法仍面臨著兩個(gè)問(wèn)題: (1)諧波平衡法的本質(zhì)依賴(lài)于非線(xiàn)性項(xiàng)的傅里葉級(jí)數(shù)展開(kāi),因此,受限于多項(xiàng)式型非線(xiàn)性系統(tǒng)求解,對(duì)于非多項(xiàng)式型復(fù)雜非線(xiàn)性問(wèn)題,諧波平衡法難以適用;(2)已有諧波平衡類(lèi)方法建立在單基頻的假設(shè)上進(jìn)行級(jí)數(shù)展開(kāi),由于擬周期響應(yīng)存在多基頻的特征,因此已有方法難以直接求解高精度擬周期解.

針對(duì)非多項(xiàng)式型非線(xiàn)性系統(tǒng)的求解問(wèn)題,目前有兩類(lèi)處理方法,分別為直接法和間接法.直接法包括HB-Taylor 法[9]與HB-AFT 法[10-11].HB-Taylor 法通過(guò)泰勒級(jí)數(shù)展開(kāi)將非線(xiàn)性函數(shù)用有限階近似多項(xiàng)式描述;HB-AFT 法通過(guò)對(duì)非線(xiàn)性項(xiàng)的時(shí)域值采樣以離散原問(wèn)題.由于直接法對(duì)原系統(tǒng)進(jìn)行了近似,導(dǎo)致求解精度低,且計(jì)算性能分別受制于高階級(jí)數(shù)描述和過(guò)采樣等問(wèn)題.Cochelin 等[12]提出的諧波平衡?重鑄法(HB-recast)是一種間接方法,通過(guò)重鑄(recast)技術(shù),成功將復(fù)雜非線(xiàn)性微分動(dòng)力學(xué)系統(tǒng)無(wú)損變換為多項(xiàng)式型微分代數(shù)方程,然后用HB 法加以求解[13-14].但是,受限于HB 法的高階計(jì)算(使用重鑄法會(huì)增加系統(tǒng)的維度,進(jìn)一步增加了高階計(jì)算的難度)與原重鑄形式的二次型限制(原方法要求新系統(tǒng)中非線(xiàn)性至多為二次多項(xiàng)式),至今難以對(duì)復(fù)雜非線(xiàn)性系統(tǒng)周期響應(yīng)進(jìn)行高效高精求解.

第2 個(gè)難題是擬周期響應(yīng)求解問(wèn)題.擬周期響應(yīng)大量出現(xiàn)在非線(xiàn)性動(dòng)力學(xué)系統(tǒng)中[15-17],其頻率響應(yīng)由多個(gè)不可約基頻及其線(xiàn)性組合描述[18].由于傳統(tǒng)HB 類(lèi)方法基于單個(gè)基頻進(jìn)行近似解逼近,不能簡(jiǎn)單通過(guò)基頻及其整數(shù)倍頻率分量對(duì)擬周期響應(yīng)描述,Chua 等[19]提出了結(jié)合廣義傅里葉級(jí)數(shù)的改進(jìn)多頻HB 法,實(shí)現(xiàn)擬周期響應(yīng)的準(zhǔn)確求解.然而,使用多頻HB 類(lèi)方法對(duì)強(qiáng)非線(xiàn)性項(xiàng)進(jìn)行頻域內(nèi)高階描述時(shí),計(jì)算效率嚴(yán)重受限于高維傅里葉分析(頻域分量由多重積分[20]、求和[21]計(jì)算得到).Liu 等[22]使用多頻HDHB 法求解受迫范德波振子的穩(wěn)態(tài)響應(yīng),避免非線(xiàn)性項(xiàng)諧波系數(shù)的直接表示以提高求解效率,但是多頻HDHB 法的精度受到非物理解的破壞[23].總之,當(dāng)前關(guān)于HB 類(lèi)方法的研究已拓寬到對(duì)擬周期響應(yīng)求解領(lǐng)域,然而由于多基頻計(jì)算中的高階頻域描述困難,對(duì)此類(lèi)復(fù)雜響應(yīng)的高性能求解尚待解決.

針對(duì)上述問(wèn)題,本文提出了RHB 法與重鑄法、多頻諧波平衡計(jì)算相結(jié)合的兩種方法: (1) RHB-重鑄法;(2)重構(gòu)多諧波平衡法(reconstruction multiple harmonic balance,RMHB).一方面RHB-重鑄法通過(guò)將一般非線(xiàn)性等價(jià)轉(zhuǎn)化為多項(xiàng)式型非線(xiàn)性系統(tǒng),再采用RHB 法以確定時(shí)域最優(yōu)采樣點(diǎn)數(shù),實(shí)現(xiàn)對(duì)復(fù)雜非線(xiàn)性動(dòng)力學(xué)系統(tǒng)的高階高精求解,仿真誤差達(dá)計(jì)算機(jī)精度.另一方面,RMHB 法通過(guò)篩選和補(bǔ)充多個(gè)基頻,借鑒RHB 法的時(shí)域等價(jià)重構(gòu)思想,完善了學(xué)術(shù)界在使用多頻HB 類(lèi)方法求解擬周期響應(yīng)時(shí)消除混淆假解的理論研究.

1 諧波平衡法及重構(gòu)諧波平衡法

1.1 諧波平衡法

對(duì)非線(xiàn)性動(dòng)力學(xué)系統(tǒng)

HB 法用有限階傅里葉級(jí)數(shù)來(lái)構(gòu)造近似解及其導(dǎo)數(shù)

其中N是HB 法的階數(shù).x0,x1,···,x2n為未知傅里葉系數(shù),也被稱(chēng)為頻域變量,將頻域變量組成待求解向量=[x0x1···x2n]T.假設(shè)多項(xiàng)式函數(shù)f(x)=x?,忽略由非線(xiàn)性項(xiàng)而出現(xiàn)的高次諧波,只需展開(kāi)前N次的傅里葉分量

其中各次分量r0,r1,···,r2n為

其中n=1,2,···,N;θ=ωt.將表示非線(xiàn)性函數(shù)的傅里葉分量記為,是待求解向量x? 的非線(xiàn)性函數(shù).

將式(2)～式(4)代入微分方程(1),令常數(shù)項(xiàng)及前n次諧波 cosnωt和 sinnωt(n=1,2,···,N) 系數(shù)配平,得到非線(xiàn)性代數(shù)方程組

其中,分塊矩陣A為

由三角函數(shù)的求導(dǎo)關(guān)系得到.

非線(xiàn)性項(xiàng)的近似表示(4)中,符號(hào)運(yùn)算的復(fù)雜度會(huì)隨算法階次的提高呈指數(shù)增長(zhǎng)[6].當(dāng)HB 法的階數(shù)超過(guò)20,即便有計(jì)算機(jī)數(shù)學(xué)軟件的輔助,非線(xiàn)性項(xiàng)的推導(dǎo)整理工作量仍難以接受.

1.2 高維諧波平衡法

為簡(jiǎn)化HB 法的符號(hào)運(yùn)算量,Hall 等[3]使用時(shí)域值替代頻域分量,即將N階HB 法中的傅里葉系數(shù)作用轉(zhuǎn)換矩陣EHDHB,建立與一個(gè)周期 2N+1 個(gè)等距配點(diǎn)上時(shí)域量間的聯(lián)系,定義

將式(7)代入HB 法代數(shù)方程組(6)得到HDHB法求解微分方程(1)對(duì)應(yīng)的代數(shù)方程組

HDHB 法顯著提高計(jì)算效率,并被認(rèn)為是HB法的一種改進(jìn),在計(jì)算流體力學(xué)領(lǐng)域得到應(yīng)用[2,24].但該算法求解強(qiáng)非線(xiàn)性動(dòng)力學(xué)問(wèn)題時(shí)產(chǎn)生非物理解(也稱(chēng)假解),嚴(yán)重影響求解精度.如圖1 所示,HDHB 法計(jì)算結(jié)果雖然是方程組(8)的數(shù)學(xué)解(使求解算法收斂),但已偏離了真實(shí)的物理響應(yīng).

圖1 HB 法和HDHB 法計(jì)算杜芬方程對(duì)比(N=2)Fig.1 Comparison of the Duffing equation computed by the HB and HDHB methods (N=2)

Liu 等[4]指出,假解現(xiàn)象產(chǎn)生于對(duì)非線(xiàn)性項(xiàng)近似的過(guò)程中,存在高階諧波對(duì)低次的混淆(污染).以立方項(xiàng)非線(xiàn)性x3為例,二階HB 法中() 表達(dá)式

HDHB 法計(jì)算中,非線(xiàn)性項(xiàng)的頻域分量(9)中不僅包含了原HB 法中的所有項(xiàng),還包括附加雜項(xiàng),這些雜項(xiàng)的混入導(dǎo)致求解中出現(xiàn)非物理解

1.3 重構(gòu)諧波平衡法

Dai 等[25]證明,HDHB 法與時(shí)域配點(diǎn)法等價(jià),時(shí)域配點(diǎn)法通過(guò)使用時(shí)域配點(diǎn)上函數(shù)值對(duì)微分方程進(jìn)行離散,從而聯(lián)立代數(shù)方程組

其中配點(diǎn)矩陣

混淆是造成計(jì)算出現(xiàn)非物理解的原因.至于高次諧波在時(shí)域配點(diǎn)計(jì)算中影響低次諧波的機(jī)理,則遵循“混淆規(guī)則”[25].

混淆規(guī)則:假設(shè) α ∈[0,π] 被以h為間隔均分為離散時(shí)間點(diǎn),那么配點(diǎn)法中可區(qū)分的最大諧波次數(shù)n∈[?L,L],其中L=π/h,L為“極限波次”.高次諧波n(|n|＞L) 被誤認(rèn)為是相應(yīng)低次諧波na

其中na∈[?L,L],m是整數(shù).

混淆規(guī)則指出,時(shí)域配點(diǎn)法可區(qū)分的最大諧波次數(shù)由配點(diǎn)數(shù)(離散時(shí)間點(diǎn)的間隔)決定.增加時(shí)域配點(diǎn)法中的配點(diǎn)數(shù)量,可區(qū)分的諧波次數(shù)越高,此時(shí)方程(10)個(gè)數(shù)多于待求解變量數(shù),配點(diǎn)矩陣E為列滿(mǎn)秩矩陣,存在Moore-Penrose 廣義逆矩陣E?

其中I2N+1為單位矩陣.在配點(diǎn)法方程同時(shí)作用矩陣E?,對(duì)方程組降維,得到RHB 法代數(shù)方程組

RHB 法在保證算法效率的同時(shí),消除非物理解,從而實(shí)現(xiàn)對(duì)原HB 法的最佳重構(gòu).定理1 圍繞如何選擇合適的配點(diǎn)數(shù),給出消除混淆確定條件[6].

定理1(條件等價(jià)性):如果配點(diǎn)數(shù)M,方法階數(shù)N,多項(xiàng)式非線(xiàn)性的冪次 ? 滿(mǎn)足

則RHB 法與HB 法等價(jià).

為說(shuō)明RHB 法消除混淆的效果.分別使用HDHB法和RHB 法計(jì)算杜芬方程,在分岔處(頻率 ω=2)進(jìn)行蒙特卡洛法模擬,選取104組隨機(jī)初值(各頻域分量在區(qū)間[?2,2]中選取),得到如表1 所示的統(tǒng)計(jì)結(jié)果.統(tǒng)計(jì)結(jié)果表明,HDHB 法計(jì)算產(chǎn)生58 組解,其中55 組都是非物理解;而RHB 法只計(jì)算得到3 個(gè)具有物理含義的解[6].

表1 RHB 法與HDHB 法在分岔處的解分布Table 1 Statistical distribution of solution by the RHB and the HDHB method at bifurcation point

此外,HB-AFT 法與HDHB 法計(jì)算流程一致,但HB-AFT 法的原理是,選用配點(diǎn)數(shù) 2?N+1 來(lái)消除混淆[11],但過(guò)采樣會(huì)占用計(jì)算機(jī)資源,在實(shí)際計(jì)算時(shí)會(huì)導(dǎo)致更大的CPU 和RAM 計(jì)算負(fù)擔(dān)[2].RHB 法基于時(shí)域配點(diǎn)法的統(tǒng)一框架,根據(jù)配點(diǎn)數(shù)的差異將所有HB 類(lèi)方法(HDHB 法、HB-AFT 法等)建立起聯(lián)系.

2 重構(gòu)諧波平衡法改進(jìn)策略

2.1 微分方程重鑄技術(shù)

針對(duì)非多項(xiàng)式型非線(xiàn)性系統(tǒng)的HB 法求解,Cochelin等[12]提出將原系統(tǒng)改寫(xiě)為二次型系統(tǒng)

其中未知向量z包含微分方程的原始變量x及一些新的變量(引入的新變量用以改寫(xiě)方程).c是關(guān)于未知量z的常數(shù)向量;l是關(guān)于z的線(xiàn)性向量值運(yùn)算符;q則是二次向量值運(yùn)算符.

因?yàn)閷?duì)任意x?次非線(xiàn)性,RHB 法都能實(shí)現(xiàn)高效計(jì)算,改寫(xiě)后方程可以是更高次多項(xiàng)式,將式(13)進(jìn)一步寫(xiě)成適配于RHB 法的微分方程重鑄形式

其中n可以是關(guān)于變量z的任意 ? 次多項(xiàng)式函數(shù).重鑄格式(14)涵蓋多類(lèi)型非多項(xiàng)式型非線(xiàn)性問(wèn)題,下面將分類(lèi)加以介紹.

(1)微分項(xiàng)耦合型

通過(guò)將方程重鑄為標(biāo)準(zhǔn)形式(14),得到

原方程轉(zhuǎn)化為典型的非線(xiàn)性度 ?=3 的動(dòng)力學(xué)系統(tǒng).得出結(jié)論: 導(dǎo)數(shù)項(xiàng)的耦合不影響非線(xiàn)性度計(jì)入,在實(shí)際HB 計(jì)算中,任意階導(dǎo)數(shù)只計(jì)入一個(gè)非線(xiàn)性度.

(2)有理分式型

以Rayleigh-Plesset 方程為例

式中,A,B,C,D均為常系數(shù).將方程兩端同除以R,并運(yùn)用倒數(shù)關(guān)系,令v=,x=1/R得到[6,12]

處理有理分式型非線(xiàn)性項(xiàng),通過(guò)額外增加方程Rx=1,將方程改寫(xiě)為非線(xiàn)性度 ?=4 的多項(xiàng)式型非線(xiàn)性系統(tǒng).

(3)根式型

相對(duì)論諧振子方程(17)中非線(xiàn)性項(xiàng)為根式

對(duì)任意形如 αq/p的根式,令 βp=α,使原根式型非線(xiàn)性轉(zhuǎn)化得到多項(xiàng)式型: αq/p=(βp)q/p=βq.

(4)初等超越函數(shù)

對(duì)初等超越函數(shù)的處理,需要導(dǎo)數(shù)信息來(lái)實(shí)現(xiàn)對(duì)原微分方程的重鑄.以非線(xiàn)性單擺方程為例

由于三角函數(shù)是超越函數(shù),不能通過(guò)簡(jiǎn)單的代數(shù)關(guān)系式改寫(xiě)轉(zhuǎn)化為多項(xiàng)式型.首先額外引入兩個(gè)變量s,c來(lái)分別表示s(t)=sinθ(t),c(t)=cosθ(t).利用三角函數(shù)的導(dǎo)數(shù)性質(zhì)

建立補(bǔ)充方程,從而實(shí)現(xiàn)對(duì)微分方程組的重鑄[13]

初等超越函數(shù)非線(xiàn)性的重鑄,以導(dǎo)數(shù)關(guān)系作為方程組改寫(xiě)的依據(jù)(部分需要用到有理分式、根式型非線(xiàn)性重鑄技巧).附錄表格中羅列了幾類(lèi)常見(jiàn)初等超越函數(shù)改寫(xiě)思路與重鑄形式,可供參考.

總之,本文主要針對(duì)非線(xiàn)性項(xiàng)為連續(xù)函數(shù)情形(包括微分項(xiàng)耦合、有理分式型、根式型以及初等超越函數(shù)4 類(lèi)),重鑄技術(shù)通過(guò)變量替換、利用特殊函數(shù)的導(dǎo)數(shù)關(guān)系(見(jiàn)附表1)等關(guān)鍵步驟將一般連續(xù)函數(shù)非線(xiàn)性項(xiàng)改寫(xiě)為更利于重構(gòu)諧波平衡法求解的多項(xiàng)式形式(14).

2.2 多諧波平衡計(jì)算

“補(bǔ)頻”(supplemental frequency)思想[17]是在原來(lái)單頻假設(shè)的基礎(chǔ)上,補(bǔ)充多個(gè)與響應(yīng)相關(guān)的頻率.假設(shè)考慮兩個(gè)基頻,將假設(shè)函數(shù)改寫(xiě)為

參數(shù)m和n滿(mǎn)足不等式

其中p代表二維傅里葉級(jí)數(shù)的截?cái)郲19-20],類(lèi)似于RHB法中階數(shù)N.

以RHB 法理論為基礎(chǔ),提出的RMHB 法利用時(shí)域信息,將動(dòng)力學(xué)問(wèn)題(1)轉(zhuǎn)化為非線(xiàn)性代數(shù)方程(11)求解.多基頻計(jì)算中,配點(diǎn)矩陣E以及轉(zhuǎn)換矩陣E?形式,記cm,n=cos(mω1+nω2)t,sm,n=sin(mω1+nω2)t.在RMHB 法計(jì)算中,存在選取合適的采樣周期和時(shí)域配點(diǎn)數(shù)用以消除混淆的結(jié)論[23].

定理2(多頻計(jì)算中條件等價(jià)性):假設(shè)多項(xiàng)式非線(xiàn)性的冪次 ? 的非線(xiàn)性系統(tǒng)響應(yīng)中包含兩個(gè)基頻,基頻之比 ω1/ω2為有理數(shù).則RMHB 法消除混淆需滿(mǎn)足采樣時(shí)間T=2π/GCD(ω1,ω2),配點(diǎn)數(shù)

其中 GCD 為兩數(shù)最大公因數(shù).

3 數(shù)值仿真結(jié)果

3.1 相對(duì)論諧振子

作為物理學(xué)中經(jīng)典問(wèn)題,有必要對(duì)諧振子模型進(jìn)行完整而嚴(yán)格的相對(duì)論處理[26].考慮一個(gè)靜質(zhì)量為m的質(zhì)點(diǎn)在一維諧振力F=?的作用下做相對(duì)論運(yùn)動(dòng).其中k為彈性常數(shù),為位移量.結(jié)合牛頓運(yùn)動(dòng)運(yùn)動(dòng)學(xué)方程以及動(dòng)量定理,可以推導(dǎo)得到相對(duì)論振子方程[27-28]

圖2 Mickens 變換解與真實(shí)物理解對(duì)比Fig.2 Comparison of Mickens transformation result and real physical solution

RHB-重鑄法按照根式型非線(xiàn)性重鑄原則,將相對(duì)論諧振子改寫(xiě)為方程(18),再使用高階RHB 法求解.結(jié)合表2 和圖3,RHB-重鑄法可對(duì)高速運(yùn)動(dòng)的諧振子(β=0.85)直接進(jìn)行高階計(jì)算.55 階RHB-重鑄法能將計(jì)算誤差控制在 10?12量級(jí)(以數(shù)值積分為參照,全文采用MATLAB 內(nèi)置函數(shù)ode15 s,相對(duì)精度容差 10?13,絕對(duì)誤差容限10?20),總計(jì)算耗時(shí)在1 s 內(nèi).

表2 各階RHB-重鑄法求解相對(duì)論諧振子Table 2 Solving relativistic harmonic oscillator by the RHBrecast method with different orders

圖3 55 階RHB-重鑄法求解相對(duì)論諧振子的計(jì)算結(jié)果Fig.3 The computing result of the RHB-recast method (N=55) for solving relativistic harmonic oscillator

除通過(guò)方法階數(shù)對(duì)求解精度的提高,還能通過(guò)兩種方式實(shí)現(xiàn)對(duì)計(jì)算性能的改善.

(1) 非線(xiàn)性代數(shù)方程組求解算法

非線(xiàn)性代數(shù)方程(11)求解算法能提高HB 類(lèi)方法計(jì)算性能,Zheng 等[31]結(jié)合Tikhonov 正則化求解,黃建亮等[32]引入狗腿法思想結(jié)合回溯線(xiàn)搜索算法求解,Thomas 等[33]使用Broyden 法以提高HB 法計(jì)算性能.本文根據(jù)算例來(lái)說(shuō)明L-M (Levenberg-Marquardt)法較之傳統(tǒng)迭代方法在求解RHB 方程組時(shí)體現(xiàn)出的優(yōu)勢(shì).圖4 展示當(dāng) β=0.85,頻域量初值估計(jì)x2=1,u0=0.7,使用兩種不同的方程組求解算法: 牛頓迭代 (Newton-Raphson)法與L-M 法得到的一個(gè)周期內(nèi)誤差曲線(xiàn).不同于牛頓迭代法,L-M 法迭代格式為[34]

圖4 非線(xiàn)性方程組求解算法對(duì)計(jì)算精度的影響(相對(duì)論諧振子)Fig.4 Influence of nonlinear equation algorithm on calculation accuracy (relativistic harmonic oscillator)

其中xk為當(dāng)前迭代未知變量值,Fk為函數(shù)值,Jk為雅可比矩陣.L-M 法通過(guò)衡量每步迭代的誤差是否發(fā)散,來(lái)決定撤回迭代并將參數(shù) λk按十倍放縮.LM 法在求解非線(xiàn)性方程組時(shí)顯示出更高的計(jì)算精度,同比牛頓迭代法,精度提高了 104以上.

(2) 合理選擇代數(shù)方程

相對(duì)論諧振子(17)為保守系統(tǒng),振動(dòng)頻率 ω 是狀態(tài)變量[35].使用RHB-重構(gòu)法對(duì)n自由度保守系統(tǒng)求解時(shí),共需考慮n(2N+1)+1 個(gè)變量,需額外增加一個(gè)代數(shù)方程使RHB 方程組適定.本文以初值約束來(lái)探討方程對(duì)計(jì)算性能的影響.圖5 為采用55 階RHB-重構(gòu)法與L-M 法求解器,使用不同代數(shù)方程(約束條件)得到的誤差曲線(xiàn).分別對(duì)應(yīng): 初始位移約束x(0)=0;初始速度約束(0)=β;同時(shí)考慮初始位移與速度約束.同時(shí)考慮位移與速度約束時(shí),RHB-重構(gòu)法計(jì)算精度更高.

圖5 改變代數(shù)方程對(duì)計(jì)算精度影響Fig.5 Influence of changed algebraic equation on calculation accuracy

相對(duì)論諧振子在高速運(yùn)動(dòng)時(shí)顯示出復(fù)雜的動(dòng)力學(xué)特性,需要高階(＞20 階)方法來(lái)進(jìn)行分析.使用重鑄技術(shù),將根式非線(xiàn)性轉(zhuǎn)化為多項(xiàng)式型,再使用RHB法實(shí)現(xiàn)快速解算,克服高階估計(jì)的限制.如圖6 所示,RHB-重鑄法適用于不同初速度條件下相對(duì)論振子的計(jì)算.

3.2 非線(xiàn)性單擺

非線(xiàn)性單擺(19)是物理學(xué)中經(jīng)典問(wèn)題,且時(shí)域響應(yīng)具有典型的周期性.但現(xiàn)有HB 類(lèi)方法不能做到高性能求解: HB-Taylor 法需要采用高的截?cái)嚯A來(lái)保證計(jì)算精度;HB-AFT 法則需進(jìn)行過(guò)采樣以抑制混淆誤差(如圖7 所示).

對(duì)非線(xiàn)性單擺重鑄,得到方程(20).使用25 階RHB-重鑄法求解(大角度擺動(dòng),θ(0)=1.5),初值估計(jì)x1=1.423,s1=1.065,c0=1.028,輔助L-M 法求解器,計(jì)算結(jié)果如圖8 所示,與數(shù)值積分參考解比較,誤差控制在 10?12數(shù)量級(jí)以下.

圖8 25 階RHB-重鑄法求解非線(xiàn)性單擺Fig.8 The computing result of the RHB-recast method (N=25) for solving nonlinear pendulum

分別采用牛頓迭代法、Tikhonov 正則法與L-M法得到圖9 計(jì)算結(jié)果.牛頓迭代法中雅可比陣奇異,求解失效;Tikhonov 法與L-M 法計(jì)算所得的誤差數(shù)量級(jí)相近,但是由于Tikhonov 法對(duì)正則參數(shù)的優(yōu)化使求解更耗時(shí)[31].本例中,使用L-M 法計(jì)算RHB 法方程組僅需0.72 s,而Tikhonov 法耗時(shí)達(dá)1.17 s,L-M法較之Tikhonov 法效率提高了62.5%,達(dá)到了相近的求解精度.對(duì)比傳統(tǒng)方法(牛頓迭代法)與優(yōu)化方法(Tikhonov 法),L-M 法都更適合于RHB 法方程組的求解.

圖9 非線(xiàn)性方程組求解算法對(duì)計(jì)算精度的影響(非線(xiàn)性單擺)Fig.9 Influence of nonlinear equation algorithm on calculation accuracy (nonlinear pendulum)

非線(xiàn)性單擺問(wèn)題是保守系統(tǒng),其周期解 θ(t) 與響應(yīng)頻率都依賴(lài)于初始振幅 θ(0).使用RHB 法求解重鑄方程組(20)時(shí),本文給出3 種代數(shù)方程組合方案,以助于對(duì)比分析.

方案1: 對(duì)未知量 θ,v,s和c全諧波平衡(從常數(shù)項(xiàng)到N次諧波項(xiàng)),計(jì) 4(2N+1) 多個(gè)方程.增加初始振幅約束 θ(0)=α.

方案2: 對(duì)未知量 θ 和v全諧波平衡,計(jì)2(2N+1)多個(gè)方程.s和c從1 次諧波開(kāi)始,平衡系數(shù)到N次,計(jì) 2N多個(gè)方程.增加初始振幅約束 θ (0)=α 與兩個(gè)對(duì)附加變量s,c的初值約束

方案3: 對(duì)未知量 θ 和v全諧波展開(kāi),計(jì)2(2N+1)多個(gè)方程.s和c從1 次諧波展開(kāi)到N次,計(jì) 2N多個(gè)方程.增加初始速度約束v(0)=0 與兩個(gè)對(duì)附加變量s,c的初值約束式(24).

分別采用3 種方案得到的計(jì)算結(jié)果如圖10 所示.對(duì)比方案2 和方案3,方案1 的計(jì)算誤差大,主要原因是沒(méi)有對(duì)附加變量s和c的初值進(jìn)行約束.方案3 同時(shí)利用了初始位移與速度信息,比方案2 計(jì)算精度更高.

圖10 不同代數(shù)方程組合方案對(duì)計(jì)算精度影響Fig.10 Influence of different combination schemes of nonlinear algebraic equation on calculation accuracy

對(duì)比求解非線(xiàn)性單擺的3 種方法(HB-AFT法、RHB-Taylor 法和RHB-重鑄法)計(jì)算結(jié)果.考慮同階截?cái)郚=25,綜合圖11 誤差曲線(xiàn)和表3 的計(jì)算結(jié)果,RHB-重鑄法確定最優(yōu)時(shí)域配點(diǎn)數(shù)M=76,降低采樣成本.RHB-Taylor 法雖可采用高階截?cái)?15 次多項(xiàng)式)保證計(jì)算精度,但超越函數(shù)的級(jí)數(shù)表示降低計(jì)算效率,計(jì)算用時(shí)達(dá)1.80 s.RHB-重鑄法與HB-AFT 法計(jì)算效率相當(dāng),但將計(jì)算精度提高了兩個(gè)數(shù)量級(jí)以上.

表3 3 種方法計(jì)算非線(xiàn)性單擺結(jié)果對(duì)比Table 3 Comparison of corresponding results of three methods for solving nonlinear pendulum

圖11 3 種方法求解非線(xiàn)性單擺對(duì)應(yīng)誤差曲線(xiàn)對(duì)比Fig.11 Comparison of corresponding error curves of three methods for solving nonlinear pendulum

RHB-重鑄法適用于求解初等超越函數(shù)非線(xiàn)性問(wèn)題.尤其是當(dāng)單擺做大幅度振動(dòng)時(shí),傳統(tǒng)的線(xiàn)性化手段無(wú)法很好地捕捉動(dòng)力學(xué)響應(yīng),而RHB-重鑄法則可以提供高精度解析解(圖12).

圖12 RHB-重鑄法求解非線(xiàn)性單擺問(wèn)題: 相平面圖Fig.12 The RHB-recast method for solving nonlinear pendulum: phase plot

3.3 非線(xiàn)性耦合的非對(duì)稱(chēng)擺

用繩將質(zhì)點(diǎn)小球懸掛于一根固定在鐵架上的彈性桿末端,而彈性桿能在水平面與質(zhì)點(diǎn)小球同步振動(dòng).此物理模型在許多其他二維擺系統(tǒng)中普遍存在,如球面擺、傅科擺等.特別地,當(dāng)傅科擺由于不完美的懸掛或由于側(cè)向運(yùn)動(dòng)引起非線(xiàn)性耦合而導(dǎo)致不對(duì)稱(chēng),可能會(huì)導(dǎo)致額外的旋轉(zhuǎn),從而掩蓋地面效應(yīng).非線(xiàn)性耦合的非對(duì)稱(chēng)單擺微分形式可以寫(xiě)作[36]

此方程中2 階導(dǎo)數(shù)耦合入非線(xiàn)性項(xiàng),導(dǎo)致使用傳統(tǒng)方法進(jìn)行直接求解變得棘手.以傳統(tǒng)的有限差分類(lèi)方法(歐拉法、龍格庫(kù)塔法、MATLAB 內(nèi)置ode算法等)而言,需額外計(jì)算代數(shù)方程組來(lái)輔助求解.而部分解析求解法的計(jì)算效果同樣不佳,Jia 等[36]曾采用多時(shí)間尺度法推導(dǎo)方程的解為

令參數(shù) κ=0.01,初始條件x(0)=0.1,y(0)=0.2,(0)=(0)=0.

以修正Chebyshev-Picard 迭代 (modified Chebyshev-Picard iteration,MCPI)法計(jì)算結(jié)果為標(biāo)準(zhǔn)解[37-38].仿真得到如圖13 所示,解(26)不僅推導(dǎo)復(fù)雜,且與真實(shí)的物理解間有較大誤差.

圖13 多時(shí)間尺度法與修正Chebyshev-Picard 迭代法(參考解)求解軌跡對(duì)比Fig.13 Comparison of the method of multiple scales and the MCPI method (reference solution) for solving trajectory diagram

對(duì)多時(shí)間尺度解(26)分析得知,動(dòng)力學(xué)響應(yīng)包含兩個(gè)基頻(所有頻率成分都是基頻的線(xiàn)性組合),通過(guò)快速傅里葉變換(FFT)得到: ω1=0.985 7,ω2=0.993 5.又從微分方程重鑄的角度分析非對(duì)稱(chēng)擺微分方程(25),因?yàn)楦唠A導(dǎo)數(shù)只算作一個(gè)非線(xiàn)性度,該方程是具有3 次多項(xiàng)式非線(xiàn)性的動(dòng)力學(xué)系統(tǒng).使用5 階RMHB 法進(jìn)行計(jì)算得到圖14 所示仿真結(jié)果,計(jì)算用時(shí)8.24 s,x方向振動(dòng)幅值誤差 9.21×10?3,y方向振動(dòng)幅值誤差 3.90×10?7.

圖14 5 階RMHB 法求解非線(xiàn)性耦合非對(duì)稱(chēng)擺Fig.14 Solving nonlinear coupling asymmetric pendulum by the RMHB5

對(duì)2 階導(dǎo)耦合型非線(xiàn)性系統(tǒng),有限差分方法不能直接求解,一些解析求解方法的計(jì)算精度低.考慮兩個(gè)基頻的RMHB 法不僅保證計(jì)算效率,還能對(duì)擬周期運(yùn)動(dòng)進(jìn)行準(zhǔn)確的預(yù)測(cè).

4 結(jié)論

圍繞復(fù)雜非線(xiàn)性動(dòng)力學(xué)系統(tǒng)高效高精度周期解算需求,本文基于微分方程重鑄、“補(bǔ)頻”等思想,分別提出RHB-重鑄法與RMHB 法,主要的工作與結(jié)論總結(jié)如下.

(1)提出RHB-重鑄法,將一般非線(xiàn)性問(wèn)題改寫(xiě)為多項(xiàng)式型非線(xiàn)性方程,配合RHB 法確定消除非物理解所需的最優(yōu)配點(diǎn)閾值,實(shí)現(xiàn)對(duì)非多項(xiàng)式型非線(xiàn)性動(dòng)力學(xué)系統(tǒng)周期響應(yīng)進(jìn)行高階預(yù)測(cè).數(shù)值實(shí)驗(yàn)說(shuō)明,對(duì)比HB-AFT 法的時(shí)域過(guò)采樣和HB-Taylor 法對(duì)非線(xiàn)性函數(shù)進(jìn)行泰勒級(jí)數(shù)截?cái)鄡煞N處理方式,RHB-重鑄法的計(jì)算精度高達(dá) 10?12,且計(jì)算時(shí)間不超過(guò)1 s,同時(shí)保證了求解的高精度與高效率.

(2)提出結(jié)合“補(bǔ)頻”思想的RMHB 法,在原來(lái)單頻假設(shè)的基礎(chǔ)上,補(bǔ)充多個(gè)與物理響應(yīng)相關(guān)的頻率.借鑒RHB 法中時(shí)域重構(gòu)等價(jià)性,給出多頻諧波平衡計(jì)算的最優(yōu)配點(diǎn)數(shù)結(jié)論,并充分利用非線(xiàn)性量的時(shí)域采樣代替復(fù)雜的高維傅里葉分析,最終實(shí)現(xiàn)對(duì)擬周期響應(yīng)的準(zhǔn)確捕捉.

(3)通過(guò)解算相對(duì)論諧振子、非線(xiàn)性單擺以及非線(xiàn)性耦合的非對(duì)稱(chēng)擺問(wèn)題驗(yàn)證了本文算法的有效性.針對(duì)傳統(tǒng)有限差分類(lèi)方法求解困難的多自由度耦合系統(tǒng),RHB 及本文提出的兩種方法都不受方程形式的限制.此外,在相同的系統(tǒng)參數(shù)設(shè)置下,合理地選擇代數(shù)方程求解器、代數(shù)約束方案,可有效提高非線(xiàn)性問(wèn)題的求解精度.

本文提出的RHB-重鑄法與RMHB 法對(duì)復(fù)雜非線(xiàn)性動(dòng)力學(xué)系統(tǒng)周期性及擬周期響應(yīng)的計(jì)算效率與精度相較于現(xiàn)行的方法與處理方式都具有優(yōu)勢(shì).在后續(xù)的研究中,將探究RHB 法在高精周期響應(yīng)求解方面的工程應(yīng)用.

數(shù)據(jù)可用性聲明

支撐本研究的科學(xué)數(shù)據(jù)已在中國(guó)科學(xué)院科學(xué)數(shù)據(jù)銀行 (science data bank) ScienceDB 平臺(tái)公開(kāi)發(fā)布,訪(fǎng)問(wèn)地址為: https://cstr.cn/31253.11.sciencedb.j00140.00022 或https://doi.org/10.57760/sciencedb.j00140.00022.

附錄A

附表A1 常見(jiàn)初等超越函數(shù)的重鑄Table A1 Recast form of the most common elementary transcendental functions