基于时间有限元的弹簧摆动力学新解法

陈德坤，丁幸波，温新婕，刘曼兰，兰 朋

(1. 哈尔滨工业大学机电工程学院，哈尔滨 150001；2. 军事科学院国防工程研究院，洛阳 471023；3. 北京起重运输机械设计研究院有限公司，北京 100007；4. 西安建筑科技大学机电工程学院，西安 710055)

弹簧摆问题是一种刚柔耦合的非线性动力学问题[1]，具有强耦合、强非线性等诸多特征。近年来随着国家电力建设的进一步推进，输电塔线系统的在恶劣环境下[2]，特别是在地震[3]，下击暴流[4]情况下的减振要求成为重点。为了研制较悬挂摆更好的减震器，催生了针对弹簧摆问题的深入研究。张鹏等[5]发现弹簧摆的减振能力优于悬挂摆，并对利用弹簧摆内共振特性设计减震器进行了探讨。霍林生等[6]、王奇等[7]、黄正等[8]在弹簧摆内共振特性上也开展了一定程度的研究。此外Wu等[9]利用弹簧摆的能量耦合特性研制了低频能量收集器。随着弹簧摆应用的深入，其模型自身的特性也受到了研究者多方面的关注。

鉴于传统方法存在一定缺陷，针对于哈密顿系统自身特性求解的时间有限元方法应运而生。由于有限元方法是基于变分原理获得的，故保辛。时间有限元最早是由Zienkiewicz等[23]提出的对于时间坐标的有限元，但并未要求保辛[24]。由于保守体系可用Hamilton体系描述，而保守体系的差分格式应当保辛[25]。故钟万勰[24]提出了保辛时间有限元的计算方法对于线性问题进行处理，体现了该方法的可靠性。吴峰等[26]利用了保辛摄动算法求解了刚-柔性杆的耦合问题，显示了保辛算法在处理小变形几何非线性问题的优势。

本文主要运用时间有限元以及矩形法作用量积分对小摆角弹簧摆进行处理，采用保辛递推，通过部分降低求解精度提高了计算速度，并能够保证积分过程的长时间稳定，保辛的特性能够避免系统误差随时间积分产生积累。

1 弹簧摆模型

时间有限元与保辛递推求解不同于传统的微分方程求解，其从动力方程的Hamilton变分原理出发，对系统作用量进行时间离散。本文以绕支点无阻尼自由摆动的弹簧摆为对象进行研究，研究模型如图1所示。采用直角坐标系描述弹簧摆的空间位置，设定二自由度坐标分别为x与y，质点质量m，弹簧在质点自重作用下的长度为l0，弹簧自然长度为r0，弹簧劲度系数为k，重力加速度为g，单位全部采用国际单位制。

图1 弹簧摆示意图Fig.1 Schematic diagram of spring pendulum

取弹簧摆的端点为原点，在此情况下建立系统动能和势能的表达式为：

作用量积分，即拉格朗日函数对于时间的积分，表示为：

这时可以采用时间有限元对于该问题进行描述，对于两个自由度x,y分别在时间区段ta,tb上进行线性插值可以构造一阶的时间有限元函数：

进行有限元变换的目的是，时间有限元变换是保辛的，能够保证系统的能量守恒。离散化后的作用量积分可以写成：

由于存在平方积分项，求解非线性方程组会耗费大量的时间。为了提高求解速度，本文在构造递推式时对平方积分进行线性化处理，采用矩形积分近似以求得近似解。单元作用量函数为：

其中，主要是对于平方积分项中的部分进行处理，用于保证后文的求解过程是前项递推的，因而只保留ta时刻的相关项，采用矩形积分方法对于该方程进行部分离散，有：

其余部分进行常规的有限元积分离散，有离散后的作用量格式：

其中：

采用保辛迭代格式求解，首先构造对偶变量：

引入状态向量并写成递推形式：

其中：

经过前面的矩形积分近似最终可以处理为S(xa,ya)的形式，可以代入初始状态向量进行递推求解。

一阶时间有限元与差分格式类似，不能显著提高有限元方法的精确度。可以通过改进时间有限元格式，增加时间有限元节点个数来提高计算精度，故而产生了二阶或者阶数更高的时间有限元。在维度更高的时间有限元推导过程中，采用矩阵形式和张量形式更为便捷，本方法中采用矩阵方法对于二阶时间有限元进行推导。

采用基函数N=[t2t1]进行时间有限元插值，即：

采用时间有限元法，即Ak1采用二阶时间有限元方法离散，Ak2、Ak3采用其他数值方法离散，这里要求的数值离散方法需要能够满足前文提到的递推要求，因而仅能通过ta端点进行构造。在本解法中采用矩形积分方法进行离散。

