基于多谐波平衡法的滞回分数阶系统稳态动力响应

孔凡 侯召旭 徐军 李书进

摘要： 结构隔震和耗能减振装置有时会同时体现出动力特性的频率依赖和滞回关系，分数阶导数模型能以较少的参数模拟力?位移关系的频率依赖性。采用谐波平衡法研究简谐激励下同时具有滞回特性和分数阶阻尼单元的系统稳态响应。利用激励和响应过程的Fourier级数展开并取谐波平衡后，求得滞回响应和位移响应Fourier级数之间的关系;将滞回微分方程写成余量的形式并结合Galerkin方法和Levenberg?Marquardt算法求得响应的Fourier级数;分别对硬化和软化Bouc?Wen滞回系统的数值模拟显示所建议方法的精度。研究结果表明，分数阶数和稳态位移幅值之间的关系依赖于系统和滞回参数以及激励频率。

关键词： 多谐波平衡; Bouc?Wen滞回; Levenberg?Marquardt算法; 分数阶导数

中圖分类号： TU311.3; O322    文献标志码： A    文章编号： 1004-4523（2021）03-0552-07

DOI：10.16385/j.cnki.issn.1004-4523.2021.03.012

引  言

极端灾害性荷载作用下的工程结构构件和节点常表现出明显的力?位移滞回关系，该关系表明恢复力不仅依赖于运动当前状态，而且依赖于其复杂的历史状态[1]。实际上，人们通常利用合理的结构设计，期望在大震作用下结构的特殊部位（如框架结构的梁端）首先出现塑性铰从而耗散地震输入能量，从而达到结构大震不倒的目的[1]。此外，一些用于结构隔震和耗能减振的控制装置，也通常表现出明显滞回特性，如约束屈曲支撑、铅芯橡胶隔震垫和黏滞阻尼器等装置[2]。为了准确合理描述结构构件和控制装置中出现的滞回特性，人们发展了多种滞回模型，包括Ramberg?Osgood模型、多线性模型以及Bouc?Wen模型等[3]。Bouc?Wen模型具有接近实际情况的光滑曲线、仅用一个微分方程就可表现刚度退化和卸载拐点等复杂滞回特性，以及与运动微分方程联合后形成易于求解的扩阶运动方程等优点。

另外，许多应用于工程结构或控制装置中的材料都存在本构关系的频率相关效应，即本构关系依赖于加载频率[4]。传统上，人们使用简单的Kelvin或Maxwell以及它们复杂的串并联形式[5]（标准线性固体模型，Standard Linear Solid Models， SLS）部分描述这种频率依赖特性。注意到，SLS模型需要大量参数才能精确地同时拟合试验得到的储能和耗能模量数据。在此过程中，SLS模型的某些拟合参数易出现负值，从而失去物理意义[6]。反之，分数阶导数[7]模型以少量参数就能精确地描述本构关系的频率依赖性[6]，引起了越来越多的关注。该模型的控制装置应用包括：黏弹性阻尼[8]、隔震装置[9]、液体黏滞阻尼器[10]、调谐液柱阻尼[11]，调谐质量阻尼器[12]等。

可见，工程结构、控制装置和受控结构可能会同时出现力?位移关系的滞回和频率依赖特性。以一种新型的广泛应用于半主动控制装置的磁敏弹性体为例，本质上作为一种橡胶黏弹性材料，具有明显的本构频率依赖性[13]。为此，Behrooz等[14]利用三参数SLS模型和Bouc?Wen滞回模型的并联形式、Eem等[15]采用Maxwell模型和Ramberg?Osgood滞回模型的并联形式模拟磁敏弹性体的力?位移关系。Behrooz和Eem的研究均关注装置的滞回力学性能，对黏弹性材料的频率依赖性却只以简单的SLS模型模拟。相反，在分数阶导数非线性动力学研究方面，人们往往较多关注分数阶导数线性[16?17]或无记忆非线性[18?19]动力系统，对具有滞回特性的分数阶系统研究较少。文献调查显示只有少量研究[20?21]关注滞回分数动力系统的随机和稳态动力响应。其中，文献[20]通过等效线性化方法将滞回分数阶系统化为等效线性参数依赖于幅值的线性系统，得到了滞回分数阶系统频率?响应幅值曲线。该曲线只具有主共振峰值，而无法反应超谐共振峰值。

