一种求非线性振动系统周期解的切比雪夫级数方法

周 薇，韩景龙，陈全龙

近十几年，国内外学者在求非线性振动系统周期解方面进行了大量研究，不仅改进了传统算法［1－5］，而且还提出了一些新方法［6－11］。周期解的稳定性问题是非线性动力学中又一重要研究内容，通常用Floquet理论研究［12］。目前，针对非线性振动系统求解方法的研究虽然已经取得了重大进展，但仍面临不断改进方法适用范围、提高求解精度以及降低计算复杂度等方面的问题［13］;已有计算近似Floquet转移矩阵(FTM)的数值方法也普遍存在精度较低的局限［6］。

切比雪夫多项式是数值逼近领域中一类重要的基函数，它不但具有一致收敛性［14］，而且在精确逼近函数的同时能有效降低高阶差值的Runge现象［15］。虽然利用切比雪夫多项式求解常微分方程的研究起步较早［16］，但在最近十几年它才被应用于非线性周期系数微分系统的各项研究［10，17－24］。

鉴于利用15～18项移位的第一类切比雪夫级数可以非常精确地逼近三角函数［25］和高维系统的FTM(如20×20维)［24］，本文将移位的第一类切比雪夫级数理论和非线性优化算法结合，提出了一种求非线性振动系统周期解的方法。本方法将状态矢量中未知切比雪夫系数的求解，转化为对主周期上系统残差求最小值的无约束最优化问题，计算出了具有较高精度的切比雪夫级数周期解。以Duffing系统和直升机旋翼运动系统为例，通过与谐波平衡法(HBM)比较，验证了本方法正确有效，算例同时表明，将切比雪夫级数理论引入到直升机气弹响应与稳定性分析领域正确、可行。

1 切比雪夫级数周期解及其稳定性分析

考察非线性振动系统:

其中:X(t)={x1(t)，x2(t)，…xn(t)}T是 n×1维列向量，f(X(t)，t)是周期为T的函数。

在非线性动力学中，与平衡点的确定相比，确定周期运动的存在性和数目是更加困难的问题。通常仅求解所关心的一部分周期运动。设系统(1)的周期解为:

将其解展开为m项移位的第一类切比雪夫级数，具有如下形式

为了在今后使用局部优化方法时能够合理估计优化初值，同时将周期解(2)展开为N阶谐波级数:

其中:ai，k和 bi，k是周期解中(2N+1)·n 个未知系数，i=1，2，…，n。对于单周期轨道，取j=1;对于倍周期轨道，取j=2n，n为周期倍化分岔的次数。以单周期轨道为例，此时(3)～(5)满足如下关系:

因为移位的第一类切比雪夫多项式的定义区间为［0，1］，所以必须对系统进行周期归一化，即做变换t=T×s且s∈［0，1］。在实际计算时，只需将系统方程展开为切比雪夫级数，并将所有ω代换为2π，t代换为s，取为关于s的多项式即可。将式(1)右端移项后，求得系统残差R(s)，用符号表示如下:

其中:［F(ai，k，bi，k)］是由系统方程决定的 n m × 1 维未知列向量，且每个行分量是关于 ai，k和 bi，k的多项式。(s)］T，⊗代表Kronecker乘法，I为单位矩阵。对精确解而言，在一个周期上任意时刻s0∈［0，1］，均有残差R(s0)=0。设Ri(s)为残差R(s)每一行分量，则所设周期解与精确解误差越小，等价于在周期［0，1］上任意时间点处Ri(s)的绝对值越小，表达成非线性最优化问题为:

利用优化算法可以求得未知系数 ai，k和 bi，k，由式(6)求得pi，j即为周期解的切比雪夫系数。

不妨采用局部优化方法中的拟牛顿法［26］来求解式(9)。由于初值选取对优化结果有影响，通过调整式(6)中谐波展开项数(也可直接参照谐波平衡法选取谐波项数)，使优化初值处于合理范围，从而避免盲目试值。若能够利用其它方法合理估计优化初值，则可以不用式(6)和(7)中的谐波展开，直接求解切比雪夫系数pi，j。特别地，一些工程问题的周期解(如算例2.2.1)虽然能够估计出一个合理的可行域闭区间，从而利用确定性方法求全局最优解，但却会大大增加时间复杂度和计算复杂度，有时并不必要。

下面分析周期解的稳定性。设X0(t)是已求得的切比雪夫级数解，ΔX(t)为小扰动，令

把式(10)代入原系统方程，略去高阶小量，最后得到以Δx(t)为未知量的线化方程

由线性周期系统理论和切比雪夫级数运算性质，只需将单位矩阵的每一列分别作为初始条件，通过积分求得一个周期末点ΔX(T)的值，即为FTM相应列的向量［25］。根据Floquet理论，若FTM特征值的模均小于1，则周期解渐近稳定性;否则不稳定。

计算系统残差、目标函数和FTM时需用到移位的第一类切比雪夫多项式及其乘法、积分算子矩阵，详见文献［10，25］。本文在计算时均取前18项移位的第一类切比雪夫多项式。

2 算例分析

2.1 Duffing系统

以具有立方非线性的Duffing系统方程为例

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

在求切比雪夫级数周期解时，先对变量t做变换t=2πs/ω (s∈［0，1］)，将时间周期归一化，再计算系统残差R。当取3阶谐波时，利用谐波平衡法和优化方法分别求得系统周期解中的未知系数如表1所示。

表1 Duffing系统周期解系数Tab.1 Periodic solution coefficients of Duffing system

