周 薇，韓景龍，陳全龍

近十幾年，國內(nèi)外學(xué)者在求非線性振動(dòng)系統(tǒng)周期解方面進(jìn)行了大量研究，不僅改進(jìn)了傳統(tǒng)算法［1－5］，而且還提出了一些新方法［6－11］。周期解的穩(wěn)定性問題是非線性動(dòng)力學(xué)中又一重要研究內(nèi)容，通常用Floquet理論研究［12］。目前，針對非線性振動(dòng)系統(tǒng)求解方法的研究雖然已經(jīng)取得了重大進(jìn)展，但仍面臨不斷改進(jìn)方法適用范圍、提高求解精度以及降低計(jì)算復(fù)雜度等方面的問題［13］;已有計(jì)算近似Floquet轉(zhuǎn)移矩陣(FTM)的數(shù)值方法也普遍存在精度較低的局限［6］。

切比雪夫多項(xiàng)式是數(shù)值逼近領(lǐng)域中一類重要的基函數(shù)，它不但具有一致收斂性［14］，而且在精確逼近函數(shù)的同時(shí)能有效降低高階差值的Runge現(xiàn)象［15］。雖然利用切比雪夫多項(xiàng)式求解常微分方程的研究起步較早［16］，但在最近十幾年它才被應(yīng)用于非線性周期系數(shù)微分系統(tǒng)的各項(xiàng)研究［10，17－24］。

鑒于利用15～18項(xiàng)移位的第一類切比雪夫級數(shù)可以非常精確地逼近三角函數(shù)［25］和高維系統(tǒng)的FTM(如20×20維)［24］，本文將移位的第一類切比雪夫級數(shù)理論和非線性優(yōu)化算法結(jié)合，提出了一種求非線性振動(dòng)系統(tǒng)周期解的方法。本方法將狀態(tài)矢量中未知切比雪夫系數(shù)的求解，轉(zhuǎn)化為對主周期上系統(tǒng)殘差求最小值的無約束最優(yōu)化問題，計(jì)算出了具有較高精度的切比雪夫級數(shù)周期解。以Duffing系統(tǒng)和直升機(jī)旋翼運(yùn)動(dòng)系統(tǒng)為例，通過與諧波平衡法(HBM)比較，驗(yàn)證了本方法正確有效，算例同時(shí)表明，將切比雪夫級數(shù)理論引入到直升機(jī)氣彈響應(yīng)與穩(wěn)定性分析領(lǐng)域正確、可行。

1 切比雪夫級數(shù)周期解及其穩(wěn)定性分析

考察非線性振動(dòng)系統(tǒng):

其中:X(t)={x1(t)，x2(t)，…xn(t)}T是 n×1維列向量，f(X(t)，t)是周期為T的函數(shù)。

在非線性動(dòng)力學(xué)中，與平衡點(diǎn)的確定相比，確定周期運(yùn)動(dòng)的存在性和數(shù)目是更加困難的問題。通常僅求解所關(guān)心的一部分周期運(yùn)動(dòng)。設(shè)系統(tǒng)(1)的周期解為:

將其解展開為m項(xiàng)移位的第一類切比雪夫級數(shù)，具有如下形式

為了在今后使用局部優(yōu)化方法時(shí)能夠合理估計(jì)優(yōu)化初值，同時(shí)將周期解(2)展開為N階諧波級數(shù):

其中:ai，k和 bi，k是周期解中(2N+1)·n 個(gè)未知系數(shù)，i=1，2，…，n。對于單周期軌道，取j=1;對于倍周期軌道，取j=2n，n為周期倍化分岔的次數(shù)。以單周期軌道為例，此時(shí)(3)～(5)滿足如下關(guān)系:

因?yàn)橐莆坏牡谝活惽斜妊┓蚨囗?xiàng)式的定義區(qū)間為［0，1］，所以必須對系統(tǒng)進(jìn)行周期歸一化，即做變換t=T×s且s∈［0，1］。在實(shí)際計(jì)算時(shí)，只需將系統(tǒng)方程展開為切比雪夫級數(shù)，并將所有ω代換為2π，t代換為s，取為關(guān)于s的多項(xiàng)式即可。將式(1)右端移項(xiàng)后，求得系統(tǒng)殘差R(s)，用符號表示如下:

其中:［F(ai，k，bi，k)］是由系統(tǒng)方程決定的 n m × 1 維未知列向量，且每個(gè)行分量是關(guān)于 ai，k和 bi，k的多項(xiàng)式。(s)］T，?代表Kronecker乘法，I為單位矩陣。對精確解而言，在一個(gè)周期上任意時(shí)刻s0∈［0，1］，均有殘差R(s0)=0。設(shè)Ri(s)為殘差R(s)每一行分量，則所設(shè)周期解與精確解誤差越小，等價(jià)于在周期［0，1］上任意時(shí)間點(diǎn)處Ri(s)的絕對值越小，表達(dá)成非線性最優(yōu)化問題為:

利用優(yōu)化算法可以求得未知系數(shù) ai，k和 bi，k，由式(6)求得pi，j即為周期解的切比雪夫系數(shù)。

不妨采用局部優(yōu)化方法中的擬牛頓法［26］來求解式(9)。由于初值選取對優(yōu)化結(jié)果有影響，通過調(diào)整式(6)中諧波展開項(xiàng)數(shù)(也可直接參照諧波平衡法選取諧波項(xiàng)數(shù))，使優(yōu)化初值處于合理范圍，從而避免盲目試值。若能夠利用其它方法合理估計(jì)優(yōu)化初值，則可以不用式(6)和(7)中的諧波展開，直接求解切比雪夫系數(shù)pi，j。特別地，一些工程問題的周期解(如算例2.2.1)雖然能夠估計(jì)出一個(gè)合理的可行域閉區(qū)間，從而利用確定性方法求全局最優(yōu)解，但卻會(huì)大大增加時(shí)間復(fù)雜度和計(jì)算復(fù)雜度，有時(shí)并不必要。

下面分析周期解的穩(wěn)定性。設(shè)X0(t)是已求得的切比雪夫級數(shù)解，ΔX(t)為小擾動(dòng)，令

把式(10)代入原系統(tǒng)方程，略去高階小量，最后得到以Δx(t)為未知量的線化方程

由線性周期系統(tǒng)理論和切比雪夫級數(shù)運(yùn)算性質(zhì)，只需將單位矩陣的每一列分別作為初始條件，通過積分求得一個(gè)周期末點(diǎn)ΔX(T)的值，即為FTM相應(yīng)列的向量［25］。根據(jù)Floquet理論，若FTM特征值的模均小于1，則周期解漸近穩(wěn)定性;否則不穩(wěn)定。

計(jì)算系統(tǒng)殘差、目標(biāo)函數(shù)和FTM時(shí)需用到移位的第一類切比雪夫多項(xiàng)式及其乘法、積分算子矩陣，詳見文獻(xiàn)［10，25］。本文在計(jì)算時(shí)均取前18項(xiàng)移位的第一類切比雪夫多項(xiàng)式。

2 算例分析

2.1 Duffing系統(tǒng)

以具有立方非線性的Duffing系統(tǒng)方程為例

其中:a=5，b=2，c=18，d=8，ω =2。

在求切比雪夫級數(shù)周期解時(shí)，先對變量t做變換t=2πs/ω (s∈［0，1］)，將時(shí)間周期歸一化，再計(jì)算系統(tǒng)殘差R。當(dāng)取3階諧波時(shí)，利用諧波平衡法和優(yōu)化方法分別求得系統(tǒng)周期解中的未知系數(shù)如表1所示。

表1 Duffing系統(tǒng)周期解系數(shù)Tab.1 Periodic solution coefficients of Duffing system

圖1 Duffing系統(tǒng)相圖(取3階諧波)Fig.1 Phase diagram of Duffing system(take 3－order harmonics)

圖2 殘差曲線(取3階諧波)Fig.2 Residual curves(take 3－order harmonics)

圖3 諧波平衡法殘差曲線Fig.3 Residual curve of HBM

圖1 中給出了同取3階諧波時(shí)，優(yōu)化方法和諧波平衡法所得到的系統(tǒng)相軌跡圖，二者吻合良好。將系統(tǒng)代入諧波平衡法求得的周期解后，系統(tǒng)絕對殘差的均值為0.516;而代入切比雪夫級數(shù)周期解，系統(tǒng)絕對殘差的均值為0.454，精度提高13.6%。圖2為此時(shí)系統(tǒng)在一個(gè)周期上的殘差曲線。同取6階諧波展開時(shí)，利用諧波平衡法所得系統(tǒng)絕對殘差的均值為0.022，而優(yōu)化方法所得系統(tǒng)絕對殘差的均值為1.22×10－5。二者在一個(gè)周期上的殘差曲線分別如圖3、圖4所示。可見，通過調(diào)整諧波展開項(xiàng)數(shù)來改變優(yōu)化算法初值可以降低殘差量級，有效保證求解精度。

