随机与谐和联合激励下分数阶非线性系统的统计线性化方法

孔凡 晁盼盼 徐军

摘要： 提出了一种计算随机与諧和联合激励下非线性分数阶系统响应二阶矩的统计线性化方法。假定位移响应可写为确定性均值和零均值随机分量之和的形式，原运动微分方程可化为关于均值分量的确定性微分方程和关于随机分量的随机微分方程组合。分别利用谐波平衡法和统计线性化方法对上述两类方程求解后，可得响应的确定性均值与随机分量。Monte Carlo模拟证实了该方法的有效性。

关键词： 非线性系统; 随机与谐和联合激励; 谐波平衡法; 统计线性化; 分数阶导数

引  言

分数阶导数模型在工程和科学问题中得到了广泛应用，如流变学[1]、扩散传输[2]和黏弹性材料力学模型等。对于黏弹性材料力学模型，Nutting等[3]在本构方程中首次使用分数阶导数的概念;Caputo[4]提出利用分数阶导数模拟地质地层的黏弹性行为;Slonimsky[5]，Smit和Vries[6]提出黏弹性介质应力与应变之间存在类似的分数阶微积分关系;Bagley和Torvik[7?8]不仅将分数阶导数模型应用于某些黏弹性材料（聚合物溶液和无交联聚合物固体），而且还证明分数阶微积分模型与材料物理原理是一致的。研究表明，黏弹性材料的分数阶导数模型具有简洁、紧凑的特点，只需少量参数就能在很宽的频率范围内描述黏弹性材料力学特性的频率依赖行为。

动力系统包含分数阶导数单元时，虽然简洁性使分数阶动力模型具有一定优势，但在解析或数值求解动力响应时却往往涉及比较复杂的计算。目前，分数阶确定动力系统分析可采用的方法比较丰富，如Laplace变换[9]、傅里叶变换[10]、特征向量展开[11]、平均方法[12]等。然而，分数阶随机动力系统研究起步较晚：Mainardi[13]建立描述布朗运动的分数阶Langevin方程，并提出用Laplace变换对其进行求解;Spanos和Zeldin[14]提出分数阶阻尼系统随机振动分析的频域方法;Agrawal[15]利用Suarez和Shokooh[11]的特征向量展开法，得到具有1/2阶阻尼的随机动力系统解析解;Ye等[16]提出附加黏弹性阻尼器的单自由度结构随机地震反应分析的傅里叶方法，得到分数阶系统的单位脉冲响应函数和Duhamel积分表达式;黄志龙等[17]提出强非线性单自由度分数阶随机系统的随机平均法，考察系统响应及其稳定性;孙春燕[18]利用几种典型的动力学方法分析分数阶随机时滞系统的动力响应。

然而，在实际工程中还存在同时受到随机激励和确定性周期荷载联合作用的一类结构（或装置）。例如，风浪荷载作用下的风力发电塔、地震中支撑旋转装置的结构[19]、受到周期性水流冲击的大坝[20]、飞行过程中直升机的机翼叶片[21]等。因此，考察随机与谐和联合激励下系统的动力响应具有重要的工程意义。目前，大量研究致力于联合激励下的整数阶非线性系统：Iyengar[22]利用高斯矩截断方法研究正弦噪声和白噪声联合激励下的Duffing振子，得到系统的多个稳态解并进行解的稳定性分析;Nayfeh和Serhan[23]结合多尺度和二阶矩截断方法，考察Duffing?Rayleigh振子在联合激励下响应的均值和均方值;Rong等[24]采用谐波平衡法和随机平均法研究联合激励下Duffing振子的响应、多稳态解和跳跃现象;Rong等[25]利用多尺度法确定van der Pol?Duffing振子在联合激励下响应幅值和相位的调制方程，得到响应稳态解并考察解的稳定性;Anh和Hieu[26]考察联合激励下的整数阶Duffing振子，假定响应可分解为确定性和随机分量，继而将原运动微分方程分解为耦合的确定性和随机微分方程，随后通过确定平均法和等效线性化分别求解确定和随机方程。Spanos等[27]采用类似的假定，利用谐波平衡法和统计线性化方法分别求解确定性和随机响应。类似的研究亦可见文献[28]。分数阶动力系统在确定性谐波和随机联合激励下的随机动力响应研究较少，目前，只有Chen和Zhu[29?33]进行相关研究：他们将随机平均方法推广到联合激励下分数阶导数阻尼的Duffing振子中，利用FPK方程求解系统响应幅值和相位的概率密度函数，研究该系统的随机跳跃、分岔和稳定性。