考察动力系统在谐波激励下的稳态响应是理解系统动力特性重要的一环，因此，本文与文献[20]一样，关注谐波激励下滞回分数阶系统稳态响应。但是，本文与文献[20]采用的方法不同，本文采用谐波平衡法研究具有不同分数阶导数的滞回系统的稳态响应，以期得到频率?响应幅值的完整信息。首先，发展分数阶滞回系统的谐波平衡法，将滞回运动微分转化为一组代数方程;随后，利用Levenberg?Marquardt算法对与该代数方程对应的最小二乘问题求解;最后，考察硬化和软化的Bouc?Wen滞回系统，得到它们在不同分数阶数情况下的频率?响应幅值曲线。与时域逐步积分法的对比显示了所建议方法的精度。

1 分数阶滞回系统的谐波平衡法

考虑具有分数阶阻尼项的单自由度滞回动力系统

式中  分别为质量和线弹性刚度系数;分别为系统位移和加速度响应;为滞回位移;为屈服前后刚度比;为分数阶阻尼系数;为阶Riemann?Liouville意义上的分数阶导数，即

下文为方便计，将Caputo分数阶导数简写为并简称为分数阶导数。

另一方面，将位移响应和激励在区间上进行周期延拓后，可将它们写为Fourier级数的形式：

式（18）表明了滞回位移的Fourier级数向量可表示为位移的Fourier级数向量的形式。然而，这两个级数向量都未知，需引进它们之间的第二个关系式。为此，引入和之间的第二个关系式，即形如

利用Galerkin法将余量式（20）在Fourier基上投影，可得到第二个与之间的关系式，随后将此关系式代入式（18）后可将向量消除，即可得到以为未知量的方程

式?中  为余量的Fouier级数。求解式（21）可得位移响应的Fourier级数，对响应的Fourier级数作Fourier逆变换后可得位移响应时程。

2 代数方程数值求解的Levenberg-Marquardt算法

式（21）的求解可通过许多数值算法实现，本文采用一种高效的Levenberg?Marquardt（LM）算法。LM算法最早由Marquardt[23]于1963年提出，是一种应用广泛的非线性最小二乘算法，可视为最速下降法和线性化方法（泰勒级数）的一种组合;可参阅相关数值计算参考书目。它将代数方程转化为求解最小二乘的形式，即

式中  均为控制滞回环形状的模型参数。

3 数值算例

利用与逐步积分法得到的相关结果对比验证所建议方法的适用性和精度。作为演示，本算例采用具有微分形式的Bouc?Wen滞回模型，但本文发展的方法同时也适用其他滞回关系。注意到，具有Bouc?Wen滞回的整数阶动力系统响应的数值解可通过增量形式的Newmark???Newton???Raphson （NNR）迭代法实现，其本质是将动力方程转化为增量静力方程，难点在于切线性变刚度的求解和滞回模型拐点的处理。在此基础上，为求解Bouc?Wen滞回的分数阶动力系统响应，作者发展了一种基于Riemann?Liouville离散分数阶导数的NNR方法：其中，与Newmark?1/6法处理加速度项对应，在对离散式（5）右边积分时，须假定离散区间内的速度呈线性变化。限于篇幅和行文重点，该增量形式的数值方法兹不赘述，拟另文发表。

首先，考察所建议方法在特定参数设置下计算得到的动力响应时程和滞回环。选取动力系统参数，;滞回参数，;激励参数，，。系统质量和刚度参数选取时体现了振子的“归一化”;阻尼参数选取时体现了工程系统的小阻尼特性;滞回参数选取时体现了明显的滞回效应。图1（a）?（c）所示为在上述参数设置情况下所建议方法和逐步积分法得到的，和滞回环对比。由图1（a）和（b）可知：除了位移时程的开始阶段，多谐波平衡法得到的平稳响应能很好地吻合逐步积分法得到的结果;时程开始的误差由多谐波平衡法忽略初始条件导致。

随后，考察分数阶数对软化Bouc?Wen系统响应幅值的影响。除分数阶数外，其他模型及激励参数与上例相同。图2所示为多谐波平衡法和逐步积分法计算得到的响应和幅值与分数阶数之间的关系曲线。由图2可见：所建议方法与时域逐步积分法吻合较好;的幅值随着的增大而减小，的幅值随着的增大稍有减小但几乎保持不变。图3给出了时的滞回环。由图3可知，滞回环的饱满程度随着的增大而减小。

接着，考察分数阶数对硬化Bouc?Wen系统响应幅值的影响。此时，取滞回参数，，其他系统、滞回和激励参数与上例同。图4所示为多谐波平衡法和逐步积分法计算得到的响应（和）幅值与分数阶数之间的关系曲线。由图可知，在该情况下，和幅值均随着的增大而增大，前者增大幅值更加明显。图5所示为时系统的滞回环。可见，所有滞回环均呈现明显的硬化特性，分数阶数对硬化滞回环的影响较软化滞回环时小。