图1 Duffing系统相图(取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阶谐波时，优化方法和谐波平衡法所得到的系统相轨迹图，二者吻合良好。将系统代入谐波平衡法求得的周期解后，系统绝对残差的均值为0.516;而代入切比雪夫级数周期解，系统绝对残差的均值为0.454，精度提高13.6%。图2为此时系统在一个周期上的残差曲线。同取6阶谐波展开时，利用谐波平衡法所得系统绝对残差的均值为0.022，而优化方法所得系统绝对残差的均值为1.22×10－5。二者在一个周期上的残差曲线分别如图3、图4所示。可见，通过调整谐波展开项数来改变优化算法初值可以降低残差量级，有效保证求解精度。

图4 优化方法残差曲线Fig.4 Residual curve of optimization method

下面进行稳定性分析。由线性周期系统理论和切比雪夫级数运算性质，经过计算，存在一个Floquet乘子的模为27.271，周期解不稳定。

2.2 直升机旋翼系统

旋翼响应与稳定性问题，是直升机动力学中一个重要研究课题。旋翼动力学模型是包含了结构、惯性和气动载荷非线性的时变常微分方程组。通常采用谐波平衡法、时间有限元法等求解。

2.2.1 铰接式旋翼系统

以3片桨叶的铰接式直升机旋翼系统运动为例，考虑旋转平面内如图5、图6所示的挥舞、摆振运动。假设旋翼每片桨叶性质相同，建立单片桨叶挥舞、摆振耦合运动方程:

图5 挥舞运动受力示意图Fig.5 Force diagram of flapping motion

图6 摆振运动受力示意图Fig.6 Force diagram of lagging motion

其中:Fx、Fz分别为平行和垂直于桨盘平面的剖面气动力，其详细表达式以及式中符号含义参见文献［27］。计算中用到的旋翼参数如表2所示。

表2 主要参数Tab.2 Main parameter

当前进比为0.2时，采用前飞条件下的Dress入流模型，总距 θ0=0.116 rad，周期变距 θC=0.035 rad，θS=－0.064 rad。将无量纲化后的参数代入系统方程，再进行周期归一化，最后得到各方程的残差Ri。在展开为2阶谐波时，谐波平衡法所得挥舞、摆振运动方程绝对残差的均值分别为1.24和0.241，由于略去谐波的量级较大，产生了错误结果。而代入优化方法所得周期解，系统绝对残差的均值分别为0.53和0.227，结果较准确。这是由于本方法所取谐波仅用于修正优化初值，而周期解仍由18项移位的第一类切比雪夫级数逼近。图7，图8显示了取3阶谐波时，两种方法所得桨叶运动的残差曲线，本方法精度提高5.2%。在同展开为7阶谐波时，谐波平衡法和优化方法所得的系统绝对残差的均值均达到10－6量级，具有较高精度，此时系统相图如图9，图10所示。

最后分析周期解的稳定性。用切比雪夫级数的积分运算来求得 FTM，Floquet乘子的模分别为0.096 3和0.680 2。此直升机旋翼系统在前进比为0.2时的挥舞、摆振周期运动均渐近稳定。

2.2.2 无铰式旋翼系统

在旋转坐标系下，以有限元法［28］建立桨叶运动方程，如图11所示。

所得桨叶方程用符号表示为［28］:

其中:q(t)是单片桨叶上总节点位移向量，q(t)={u1，v1，v1'，w1，w1'，φ1，…u14，φ10，u15}，u、v、w、φ 分别为桨叶弹性轴上各节点拉伸、摆振、挥舞、扭转自由度的弹性位移。算例采用文献［28］中BO105旋翼参数、准定常气动力和Drees入流模型。

图7 挥舞方程残差曲线Fig.7 Flapping equation residual curve

图8 摆振方程残差曲线Fig.8 Lagging equation residual curve

图9 挥舞相图Fig.9 Flapping phase graph

当前进比为0.1时，将各参数无量纲化后代入方程(14)。为降低方程维数、节省求解时间，计算时取桨叶前六阶固有模态，将式(14)转化为模态坐标方程，整理成状态方程形式如下:

图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 桨尖扭转相图Fig.14 Blade － tip torsion phase graph

将系统(15)在周期解p(t)处施加小扰动，求得线化系统 Floquet乘子的模分别为 0.441，0.093，0.192，0.576，0.107 和0.153，均小于 1，所以该旋翼系统在前进比为0.1时的周期响应渐近稳定。

3 结论

本文提出了一种求解非线性振动系统周期解的切比雪夫级数方法，将未知状态矢量用切比雪夫级数表示，相当于考虑了系统中高阶谐波的影响，因而求解时能够获得较高精度。方法中的优化问题还可等价于直接求解切比雪夫系数的含约束非线性最优化问题，但在相同条件下会增加计算复杂度。通过计算刚性铰接式与弹性无铰式旋翼的周期运动，表明将切比雪夫级数理论引入直升机气弹响应和稳定性研究正确、可行，在求得的周期解误差较小时更利于获得系统较准确的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］.振动与冲击，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］张琪昌，赵蔷薇，王炜.强非线性振动系统的通用化求解程序及应用［J］.振动与冲击，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］李银山，李树杰.构造一类非线性振子解析逼近周期解的初值变换法［J］.振动与冲击，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］胡海岩.应用非线性动力学［M］.北京:航空工业出版社，2000.

［13］陈洋洋，燕乐纬，佘锦炎，等.非线性自治振动系统同宿解的广义双曲函数摄动法［J］.应用数学和力学，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］袁亚湘.非线性优化计算方法［M］.北京:科学出版社，2008.

［27］约翰逊W.直升机理论［M］.北京:航空工业出版社，1991.

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