本文提出一种用于计算联合激励下的分数阶非线系统响应的方法。首先，采用与文献[22，26?28]类似的假定，即系统响应可写为谐和均值过程与随机零均值过程之和的形式，将原系统运动方程化为耦合的确定性和随机微分方程;其次，利用谐波平衡法和统计线性化对上述两类方程耦合求解，可得到响应的确定与随机分量;最后，与时域数值模拟的对比证实了该方法的精度和计算效率。由于统计线性化方法的广泛适用性，本文所建议方法可方便地推广至多自由度、甚至滞回非线性系统。

1 联合激励下非线性分数阶系统的近似解

考虑随机与谐和联合激励下的具有分数阶导数阻尼的单自由度非线性系统

可见，联合激励下的分数阶导数系统能近似为谐和激励下的确定性动力系统（6）和随机激励下的动力系统（7）。求解这两类系统可得到式（1）的近似解。注意到，式（6）与（7）是相互耦合的，即式（6）包含随机响应分量特征值，式（7）中包含确定性响应分量，必须同时考虑它们才能实现对二者的求解。

2 随机和确定性分量求解

为此，本文采用迭代法求解上述方程组，具体步骤如下：

1） 确定随机响应方差初值。忽略非线性项，即令，可得分数阶线性系统在随机激励下的均方初值;将代入到方程（11）和（12）中，联立线性方程组求解谐和项初值和。

2） 将代入式（19）求解，并根据式（20）求解更新后的;

3） 将更新后的代入式（11）?（12）求解更新后的谐和分量幅值;

4） 重复步骤2）和3）直至达到相关收敛准则。

3 数值算例

3.1 简谐和白噪声联合激励下

考察隨机激励为白噪声的情况。归一化后的激励功率谱密度可以由式（35）和（36）给出。首先考虑非共振情形。作为演示算例，选择系统参数，，，;确定性激励参数，。本文所提方法计算得到响应功率谱密度和均方值与数值模拟方法所得结果对比验证了该方法的适用性。如图1和2所示，Monte Carlo模拟方法与本文所提方法得到的功率谱密度和稳态响应均方值吻合良好。其中，图1中的箭头表示确定性分量的功率谱密度。此外，Monte Carlo模拟得到的响应均方值为，统计线性化方法得到的响应均方值为，误差仅为1.1%。

为进一步验证所提方法的有效性，考虑其他几组系统和激励参数的情况。简谐激励具有不同幅值时，其频率与响应均方值之间的关系（其他参数均与图2所示的系统参数相同）如图3所示。图中，虚线代表统计线性化法的计算结果，标记线为时域数值模拟方法计算结果。统计线性化方法中，使用了线性初值启动了Newton迭代计算;时域数值模拟方法中假定系统初始状态为零。

由图3可知，多数情况下，统计线性化方法和数值模拟的结果吻合较好。在大幅值确定性激励作用下，曲线产生跳跃现象。确定性激励幅值越大，曲线的峰值越大且越呈明显不对称状态，左坡缓而右坡陡，产生的跳跃越剧烈。图3所示的均方值跳跃现象可能与谐波激励下分数阶Duffing振子的幅频曲线的跳跃现象有关，分析如下。