分数阶数对硬化和软化系统响应幅值的影响不同。为深入考察这一点，本文利用谐波平衡法和逐步积分法得到了分数阶导数的情况下，上述软化和硬化系统的激励频率?响应幅值曲线，如图6和7所示。由图6?7可见，在所有情况下，谐波平衡法均能很好地吻合逐步积分法给出的结果。此外，应注意到，文[20]利用线性化方法研究滞回分数阶系统的频率响应曲线，只能得到主共振峰值，无法得到超谐共振峰值;而本文采用的谐波平衡法能同时得到主共振和超谐共振峰值。文[24]采用增量谐波平衡法和平均法研究了分数阶硬Duffing振子的频率响应曲线，也得到类似的结论。由图6?7可见：软特性Bouc?Wen分数阶系统的激励频率?响应幅值曲线左侧较陡;反之，硬特性Bouc?Wen分数阶系统的激励频率?响应幅值曲线右侧较陡。本例采用的激励频率（）处于图6所示软特性峰值右侧，分数阶数增大时响应幅值呈明显递减趋势;同时，激励频率处于图7所示硬特性峰值左侧，分数阶数增大时响应幅值呈小幅增加趋势。由以上分析不难得出：对于特定滞回系统，分数阶数对响应幅值的影响与激励频率大小有密切关系。例如，对于如图6所示的软特性分数阶系统，与时的情况不同，当激励频率时，响应幅值随分数阶数的增加而增加。另外，分数阶数对响应幅值的影响还与滞回参数有关，如本例所示软或硬特性系统参数。值得注意的是，当系统阻尼和滞回耗能较小时，非线性系统的激励频率?响应幅值曲线易出现分岔现象，进而导致更加复杂的分数阶数和响应幅值之间的关系。

现对所提出方法在参数设置和计算效率方面作如下评述。首先，在计算式（24）所示的Jacobian矩阵时，每个迭代步均涉及时域和频域之间的相互转化，必须进行Fourier变换[25]。因此，使用快速Fourier变换（Fast Fourier Transform， FFT）及其逆变换能加快速度。其次，系统非线性的存在，可能会产生高频响应;因此，在设置响应上限截止频率时，须覆盖响应可能的频率区间。本文采用的上限截止频率大于10倍激励频率。最后，频率泄露可能影响收敛速度并降低计算精度。为此，在设置采样频率时，建议在激励频率上取离散频率点，即激励频率为离散频率步长的整数倍。计算每个激励频率点的响应幅值时，均须综合考虑上述要求，合理指定计算参数。例如，采用与图1相同的系统参数设置，选取上限截止频率，频域内采样间隔，共100个采样点。此时，第10个频率离散点等于激励频率。试算表明，仅5次迭代就可达到预定收敛目标，耗时0.07 s。

本文结合Riemman?Liouville定义的离散格式和Newmark?Newton?Raphson迭代方法在时域中求解Bouc?Wen分数阶导数系统的动力响应。Riemman?Liouville离散格式中，当前积分步的运动状态须表示为之前所有离散时间点运动状态的加权和。其中，每个历史时间点的权重须及时更新，大大降低了计算效率。此外，与整数阶非线性动力系统的Newton?Raphson迭代方法一样，更新切线刚度与寻找滞回位移拐点也增加了时域方法的计算量。计算经验表明：时域方法的计算时长比时间离散点数增长更快;同时，须采用非常小的时间步长以保证计算精度。例如，采用与上述频域方法相同的系統参数设置，使时域方法达到与频域方法相同的计算精度时，计算耗时0.71 s。

4 結论与展望

本文发展了谐和激励下滞回分数阶动力系统平稳解的谐波平衡法。首先，利用谐波平衡法将滞回分数阶微分方程转化为了一组以响应Fourier级数为未知量的非线性代数方程;利用Levenberg?Marquardt算法求解了与非线性代数方程对应的最小二乘问题，得到了响应Fourier级数。研究表明，谐波平衡法得到的响应时程、滞回环、频率?幅值曲线均与逐步积分法得到的相关结果吻合较好。该方法能同时得到滞回分数阶系统频率?幅值曲线的主共振和超谐共振峰值。

对具有不同滞回特性的系统，考察了分数阶数对响应幅值大小的影响。文中给出了在特定软特性和硬特性参数的情况下，激励频率等于线性系统自振频率时，分数阶数对响应幅值影响的规律。研究发现，分数阶数对响应幅值影响较复杂，与系统特性、滞回特性、激励频率均有一定关系;拟另文发表。

进一步的工作可在其他滞回模型、多自由度系统以及滞回分数阶系统的混沌和分岔等方面展开。

参考文献：

[1] 李 杰，李国强.地震工程学导论[M].北京：地震出版社，1992.

