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

由于弹簧摆兼具耦合性，强非线性等诸多特性，弹簧摆的运动特征受其初始位置、质点质量、弹簧刚度以及弹簧在自重作用下的长度影响很大。一般对于内共振特性的研究主要采用振动模态圆频率${\omega _{\rm{s}}} = \sqrt {\dfrac{k}{m}}$与摆动模态圆频率${\omega _{\rm{p}}} = \sqrt {\dfrac{g}{{{l_0}}}}$两个重要参数的比值进行描述。在$\lambda = {\omega _{\rm s}}{\text{∶}}{\omega _{\rm p}} =$$1{\text{∶}}1$时系统呈现准周期运动，而在$\lambda = {\omega _{\rm s}}{\text{∶}}{\omega _{\rm p}} = 1{\text{∶}}2$时系统则出现明显的内共振现象，在比值更大时会出现混沌现象。于洪洁等^[10]、李银山等^[11]、赵聪等^[12]对内共振现象与混沌现象的产生原因进行了深入研究。国际上van de Weele等^[13]通过频率比和能量比研究了弹簧摆的有序与混沌的动力学过程，de Sousa等^[14]研究了参数空间下的能量耦合传递情况，提出了一种应用于非线性耦合系统，揭示耦合如何介导内部能量交换，以及能量分布如何随系统参数而变化的方法。司丽荣等^[15]、Rusbridge等^[16]对于弹簧摆的内共振问题进行了实验，为弹簧摆在不同λ下的内共振规律提供了实验支撑。Grażyna Sypniewska-Kamińska等^[17]采用给定区间逼近三角函数，并提出了基于此的多尺度法，提高了近似效果的精度。Amer等^[18-19]采用多尺度法，对于含阻尼的弹簧摆运动形式进行了持续研究，分析了相关物理参数对解的影响。

弹簧摆的特性决定了其求解属于非线性方程的求解。现有弹簧摆问题主要有解析法和数值两种方法，其中对于混沌弹簧摆问题的求解的解析方法采用的是基于Hamilton系统的Holmes方法^[20]，该方法是少数能够精确求解混沌系统的方法。李银山等^[11]通过此方法求解了不同λ下弹簧摆的相图。在2000年后随着计算机技术的成熟，国内对于弹簧摆的研究多采用非线性微分方程组的数值方法^[9-10]与谐波平衡法^[21-22]进行求解。于洪洁等^[10]详细地研究了现有不同类型非线性常微分方程组解法并给出了计算效率与计算精度。对于内共振问题计算则多采用谐波平衡法处理。这两类方法各有千秋，非线性方程组的数值解法多以Runge-Kutta法，Gear法与Adams法为基础^[10]，由于这类求解方法没有体现哈密顿系统本身的保守特性，因而容易产生累计误差，还容易出现系统的不稳定甚至发散的情况。而针对于谐波平衡法，其对于弹簧摆的初始位置有一定要求，一般要求$\left| {\dfrac{x}{{{l_0}}}} \right| \ll 1$，$\left| {\dfrac{{y + {l_0}}}{{{l_0}}}} \right| \ll 1$。这限制了其对于在大摆角情况下的分析能力。

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

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

1   弹簧摆模型

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

图  1  弹簧摆示意图

Figure  1.  Schematic diagram of spring pendulum

下载: 全尺寸图片 幻灯片

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

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

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

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

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

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

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

其中：

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

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

其中：

经过前面的矩形积分近似最终可以处理为${\boldsymbol{S}}\left( {{x_{\rm{a}}},{y_{\rm{a}}}} \right)$的形式，可以代入初始状态向量进行递推求解。

以上是针对无阻尼系统的推导，对于有阻尼系统，时间有限元也能展现其适应性。对文献[27]中提到的阻尼力与时间成正比阻尼形式，若阻尼系数不是时变参数，可采用变量代换方法求解。即令阻尼系数为γ，${{\boldsymbol{q}}_1}\left( t \right) = {\boldsymbol{q}}\left( t \right){{\rm e}^{ - \tfrac{\gamma }{m}t}}$。利用该变换可以消除耦合项并将作用量${A_k}$转化为：

转化后的${A_{k1}}$与${A_k}$形式类似，在${A_{k11}}$部分可以直接利用时间有限元方法求解；在${A_{k21}}$、${A_{k31}}$两部分可以继续通过构造梯形积分公式的形式进行积分。只需要再对求得的${{\boldsymbol{q}}_1}\left( t \right)$变量乘系数${{\rm e}^{\tfrac{\gamma }{m}t}}$即可获得原变量的解。

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

采用基函数${\boldsymbol{N}} = [ {{t^2}}\;\;t\;\;1]$进行时间有限元插值，即：

采用时间有限元法，即${A_{k1}}$采用二阶时间有限元方法离散，${A_{k2}}$、${A_{k3}}$采用其他数值方法离散，这里要求的数值离散方法需要能够满足前文提到的递推要求，因而仅能通过${t_{\rm{a}}}$端点进行构造。在本解法中采用矩形积分方法进行离散。

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

