一种基于Hamel 形式的无条件稳定动力学积分算法1)

顾崴 刘 铖†, 安志朋 史东华,†††

* (北京理工大学信息与电子学院,北京 100081)

† (中国工程物理研究院高性能数值模拟软件中心,北京 100088)

** (北京应用物理与计算数学研究所,北京 100088)

†† (北京理工大学宇航学院,北京 100081)

*** (北京理工大学数学与统计学院,北京 100081)

††† (复杂信息数学表征分析与应用北京市重点实验室,北京 100081)

引言

时间积分算法被广泛应用于动力学系统的求解中.动力学方程的时间积分经常会出现数值不稳定现象,有限元空间离散也通常会造成伪高频振荡,因此,针对上述问题的数值积分算法的研究一直备受关注.本文将基于活动标架法和离散力学,发展对动力学系统无条件稳定的数值积分格式.

Newmark 方法[1]是广泛使用的时间积分方法,它包括许多格式,比如中心差分法(CDM)和梯形规则(TR)等,其中隐式TR 具有二阶精度和无条件稳定性.此外,对于线性保守系统,TR 可以长时间保持能量.但对于波传播和刚性系统,TR 不能过滤掉不需要的高频.为解决此问题,研究人员开发了更理想的方法,如广义 α 方法[2]、WBZ-α 方法[3]和HHT-α方法[4].

广义 α 方法[2]是求解结构动力学问题的一种隐式方法,可以通过单个参数来控制数值耗散,使得在耗散系统高频响应的同时较好地保持系统的低频响应.广义 α 方法在时间上具有二阶精度且无条件稳定,已成功用于广泛的工程应用[5-8].为应用于多体动力学,Yen 等[9]将其扩展到约束系统.文献[10]分析了广义 α 格式和Lagrange 乘子的二阶收敛性.张雄和王天书[11]和田强[12]比较了上述算法的精度和耗散特性,研究了力学系统守恒量的变化特性.此外,文献[13]将二阶广义 α 方法推广到三阶,并且可以通过单个参数控制数值耗散,同时也证明了该方法是无条件稳定的.季奕等[14-15]提出了非线性动力问题的单步以及两步可控耗散时间积分法,这些方法是二阶精度和无条件稳定的,并且不需要初始加速度矢量.

为了研究具有李群结构的柔性多体系统,Brüls等[16]提出了李群广义 α 时间积分算法,并证明了该方法具有全局二阶精度.Arnold 等[17]详细分析了李群广义 α 方法的收敛性.文献[18]给出了SE(3)群的广义 α 方法.文献[19] 基于BDF(backward differentiation formulae)方法提出了更高精度的李群积分BLieDF 方法.

在几何力学领域,变分积分子[20-21]也是一种数值积分算法.它源于离散力学中离散Hamilton 原理,能如实反映动力学系统性质,目前已被广泛应用于Lagrange 系统[22]和约束系统[23]等各种问题,并取得良好的效果.变分积分子的核心思想是先对系统相应的变分原理进行离散,然后通过离散变分运算,得到连续Euler-Lagrange 方程的离散形式.Lew 等[21]指出,在有限维情形下,由于直接从离散的变分原理出发,变分积分子能够如实反映连续系统的动力学性质,如保守系统的保动量和保辛结构,保守和非保守系统的离散能量守恒,因而避免了传统数值方法存在的数值耗散问题,在数值上体现出更好的性质.随着离散力学的发展,研究人员提出了辛能量–动量变分积分子[24]和李群变分积分子[25]等.文献[26]结合变分积分子,提出了二阶数值时间积分方法RATTLie,并将它们应用于构型空间为李群的约束力学系统.同时,研究人员又发展了高阶变分积分子,如Galerkin 变分积分子[27]和谱变分积分子[28]等,这些高阶变分积分子既能保辛和保动量,又能实现任意阶精度或具有指数收敛性.