[2] 欧进萍.结构振动控制[M].北京：科学出版社，2003.

[3] 潘鹏，张耀庭.建筑抗震设计理论与方法[M].北京：科学出版社，2017.

[4] Xu Z-D， Wang D-X， Shi C-F. Model， tests and application design for viscoelastic dampers[J]. Journal of Vibration and Control， 2011， 17： 1359-1370.

[5] Xu Z-D. Earthquake mitigation study on viscoelastic dampers for reinforced concrete structures[J]. Journal of Vibration and Control， 2007， 13： 29-43.

[6] Di Paola M， Pirrotta A， Valenza A. Visco-elastic behavior through fractional calculus： An easier method for best fitting experimental results[J]. Mechanics of Materials， 2011， 43（12）： 799-806.

[7] Oldham K， Spanier J. The Fractional Calculus Theory and Applications of Differentiation and Integration to Arbitrary Order[M]. Elsevier， 1974.

[8] Xu Z-D， Xu C， Hu J. Equivalent fractional Kelvin model and experimental study on viscoelastic damper[J]. Journal of Vibration and Control， 2015， 21（13）： 2536-2552.

[9] Hwang J， Wang J. Seismic response prediction of HDR bearings using fractional derivative Maxwell model[J]. Engineering Structures， 1998， 20： 849-856.

[10] Makris N， Dargush G， Constantinou M. Dynamic analysis of generalized viscoelastic fluids[J]. Journal of Engineering Mechanics， 1993， 119： 1663-1679.

[11] Di Matteo A， Iacono F L， Navarra G， et al. Innovative modeling of tuned liquid column damper motion?[J]. Communications in Nonlinear Science and Numerical Simulation， 2015， 23： 229-244.

[12] Rüdinger F. Tuned mass damper with fractional derivative damping?[J]. Engineering Structures， 2006， 28： 1774-1779.

[13] Li Y， Li J. On rate-dependent mechanical model for adaptive magnetorheological elastomer base isolator [J]. Smart Materials and Structures， 2017， 26： 045001.

[14] Behrooz M， Wang X， Gordaninejad F. Modeling of a new semi-active/passive magnetorheological elastomer isolator[J]. Smart Materials and Structures， 2014， 23： 045013.

[15] Eem S H， Jung H J， Koo J H. Modeling of magneto-rheological elastomers for harmonic shear deformation[J]. IEEE Transactions on Magnetics， 2012， 48： 3080-3083.

[16] Chang T S， Singh M P. Seismic analysis of structures with a fractional derivative model of viscoelastic dampers[J]. Earthquake Engineering and Engineering Vibration， 2002， 1： 251-260.

[17] Singh M P， Chang T S， Nandan H. Algorithms for seismic analysis of MDOF systems with fractional derivatives?[J]. Engineering Structures， 2011， 33： 2371-2381.

[18] Failla G， Pirrotta A. On the stochastic response of a fractionally-damped Duffing oscillator[J]. Communications in Nonlinear Science and Numerical Simulation， 2012， 17： 5131-5142.

[19] Chen L， Wang W， Li Z， et al. Stationary response of Duffing oscillator with hardening stiffness and fractional derivative[J]. International Journal of Non-Linear Mechanics， 2013， 48： 44-50.

[20] Spanos P D， Di Matteo A， Pirrotta A. Steady-state dynamic response of various hysteretic systems endowed with fractional derivative elements[J]. Nonlinear Dynamics， 2019，98（4）：3113-3124.

[21] Di Matteo A， Spanos P D， Pirrotta A. Approximate survival probability determination of hysteretic systems with fractional derivative elements[J]. Probabilistic Engineering Mechanics， 2018， 54： 138-146.

[22] Tren?evski K， Tomovski ?. On fractional derivatives of some functions of exponential type[J]. Publikacije Elektrotehni?kog fakulteta-serija： Matematika， 2002，（13）： 77-84.

[23] Marquardt D W. An algorithm for least-squares estimation of nonlinear parameters[J]. Journal of the Society for Industrial and Applied Mathematics， 1963， 11（2）： 431-441.

[24] Shen Y J， Wen S F， Li X H， et al. Dynamical analysis of fractional-order nonlinear oscillator by incremental harmonic balance method[J]. Nonlinear Dynamics， 2016， 85： 1457-1467.

[25] Wong C W， Ni Y Q， Lau S L. Steady-state oscillation of hysteretic differential model. I： Response analysis[J]. Journal of Engineering Mechanics， 1994， 120（11）： 2271-2298.

作者簡介： 孔  凡（1984-），男，博士，副教授，硕士生导师。E-mail： kongfan@whut.edu.cn