图4（a）?（d）给出了谐和激励单独作用下或联合激励共同作用下系统响应与谐和激励频率之间的关系。为对比方便，纵坐标为谐和激励单独作用下的确定性幅值或联合激励作用下的随机响应均方值。与图3不同的是，为得到图4（a）?（d）中的数值解，使用了前（从小到大扫频）或后（从大到小扫频）一个激励频率点的稳态响应作为下一个频率点的Newton迭代初始值（对于所建议方法的式（11）?（12））或时域初始条件（对于时域逐步积分法）以启动计算;以下将其称为扫频计算方式。图4（a）?（d）中，实线（红色）为谐和激励单独作用下响应的频域解析解（式（37）），菱形（绿色）标记线为利用Newton迭代得到的式（37）的频域数值解，二者均为确定性响应幅值。此外，三角标记线（黑色）为联合激励下所建议方法得到的联合响应均方值;圆形标记线（粉色）为所建议方法得到的确定性分量幅值;星形标记线（蓝色）为时域数值模拟得到的联合响应均方值。

由图4（a）?（d）可得如下结论。首先，对比实线和菱形标记线可知，谐和激励幅值越大，幅频曲线产生跳跃的三值区间越大。其中，左支和右支均为稳定解，可通过扫频数值计算方式得到，与解析解吻合非常好;中支为非稳定解，无法通过一般数值算法得到。第二，对比圆形标记线和实线（菱形标记线）可知，随机激励使确定性响应幅频曲线产生偏移，且随机激励相对谐和激励越大，向高频偏移量越大。文献[26]也得到了类似的结论。第三，谐和激励相对随机激励越大，确定性谐和响应在联合响应中占比越大，随机响应占比越小。这与直观结论是一致的。例如，在的非共振区间，随机响应占优，在其共振区间，确定性响应占优;在的所有频率区间，谐和响应占优。第四，采用扫频计算的方式，可得到频率?响应均方值曲线的跳跃区间稳态解。虽然时域扫频和本文所建议方法在非跳跃区间吻合较好，但在跳跃区间吻合欠佳，值得进一步深入考察。这可能与3个因素相关：首先，将总体响应分解为确定性谐和与随机响应之和的假定（式（3））在该区间的合理性;其次，求解随机响应二阶矩的统计线性化方法在该区间的精确性;最后，时域积分方法的累积误差。现就以上3个因素作以下评述。

将总体响应分解为谐波响应和随机响应之和是一种启发式的假定，也是发展本文所建议方法的基础。然而，目前尚没有相关报道对于这种假定的合理性与适用范围作出严格数学证明或定性的物理解释。显然，对于线性系统，这种分解方式是完全精确的。有理由推定，该假定同样也适用于弱非线性系统和小强度随机激励的情况。例如，Rong等[24]认为随机响应对于确定性响应是微小摄动时可采用这种假定;Cai和Lin[34]认为该方法不适用于强非线性和乘性激励的情况。其他采用这种假定的文献可见[22，26，28，35?37]。其次，本文应用了统计线性化方法，这种方法须假定响应是高斯或近似高斯分布的随机过程。然而，在发生随机跳跃的频率区间内，系统响应的概率密度函数具有双峰值，明显偏离高斯分布。最后，由图4（d）可知，本文所建议方法计算得到的联合响应均方值与只有谐和激励时确定性响应幅值的解析解吻合较好，证明了本文所建议方法在谐和激励占优的情况下的合理性。然而，时域模拟方法与本文建议方法得到的分岔点并不吻合（图4（c）?（d）中箭头）。因此，建议对采用的时域数值模拟差分算法[38]进一步进行误差分析。

联合激励下非线性整数阶与分数阶Duffing振子的跳跃与分岔的物理机制可分别参阅文献[30，39]，它们与确定性跳跃有区别和联系。不同于只发生在确定性频率点上的确实性单向跳跃，随机双向跳跃发生在可跳跃频率区间上的任意频率点。发生随机跳跃的Duffing振子的响应概率密度具有双峰值。此时，发生概率最高的两个幅值与幅频响应曲线三值区间的两个稳定解对应。然而，由本文建议方法得到的跳跃区间内两个稳定响应均方值与多峰响应概率密度之间的关系，仍有待深入考察。