Bloch 等[29]将活动标架法融入变分原理,推导了有限维力学系统的Hamel 动力学方程,并提出了Hamel 变分积分子[30],其框架可统一描述李群和李代数变分积分子.Hamel 变分积分子产生了一种更简单的算法,目前已应用于机械臂最优控制问题[31],展现了良好性能.Shi 等[32]提出了无穷维力学系统的Hamel 形式,进而研究了经典场论的Hamel 形式[33-34],并在此基础上构建了Hamel场变分积分子.王亮等[35]利用Hamel 场变分积分子,给出了几何精确梁的动力学建模和数值仿真,该算法能长时间保持系统的能量、动量和几何结构,高山等[36]对其保持动量给出了理论证明.因此,Hamel场变分积分子可以对一类无穷维力学系统进行准确刻画,是力学系统长时间数值模拟的理想选择.

Newmark 方法和广义 α 方法等一系列时间积分算法最初均是通过经验推导或总结出来的,后来赋予了各种理论基础,如泰勒级数和Galerkin 弱形式等,最后利用数值算法的分析手段确定其收敛性和稳定性.这些方法都不是为了保持能量和动量而设计的,在一些特殊参数情况下,它们可以保持动量[37].文献[38]从Hamilton 方法出发,指出Newmark 族以及相关的积分算法,在离散力学的Veselov 形式下才是变分的[38].迄今为止,在变分积分子的框架下一直未出现一种能耗散系统高频响应同时又能较好保持系统低频响应的数值积分算法.因此,在该领域中,构造一种无条件稳定的时间积分算法是必要的,这样既可以拓展变分积分子的理论体系,同时又可以丰富时间积分算法的理论基础.在这种背景下,本文基于离散力学框架和活动标架法,采用特殊的变分形式[39-40],提出了一种Hamel 广义 α 方法,并推广到了一般李群空间,理论分析和数值测试均显示出了本文方法在精度、耗散和稳定性方面的优势.从而为进一步研究无条件稳定的高阶数值积分算法提供参考依据.

1 Hamel 形式与Hamel 变分积分子

1.1 Hamel 场形式

令Q为向量空间W上的无穷维光滑流形.Lagrange函数为L:TQ→R,其中TQ是位形空间Q的切丛.系统的动力学由Euler-Lagrange 方程组决定.若选择合适的标架,可获得系统动力学的更简洁形式——Hamel 方程[34].令U是q∈Q的邻域,定义一族可逆的有界线性算子 Ψq:W→TqQ为活动标架算子,∀q∈U,同时也是可逆的有界线性算子.因而,对于 ∀ξ ∈W定义向量场

给定两个向量 ξ,η ∈W,定义一个反对称括号算子[·,·]q:W×W→W

其中 [·,·] 是流形上两个向量场的Jacobi-Lie 括号.对于任意 ξ,η ∈W和 α ∈W*,括号算子 [·,·]q的对偶算子通过下式给出

给定非保守力F∈Γ(T*Q),其中Γ (T*Q) 表示T*Q的所有截面,则虚功记为

定理1 Hamilton-Pontryagin 原理的Hamel 形式[34]

令力学系统位形空间为Q,TQ是切丛,T*Q为余切丛.令q∈Q,(q,v)∈TQ且 (q,p)∈T*Q.令L:TQ→R为Lagrange 函数,l是关于活动标架下局部坐标(q,ξ)的Lagrange 函数,则下列陈述等价.

(1) 令 (q,v,p) 为Pontryagin 丛TQ⊕T*Q的局部坐标,则曲线 (q(t),v(t),p(t)),0 ≤t≤T为其相应的作用积分

关于独立变分 δq,δv,δp的临界点,即

且δq(0)=δq(T)=0.

(2) 曲线 (q(t),v(t),p(t)) 满足隐式Euler-Lagrange 方程组

(3) 曲线 (q,ξ,µ)∈U×(W⊕W*) 是作用泛函

关于独立变分 δq=Ψqη ,δ ξ 和 δµ 的临界点

(4) 曲线 (q,ξ,µ) 满足隐式Hamel 方程组的弱形式

对于任意点q∈Q,曲线 (q,ξ,µ) 满足隐式Hamel 方程组的强形式

1.2 Hamel 变分积分子

假设Q是向量空间,连续曲线q:[0,T]→Q用离散序列近似表示.定义

式中参数 α ∈[0,1],Δt是时间步长,令

为关于活动标架 Ψqk+α的速度分量.

通过连续Lagrange 函数,定义其离散形式

相应离散作用和为

同理,变量 η 的离散形式 ηk+α可以线性展开