圖4 優(yōu)化方法殘差曲線Fig.4 Residual curve of optimization method

下面進(jìn)行穩(wěn)定性分析。由線性周期系統(tǒng)理論和切比雪夫級數(shù)運(yùn)算性質(zhì)，經(jīng)過計(jì)算，存在一個(gè)Floquet乘子的模為27.271，周期解不穩(wěn)定。

2.2 直升機(jī)旋翼系統(tǒng)

旋翼響應(yīng)與穩(wěn)定性問題，是直升機(jī)動(dòng)力學(xué)中一個(gè)重要研究課題。旋翼動(dòng)力學(xué)模型是包含了結(jié)構(gòu)、慣性和氣動(dòng)載荷非線性的時(shí)變常微分方程組。通常采用諧波平衡法、時(shí)間有限元法等求解。

2.2.1 鉸接式旋翼系統(tǒng)

以3片槳葉的鉸接式直升機(jī)旋翼系統(tǒng)運(yùn)動(dòng)為例，考慮旋轉(zhuǎn)平面內(nèi)如圖5、圖6所示的揮舞、擺振運(yùn)動(dòng)。假設(shè)旋翼每片槳葉性質(zhì)相同，建立單片槳葉揮舞、擺振耦合運(yùn)動(dòng)方程:

圖5 揮舞運(yùn)動(dòng)受力示意圖Fig.5 Force diagram of flapping motion

圖6 擺振運(yùn)動(dòng)受力示意圖Fig.6 Force diagram of lagging motion

其中:Fx、Fz分別為平行和垂直于槳盤平面的剖面氣動(dòng)力，其詳細(xì)表達(dá)式以及式中符號含義參見文獻(xiàn)［27］。計(jì)算中用到的旋翼參數(shù)如表2所示。

表2 主要參數(shù)Tab.2 Main parameter

當(dāng)前進(jìn)比為0.2時(shí)，采用前飛條件下的Dress入流模型，總距 θ0=0.116 rad，周期變距 θC=0.035 rad，θS=－0.064 rad。將無量綱化后的參數(shù)代入系統(tǒng)方程，再進(jìn)行周期歸一化，最后得到各方程的殘差Ri。在展開為2階諧波時(shí)，諧波平衡法所得揮舞、擺振運(yùn)動(dòng)方程絕對殘差的均值分別為1.24和0.241，由于略去諧波的量級較大，產(chǎn)生了錯(cuò)誤結(jié)果。而代入優(yōu)化方法所得周期解，系統(tǒng)絕對殘差的均值分別為0.53和0.227，結(jié)果較準(zhǔn)確。這是由于本方法所取諧波僅用于修正優(yōu)化初值，而周期解仍由18項(xiàng)移位的第一類切比雪夫級數(shù)逼近。圖7，圖8顯示了取3階諧波時(shí)，兩種方法所得槳葉運(yùn)動(dòng)的殘差曲線，本方法精度提高5.2%。在同展開為7階諧波時(shí)，諧波平衡法和優(yōu)化方法所得的系統(tǒng)絕對殘差的均值均達(dá)到10－6量級，具有較高精度，此時(shí)系統(tǒng)相圖如圖9，圖10所示。

最后分析周期解的穩(wěn)定性。用切比雪夫級數(shù)的積分運(yùn)算來求得 FTM，F(xiàn)loquet乘子的模分別為0.096 3和0.680 2。此直升機(jī)旋翼系統(tǒng)在前進(jìn)比為0.2時(shí)的揮舞、擺振周期運(yùn)動(dòng)均漸近穩(wěn)定。

2.2.2 無鉸式旋翼系統(tǒng)

在旋轉(zhuǎn)坐標(biāo)系下，以有限元法［28］建立槳葉運(yùn)動(dòng)方程，如圖11所示。

所得槳葉方程用符號表示為［28］:

其中:q(t)是單片槳葉上總節(jié)點(diǎn)位移向量，q(t)={u1，v1，v1'，w1，w1'，φ1，…u14，φ10，u15}，u、v、w、φ 分別為槳葉彈性軸上各節(jié)點(diǎn)拉伸、擺振、揮舞、扭轉(zhuǎn)自由度的彈性位移。算例采用文獻(xiàn)［28］中BO105旋翼參數(shù)、準(zhǔn)定常氣動(dòng)力和Drees入流模型。

圖7 揮舞方程殘差曲線Fig.7 Flapping equation residual curve

圖8 擺振方程殘差曲線Fig.8 Lagging equation residual curve

圖9 揮舞相圖Fig.9 Flapping phase graph

當(dāng)前進(jìn)比為0.1時(shí)，將各參數(shù)無量綱化后代入方程(14)。為降低方程維數(shù)、節(jié)省求解時(shí)間，計(jì)算時(shí)取槳葉前六階固有模態(tài)，將式(14)轉(zhuǎn)化為模態(tài)坐標(biāo)方程，整理成狀態(tài)方程形式如下:

圖10 擺振相圖Fig.10 Lagging phase graph

圖11 槳葉的空間有限元Fig.11 Blade space finite element

圖12 槳尖揮舞相圖Fig.12 Blade-tip flapping phase graph

圖13 槳尖擺振相圖Fig.13 Blade-tip lagging phase graph

圖14 槳尖扭轉(zhuǎn)相圖Fig.14 Blade － tip torsion phase graph

將系統(tǒng)(15)在周期解p(t)處施加小擾動(dòng)，求得線化系統(tǒng) Floquet乘子的模分別為 0.441，0.093，0.192，0.576，0.107 和0.153，均小于 1，所以該旋翼系統(tǒng)在前進(jìn)比為0.1時(shí)的周期響應(yīng)漸近穩(wěn)定。

3 結(jié)論

本文提出了一種求解非線性振動(dòng)系統(tǒng)周期解的切比雪夫級數(shù)方法，將未知狀態(tài)矢量用切比雪夫級數(shù)表示，相當(dāng)于考慮了系統(tǒng)中高階諧波的影響，因而求解時(shí)能夠獲得較高精度。方法中的優(yōu)化問題還可等價(jià)于直接求解切比雪夫系數(shù)的含約束非線性最優(yōu)化問題，但在相同條件下會(huì)增加計(jì)算復(fù)雜度。通過計(jì)算剛性鉸接式與彈性無鉸式旋翼的周期運(yùn)動(dòng)，表明將切比雪夫級數(shù)理論引入直升機(jī)氣彈響應(yīng)和穩(wěn)定性研究正確、可行，在求得的周期解誤差較小時(shí)更利于獲得系統(tǒng)較準(zhǔn)確的FTM。

［1］Dunne J F，Hayward P.A split-frequency harmonic balance method for nonlinear oscillators with multi-harmonic forcing［J］.Journal of Sound and Vibration，2006，295:939－963.

［2］Chen Y M，Liu J K.A new method based on the harmonic balance method for nonlinear oscillators［J］.Physics Letters A，2007，368:371－378.

［3］ Cochelin B，Vergez C.A high order purely frequency-based harmonic balance formulation for continuation of periodic solutions［J］.Journal of Sound and Vibration，2009，324:243－262.

［4］ Grolet A，Thouverez F.On a new harmonic selection technique for harmonic balance method［J］.Mechanical Systems and Signal Processing，2012，30:43－60.

［5］Lu C J，Lin Y M.A modified incremental harmonic balance method for rotary periodic motions［J］.Nonlinear Dynamic，2011，66:781－788.

［6］唐進(jìn)元，陳思雨.阻尼和剛度項(xiàng)含時(shí)變參數(shù)的強(qiáng)非線性振動(dòng)系統(tǒng)周期解研究［J］.振動(dòng)與沖擊，2007，26(10):96－100.TANG Jin-yuan，CHEN Si-yu.Study on periodic solutions of strong nonlinear systems with time-varying damping and stiffness coefficient［J］.Journal of Vibration and Shock，2007，26(10):96－100.

［7］ Hu H，Tang J S.A convolution integral method for certain strongly nonlinear oscillations［J］.Journal of Sound and Vibration，2005，285(45):1235－1241.