进一步考察其他因素对响应均方值的影响。以下各图为使用本文建议方法和时域数值模拟的非扫频计算方式获得数值结果。图5所示为简谐激励幅值不同时，系统阻尼比与响应均方值之间的关系（其他参数均与图2所示的系统参数相同）。可见，多数情况下，所提方法得到的结果与数值模拟结果吻合较好;谐和激励幅值较大（）且较小时，两种方法所得结果有较大差距。图6所示为简谐激励具有不同幅值时，与响应均方值之间的关系（其他参数均与图2采用的系统参数相同）。可见，每一种确定性激励幅值都具有两种方法无法吻合的非线性强度区间。在此区间之外，二者吻合较好。激励幅值越大，无法吻合的幅度也越大，且与之对应的区间越小。图7显示了分数阶数不同时两种方法所得结果的对比。可见，多数情况下二者吻合较好，但是在谐和激励和均较大的情况下二者吻合欠佳。由图4（d）可知，激励幅值时，激励频率对应着系统响应均方值产生跳跃的多值区间。此时，采用非扫频计算方式得到的响应均方值可能存在两种方法无法吻合的情况。上文已初步分析了导致这种情况的三种原因。图5?7中所示的两种方法吻合欠佳的情况，也均对应着产生随机跳跃的参数区间。例如，图6中，时会产生随机跳跃。

3.2 簡谐和色噪声联合激励下

考虑功率谱密度为（38）的色噪声，则归一化后的激励的功率谱密度可以由式（35）?（36）和（38）得出。作为演示算例，系统参数取，，;确定性激励参数取;随机激励取和，无量纲截止频率取。

两种方法的均方值比较如图8?11所示（图中虚线代表本文建议方法的计算结果，标记实线为时域数值模拟结果），均采用非扫频计算方式。图8所示为简谐激励具有不同幅值时，与响应均方值之间的关系。可见，多数情况下，建议方法和时域数值模拟吻合较好。与白噪声情况类似，在产生随机跳跃的简谐激励频率附近两种方法吻合欠佳。

由图8可知，激励频率处于外主共振区间。使用该激励频率用于考察其他参数与响应均方值之间的关系，如图9?11所示。其中，图9所示为简谐激励幅值不同时，系统阻尼比与响应均方值之间的关系;图10所示为简谐激励具有不同幅值时，与响应均方值之间的关系;图11所示为分数阶数不同时所建议方法和时域数值模拟的结果对比。可见，简谐激励频率处于外主共振区间时，所提方法得到的结果与时域数值模拟结果均吻合较好。

4 结论与展望

本文结合谐波平衡和统计线性化方法考察了随机和确定性谐和激励联合作用下分数阶非线性振子的响应。该方法的重点在于将联合激励下的稳态响应分解为确定性简谐和零均值随机分量之和，并将原运动方程化为关于确定性与随机响应分量的子运动方程。分别利用谐波平衡法和分数阶非线性系统的统计线性化方法求解了耦合的子运动方程，得到了系统响应的二阶矩。数值模拟表明，多数参数设置情况下，所建议方法有较好的精度。甚至对于强非系统，所建议方法都可获得理想结果。然而，该方法在产生随机跳跃的参数空间内的适用性和计算精度有待进一步考察。显然，本文所建议方法继承了统计线性化方法的广泛适用性。因此，发展该方法的意义是为联合激励下的具有分数阶导数模型的工程结构提供高效且具有理想精度的随机响应计算方法。

本文所建议方法可顺利拓展到联合激励下具有其他非线性形式和（或）非平稳随机激励的多自由度分数阶导数系统。

参考文献：

[1] Scott Blair G W， Reiner M. The rheological law underlying the nutting equation[J]. Applied Scientific Research， 1951，2（1）：225-234.

[2] Oldham K B， Spanier J. The replacement of Fick's law by a formulation involving semidifferentiation[J]. Journal of Electroanalytical Chemistry， 1970，26（2）：331-341.

[3] Nutting P G. A new general law of deformation[J]. Journal of the Franklin Institute， 1921，191（5）：679-685.