通过离散变分原理,可得隐式Hamel 方程组强形式(12)的离散形式—–Hamel 变分积分子

其中Dil为Lagrange 函数对第i个参数的偏导数.

2 向量空间上的Hamel 广义 α 方法

本节考虑向量空间中力学系统,基于活动标架法,构造力学系统的Hamel 广义 α 方法.

2.1 Hamel 方程

考虑一个受控系统,其位形空间为向量空间Q,状态空间为TQ,T*Q表示相空间.对 ∀q∈Q取标架算子 Ψq为恒算子,定义系统的Lagrange 函数L:TQ→R

其中fτ,fext分别为反馈控制力和外力.

由前文定理1,可知下面的表述是等价的.

(1) 定义在TQ⊕T*Q上的曲线(q(t),v(t),p(t)) 是作用泛函

关于独立变分 δq,δv,δp的临界点,即

(2) 曲线 (q(t),v(t),p(t)) 满足隐式Hamel 方程组的弱形式

曲线 (q(t),v(t),p(t)) 满足隐式Hamel 方程组的强形式

此时Hamel 方程组为隐式Euler-Lagrange 方程组.所以由式(24)可定义加速度函数

2.2 Hamel-广义 α 方法

给定时间步长 Δt,区间[ 0,T] 等分为N个子区间,为简单起见,令区间 [tn,tn+1] 为 [0,Δt],为构建数值积分算法,定义离散形式

并且满足Mv0=p0,最后根据加速度函数定义相应的离散加速度

此外,令外力和反馈控制律离散形式为

因此,可知力f所对应的离散形式为

通过控制律fτ可获得无条件稳定的数值积分算法.定义辅助未知量 δq,δp和 δv的离散形式

其中,A1,A2,A3,A4为四个参数.将离散形式(26)～式(32)带入隐式Hamel 方程组的弱形式(23),求解关于的驻点方程,可得到比广义 α 方法更一般的数值积分算法

式中参数分别为

其次,tn+1时刻离散加速度函数为

为了得到更稳定的算法,令

综上,方程(33)～(35)称之为Hamel 广义 α 方法.

注1.Hamel 广义 α 方法可以退化为一些经典算法,如Newmark、HHT-α 和广义 α 方法等.如表1 所示,令A3=2(3β-1),A4=-3(4β-1),表格中参数α,αm,αf的含义参见文献[1-2,4].

表1 Hamel 广义 α 方法与经典算法比较Table 1 Comparison between Hamel generalized-α method and classical algorithm

注2.对离散形式(26)进行改进如下

可得到更高阶数值积分算法(33),其中参数µ0,µ1,γ0,γ1变为

其余参数不变,其中I为单位阵,K=∂2V(∂qi∂qj) 为弹性矩阵.选取特定的参数,该数值格式会得到更高的收敛速率,具体数值算例及其收敛速率见3.3 节.

3 算法性质

3.1 算法分析

考虑一个力学系统

相应的Euler-Lagrange 方程组为

Hamel 广义 α 方法应用于上述方程,令f=0,则离散方程可写成矩阵形式

故本文所提方法为二阶精度算法.

算法的稳定性和数值耗散取决于其辅助矩阵的特征值.谱半径是数值耗散的一种度量,谱半径越小,数值耗散越大.算法的谱半径由矩阵的特征值决定,定义为

其中 λi为矩阵的第i个特征值.若 ρ ≤1,则对于线性问题算法是无条件稳定的.

理想算法应该是严格控制高频,但是低频区域不引起过度耗散,在低频域的谱半径值接近于1.随着 Ω 的增加,谱半径值平滑地减小.若辅助矩阵的主根形式为

其中a和b是实数,只要

成立,则高频耗散最大化.由(41)可知

满足这个条件,Hamel 广义 α 方法便可实现高频耗散最大化.此外,当A=0.5 时,可得A4=0,无法确定A3的值,故而设定A3=-0.5.

令ρ∞表示高频极限处的谱半径,变量可表示为

为方便描述方程(44),借助 ρ∞,定义B1和B2

综上所述,同时满足方程式(40)、式(42)和式(45),则Hamel 广义 α 方法是无条件稳定、二阶精度的算法,并且可实现高频耗散最大,低频耗散最小.

3.2 算法稳定性

