1. 引言
梁结构作为一种常见的工程结构,被广泛应用于房屋建设、桥梁工程、航空航天等领域。在实际生活中,梁自身容易受到外界因素的干扰而引起结构的振动,影响整体的稳定性,所以对梁振动问题的研究一直是人们关注的焦点。18世纪,文献 [1] 就简化线性弹性理论研究了梁的受力和变形问题,得到了Euler-Bernoulli梁方程。文献 [2] 研究了一维区间上非线性屈曲梁方程的初边值问题;文献 [3] 考虑刚体和柔性体运动相互作用的平移梁模型,推导出梁横向振动下的解析解,并在不同边界条件下验证了解析解的有效性;文献 [4] 针对轴向运动系统提出了一种简单的粘滞阻尼机制,研究了超临界转速范围内简支梁在横向激励下的动力学响应;文献 [5] 研究了两端固定轴向运动梁的横向振动,得到了固有频率和模态函数;文献 [6] 就弯曲和扭转联合作用下薄壁梁振动问题,研究了Sobolve空间中的一类非线性梁方程组的初边值问题;文献 [7] 考虑材料的不同拉压弹性模量,研究了铁木辛柯梁自由振动问题。上述研究中的模型方程多是一些高阶、非线性偏微分方程,在数值求解的精确性上研究较多,但在稳定性方面的研究还比较少。
L-稳定方法是一种基于Taylor展开研究数值稳定性的方法。文献 [8] 通过具有A-稳定的数值方法求解了刚性微分方程;文献 [9] 在此基础上进一步发展并将其应用到非线性模型方程;文献 [10] 构造了一类求解Stiff方程组的L-稳定高精度显示单步法;文献 [11] 针对刚性常微分方程组构造了L-稳定方法,该方法是二阶显式单步法;文献 [12] 通过W-变换构造了一类具有L-稳定和B-稳定的数值方法;文献 [13] 基于L-稳定的数值方法有效地抑制了求解过程中出现的高频振荡问题;文献 [14] 针对结构力学和多体系统动力学方程构造了L-稳定的块格式,分析了L-稳定块格式的精度;文献 [15] 提出了L-稳定实时协调算法的等效力控制方法并验证了其可行性。上述研究中,L-稳定方法多用于常微分方程与微分–代数方程稳定性问题的求解,尚未应用于求解梁振动这类高阶非线性偏微分方程。
文献 [16] [17] [18] 是求解梁振动方程的三种常见数值方法:有限差分法、有限元法、微分求积法(Differential Quadrature, DQ),具有求解格式简单或精度较高等优点。本文基于微分求积法进行改进,将偏微分方程的初边值问题转化为微分–代数方程求解,并应用L-稳定求解格式,以简支梁振动方程为例进行数值仿真。
2. 数学模型
高阶非线性梁振动偏微分方程的一般形式可以表示为
(1)
式(1)中,
是2维的未知函数,其关于时间t的偏导数的最高阶为二阶,G是关于函数
及其偏导数
的一个已知非线性函数,
是
的k阶偏导数:
,
,
是整数,
。
由文献 [19],根据梁两端连接情况,可以将边界条件大致分为四类,分别为
两端固支:
(2)
两端简支:
(3)
两端滑支:
(4)
两端自由:
(5)
将边界条件看作约束函数
,引入拉格朗日乘子
,可以将带有简支边界条件的梁振动偏微分方程转化为指标-3微分–代数方程组:
(6)
其中,
是
关于
的偏导数的转置,对
分别关于时间t求一、二阶导数得速度级约束
,加速度级约束
,令
,
,可以得到指标-2、指标-1的微分–代数方程组:
(7)
(8)
3. L-稳定方法
将区间
离散为:
。设函数
在
上k阶可导,选取三角插值基函数,根据插值定理,通过求导运算,则函数
的k阶偏导数
可以表示为
(9)
这里,
是
关于
的
阶权系数矩阵,不同阶权系数矩阵之间满足递推公式:
(10)
在区间
上取等分节点:
,步长设为
,
。将式(9)、(10)代入式(7),再进行降阶可得
(11)
在区间
上取等分节点:
,则有
,
。设
,
,
,
。由文献 [20],根据Ehle定理及猜想,通过有理逼近构造L-稳定求解格式:
(12)
其中,
,
,
是全1矩阵,
是单位矩阵,
,
,
,
为Kronecker积。将式(12)代入式(11),通过牛顿迭代得到
,代入式(12)可得
。
特别地,r = 4时,可得
(13)
4. 数值算例
无轴向运动简支梁在外力作用下的振动方程为
(14)
其中,
为弯曲刚度,
为梁的质量密度,A为梁的横截面积,
为梁外载荷,l为梁的横长度。假定相关参数为
初始条件为
4.1. 梁振动方程时间历程图分析
图1是仿真时间为30 s,时间步长h = 0.01,空间节点选取7个切比雪夫节点,对指标-2方程求解得到的梁中点位移时程图;图2是对应的约束函数、速度级约束、加速度级约束误差时程图。可以看出,图1展示了简支梁在外激励作用下的来回往复运动轨迹,位移图像的不规则变化也体现了梁振动方程的非线性特征;图2的约束函数中的四个约束分量的误差量级都在10~12,这表示边界条件得到了较好地满足,其余约束的误差也很小,说明L-稳定方法在仿真过程中保持了良好的稳定性。
![](//html.hanspub.org/file/6-2622006x78_hanspub.png?20220111080647714)
Figure 1. Time history diagram of beam midpoint displacement
图1. 梁中点位移时程图
![](//html.hanspub.org/file/6-2622006x79_hanspub.png?20220111080647714)
Figure 2. Time history diagrams of constraint function and velocity level constraint and acceleration level constraint error
图2. 约束函数、速度级约束、加速度级约束误差时程图
4.2. L-稳定方法在不同参数下的结果比较
表1是仿真时间为30 s,空间节点选取7个切比雪夫节点,时间步长h = 0.01,不同指标方程下的对比结果。可以看出,指标-1方程在加速度级约束上误差最小,指标-2方程在速度级约束上误差最小,指标-3在约束函数上误差最小,这体现了不同代数约束对指标方程约束误差的影响;综合三种约束误差,指标-2方程对约束的满足程度更高。
![](Images/Table_Tmp.jpg)
Table 1. Comparison of indexes-1, -2, -3 solutions to differential algebraic equations with L-stable method
表1. L-稳定方法求解指标-1、-2、-3的微分–代数方程组的结果比较
表2是仿真时间30 s,空间节点为切比雪夫节点,时间步长h = 0.01,针对指标-2方程关于空间节点数目的结果比较:节点数取5时,无法求解,节点数从7开始,随着节点数增多,约束误差逐渐增大,运行时间延长,故而节点选取7时,求解时间最短,约束保持更好。
![](Images/Table_Tmp.jpg)
Table 2. Comparison of results of L-stable method with different space node numbers
表2. L-稳定方法在不同的空间节点数下的结果比较
表3是仿真时间30 s,空间节点选取7个切比雪夫节点,针对指标-2方程在不同时间步长下的结果比较:h = 0.1时,稳定性最好,运行时间最短,随着时间步长变小,各级约束误差虽然增大,但变化幅度较小,仍保持良好的稳定性。
![](Images/Table_Tmp.jpg)
Table 3. Comparison of results of L-stable method with different time steps
表3. L-稳定方法在不同时间步长下的结果比较
4.3. L-稳定方法与龙格–库塔法、微分求积法的结果比较
图3是仿真时间30 s,空间节点选取7个切比雪夫节点,时间步长h = 0.1,L-稳定方法与四阶龙格–库塔法(RK4)、DQ法的梁中点位移时程图,DQ法采用权系数修正处理边界条件。可以看出h = 0.1时,L-稳定方法仍能求得梁中点位移近似解,而另两种方法求解的结果发散,这体现了L-稳定方法在大步长下仍具有很好的求解精度,求解效果更好。
图4是时间步长为h = 0.01,空间节点选取7个切比雪夫节点,L-稳定方法与RK4法、DQ法(T = 1000 s)时仿真结果对比图。图5是相应条件下L-稳定方法的约束时间历程图。可以看出,L-稳定方法与DQ法求解结果较为接近,而RK4法求解结果差异稍大;在1000 s时,其约束误差较小,长时间仿真精度高,稳定性好。
![](//html.hanspub.org/file/6-2622006x80_hanspub.png?20220111080647714)
Figure 3. Time history diagram of beam midpoint displacement of L-stable method, RK4 method and DQ method (h = 0.1)
图3. L-稳定方法与RK4法、DQ法的梁中点位移时间历程图
![](//html.hanspub.org/file/6-2622006x81_hanspub.png?20220111080647714)
Figure 4. Comparison time history diagrams of beam midpoint displacement of L-stable method, RK4 method and DQ method (T = 1000 s)
图4. L-稳定方法与RK4法、DQ法的梁中点位移时间历程图(T = 1000 s)
![](//html.hanspub.org/file/6-2622006x82_hanspub.png?20220111080647714)
Figure 5. L-stable method constraint time history diagram (T = 1000 s)
图5. L-稳定方法约束时间历程图(T = 1000 s)
4.4. L-稳定方法与RK4法能量图分析
考虑去除外部激励,给定一个初始扰动:
图6是仿真时间30 s,空间节点选取7个切比雪夫节点,时间步长h = 0.01,L-稳定方法与RK4法总能量对比图。可以看出,L-稳定方法能量总体保持在一个较小的范围内,而RK4法的总能量随时间延长而逐渐增大,误差累积较大。
![](//html.hanspub.org/file/6-2622006x84_hanspub.png?20220111080647714)
Figure 6. Comparison diagrams of total energy between L-stable method and RK4 method
图6. L-稳定方法与RK4法总能量对比图
5. 结论
本文针对高阶非线性梁振动偏微分方程的一般形式,基于插值定理进行空间离散,给出了偏微分方程求解的L-稳定格式。数值结果表明,L-稳定方法求解结果较好地展示了梁中点振动的非线性特征;针对指标-2方程,L-稳定方法的约束保持整体最优;时间步长对约束误差的影响较小,可以在大步长下满足各级约束;空间节点数目取7时稳定性最好;与RK4法、DQ法相比,在大步长下L-稳定方法求解精度更高,约束误差更小,同时长时间求解也能保持较高的计算精度和数值稳定性;忽略外部激励,给定初始扰动,L-稳定方法求解的系统总能量保持在一个较小的范围内,而RK4法的能量无法保持,误差累计较大。此方法构造原理简单,特别在边界处理上能更好地满足约束条件,后续可进一步加以改进,提高求解效率。
基金项目
国家自然科学基金(11772166, 12172186)。
NOTES
*第一作者。
#通讯作者。