[4] Caputo M. Vibrations of an infinite plate with a frequency independent Q[J]. Journal of the Acoustical Society of America， 1976，60（3）：634-639.

[5] Slonimsky G L. Laws of mechanical relaxation processes in polymers[J]. Journal of Polymer Science Polymer Symposia， 1967，16（3）：1667-1672.

[6] Smit W， Vries H D. Rheological models containing fractional derivatives[J]. Rheologica Acta， 1970，9（4）：525-534.

[7] Bagley R L， Torvik P J. A theoretical basis for the application of fractional calculus to viscoelasticity[J]. Journal of Rheology， 1983，27（3）.

[8] Bagley R L， Torvik P J. On the fractional calculus model of viscoelastic behavior[J]. Journal of Rheology， 1986，30（1）：133-155.

[9] Bagley R L， Torvik P J. Fractional calculus in the transient analysis of viscoelastically damped structures[J]. AIAA Journal， 1985，23（6）：918-925.

[10] Gaul L， Klein P， Kempfle S. Impulse response function of an oscillator with fractional derivative in damping description[J]. Mechanics Research Communications， 1989，16（5）：297-305.

[11] Suarez L E， Shokooh A. An eigenvector expansion method for the solution of motion containing fractional derivatives[J]. Journal of Applied Mechanics， 1997，64（3）：629-635.

[12] Wahi P， Chatterjee A. Averaging oscillations with small fractional damping and delayed terms[J]. Nonlinear Dynamics， 2004，38（1）：3-22.

[13] Mainardi F. Fractional Calculus： Some Basic Problems in Continuum and Statistical Mechanics[M]. New York： Springer-Veriag， 1997.

[14] Spanos P D， Zeldin B A. Random vibration of systems with frequency-dependent parameters or fractional derivatives[J]. ASCE Journal of Engineering Mechanics， 1997，123（3）：290-292.

[15] Agrawal O P. An analyticai scheme for stochastic dynamic systems containing fractional derivatives[J]. Proceedings of the American Society of Mechanical Engineers Design Engineering Technical Conferences， Las Vegas， 1999： 243-249.

[16] Ye K， Li Li， Tang Jiaxiang. Stochastic seismic response of structures with added viscoelastic dampers modeled by fractional derivative[J]. Earthquake Engineering and Engineering Vibration， 2003，2（1）：133-139.

[17] 黄志龙， 金肖玲. 具有分數阶导数阻尼的单自由度强非线性随机系统的响应与稳定性[C]. 第十一届全国非线性振动学术会议暨第八届全国非线性动力学和运动稳定性学术会议论文集， 2007.

HUANG Zhilong， JIN Xiaoling. Response and stability of SDOF strongly nonlinear stochastic system with light fractional derivative damping [C]. Proceedings of the Eleventh National Conference on Nonlinear Vibration and the Eighth National Conference on Nonlinear Dynamics and Stability of Motion， 2007.

[18] 孙春艳. 几类典型随机系统的动力学响应分析[D]. 西安：西北工业大学， 2015.

SUN Chunyan. Dynamic response analysis for several classes of typical stochastic systems[D]. Xi'an： Northwestern Polytechnical University， 2015.

[19] Srinivasan V， Soni A H. Seismic analysis of rotating mechanical systems—A review[J]. The Shock and Vibration Digest， 1982，14（6）：13-19.

[20] Spencer B F， Tang J， Hilal C G. Reliability of non-linear oscillators subjected to combined periodic and random loading[J]. Journal of Sound and Vibration， 1990，140（1）：163-169.

[21] Namachchivaya N S. Almost sure stability of dynamical systems under combined harmonic and stochastic excitations[J]. Journal of Sound and Vibration， 1991，151（1）：77-90.

[22] Iyengar R N. A nonlinear system under combined periodic and random excitation[J]. Journal of Statistical Physics， 1986，44（5-6）：907-920.

[23] Nayfeh A H， Serhan S J. Response statistics of non-linear systems to combined deterministic and random excitations[J]. International Journal of Nonlinear， 1990，25（5）：493-509.