谱分析是分析时间积分算法稳定性与精度的常用方法.下文通过力学系统(37)来探讨算法性质.令式(38)中 ω=2π.若 0 ≤ρ ≤1,则Hamel 广义 α 方法无条件稳定;在此基础上,若 ρ ＜1,该方法是数值耗散的,当ρ=1 时,该方法是非耗散的.定义数值阻尼比和周期延长率为

图1 展示了Hamel 广义 α 方法的谱半径,发现是无条件稳定的,并且在高频区域内可通过 ρ∞精确控制数值耗散程度.

图1 谱半径Fig.1 Spectral radius

在0 ≤ρ∞≤1 情况下,当 Ω ≤0.1 时,谱半径接近于1,保持了基本的低频模式,随后快速降至 ρ∞,过滤掉不需要的高频信息.图2 和图3 分别给出了该方法的数值阻尼比和周期延长率,可以发现振幅和相位精度随着 ρ∞的增加而提高.

图2 数值阻尼比Fig.2 Numerical damping ratio

图3 周期延长率Fig.3 Period elongation

3.3 算法收敛性

本节通过数值算例测试算法的收敛速度.收敛速度是特定时间内随着时间步长递减,绝对误差的下降率.令系统方程(38)式中 ω=π.初始条件为

该系统的精确解为

通过Hamel-广义 α 方法t=1 时刻的绝对误差,如图4可以发现本文方法关于位移是二阶精度,而且随着ρ∞增大,算法准确度也在提升.

图4 Hamel 广义 α 方法t=1 时刻的误差Fig.4 Errors of Hamel generalized-α method at time t=1

若A=0.5,A3=-1,A4=1,B1=B2=0,通过Hamel广义 α 方法t=1 时刻绝对误差,如图5,发现注2 的数值积分格式关于位移接近于四阶精度.

图5 Hamel 广义 α 方法t=1 时刻的误差Fig.5 Errors of Hamel generalized-α method at time t=1

4 数值算例

本节通过数值算例,在 ρ∞=0.8 的情况下将Hamel 广义 α 方法与广义α 方法,Newmark 方法(β=0.25,γ=0.5)的结果做对比,来验证其性质.

考虑一个两自由度系统,其Lagrange 函数为

系数矩阵为

其中m1=m2=1 和k1=104,k2=1.

相应Euler-Lagrange 方程组为

采用如下初值

图6 和图7 分别给出了系统1 和系统2 三种方法的运行结果,通过位移、速度和加速度曲线,可以看到Hamel 广义 α 方法和广义 α 方法在稳定性方面要比Newmark 方法更具优势.图8 和图9 分别给出了系统1 和系统2 Hamel 广义 α 方法和广义 α 方法两种方法的数值误差,通过与经典算法对比发现,本文提出的方法具有良好的稳定性和精度性质.

图6 系统1 各种方法数值结果对比图Fig.6 Comparison of numerical results of various methods for system 1

图7 系统2 各种方法数值结果对比图Fig.7 Comparison of numerical results of various methods for system 2

图8 系统1 Hamel 广义 α 方法与广义 α 方法数值结果误差Fig.8 Errors of numerical results between Hamel generalized-α method and generalized-α method for system 1

图9 系统2 Hamel 广义 α 方法与广义 α 方法数值结果误差Fig.9 Errors of numerical results between Hamel generalized-α method and generalized-α method for system 2

5 李群空间上的Hamel-广义 α 方法

本节考虑李群空间中力学系统的Hamel 形式,并提出相应的运动方程—–Hamel 方程,在此基础上,构造李群空间上力学系统的时间积分算法.

5.1 Hamel 方程

考虑李群空间上的力学系统,其位形流形为李群G,状态空间为TG,T*G表示相空间,定义系统的Lagrange 函数L:TG→R.对 ∀g∈G令标架算子Ψg为左不变群TLg,因此,定义活动标架下的Lagrange函数l:G×g→R

给定非保守力

其中fτ,fext,fcon分别为反馈控制力,外力和约束力

λ为Lagrange 乘子,Φ(g,ξ)=0 为系统的约束方程.由前文定理1,可得出下面的表述是等价的.

(1) 定义在G×(g⊕g*) 上的曲线 (g(t),ξ(t),µ(t)) 是作用泛函

关于独立变分 η ,δ ξ 和 δµ 的临界点,即

