图1 多自由度分段线性系统
Fig. 1 A piecewise linear system with multiple degrees of freedom
何冬东， 高 强， 钟万勰
(工业装备结构分析国家重点实验室(大连理工大学)；大连理工大学 运载工程与力学学部 工程力学系， 辽宁 大连 116024)
摘要： 基于参变量变分原理，提出了一种求解具有大量间隙弹簧的周期性分段线性系统动态响应的高效率数值方法．通过参变量变分原理来描述间隙弹簧，将复杂的非线性动力问题转化为线性互补问题求解，避免了求解过程中的迭代和刚度阵更新，该算法能准确判断间隙弹簧的压缩和松弛状态．基于结构的周期性和能量传播速度的有限性，提出了一种求解系统动态响应的高效率精细积分方法．该算法指出周期结构的矩阵指数中存在大量的相同元素和零元素，从而不需要重复计算和存储这部分元素，节省了计算量并降低了计算机存储要求．分析了一个五自由度分段线性系统在简谐荷载作用下的动力学行为，包括稳定的周期运动、准周期运动和混沌运动．通过与Runge-Kutta方法的比较，该文方法的正确性和高效率得到了验证．
关 键 词： 分段线性系统； 动态响应； 参变量变分原理； 周期性； 精细积分方法
在自然科学和工程技术中广泛存在着各种非线性问题，为了简化非线性问题，常常采用分段线性理论来近似系统的非线性行为．在众多领域中，均能发现很多非线性动力学问题的理论模型中含有分段线性，例如电力工程[1]、生物工程[2]、航空航天工程[3]等．在工程结构系统中较为常见的是分段线性刚度或阻尼，如由压杆和绳索组成的tensegrity结构[4]、含有间隙的振动机械[5]和振动压实系统[6]等．由于分段线性系统的动力学特性有利于工程设计和应用[7]，因此其被广泛应用于工程实践中．各种理论和实践分析表明，含有分段线性特点的动力学系统，受周期荷载作用将会展现非常复杂的动力学行为[8-9]．
尽管一些分段线性动力学系统的自由振动可以被解析求解，但是要解析求解一个分段线性系统的强迫振动是非常困难的，即使是很简单的分段线性系统．一些经典的解析方法和近似解析法，如摄动法、平均法、谐波平衡法和增量谐波平衡法等[10-12]，常常被用于求解分段线性动力学系统的稳态响应．Narimani等[13]和Deshpande等[14]采用平均法求解了两段线性系统在稳态情况下的频率响应．Moussi等[15]结合谐波平衡法和渐进数值方法提出了一种求解双自由度分段线性系统正则问题的计算方法．Xu等[16]采用增量谐波平衡法给出了具有分段线性黏性阻尼特性动力系统的周期解．但是这些方法应用范围仅限于弱非线性系统，对于具有强非线性特点的分段线性系统则很难给出准确解．此外，对于具有大量自由度和非线性的分段线性系统，这些方法也很难给出准确解．
对于多自由度的强非线性系统，通常只能采用数值方法求解．若给定初值，很多数值方法如Runge-Kutta方法、直接积分法结合Newton-Raphson方法等都能给出分段线性系统的高精度瞬态响应和稳态响应．但这些算法通常需要反复迭代计算，使得计算效率非常低，特别是收敛精度设置得很高时．Zhong(钟万勰)等[17]提出的参变量变分原理对于分段线性分析具有明显优势．因此，采用参变量变分原理来描述分段线性问题，将其转化为线性互补问题[18-19]，可以避免求解过程中的反复迭代．针对结构的周期性特点，基于精细积分方法[20]，建立一种求解分段线性系统动态响应的高效率数值积分方法．该算法指出周期结构的矩阵指数中存在大量的相同元素与零元素，从而不需要重复计算和存储这部分元素，节省了计算量并降低了计算机存储要求．数值结果验证了本文方法的正确性．分析了一个五自由度不对称分段线性系统稳态运动的动力学行为，如稳定的周期运动、准周期运动和混沌运动．此外，通过求解一个含有1 000个自由度和1 000个间隙弹簧的分段线性系统在周期荷载作用下的动态响应，比较了本文方法和Runge-Kutta方法的计算时间，证明了本文方法对于求解大规模周期性分段线性系统动态响应的高效率．
图1 多自由度分段线性系统
Fig. 1 A piecewise linear system with multiple degrees of freedom
考虑一个如图1所示具有大量间隙弹簧的多自由度分段线性系统．该系统含有n个质量块，n个正常连接的弹簧和n个阻尼器，在每个质量块的右端都有一根刚度为k2的弹簧，质量块与该弹簧之间存在一个初始间隙Δi(i=1,2,…,n)．每个质量块均受一简谐荷载fi(t)=Aisin(ωit)作用．间隙弹簧只有在质量块的位移超过初始间隙Δi时才会受力，而当质量块的位移小于初始间隙时，该弹簧处于松弛状态，而刚度为k1的弹簧则始终提供弹性恢复力．于是该分段线性系统的动力平衡方程可以写作
(1)
其中，
和
分别为系统的位移、速度和加速度向量，M为质量阵，C为阻尼阵，K1是所有刚度为k1的弹簧组成的刚度阵， f为外荷载向量，P(x)为刚度为k2的间隙弹簧提供的弹性恢复力，即
(2)
式中，In是一个n×n维的单位矩阵．
对于第i根间隙弹簧，其本构关系如图2所示，当质量块的位移未超过初始间隙时，该弹簧处于松弛状态，其提供的弹性恢复力为0；当质量块的位移超过初始间隙时，该弹簧处于压缩状态，提供的弹性恢复力为一个正值．其物理方程可以表示为
(3)
从式(3)可以看出，这是一个分段线性刚度问题，其本构关系在xi=Δi处存在转折点．因此，在求解方程(1)的过程中，在每个积分步内都必须准确判断质量块的位移是否超过了初始间隙，然后更新间隙弹簧的刚度阵．这使得求解分段线性系统动态响应的效率非常低下，尤其是该系统还含有大量自由度和大量间隙弹簧．因此，有必要发展一种求解具有大量间隙弹簧的多自由度分段线性动力学系统动态响应的高效率算法．在本文中，将采用参变量变分原理[17]来描述间隙弹簧的分段线性特点，从而避免在积分过程中更新间隙弹簧的刚度阵．
图2 间隙弹簧的本构关系
Fig. 2 The constitutive relation of the gap-activated spring
显然，式(3)中间隙弹簧的弹性恢复力Pi是一个非负的值．引入一个非负的松弛变量νi，其本构关系(3)可以写作如下互补关系式：
Pi-k2(xi-Δi)-νi=0， Pi≥0, νi≥0, Piνi=0．
(4)
将运动平衡方程(1)中间隙弹簧的弹性恢复力作为外力项移到方程右边，并结合互补关系式，可以得到分段线性系统的动力学控制方程为
(5)
P-K2(x-Δ)-ν=0,
(6)
式中,Pi≥0, νi≥0, Piνi=0,Δ是初始间隙向量，K2是间隙弹簧压缩刚度组成的刚度阵，即
K2=k2In．
(7)
间隙弹簧的应力状态由其弹性恢复力Pi来表征，即当Pi>0时，间隙弹簧处于压缩状态；当Pi=0时，间隙弹簧为松弛状态．因此，在实际求解过程中，只需计算Pi的值即可，避免了不断判断间隙弹簧的状态，从而避免了其刚度矩阵的更新．
基于参变量变分原理，分段线性动力学方程(1)被转换成了由线性动力学方程(5)和互补关系式(6)组成的动力学控制方程．对于方程(5)，可以采用直接积分法求解，如Runge-Kutta、Newmark-β或Wilson-θ方法等．由于稳定性和精度的要求，这些方法的积分步长都需要取得很小，这使得问题的计算规模非常大、效率非常低．由Zhong(钟万勰)等[20]提出的精细积分方法允许采用很大的积分步长，并具有精度高和稳定性好的特点[21]．对于小规模结构，精细积分方法可以高效精确地求解其动态响应，但是这种方法在处理规模很大的系统时，矩阵指数的计算量非常大．因此有必要发展一种求解多自由度分段线性系统响应的高效率积分方法．
方程(5)在状态空间中可以写作
(8)
其中v是状态向量，且
(9)
将积分区间等分为步长为η的时间区间，即t0，t1=t0+η，…，tk=t0+kη，…．若记vk=v(tk)，那么方程(8)的解为
vk+1=Tvk+
exp[H(η-ξ)]F(tk+ξ)dξ，
(10)
其中T是系统的矩阵指数，其定义如下:
T=exp(Hη)．
(11)
假设在一个积分步长t∈[tk,tk+1]内，参数变量P(t)采用线性近似，即
(12)
则方程(10)中积分可以解析求得
(13)
其中
(14)
式中，T12和T22分别是矩阵指数T中对应的分块矩阵，即
(15)
当得到方程(13)中的积分格式后，现在的主要问题是求解矩阵指数T．对于自由度数目较少的结构，其矩阵指数可以采用精细积分方法[20-21]快速精确地求解．但是当结构的自由度数目很多时，原始精细积分方法的计算量将非常大．在本文中，考虑的结构可以被看作是一个周期结构，两个相邻质量块之间的部分可以作为一个原胞，其包含一根间隙弹簧．本文将利用结构的周期性和结构动态响应的物理特点发展一种计算周期结构矩阵指数的高效率方法．
若方程(8)中外力F(t)=0，则方程(13)可以写作如下分块形式：
(16)
方程(16)中矩阵指数的分块矩阵中元素的物理意义如下：仅第j个自由度上有单位位移，经过一个时间积分步长η后，第i个自由度的位移和速度分别为T11(j,i)和T21(j,i)．同理，仅第j个自由度上有单位速度，经过一个时间积分步长η后，第i个自由度的位移和速度分别为T12(j,i)和T22(j,i)．假设一个原胞的自由度数目为d，若仅在第d+i个自由度上给定单位位移，则经过一个时间积分步长η后，第d+j个自由度的位移为T11(d+j,d+i)和速度为T21(d+j,d+i)．根据结构的周期性，两个位移T11(j,i)和T11(d+j,d+i)相同，两个速度T21(j,i)和T21(d+j,d+i)也相同，即
T11(j,i)=T11(d+j,d+i), T21(j,i)=T21(d+j,d+i)，
(17)
同理
T12(j,i)=T12(d+j,d+i), T22(j,i)=T22(d+j,d+i)．
(18)
方程(17)和(18)表明，在周期结构对应的矩阵指数中有大量相同的元素，为了提高计算效率，这些相同的元素不需要重复计算和存储．
众所周知，在一段时间内，能量在结构中只能传递有限距离．假设在某一个自由度上施加一外力，在一个较小的时间步长内，该外力显然只能影响附近有限自由度的响应，而不能影响到所有自由度．因此，矩阵指数的结构特点一定是带状分布．总结起来，基于结构的周期性和能量在结构中传播速度的有限性，周期结构的矩阵指数的结构特点是周期性和带状分布，即
(19)
其中标量ai(i=1,2,…,d)对应于矩阵指数的对角线元素，向量bi的第一个元素不为零，向量ci的最后一个元素不为零．显然，式(15)中矩阵指数的其他分块矩阵T12，T21和T22均具有上述结构特点．因此，只需要几个原胞组成的子结构对应的矩阵指数即可得到完整结构对应的矩阵指数的全部信息．相对于完整结构，几个原胞的自由度少很多，可以快速计算出其对应的矩阵指数．因此，采用上述方法可以提高整体结构矩阵指数的计算效率并降低存储要求．
图3 周期结构
Fig. 3 A periodic structure
考虑如图3所示的结构，并假设在一个积分步长内，能量传播距离不超过3个节点．提取出子结构①，并计算其对应的矩阵指数．假设仅节点6上有单位位移，一个积分步长后，其能量不能传播到节点2，因而其边界对结构的响应不产生影响．显然，子结构①与完整结构②对应节点的位移和速度相同，而完整结构多出的左右两节点的响应都为0．因此，只需计算子结构①的矩阵指数，并从矩阵指数中恰当位置提取节点1对应的数据，即可得到式(19)中需要的ai，bi和ci．
以上给出了计算周期结构对应矩阵指数的主要思想．但是，以上讨论基于预知了能量传播的速度，但对于复杂问题，能量传播速度很难事先计算．因此，本文采用如下的处理方法．为了计算周期结构的矩阵指数，首先，给定两个参数，一个是期望采用的时间步长η1，另一个是一个整数q，它表示在一个时间步长内能量不能传递出q个单胞．假设能量传递出q个单胞所需最小时间为η2，则最后采用η=min(η1,η2)作为时间步长．先将η1作为时间步长，并形成2q+1个单胞组成的子结构对应的系统矩阵，再执行精细积分方法．在执行精细积分法的计算过程中，由于将时间步长η1划分为2N份，再计算时间区段η1/2N对应的矩阵指数，此时η1/2N是一个非常小的时间区段，能量不可能传递出q个单胞．然后，执行精细积分法的增量加法定理，每执行一次相当于时间步长比上一次加倍，这个过程中检查能量是否传递出q个单胞．如果执行完N次加法定理后，能量没有传递出q个单胞，则采用η1作为最终的时间步长，此时矩阵指数也已经计算出．如果在执行完第
次加法定理后，检测到能量传递超出q个单胞，则取
作为最终的时间步长，对应的矩阵指数在上一次的加法定理中也已经计算出．
采用这样的方法，仅需用很少的单胞组成质量、阻尼和刚度矩阵．然后计算对应的矩阵指数，提取一些矩阵元素，即可得到周期结构对应矩阵指数的全部信息．从而可以高效地计算周期结构的动态响应．
由间隙弹簧的参数变量λk+1表征的分段线性系统动态响应已经在方程(13)中给出，通过分块，即可得到系统位移的速度为
(20)
将位移代入到互补关系(6)中，得到间隙弹簧的线性互补问题：
ν-Z1Pk+1-Z2=0， Pi,k+1νi=0, νi≥0, Pi,k+1≥0，
(21)
其中
(22)
采用Lemke方法[22]，方程(21)中间隙弹簧的参数变量很容易计算得到．对于第i根弹簧，其应力状态由其参数变量Pi来表征，即Pi=0表明间隙弹簧处于松弛状态，而Pi>0表明间隙弹簧处于压缩状态．将参数变量Pk+1代入方程(13)，则可得到分段线性系统在tk+1时的响应．
考虑图1所示的分段线性系统，取n=5，研究其复杂的动力学行为．系统的初始位移和速度均为0，在第1个、第3个和第5个节点作用外荷载Asin(ωt)，在第2个和第4个节点作用外荷载-Asin(ωt)．为了方便分析分段线性系统各参数对其动力学行为的影响，将动力学控制方程(5)和(6)写作如下无量纲形式：
(23)
(24)
其中
(25)
可以发现，该五自由度分段线性系统的动力学行为由
和θ这4个无量纲参数确定．对这4个无量纲参数取4组不同的值，分别采用本文方法和Runge-Kutta方法来计算该系统的动态响应，计算时间区间为[0, 1 000T]，其中T=2π/θ是外荷载的周期．由于仅有5个自由度，本文方法计算完整结构的矩阵指数，时间步长取为η=T/200．Runge-Kutta方法采用MATLAB的Ode45函数作为求解器，并将容差设置为δ=10-10．
先积分500个周期使系统响应达到稳态，记录系统在外荷载作用500个周期以后(即500T-1 000T)的响应．当
和θ=1时，以及
和θ=1时，两种方法得到的系统第5个节点稳态运动的相轨迹和Poincaré映射如图4所示．可以发现，本文方法与Runge-Kutta方法得到的相轨迹都是一条闭合曲线，并且完全重合，表明本文方法的正确性．图4(a)和(b)中Poincaré映射分别为2个和7个不动点，这表明当和θ=1时系统出现2倍周期运动，当和θ=1时系统出现7倍周期运动．当
和θ=2时，本文方法积分得到的系统第5个节点稳态运动的相轨迹和Poincaré映射如图5所示．可以发现，相轨迹是在有界区域里的一条杂乱曲线，而Poincaré映射为分布在一条闭合曲线上的点集，这表明此时系统出现了准周期运动．当
和θ=1.5时，图6给出了本文方法积分得到的系统第5个节点稳态运动的相轨迹和Poincaré映射，可以发现，相轨迹是在有界区域里的一条杂乱曲线，而Poincaré映射为分布在一个杂乱无章的点集，这表明系统在这组参数下出现了混沌运动．比较图5和图6，仅从相轨迹上很难区分准周期运动和混沌运动，而在Poincaré映射图上很容易区分．以上结果表明本文方法的正确性，并且分段线性系统在周期荷载作用下会产生非常复杂的动力学行为．
图4 本文方法和Runge-Kutta方法给出的相轨迹和Poincaré映射
Fig. 4 Phase portraits and the Poincaré map obtained from the proposed method and the Runge-Kutta method
(a) 相轨迹 (b) Poincaré映射
(a) The phase portrait (b) The Poincaré map
图5 当
和θ=2时，本文方法给出的相轨迹和Poincaré映射
Fig. 5 When α=5, ζ=0.2, 
and θ=2, the phase portrait and the Poincaré map obtained from the proposed method
考虑图1所示的不对称分段线性系统，取n=1 000，即有1 000个自由度和1 000根间隙弹簧，取m=1 kg，k1=1 N/m，k2=1 N/m，c=0.2 (N·s)/m和Δ=0 m．初始位移和速度均为0，从左至右在奇数节点作用外荷载fi(t)=sin(3t) N(i=1,3,…,999)，在偶数节点作用外荷载fi(t)=-2sin(5t) N(i=2,4,…,1 000)．分别采用本文方法和Runge-Kutta方法计算该系统的动态响应，积分区间为[0, 200] s．对于本文方法，选取11个自由度作为基本结构，由算法确定的积分步长为η=0.02 s．Runge-Kutta方法采用MATLAB的Ode45函数作为求解器，并分别取2组容差δ=10-5和δ=10-10来控制其精度．
(a) 相轨迹 (b) Poincaré映射
(a) The phase portrait (b) The Poincaré map
图6 当
和θ=1.5时，本文方法给出的相轨迹和Poincaré映射
Fig. 6 When α=15, ζ=0.2, 
and θ=1.5, the phase portrait and Poincaré map obtained from the proposed method
图7 本文方法给出的相轨迹
Fig. 7 The phase portrait obtained from the proposed method
将Runge-Kutta方法取容差δ=10-10积分得到的结果作为参考解，分别计算本文方法和Runge-Kutta方法取容差δ=10-5的绝对误差．位移和速度的绝对误差分别定义为
δx=x(t)-xRK(t), δv=v(t)-vRK(t),
(26)
其中，xRK(t)和vRK(t)分别表示Runge-Kutta方法取容差δ=10-10得到的第一个质量块的位移和速度;x(t)和v(t)分别表示本文方法或Runge-Kutta方法取容差δ=10-5得到的第一个质量块的位移和速度．
(a) 位移的绝对误差 (b) 速度的绝对误差
(a) The absolute errors of the displacements (b) The absolute errors of the velocities
图8 本文方法的绝对误差
Fig. 8 The absolute errors of the displacements and velocities obtained from the proposed method
(a) 位移的绝对误差 (b) 速度的绝对误差
(a) The absolute errors of the displacements (b) The absolute errors of the velocities
图9 Runge-Kutta方法取容差δ=10-5的绝对误差
Fig. 9 The absolute errors of the displacements and velocities obtained from the Runge-Kutta method with δ=10-5
由本文方法积分得到的结构第一个节点在t∈[100, 200] s的相轨迹如图7所示．本文方法和容差为δ=10-5的Runge-Kutta方法的绝对误差分别如图8和9所示．可以发现， 本文方法的精度非常高， 比容差为δ=10-5的Runge-Kutta方法高约10倍， 这说明了本文方法的正确性．本文方法和Runge-Kutta方法取不同容差的计算时间如表1所示，可以发现，本文方法的计算效率分别约是Runge-Kutta方法取容差δ=10-5和δ=10-10的5倍和56倍．以上结果说明：本文方法对于求解具有大量自由度和大量间隙弹簧的分段线性系统的动态响应具有高精度和高效率．
表1 计算效率比较
Table 1 Efficiency comparison
本文提出一种求解具有大量间隙弹簧的周期性分段线性系统在简谐荷载作用下动态响应的高效率时域数值积分方法．通过采用参变量变分原理来描述间隙弹簧，分段线性问题被转化为线性互补问题，从而可以通过Lemke方法很容易求解．该算法避免了计算过程中的迭代和刚度阵更新，并能准确判断间隙弹簧的压缩和松弛状态．基于结构的周期性和能量传播速度的有限性，提出一种求解系统动态响应的高效率精细积分方法．该方法利用周期结构的周期性和能量传播有限性特点，指出周期结构的矩阵指数中存在大量的相同元素与零元素，从而不需要重复计算和存储这部分相同的元素，并避免计算和存储零元素，节省了计算量并降低了计算机存储要求．分析一个五自由度不对称分段线性系统，其在周期荷载作用下会发生非常复杂的动力学行为，如稳定的周期运动、准周期运动和混沌运动．通过求解一个大自由度的分段线性系统在周期荷载下的动态响应，比较了本文方法和Runge-Kutta方法的计算结果与计算时间，本文方法的正确性和高效率得到了验证．
参考文献(References)：
[1] INABA N, SEKIKAWA M. Chaos disappearance in a piecewise linear Bonhoeffer-Van der Pol dynamics with a bistability of stable focus and stable relaxation oscillation under weak periodic perturbation[J]. Nonlinear Dynamics, 2014, 76(3): 1711-1723.
[2] GOUZÉ J, SARI T. A class of piecewise linear differential equations arising in biological models[J]. Dynamics & Stability of Systems, 2002, 17(4): 299-316.
[3] KALM
R-NAGY T, CSIKJA R, ELGOHARY T A. Nonlinear analysis of a 2-DOF piecewise linear aeroelastic system[J]. Nonlinear Dynamics, 2016, 85(2): 739-750.
[4] MOTRO R. Tensegrity: Structural Systems for the Future[M]. London: Kogan Page Science, 2003.
[5] PUN D, LIU Y B. On the design of the piecewise linear vibration absorber[J]. Nonlinear Dynamics, 2000, 22(4): 393-413.
[6] 滕云楠, 韩明. 振动压实-土系统非线性滞回模型研究[J]. 真空, 2016, 53(4): 61-64.(TENG Yunnan, HAN Ming. Study on nonlinear hysteretic model of vibration compaction-soil system[J]. Vacuum, 2016, 53(4): 61-64.(in Chinese))
[7] ZHANG C. Theoretical design approach of four-dimensional piecewise-linear multi-wing hyperchaotic differential dynamic system[J]. Optik-International Journal for Light and Electron Optics, 2016, 127(11): 4575-4580.
[8] CHOI Y S, NOAH S T. Forced periodic vibration of unsymmetric piecewise-linear systems[J]. Journal & Sound Vibration, 1988, 121(1): 117-126.
[9] HUDSON J L, ROSSLER O E, KILLORY H C. Chaos in a four-variable piecewise-linear system of differential equations[J]. IEEE Transaction on Circuits & Systems, 1988, 35(7): 902-908.
[10] NAYFEH A H. Introduction To Perturbation Techniques[M]. New York: John Wiley & Sons, 1981.
[11] 吴志强, 雷娜. 分段线性系统的振动性能分析[J]. 振动与冲击, 2015, 34(18): 94-99.(WU Zhiqiang, LEI Na. Vibration performance analysis of piecewise linear system[J]. Journal of Vibration and Shock, 2015, 34(18): 94-99.(in Chinese))
[12] 钟顺, 陈予恕. 分段线性非线性汽车悬架系统的分岔行为[J]. 应用数学和力学, 2009, 30(6): 631-638.(ZHONG Shun, CHEN Yushu. Bifurcation of piecewise-linear nonlinear vibration system of vehicle suspension[J]. Applied Mathematics and Mechanics, 2009, 30(6): 631-638.(in Chinese))
[13] NARIMANI A, GOLNARAGHI M E, JAZAR G N. Frequency response of a piecewise linear vibration isolator[J]. Journal of Vibration Control, 2004, 10(12): 1775-1794.
[14] DESHPANDE S, MEHTA S, JAZAR G N. Optimization of secondary suspension of piecewise linear vibration isolation systems[J]. International Journal of Mechanical Sciences, 2006, 48(4): 341-377.
[15] MOUSSI E H, BELLIZZI S, COCHELIN B, et al. Nonlinear normal modes of a two degrees-of-freedom piecewise linear system[J]. Mechanical Systems & Signal Processing, 2015, 64(S1): 266-281.
[16] XU L, LU M W, CAO Q. Nonlinear vibrations of dynamical systems with a general form of piecewise-linear viscous damping by incremental harmonic balance method[J]. Physics Letters A, 2002, 301(1): 65-73.
[17] ZHONG W, ZHANG R. Parametric variational principles and their quadratic programming solutions in plasticity[J]. Computers & Structures, 1988, 30(4): 887-896.
[18] YU S D. An efficient computational method for vibration analysis of unsymmetric piecewise-linear dynamical systems with multiple degrees of freedom[J]. Nonlinear Dynamics, 2013, 71(3): 493-504.
[19] ACARY V, DE JONG H, BROGLIATO B. Numerical simulation of piecewise-linear models of gene regulatory networks using complementarity systems[J]. Physica D: Nonlinear Phenomena, 2014, 269(2): 103-119.
[20] ZHONG W, WILLIAMS F. A precise time step integration method[J]. Proceedings of the Institution of Mechanical Engineers, Part C: Journal of Mechanical, 1994, 208(63): 427-430.
[21] ZHONG W X. On precise integration method[J]. Journal of Computational and Applied Mathematics, 2004, 163(1): 59-78.
[22] SHA D, SUN H, ZHANG Z, et al. A variational inequality principle in solid mechanics and application in physically non-linear problems[J]. International Journal for Numerical Methods in Biomedical Engineering, 1990, 6(1): 35-45.
HE Dongdong, GAO Qiang, ZHONG Wanxie
(State Key Laboratory of Structural Analysis for Industrial Equipment(Dalian University of Technology); Department of Engineering Mechanics, Faculty of Vehicle Engineering and Mechanics, Dalian University of Technology, Dalian, Liaoning 116024, P.R.China) (Contributed by ZHONG Wanxie, M. AMM Editorial Board)
Abstract: An efficient method based on the parametric variational principle (PVP) was proposed for computing the dynamic responses of periodic piecewise linear systems with multiple gap-activated springs. Through description of gap-activated springs with the PVP, the complex nonlinear dynamic problem was transformed to a standard linear complementary problem. This method can avoid iterations and updating the stiffness matrix in the computing process and can accurately determine the states of the gap-activated springs. Based on the periodicity of the system and the precise integration method (PIM), an efficient numerical time-integration method was developed to obtain the dynamic responses of the system. This method indicates that there are a large number of identical elements and zero elements in the matrix exponents of a periodic structure, and saves computation load and computer storage by avoiding repeated calculation and storage of these elements. Numerical results validate the proposed method. The dynamic behaviors of a 5-DOF piecewise linear system under harmonic excitations were analyzed, including the stable periodic motion, the quasi-periodic motion and the chaotic motion. In comparison with the Runge-Kutta method, the proposed method has satisfactory correctness and efficiency.
Key words: piecewise linear system; dynamic response; parametric variational principle; periodicity; precise integration method
Foundation item: The National Natural Science Foundation of China(11572076；914748203)；The National Basic Research Program of China(973 Program) (2014CB049000)
文章编号：1000-0887(2018)07-0737-13
ⓒ 应用数学和力学编委会，ISSN 1000-0887
*收稿日期： 2018-02-04
基金项目： 国家自然科学基金(11572076；914748203)；国家重点基础研究发展计划(973计划)(2014CB049000)
作者简介： 何冬东(1988—)，男，博士生(E-mail: d.dong.he@mail.dlut.edu.cn);高强(1978—)，教授，博士，博士生导师(通讯作者. E-mail: qgao@dlut.edu.cn).
(我刊编委钟万勰来稿)
中图分类号： TB122; O322
文献标志码：A
DOI: 10.21656/1000-0887.390055
引用本文/Cite this paper:
何冬东， 高强， 钟万勰. 求解周期性分段线性系统动态响应的高效数值方法[J]. 应用数学和力学, 2018, 39(7): 737-749.
HE Dongdong, GAO Qiang, ZHONG Wanxie. An efficient numerical method for computing dynamic responses of periodic piecewise linear systems[J]. Applied Mathematics and Mechanics, 2018, 39(7): 737-749.