由于引入了中间节点qi对应，采用二阶时间有限元的推导需要进行节点凝聚以满足递推格式，而且可以推广到更高阶数的时间有限元上。需要注意，对于不采用有限元离散的非线性项部分，这部分在计算过程中不能引入中间节点进行节点凝聚的方法进行计算。因为下文给出的凝聚格式不能针对于非线性项进行处理，故而应采用两端时间节点位置进行离散。

节点凝聚格式由变分原理推导，由于Ak2、Ak3两作用量仅利用端点离散，与qi无关，故Ak的节点凝聚格式同Ak1的相一致，可以采用齐次方程下的节点凝聚法进行推导。通过可知对于形式为：

的作用量，可以转换为：

其中：

节点凝聚可以保证本方法构造成保辛递推格式，完成递推求解。在节点凝聚后仍然需要构造对偶变量并确立保辛递推格式以完成求解。

由于该计算方法的时间有限元部分Ak1保辛，误差的产生主要源于梯形积分在Ak2、Ak3处产生的误差。这里引入代数精度的概念，根据文献[28]，梯形积分仅具有一阶精度，是导致整个积分算法精度较低的一个关键因素。时间有限元又因为采用迭代算法不能在其迭代矩阵中引入非线性的vk，仅可选用vk-1相关值进行积分，故对积分精度有一定影响。本方法中构造了一种梯形算法以提高数值积分的积分精度。该积分算法采用ta时刻q(t)及其导数构造积分格式，即：

对于Ak2、Ak3采用该方程构造积分格式，由代数精度的定义可知其具有二阶代数精度，较矩形公式有进一步的精度提高。

2 模拟分析

为了检验方法的准确性，本文采用数值模拟同实验数据[15]进行对比。通过计算可知λ=2，处于内共振区间。同时与广泛使用的谐波平衡法[22]进行了对比，共进行2个状态下的数值模拟。

2.1 无阻尼情况

通过已有文献[10 − 11, 21 − 22, 29 − 31]，对于内共振现象的数值解法通常采用不考虑阻尼的方程进行求解。采用表1的实验参数进行了模拟，采用传递辛矩阵法，单步时间为0.002 s时的模拟结果图2、图3所示。

图2 摆球在两次模拟下的摆动曲线Fig.2 Swing curve of swinging ball in two simulations

图3 摆球位移随时间变化关系Fig.3 Relationship between swinging ball's displacement and time

表1 模拟参数1Table 1 Simulation parameters No.1

由图2～图4可以观察到以下几个现象：

图4 系统能量偏差随时间变化关系Fig.4 Relationship between system energy error and time

1) 对于振幅的分析：在初始状态1下出现了明显的共振现象，而初始状态2下内共振现象不明显，这与实验现象一致，与文献[22]中提到的只有在小角度摆动才出现强烈耦合的现象是一致的，而且该共振显现的周期与幅值同文献[22]中采用高阶Runge-Kutta法的结果也一致。这说明，本算法能够对内共振现象进行初步的预测。但实验中振动幅度存在明显的衰减现象，且内共振现象发生的周期与实验不一致。认为运动过程中存在阻尼对于内共振现象的发生是有影响的，而实验过程中的阻尼主要由空气阻尼或实验器材本身阻尼影响产生。需要进行进一步分析。

2) 对于频率的分析：40 s内的平均垂直振动频率为ωs=7.41Hz ， 水平振动频率ωp=3.65Hz，与实验中测定的ωs=7.27Hz ，ωp=3.63Hz分别存在1.93%与0.54%的误差。相较于其进行简谐运动近似的相对偏差4.9%与5.5%有较大提高。

3) 对于机械能守恒效果的分析：由于对保守系统进行分析，观察总能量的变化规律是评判算法好坏的重要指标。在该算法中能量偏差存在着周期性的变化，在计算步长为0.002 s时可以保证能量偏差与输入能量初值比维持在5‰以下。在计算过程中可以看出总能量偏差同初始能量之间的变化关系是周期性质的，而不像[10]中提到的对于微分方程求解方法的能量误差会出现发散的现象。这能够保证能量积分的相对稳定，对于长历程积分有很好适应性。

2.2 有阻尼情况

由于考虑引入阻尼力后，考虑阻尼系数γ=0.0013的情况对于表1中模拟次号1的初始状态进行分析，分析结果如下：