(2) 曲线 (g(t),ξ(t),µ(t)) 满足隐式Hamel 方程组的弱形式

曲线 (g,ξ,µ) 满足隐式Hamel 方程组的强形式

证明参见文献[32].

由式(53)可定义加速度函数

设计控制力fτ可获得稳定的时间积分数值算法.

5.2 李群空间上的Hamel 广义 α 方法

给定时间步长 Δt,区间[ 0,T] 等分为N个子区间,为简单起见,令区间 [tn,tn+1] 为 [0,Δt],为构建李群上数值积分算法,定义G×g*上局部坐标的离散形式

其中变量 Δgn表示时间区间 [tn,tn+1] 的平均速度场,同时令

并且满足 µ0=Mξ0,最后根据加速度函数给出相应的离散加速度

式中外力,约束力和反馈控制律离散形式为式(31).

定义辅助未知量 ηs(t) 和 µs(t) 的离散形式

将离散形式(55)～(59)代入方程组(51),求解关于的驻点方程,其中关于的驻点方程可推导出

结合方程(55)可得出

其次,对于无约束系统tn+1时刻离散加速度函数为

若为带约束系统,则

为了得到更稳定的算法,令

综上所述,方程(60)～(65)为一般李群上的Hamel广义 α 方法.

注: Hamel 广义 α 方法可以退化为李群广义 α 方法等经典算法(见表2).令A3=2(3β-1),A4=-3(4β-1),表格中参数αm,αf参见文献[16].

表2 Hamel 广义 α 方法与经典算法比较Table 2 Comparison between Hamel generalized-α method and classical algorithm

5.3 数值算例

本节通过空间双摆模型[18,40]的数值仿真来说明Hamel 广义 α 方法的性能,模型采用的是一阶SE(3)几何精确梁单元.最终使用Matlab2016 b 实现了Hamel-广义 α 方法.

空间双摆模型初始放置于XOY水平面上,如图10 所示,仅在重力作用下运动.两个摆臂的长度分别为0.3 m 和0.5 m,截面均为宽与高0.5 mm 的矩形,物理参数为:ρ=2700 kg/m3,E=2.1×1010Pa,v=0.3,重力加速度为g=9.81 m/s2.我们采用广义α方法作为参考,步长均为5×10-4s,仿真时长1 s,算法参数均为 ρ∞=0.8,并且对两种方法的数值结果做了对比,发现与传统广义 α 方法结果保持一致,如图11 和图12,说明本方法适用于李群空间上的系统.

图10 柔性双摆系统的初始构型Fig.10 Initial configuration of a flexible double pendulum

图11 柔性双摆系统两杆连接处轨迹Fig.11 Trajectory at the joint of a flexible double pendulum

图12 柔性双摆系统第二杆末端轨迹Fig.12 Trajectory at the end of a flexible double pendulum

6 结论

本文基于活动标架法和Hamel 场变分积分子,讨论如何发展新的数值积分算法的构造方法,并提出了一种无条件稳定的Hamel 广义 α 方法.该方法具有如下特点: 无条件稳定,二阶精度和可以控制数值耗散;选择特殊算法参数,可以退化成传统的广义α方法、Newmark 等;由于活动标架法的特性,该方法既可以用于处理一般线性空间上的问题,也能够用于处理李群空间上的问题;通过构造特殊的变分形式,该方法可以演化出更高精度的数值格式.本文通过理论分析,给出了算法满足精度,数值耗散等条件,指出该算法继承了传统数值积分算法的特点,可以通过单个参数来控制数值耗散,可实现高频损耗并同时最小化低频损耗的目的;甚至可以提出更高阶精度的数值格式.算法也可以从线性空间拓展到李群空间.通过和Newmark 方法、广义 α 方法比较的结果表明,本文方法具有理想的精度和稳定性,并且可以实现更高的精度.与李群广义 α 方法相比,本文方法可以实现同等效果,这些数值算例验证了该方法的上述特点和可行性.

目前,上述方法和算法的研究仍需进一步完善,进而提出具无条件稳定性,高阶精度和可以控制数值耗散以及可变时间步长的数值积分方法,同时发展出相应的李群形式.

致谢

感谢北京理工大学宇航学院季奕博士在广义α方法方面提供的帮助.