节点凝聚格式由变分原理推导，由于${A_{k2}}$、${A_{k3}}$两作用量仅利用端点离散，与${{\boldsymbol{q}}_i}$无关，故${A_k}$的节点凝聚格式同${A_{k1}}$的相一致，可以采用齐次方程下的节点凝聚法进行推导。通过可知对于形式为：

的作用量，可以转换为：

其中：

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

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

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

2   模拟分析

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

2.1   无阻尼情况

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

表  1  模拟参数1

Table  1.  Simulation parameters No.1

因素	模拟次号
	1	2
水平初位置 x0/m	0.1414	0.1960
竖直初位置 y0/m	−0.9290	−0.8774
水平初动量 px0/(kg·m/s)	0.0000	0.0000
竖直初动量 Py0/(kg·m/s)	0.0000	0.0000
初始摆长 r0/m	0.6162	0.6162
振子质量 m/kg	0.0484	0.0484
弹簧刚度 k/(N/m)	2.3000	2.3000

下载: 导出CSV 

| 显示表格

图  2  摆球在两次模拟下的摆动曲线

Figure  2.  Swing curve of swinging ball in two simulations

下载: 全尺寸图片 幻灯片

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

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

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

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

图  3  摆球位移随时间变化关系

Figure  3.  Relationship between swinging ball's displacement and time

下载: 全尺寸图片 幻灯片

图  4  系统能量偏差随时间变化关系

Figure  4.  Relationship between system energy error and time

下载: 全尺寸图片 幻灯片

2.2   有阻尼情况

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

图  5  摆动曲线与方向位移同时间关系

Figure  5.  The relationship between direction displacement swing curve and direction displacement with time

下载: 全尺寸图片 幻灯片

图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方向几乎相对独立的进行振动。

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

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

表  2  积分方式对于模拟效果的影响

Table  2.  The influence of integral method on simulation

积分方法	模拟次号1
	精度/(×10−5 J)	精度偏差/(%)	时间/s	时间偏差/(%)
一阶矩形	6.2154	0.000	0.9002	0.00
一阶梯形	6.2154	0.000	1.2032	33.66
二阶矩形	6.2173	0.030	1.4258	58.39
二阶梯形	6.2175	0.033	1.7018	89.05
积分方法	模拟次号2
	精度/(×10−5 J)	精度偏差/(%)	时间/s	时间偏差/(%)
一阶矩形	7.8835	0.000	0.9047	0.00
一阶梯形	7.8835	0.000	1.1992	32.55
二阶矩形	7.8859	0.030	1.4203	56.99
二阶梯形	7.8863	0.036	1.5658	73.07
注：偏差是相对本列中最小量取得的，计算方法为偏差量除以最小量。

下载: 导出CSV 

| 显示表格

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

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

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

图  6  总能量偏差同步长关系(模拟次号1)

Figure  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$，$\delta = 6.2 \times {10^{ - 7}}$时，同4阶R-K法的计算速度类似，显著高于Ode15 s。确实在降低求解步长的情况下提高了计算速度，并保持了积分误差的稳定性。

表  3  不同频率下数值方法比较^[10]

Table  3.  Comparison of numerical methods at different frequencies^[10]

数值方法	λ=1
	步长h/s	所需误差精度 ε	能量相对误差 δ	计算时间 t/s
Adams	10−2	10−6	1.86×10−6	0.24
4阶R-K	10−2	10−6	8.22×10−8	0.31
Ode113	−	10−10	1.99×10−10	0.42
Ode45	−	10−10	2.28×10−8	0.64
Ode15 s	−	10−10	5.44×10−7	1.50

下载: 导出CSV 

| 显示表格

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

由于不要求$\left| {\dfrac{x}{{{l_0}}}} \right| \ll 1$，$\left| {\dfrac{{y + {l_0}}}{{{l_0}}}} \right| \ll 1$。本方法除了可以解决耦合共振问题，也可以求解几何非线性问题即大角度摆动产生的非线性问题^{[16, 32]}，本算例通过限制总能量偏差的方法求解大角度摆动问题，设计条件如表4所示。

表  4  大角度摆动初始条件

Table  4.  Initial condition of wide-angle swing

因素	模拟次号
	3
水平初位置 x0/m	0.7000
竖直初位置 y0/m	−0.4000
水平初动量 px0/(kg·m/s)	0.0000
竖直初动量 Py0/(kg·m/s)	0.0000
初始摆长 r0/m	0.6162
振子质量 m/kg	0.0484
弹簧刚度 k/(N/m)	2.3000

下载: 导出CSV 

| 显示表格

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

下载: 全尺寸图片 幻灯片

图  7  大角度摆动模拟状态

Figure  7.  Wide-angle swing simulation

下载: 全尺寸图片 幻灯片

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

5   结论

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

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

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