图5反映了含有阻尼力的摆动曲线与实验曲线更为类似，也更符合实验曲线的衰减规律。此外在增加系统阻尼的情况下，模拟轨迹出现了较之前的一些不同。首先，是阻尼系数导致振动频率的降低，增加阻尼后，两个方向振动频率分别降低至3.62 Hz与7.24 Hz，与实验计算的频率偏差分别为0.27%与0.41%，相较于不加阻尼分别提高了1.66%与0.13%，效果很好。其次，是阻尼系数的引入能够导致内共振现象发生频率产生变化，对于模拟次号1实验，内共振现象峰值出现频数由5次降低至4次，说明阻尼对内共振现象确实存在较大影响。最后，通过多次改变初始条件的值可以发现，通过使初值在一定范围(x方向不变，y方向0.01 m～0.1 m)内浮动可以发现，内共振现象的发生是依赖与初始条件的，且这种依赖现象同阻尼关系不大。在计算中可以发现，在y在区间0.05 m～0.06 m时，内共振现象很难发生，x方向与y方向几乎相对独立的进行振动。

图5 摆动曲线与方向位移同时间关系Fig.5 The relationship between direction displacement swing curve and direction displacement with time

2.3 时间有限元方法之间的对比

除一阶矩形法(仅对Ak2使用)外，一阶梯形(对Ak2、Ak3使用)法，二阶矩形法(仅对Ak2使用)以及二阶梯形法(对Ak2、Ak3使用)均给出了类似的结果，并对各组结果进行比较。表2采用计算步长0.2 ms，计算总时长20 s，分别对于第一次实验和第二次实验，针对于这四种方法的运算速度以及计算精度通过下表进行比较。

从表2可看出四种方法精度相近，采用增加阶数与增加代数精度的方法并没有显著增加计算精度，反而显著增加了计算时间。综上，采用一阶矩形法求解该问题即可获得理想的结果。

表2 积分方式对于模拟效果的影响Table 2 The influence of integral method on simulation

关于本方法的计算精度，整个推导过程存在2处近似，但仅采用梯形法积分是存在能量误差的。这部分会导致能量计算不够精确，是误差的主要来源。因此对总能量偏差与步长关系进行研究。

由图6可以看出，计算步长对于总能量偏差影响较大，在采用线性时间有限元的情况下在步长高于0.002 s时总能量偏差即超过5‰，因而在计算过程中需要尽量选取低于0.002 s的步长，以保证积分的准确性。

图6 总能量偏差同步长关系(模拟次号1)Fig.6 Relationship between total energy deviation and step size (simulation No.1)

3 同其他非线性方程解法的对比

时间有限元方法在计算精度相同的情况下同其他微分方程求解算法进行对比。计算机CPU采用四核i7-4710 MQ，主频2.5G Hz，内存4 GB。在计算步长为10 µs时计算2 s计算时间为0.64 s，且收敛速度较其他数值方法有明显提高。同表3中的非线性微分方程求解方法相比较，在λ =1，δ=6.2×10-7时，同4阶R-K法的计算速度类似，显著高于Ode15 s。确实在降低求解步长的情况下提高了计算速度，并保持了积分误差的稳定性。

表3 不同频率下数值方法比较[10]Table 3 Comparison of numerical methods at different frequencies[10]

4 对强非线性问题的适应性

表4 大角度摆动初始条件Table 4 Initial condition of wide-angle swing

在文献[11]中提到，在两振动方向频率比为1∶2的情况下，方程的解是初值敏感的，在不同的初值条件下能够表现出准周期性或混沌的特性。这里通过模拟次号2与模拟次号3相对照以检验该算法对于混沌问题的适应性。采用传递辛矩阵法，步进时间为0.002 s时的模拟结果如下：

从图7(b)和图7(c)的比较中可以发现在不同初值下，相图x方向出现了一定的混沌特征，这与文献[11]的结论一致。但从总能能量偏差的角度可以看出能量偏差虽然较小角度摆动的能量偏差更大，但在合理范围内且并没有出现发散趋势，通过前文的模拟也可知通过减小步长可以大幅降低总能量偏差。该模拟表明了时间有限元的保辛递推算法具有非线性长时间积分下稳定的特点。

图7 大角度摆动模拟状态Fig.7 Wide-angle swing simulation

5 结论

本文通过模拟证明了时间有限元能够处理弹簧摆的内共振问题。具有收敛较快、精度较高，并在长尺度积分下保持稳定的特性。

在对时间有限元步长进行了误差分析与速度分析后发现，该方法不存在误差累积的现象，计算精度主要取决于步长，在步长为0.001 s～0.01 s的范围内最大能量偏差小于2.17%，且能够保证计算速度。一阶矩形时间有限元的计算精度与时间经济性更好。此外对于更一般的大角度非线性弹簧摆问题也进行了求解，得到了长时间积分下误差可控的模拟结果。

由于该方法体现的刚柔耦合性不需要通过相对坐标进行参与，在绝对坐标下即可实现。虽相较于保辛摄动迭代法[26]精度有所下降，但更适于绝对节点坐标法的求解。