［8］張琪昌，趙薔薇，王煒.強(qiáng)非線性振動(dòng)系統(tǒng)的通用化求解程序及應(yīng)用［J］.振動(dòng)與沖擊，2012，31(8):1－4.ZHANGQi-chang，ZHAO Qiang-wei，WANGWei.Universal solving program and its application in a strongly nonlinear oscillation system［J］.Journal of Vibration and Shock，2012，31(8):1－4.

［9］李銀山，李樹杰.構(gòu)造一類非線性振子解析逼近周期解的初值變換法［J］.振動(dòng)與沖擊，2010，29(8):99－102.LI Yin-shan， LI Shu-jie. Method of initial-value transformation for obtaining approximate analytic periods of a class of nonlinear oscillator［J］.Journal of Vibration and Shock，2010，29(8):99－102.

［10］ Zhou T，Xu JX.Research on the periodic orbit of nonlinear dynamic systems using Chebyshev polynomials［J］.Journal of Sound and Vibration，2001，245(2):239－250.

［11］Feng Z X，Xu X，Ji S G.Finding the periodic solution of differential equation via solving optimization problem ［J］.J Optim Theory Appl，2009，143:75 －86.

［12］胡海巖.應(yīng)用非線性動(dòng)力學(xué)［M］.北京:航空工業(yè)出版社，2000.

［13］陳洋洋，燕樂緯，佘錦炎，等.非線性自治振動(dòng)系統(tǒng)同宿解的廣義雙曲函數(shù)攝動(dòng)法［J］.應(yīng)用數(shù)學(xué)和力學(xué)，2012，33(9):1064－1077.CHEN Yang-yang， YAN Le-wei， SZE Kam-yim， et al.Generalized hyperbolic perturbation method for homoclinic solutions of strongly nonlinear autonomous systems［J］.Applied Mathematics and Mechanics，2012，33(9):1064－1077.

［14］ Johnson L W，Riess R D.On the convergence of polynomials interpolating at the zeros of Tn(x)［J］.Math Z 1970，116(4):355－358.

［15］ Berrut J，Trefethen L N.Barycentric Lagrange interpolation.SIAM Rev，2004，46(3):501 －507.

［16］ Clenshaw C.The numerical solution of linear differential equations in Chebyshev series［J］.Proc Camb Philos Soc，1957，53:134－149.

［17］ Redkar S，Sinha S C.Reduced order modeling of nonlinear time periodic systems subjected to external periodic excitations［J］.Commun Nonlinear Sci Numer Simulat，2011，16:4120－4133.

［18］ Khasawneh F A，Mann B P，Butcher E A.A multi-interval Chebyshev collocation approach for the stability of periodic delay systems with discontinuities［J］.Commun Nonlinear Sci Numer Simulat，2011，16:4408 －4421.［19］ Celik I，Gokmen G.Approximate solution of periodic Sturm-Liouville problems with Chebyshev collocation method［J］.Applied Mathematics and Computation，2005，170:285－295.

［20］ Celik V.Collocation method and residual correction using Chebyshev series［J］. Applied Mathematics and Computation，2006，174:910－920.

［21］ Butcher E A，Bobrenkov O A.On the Chebyshev spectral continuous time approximation for constant and periodic delay differential equations［J］.Commun Nonlinear Sci Numer Simulat，2011，16:1541 －1554.

［22］Sedaghat S，Ordokhani Y，Dehghan M.Numerical solution of the delay differential equations of pantograph type via Chebyshev polynomials［J］.Commun Nonlinear Sci Numer Simulat，2012，17:4815 －4830.

［23］ Sinha S C，Pandiyan R.Analysis of quasilinear dynamical systems with periodic coefficients via Liapunov-Floquet transformation［J］.Int.J.Non-Linear Mechanics，1994，29(5):687－702.

［24］ Pandiyan R，Sinha S C.Analysis of time-periodic nonlinear dynamical systems undergoing bifurcations［J］.Nonlinear Dynamics，1995，8:21 －43.

［25］Sinha SC，Wu D H.An efficient computational scheme for the analysis of periodic systems［J］.Journal of Sound and Vibration，1991，151(1):91 －117.

［26］袁亞湘.非線性優(yōu)化計(jì)算方法［M］.北京:科學(xué)出版社，2008.

［27］約翰遜W.直升機(jī)理論［M］.北京:航空工業(yè)出版社，1991.

［28］ Gunji B， Chopra T. University of Maryland advanced rotorcraft code(UMARC)theory manual［R］.UMAERO Report，1994.