[24] Rong Haiwu， Xu Wei， Meng Guang， et al. Response of a duffing oscillator to combined deterministic harmonic and random excitation[J]. Journal of Sound & Vibration， 2001，242（2）：362-368.

[25] Rong Haiwu， Meng Guang， Wang Xiangdong， et al. Response statistic of strongly non-linear oscillator to combined deterministic and random excitation[J]. International Journal of Non-Linear Mechanics， 2004，39（6）：871-878.

[26] Anh N D， Hieu N N. The Duffing oscillator under combined periodic and random excitations[J]. Probabilistic Engineering Mechanics， 2012，30（5）：27-36.

[27] Spanos P D， Zhang Y， Kong F. Formulation of statistical linearization for MDOF systems subject to combined periodic and stochastic excitations[J]. Journal of Applied Mechanics， 2019，86（10）：101003.

[28] Zhu H T， Guo S S. Periodic response of a Duffing oscillator under combined harmonic and random excitations[J]. Journal of Vibration and Acoustics， 2015， 137（4）： 041015.

[29] Chen Lincong， Zhu Weiqiu. Stochastic stability of Duffing oscillator with fractional derivative damping under combined harmonic and white noise parametric excitations[J]. Acta Mechanica， 2009，207（1-2）：109-120.

[30] Chen Lincong， Zhu Weiqiu. Stochastic jump and bifurcation of Duffing oscillator with fractional derivative damping under combined harmonic and white noise excitations[J]. International Journal of Non-Linear Mechanics， 2011，46（10）：1324-1329.

[31] Chen Lincong， Zhu Weiqiu. Stochastic averaging of strongly nonlinear oscillators with small fractional derivative damping under combined harmonic and white noise excitations[J]. Nonlinear Dynamics， 2009，56（3）：231-241.

[32] 陳林聪， 李海锋， 梅真， 等. 宽带噪声激励下含分数阶导数的van der Pol-Duffing振子的可靠性[J]. 西南交通大学学报， 2014，49（01）：45-51.

CHEN Lincong， LI Haifeng， MEI Zhen， et al. Relability of van der Pol-Duffing oscillator with fractional derivative under wide-band noise excitations[J]. Journal of Southwest Jiaotong University， 2014，49 （01）： 45-51.

[33] 陈林聪， 朱位秋. 谐和与宽带噪声联合激励下具有分数导数型阻尼的Duffing振子的平稳响应[C]. 中国力学学会学术大会， 2009.

CHEN Lincong， ZHU Weiqiu. Stationary response of Duffing oscillator with fractional derivative damping under combined harmonic and wide band noise excitations[C]. Conference of the Chinese Society of Mechanics， 2009.

[34] Cai G O， Lin Y K. Nonlinearly damped systems under simultaneous broad-band and harmonic excitations[J]. Nonlinear Dynamics， 1994，6（2）：163-177.

[35] Manohar S C， Iyengar N R. Entrainment in van der Pol's oscillator in the presence of noise[J]. International Journal of Non?Linear Mechanics， 1991，26（5）：679-686.

[36] Bulsara A R， Lindenberg K， Shuler K E. Spectral analysis of a nonlinear oscillator driven by random and periodic forces. I. Linearized theory[J]. Journal of Statistical Physics， 1982，27（4）：787-808.

[37] Budgor A B. Studies in nonlinear stochastic processes. Ⅲ. Approximate solutions of nonlinear stochastic differential equations excited by Gaussian noise and harmonic disturbances[J]. Journal of Statistical Physics， 1977，17（1）：21-44.

[38] Koh C G， Kelly J M. Application of fractional derivatives to seismic analysis of base-isolated models[J]. Earthquake Engineering and Structural Dynamics， 1990，19（2）：229-241.

[39] Huang Z L， Zhu W Q， Suzuki Y. Stochastic averaging of strongly non-linear oscillators under combined harmonic and white-noise excitations[J]. Journal of Sound & Vibration， 2000，238（2）：233-256.