AN UNCONDITIONALLY STABLE DYNAMICAL INTEGRATION ALGORITHM BASED ON HAMEL’S FORMALISM
-
摘要: 时间积分算法是求解动力学系统的一个核心问题. 动力学方程的时间积分经常会出现数值不稳定现象, 有限元空间离散也通常会造成伪高频振荡, 因而, 发展解决上述问题的数值积分算法具有重要的理论价值. 本文基于Hamel场变分积分子, 通过新的数值积分算法的构造方法, 提出了一种无条件稳定的Hamel广义
$ \alpha $ 方法, 具体内容包括: 构造特殊的变分形式, 利用变分积分子等工具, 建立无条件稳定的数值积分算法; 在相同框架下, 提出更高精度的数值格式; 结合活动标架法的特性, 将算法的一般形式推广到李群空间, 得到Hamel广义$ \alpha $ 方法李群形式; 对算法的收敛性和稳定性等性质进行了讨论, 并通过算例验证了结论. 理论分析的结果表明本文所提出的Hamel广义$ \alpha $ 方法是无条件稳定的, 具有二阶精度并且能够快速过滤掉虚假的高频振荡. 数值算例的结果显示, 本文所提方法具备了传统方法的精度、耗散和稳定性优势, 既适合一般的线性空间, 也适用于李群空间, 同时还可以发展高阶精度算法. 本文发展了构造变分积分子的新模式.-
关键词:
- 广义α方法 /
- Hamel场变分积分子 /
- Hamel广义α方法 /
- 李群 /
- 活动标架法
Abstract: Time integration algorithm is a key issue in solving dynamical system. An unconditionally stable Hamel generalized α method is proposed to solve the instability issue arising in the time integration of dynamic equations and to eliminate the pseudo high order harmonics incurred by the spatial discretization of finite element simultaneously. Therefore, the development of numerical integration algorithm to solve the above-mentioned problems has important theoretical and application value. The algorithm proposed in this paper is developed based on the moving frame method and Hamel’s field variational integrators along with the strategy to construct an unconditionally stable Hamel generalized α method. It is shown that a new numerical formalism with higher accuracy can be derived under the same framework of the unconditional stable algorithm established through a special variational formalism and variational integrators. The above-mentioned formalism can be extended from general linear space to Lie group by utilizing the moving frame method and the Lie group formalism of the Hamel generalized α method has been obtained. Both the convergence and stability of the algorithm are discussed, and some numerical examples are presented to verify the conclusion. It is demonstrated by the theoretical analysis that the Hamel generalized α method proposed in the paper is unconditionally stable, second-order accurate and can quickly filter out pseudo high-frequency harmonics. Both conventional and proposed methods have been applied to numerical examples respectively. Comparisons between results of numerical examples show that the aforementioned advantages of the proposed method in terms of accuracy, dissipation and stability are tested and verified. At the same time, it can be developed that new numerical integration algorithms with even higher order accuracy. The scheme can also be proposed, which is suitable for both general linear space and Lie group space. A new way for constructing variational integrators is also obtained in this paper. -
引 言
时间积分算法被广泛应用于动力学系统的求解中. 动力学方程的时间积分经常会出现数值不稳定现象, 有限元空间离散也通常会造成伪高频振荡, 因此, 针对上述问题的数值积分算法的研究一直备受关注. 本文将基于活动标架法和离散力学, 发展对动力学系统无条件稳定的数值积分格式.
Newmark方法[1]是广泛使用的时间积分方法, 它包括许多格式, 比如中心差分法(CDM)和梯形规则(TR)等, 其中隐式TR具有二阶精度和无条件稳定性. 此外, 对于线性保守系统, TR可以长时间保持能量. 但对于波传播和刚性系统, TR不能过滤掉不需要的高频. 为解决此问题, 研究人员开发了更理想的方法, 如广义
$ \alpha $ 方法[2]、WBZ-$ \alpha $ 方法[3]和HHT-$ \alpha $ 方法[4].广义
$ \alpha $ 方法[2]是求解结构动力学问题的一种隐式方法, 可以通过单个参数来控制数值耗散, 使得在耗散系统高频响应的同时较好地保持系统的低频响应. 广义$ \alpha $ 方法在时间上具有二阶精度且无条件稳定, 已成功用于广泛的工程应用[5-8]. 为应用于多体动力学, Yen等[9]将其扩展到约束系统. 文献[10]分析了广义$ \alpha $ 格式和Lagrange乘子的二阶收敛性. 张雄和王天书[11]和田强[12]比较了上述算法的精度和耗散特性, 研究了力学系统守恒量的变化特性. 此外, 文献[13]将二阶广义$ \alpha $ 方法推广到三阶, 并且可以通过单个参数控制数值耗散, 同时也证明了该方法是无条件稳定的. 季奕等[14-15]提出了非线性动力问题的单步以及两步可控耗散时间积分法, 这些方法是二阶精度和无条件稳定的, 并且不需要初始加速度矢量.为了研究具有李群结构的柔性多体系统, Brüls等[16]提出了李群广义
$ \alpha $ 时间积分算法, 并证明了该方法具有全局二阶精度. Arnold等[17]详细分析了李群广义$ \alpha $ 方法的收敛性. 文献[18]给出了SE(3)群的广义$ \alpha $ 方法. 文献[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方法和广义
$ \alpha $ 方法等一系列时间积分算法最初均是通过经验推导或总结出来的, 后来赋予了各种理论基础, 如泰勒级数和Galerkin弱形式等, 最后利用数值算法的分析手段确定其收敛性和稳定性. 这些方法都不是为了保持能量和动量而设计的, 在一些特殊参数情况下, 它们可以保持动量[37]. 文献[38]从Hamilton方法出发, 指出Newmark族以及相关的积分算法, 在离散力学的Veselov形式下才是变分的[38]. 迄今为止, 在变分积分子的框架下 一直未出现一种能耗散系统高频响应同时又能较好保持系统低频响应的数值积分算法. 因此, 在该领域中, 构造一种无条件稳定的时间积分算法是必要的, 这样既可以拓展变分积分子的理论体系, 同时又可以丰富时间积分算法的理论基础. 在这种背景下, 本文基于离散力学框架和活动标架法, 采用特殊的变分形式[39-40], 提出了一种Hamel广义$ \alpha $ 方法, 并推广到了一般李群空间, 理论分析和数值测试均显示出了本文方法在精度、耗散和稳定性方面的优势. 从而为进一步研究无条件稳定的高阶数值积分算法提供参考依据.1. Hamel形式与Hamel变分积分子
1.1 Hamel场形式
令
$ Q $ 为向量空间$ W $ 上的无穷维光滑流形. Lagrange函数为$ L:TQ \to \mathbb{R} $ , 其中$ TQ $ 是位形空间$ Q $ 的切丛. 系统的动力学由Euler-Lagrange方程组决定. 若选择合适的标架, 可获得系统动力学的更简洁形式——Hamel方程[34]. 令$ U $ 是${\boldsymbol{q}} \in Q$ 的邻域, 定义一族可逆的有界线性算子$ {\Psi _q}:W \to {T_q}Q $ 为活动标架算子,$\forall {\boldsymbol{q}} \in U$ , 同时,$ \Psi _q^{ - 1}:{T_q}Q \to W $ 也是可逆的有界线性算子. 因而, 对于$\forall {\boldsymbol{\xi}} \in W$ 定义向量场$$ \Psi {\boldsymbol{\xi}} : = \bigcup {_{q \in U}} \left( {{\boldsymbol{q}},{\Psi _q}{\boldsymbol{\xi}} } \right) $$ (1) 给定两个向量
${\boldsymbol{\xi}} ,{\boldsymbol{\eta}} \in W$ , 定义一个反对称括号算子$ {\left[ { \cdot , \cdot } \right]_q}:W \times W \to $ W$$ {\left[ {{\boldsymbol{\xi}} ,{\boldsymbol{\eta}} } \right]_q}: = \Psi _q^{ - 1}\left[ {\Psi {{{\boldsymbol{\xi}}}} ,\Psi {\boldsymbol{\eta}} } \right]({\boldsymbol{q}}) $$ (2) 其中
$ \left[ { \cdot , \cdot } \right] $ 是流形上两个向量场的Jacobi-Lie括号.对于任意${\boldsymbol{\xi}} ,{\boldsymbol{\eta}} \in W$ 和${\boldsymbol{\alpha }} \in {W^ * }$ , 括号算子$ {\left[ { \cdot , \cdot } \right]_q} $ 的对偶算子$ \left[ { \cdot , \cdot } \right]_q^ * :W \times W_{}^ * \to W_{}^ * $ 通过下式给出$$ \left\langle {\left[ {{\boldsymbol{\xi}} ,{\boldsymbol{\alpha}} } \right]_q^ * ,{\boldsymbol{\eta}} } \right\rangle : = \left\langle {{\boldsymbol{\alpha}} ,{{\left[ {{\boldsymbol{\xi}} ,{\boldsymbol{\eta}} } \right]}_q}} \right\rangle $$ (3) 给定非保守力
${\boldsymbol{F}} \in \varGamma ({T^ * }Q)$ , 其中$\varGamma ({T^ * }Q)$ 表示$ {T^ * }Q $ 的所有截面, 则虚功记为$$ \left\langle {{\boldsymbol{F}}({\boldsymbol{q}}),{\text{δ }}{\boldsymbol{q}}} \right\rangle = \left\langle {{\boldsymbol{f}}({\boldsymbol{q}}),{\boldsymbol{\eta}} } \right\rangle $$ (4) 其中
${\boldsymbol{f}}({\boldsymbol{q}}) = \Psi _q^*{\boldsymbol{F}}({\boldsymbol{q}})$ 和${\text{δ }}{\boldsymbol{q}} = \Psi _q^{}{\boldsymbol{\eta}}$ , 并且对偶标架算子定义如下$$ \left\langle {\Psi _q^*{\boldsymbol{\alpha}} ,{\boldsymbol{\xi}} } \right\rangle = \left\langle {{\boldsymbol{\alpha}} ,\Psi _q^{}{\boldsymbol{\xi}} } \right\rangle $$ (5) 定理1 Hamilton-Pontryagin原理的Hamel形式[34]
令力学系统位形空间为
$ Q $ ,$ TQ $ 是切丛,$ {T^ * }Q $ 为余切丛. 令${\boldsymbol{q}} \in Q$ ,$\left( {{\boldsymbol{q}},{\boldsymbol{v}}} \right) \in TQ$ 且$\left( {{\boldsymbol{q}},{\boldsymbol{p}}} \right) \in {T^ * }Q$ . 令$L:TQ \to \mathbb{R}$ 为Lagrange函数,$ l $ 是关于活动标架下局部坐标$({\boldsymbol{q}},{\boldsymbol{\xi}} )$ 的Lagrange函数, 则下列陈述等价.(1) 令
$\left( {{\boldsymbol{q}},{\boldsymbol{v}},{\boldsymbol{p}}} \right)$ 为Pontryagin丛$ TQ \oplus {T^ * }Q $ 的局部坐标, 则曲线$\left( {{\boldsymbol{q}}(t),{\boldsymbol{v}}(t),{\boldsymbol{p}}(t)} \right)$ ,$ 0 \leqslant t \leqslant T $ 为其相应的作用积分$$ S\left( {{\boldsymbol{q}},{\boldsymbol{v}},{\boldsymbol{p}}} \right) = \int_0^T {\left[ {L\left( {{\boldsymbol{q}},{\boldsymbol{v}}} \right) + \left\langle {{\boldsymbol{p}},\dot {\boldsymbol{q}} - {\boldsymbol{v}}} \right\rangle } \right]} {\text{ d}}t $$ (6) 关于独立变分
${\text{δ }}{\boldsymbol{q}},{\text{δ }}{\boldsymbol{v}},{\text{δ }}{\boldsymbol{p}}$ 的临界点, 即$$ {\text{δ }}S\left( {{\boldsymbol{q}},{\boldsymbol{v}},{\boldsymbol{p}}} \right) + \int_0^T {\left\langle {{\boldsymbol{F}}({\boldsymbol{q}}),{\text{δ }}{\boldsymbol{q}}} \right\rangle } {\text{ d}}t = 0 $$ (7) 且
${\text{δ }}{\boldsymbol{q}}(0) = {\text{δ }}{\boldsymbol{q}}(T) = 0$ .(2) 曲线
$\left( {{\boldsymbol{q}}(t),{\boldsymbol{v}}(t),{\boldsymbol{p}}(t)} \right)$ 满足隐式Euler-Lagrange方程组$$ {\boldsymbol{p}} = \frac{{\partial L}}{{\partial {\boldsymbol{v}}}},{\text{ }}\dot {\boldsymbol{q}} = {\boldsymbol{v}},{\text{ }}{\mathbf{ }}\dot {\boldsymbol{p}} = \frac{{\partial L}}{{\partial q}} + {\boldsymbol{F}} $$ (8) (3) 曲线
$({\boldsymbol{q}},{\boldsymbol{\xi}} ,{\boldsymbol{\mu}} ) \in U \times (W \oplus {W^ * })$ 是作用泛函$$ s({\boldsymbol{q}},{\boldsymbol{\xi}} ,{\boldsymbol{\mu}} ) = \int_0^T {\left[ {l\left( {{\boldsymbol{q}},{\boldsymbol{\xi}} } \right) + \left\langle {{\boldsymbol{\mu}} ,\Psi _q^{ - 1}\dot {\boldsymbol{q}} - {\boldsymbol{\xi}} } \right\rangle } \right]} {\text{ d}}t $$ (9) 关于独立变分
${\text{δ }}{\boldsymbol{q}} = {\Psi _q}{\boldsymbol{\eta}}$ ,${\text{δ }}{\boldsymbol{\xi}}$ 和${\text{δ }}{\boldsymbol{\mu}}$ 的临界点$$ {\text{δ }}s\left( {{\boldsymbol{q}},{\boldsymbol{\xi}} ,{\boldsymbol{\mu}} } \right) + \int_0^T {\left\langle {{\boldsymbol{f}}({\boldsymbol{q}}),{\boldsymbol{\eta}} } \right\rangle } {\text{ d}}t = 0 $$ (10) 其中
${\boldsymbol{ \eta}} \in \Psi _q^{ - 1}({T_q}Q)$ 且${\boldsymbol{\eta}} (0) = {\boldsymbol{\eta}} (T) = 0$ .(4) 曲线
$({\boldsymbol{q}},{\boldsymbol{\xi}} ,{\boldsymbol{\mu}} )$ 满足隐式Hamel方程组的弱形式math-ident="5em" $ \left.{ \begin{aligned} &\int_0^T {\left\langle {\Psi _q^*\frac{{{\text{δ }}l}}{{{\text{δ }}{\boldsymbol{q}}}} - \dot {\boldsymbol{\mu}} + \left[ {{\boldsymbol{\xi}} ,{\boldsymbol{\mu}} } \right]_q^* + {\boldsymbol{f}},{\boldsymbol{\eta }}} \right\rangle } {\text{ d}}t = 0 \\ & {\boldsymbol{\mu}} = \frac{{{\text{δ }}l}}{{{\text{δ }}{\boldsymbol{\xi}} }},\;\;{\text{ }}{\boldsymbol{\xi}} = \Psi _q^{ - 1}\dot {\boldsymbol{q}}\end{aligned} } \right\} $ (11) 对于任意点
${\boldsymbol{q}} \in Q$ ,曲线$({\boldsymbol{q}},{\boldsymbol{\xi}} ,{\boldsymbol{\mu }})$ 满足隐式Hamel方程组的强形式$$ \Psi _q^*\frac{{{\text{δ }}l}}{{{\text{δ }}{\boldsymbol{q}}}} - \dot {\boldsymbol{\mu}} + \left[ {{\boldsymbol{\xi}} ,{\boldsymbol{\mu}} } \right]_q^* = {\mathbf{0}},\;\;{\text{ }}{\boldsymbol{\mu}} = \frac{{{\text{δ }}l}}{{{\text{δ }}{\boldsymbol{\xi}} }} $$ (12) 且满足约束
${\boldsymbol{\xi}} = \Psi _q^{ - 1}\dot {\boldsymbol{q}}$ .1.2 Hamel变分积分子
假设
$ Q $ 是向量空间, 连续曲线${\boldsymbol{q}}:\left[ {0,T} \right] \to Q$ 用离散序列$\left\{ {{{\boldsymbol{q}}_k}} \right\}_{k = 0}^N$ 近似表示. 定义$$ {{\boldsymbol{q}}_{k + \alpha }} = \left( {q_{k + \alpha }^1,q _{k + \alpha }^2,\cdots,q_{k + \alpha }^n} \right) = \left( {1 - \alpha } \right){{\boldsymbol{q}}_k} + \alpha {{\boldsymbol{q}}_{k + 1}} $$ (13) 式中参数
$ \alpha \in \left[ {0,1} \right] $ ,$ \Delta t $ 是时间步长, 令$$ {{\boldsymbol{\xi}} _{k + \alpha }} = \left( {\xi _{k + \alpha }^1,\xi _{k + \alpha }^2,\cdots,\xi _{k + \alpha }^n} \right) $$ (14) 为关于活动标架
$ \Psi _{{q_{k + \alpha }}}^{} $ 的速度分量.通过连续Lagrange函数, 定义其离散形式
$$ {l^d} = {l^d}\left( {{{\boldsymbol{q}}_{k + \alpha }},{{\boldsymbol{\xi}} _{k + \alpha }}} \right): = \Delta tl\left( {{{\boldsymbol{q}}_{k + \alpha }},{{\boldsymbol{\xi}} _{k + \alpha }}} \right) $$ (15) 相应离散作用和为
$$ {s^d} = \sum\limits_{k = 0}^{N - 1} {{l^d}\left( {{{\boldsymbol{q}}_{k + {\boldsymbol{\alpha}} }},{{\boldsymbol{\xi}} _{k + \alpha }}} \right)} $$ (16) 同理, 变量
${\boldsymbol{\eta }}$ 的离散形式${{\boldsymbol{\eta}} _{k + \alpha }}$ 可以线性展开$$ {{\boldsymbol{\eta}} _{k + \alpha }} = \left( {1 - \alpha } \right){{\boldsymbol{\eta}} _k} + \alpha {{\boldsymbol{\eta}} _{k + 1}} $$ (17) 通过离散变分原理, 可得隐式Hamel方程组强形式(12)的离散形式—–Hamel变分积分子
math-ident="5em" $ \left.\begin{aligned} & \frac{1}{{\Delta t}}\left( {{{\boldsymbol{\mu}} _{k + 1}} - {{\boldsymbol{\mu}} _k}} \right) = \Psi _{{q_{k + \alpha }}}^*{D_1}{l^d}\left( {{{\boldsymbol{q}}_{k + \alpha }},{{\boldsymbol{\xi}} _{k + \alpha }}} \right) + \\ &\qquad \left[ {{{\boldsymbol{\xi}} _{k + \alpha }},{{\boldsymbol{\mu}} _{k + \alpha }}} \right]_{{q_{k + \alpha }}}^* \\ & {{\boldsymbol{\mu}} _{k + \alpha }} = {D_2}{l^d}\left( {{{\boldsymbol{q}}_{k + \alpha }},{{\boldsymbol{\xi}} _{k + \alpha }}} \right) \\ & \frac{1}{{\Delta t}}\left( {{{\boldsymbol{q}}_{k + 1}} - {{\boldsymbol{q}}_k}} \right) = {\Psi _{{q_{k + \alpha }}}}{{\boldsymbol{\xi}} _{k + \alpha }}\end{aligned}\right\} $ (18) 其中
$ {D_i}l $ 为Lagrange函数对第i个参数的偏导数.2. 向量空间上的Hamel广义
$ \alpha $ 方法本节考虑向量空间中力学系统, 基于活动标架法, 构造力学系统的Hamel广义
$ \alpha $ 方法.2.1 Hamel方程
考虑一个受控系统, 其位形空间为向量空间
$Q$ , 状态空间为$ TQ $ ,$ {T^*}Q $ 表示相空间. 对$\forall {\boldsymbol{q}} \in Q$ 取标架算子$ {\Psi _q} $ 为恒算子, 定义系统的Lagrange函数$L: $ TQ$\to \mathbb{R} $ $$ L\left( {{\boldsymbol{q}},{\boldsymbol{v}}} \right) = \mathcal{T}\left( {\boldsymbol{v}} \right) - \mathcal{V}\left( {\boldsymbol{q}} \right) $$ (19) 其中,
$\mathcal{T}\left( {\boldsymbol{v}} \right) = \dfrac{1}{2}\left\langle {{\boldsymbol{Mv}},{\boldsymbol{v}}} \right\rangle$ 表示系统的动能,${\boldsymbol{v }}= \dot {\boldsymbol{q}}$ 为速度,${\boldsymbol{M}}$ 为质量矩阵,$\mathcal{V}\left( {\boldsymbol{q}} \right)$ 表示势能. 给定力$$ {\boldsymbol{f}} = {{\boldsymbol{f}}_\tau } + {{\boldsymbol{f}}_{{\text{ext}}}} $$ (20) 其中
${{\boldsymbol{f}}_\tau }$ ,${{\boldsymbol{f}}_{{\text{ext}}}}$ 分别为反馈控制力和外力.由前文定理1, 可知下面的表述是等价的.
(1) 定义在
$ TQ \oplus {T^ * }Q $ 上的曲线$\left( {{\boldsymbol{q}}(t),{\boldsymbol{v}}(t),{\boldsymbol{p}}(t)} \right)$ 是作用泛函$$ S\left( {{\boldsymbol{q}},{\boldsymbol{v}},{\boldsymbol{p}}} \right) = \int_0^T {\left[ {L\left( {{\boldsymbol{q}},{\boldsymbol{v}}} \right) + \left\langle {{\boldsymbol{p}},\dot {\boldsymbol{q}} - {\boldsymbol{v}}} \right\rangle } \right]} {\text{ d}}t $$ (21) 关于独立变分
${\text{δ }}{\boldsymbol{q}},{\text{δ }}{\boldsymbol{v}},{\text{δ }}{\boldsymbol{p}}$ 的临界点, 即$$ {\text{δ }}S\left( {{\boldsymbol{q}},{\boldsymbol{v}},{\boldsymbol{p}}} \right) + \int_0^T {\left\langle {{\boldsymbol{f}}({\boldsymbol{q}}),{\text{δ }}{\boldsymbol{q}}} \right\rangle } {\text{ d}}t = 0 $$ (22) 且
${\text{δ }}{\boldsymbol{q}}(0) = {\text{δ }}{\boldsymbol{q}}(T) = 0$ .(2) 曲线
$\left( {{\boldsymbol{q}}(t),{\boldsymbol{v}}(t),{\boldsymbol{p}}(t)} \right)$ 满足隐式Hamel方程组的弱形式math-ident="5em" $ \begin{split} & \int_0^T {\left\langle {\frac{{\partial L}}{{\partial {\boldsymbol{q}}}} - \dot {\boldsymbol{p}} + {\boldsymbol{f}},{\text{δ }}{\boldsymbol{q}}} \right\rangle } {\text{ d}}t+ \\ & \qquad\int_0^T {\left\langle {\frac{{\partial L}}{{\partial {\boldsymbol{v}}}} - {\boldsymbol{p}},{\text{δ }}{\boldsymbol{v}}} \right\rangle } {\text{ d}}t + \int_0^T {\left\langle {\dot {\boldsymbol{q}} - {\boldsymbol{v}},{\text{δ }}{\boldsymbol{p}}} \right\rangle } {\text{ d}}t = 0 \\ \end{split} $ (23) 曲线
$ \left( {{\boldsymbol{q}}(t),{\boldsymbol{v}}(t),{\boldsymbol{p}}(t)} \right) $ 满足隐式Hamel方程组的强形式$$ {\boldsymbol{p}} = \frac{{\partial L}}{{\partial {\boldsymbol{v}}}},{\text{ }}\dot {\boldsymbol{q}} = {\boldsymbol{v}},{\text{ }}{\mathbf{ }}\dot {\boldsymbol{p}} = \frac{{\partial L}}{{\partial {\boldsymbol{q}}}} + {\boldsymbol{f}} $$ (24) 此时Hamel方程组为隐式Euler-Lagrange方程组.
所以由式(24)可定义加速度函数
$$ {\boldsymbol{a}}(t) = {{\boldsymbol{M}}^{ - 1}}\left[ {\frac{{\partial L}}{{\partial {\boldsymbol{q}}}} + {\boldsymbol{f}}(t)} \right] $$ (25) 2.2 Hamel-广义
$ \alpha $ 方法给定时间步长
$ \Delta t $ , 区间$\left[ {0,T} \right]$ 等分为$N$ 个子区间, 为简单起见, 令区间$\left[ {{t_n},{t_{n + 1}}} \right]$ 为$ \left[ {0,\Delta t} \right] $ , 为构建数值积分算法, 定义离散形式$$ {\boldsymbol{q}}(t) = {{\boldsymbol{q}}_n} + \frac{{{{\boldsymbol{q}}_{n + 1}} - {{\boldsymbol{q}}_n}}}{{\Delta t}}t $$ (26) $$ {\boldsymbol{p}}(t) = {{\boldsymbol{p}}_n} + \frac{{{{\boldsymbol{p}}_{n + 1}} - {{\boldsymbol{p}}_n}}}{{\Delta t}}t $$ (27) $$ {\boldsymbol{v}}(t) = {{\boldsymbol{v}}_n} + \frac{{{{\boldsymbol{v}}_{n + 1}} - {{\boldsymbol{v}}_n}}}{{\Delta t}}t $$ (28) 并且满足
${\boldsymbol{M}}{{\boldsymbol{v}}_0} = {{\boldsymbol{p}}_0}$ , 最后根据加速度函数定义相应的离散加速度math-ident="5em" $ \left.\begin{aligned} & {\boldsymbol{a}}(t) = {{\boldsymbol{a}}_n} + \frac{{{{\boldsymbol{a}}_{n + 1}} - {{\boldsymbol{a}}_n}}}{{\Delta t}}t \\ & {{\boldsymbol{a}}_n} = {{\dot {\boldsymbol{v}}}_n} + {{\boldsymbol{M}}^{ - 1}}{\left[ {{{\boldsymbol{f}}_\tau }} \right]_n} \\ &{{\dot {\boldsymbol{v}}}_n} = {{\boldsymbol{M}}^{ - 1}}{\left( {\frac{{\partial L}}{{\partial {\boldsymbol{q}}}} + {{\boldsymbol{f}}_{{\text{ext}}}}} \right)_n} \end{aligned} \right\} $ (29) 此外, 令外力和反馈控制律离散形式为
math-ident="5em" $ \left.\begin{aligned} & {{\boldsymbol{f}}_{{\text{ext}}}}(t) = {\left[ {{{\boldsymbol{f}}_{{\text{ext}}}}} \right]_n} + \frac{{{{\left[ {{{\boldsymbol{f}}_{{\text{ext}}}}} \right]}_{n + 1}} - {{\left[ {{{\boldsymbol{f}}_{{\text{ext}}}}} \right]}_n}}}{{\Delta t}}t \\ & {{\boldsymbol{f}}_\tau }(t) = {\left[ {{{\boldsymbol{f}}_\tau }} \right]_n} + \frac{{{{\left[ {{{\boldsymbol{f}}_\tau }} \right]}_{n + 1}} - {{\left[ {{{\boldsymbol{f}}_\tau }} \right]}_n}}}{{\Delta t}}t \end{aligned}\right\}$ (30) 因此, 可知力
${\boldsymbol{f}}$ 所对应的离散形式为$$ {\boldsymbol{f}}(t) = {{\boldsymbol{f}}_n} + \frac{{{{\boldsymbol{f}}_{n + 1}} - {{\boldsymbol{f}}_n}}}{{\Delta t}}t $$ (31) 通过控制律
$ {{\boldsymbol{f}}_\tau } $ 可获得无条件稳定的数值积分算法.定义辅助未知量
$ {\text{δ }}{\boldsymbol{q}},{\text{δ }}{\boldsymbol{p}} $ 和$ {\text{δ }}{\boldsymbol{v}} $ 的离散形式math-ident="5em" $ \left.\begin{aligned} & {\text{δ }}{\boldsymbol{q}}(t) = \left( {{A_1} + {A_2}\frac{t}{{\Delta t}}} \right){\text{δ }}\hat {\boldsymbol{q}} + \left( {A_3}\Delta t +\right. \\ &\qquad\qquad\left.{A_4}t \right){{\boldsymbol{M}}^{ - 1}}{\text{δ }}\hat {\boldsymbol{v}} \\ & {\text{δ }}{\boldsymbol{p}}(t) = {\text{δ }}\hat {\boldsymbol{p}} \\ & {\text{δ }}{\boldsymbol{v}}(t) = {\text{δ }}\hat {\boldsymbol{v}} \end{aligned}\right\} $ (32) 其中,
$ {A_1},{A_2},{A_3},{A_4} $ 为四个参数. 将离散形式(26) ~ 式(32)带入隐式Hamel方程组的弱形式(23), 求解关于${\text{δ }}\hat {\boldsymbol{q}},{\text{δ }}\hat {\boldsymbol{p}}$ 和${\text{δ }}\hat {\boldsymbol{v}}$ 的驻点方程, 可得到比广义$ \alpha $ 方法更一般的数值积分算法math-ident="5em" $ \left.\begin{aligned} & \sum\limits_{k = 0}^1 {{\alpha _k}{{\boldsymbol{q}}_{n + k}}} + \sum\limits_{k = 0}^1 {\Delta t{\mu _k}{{\boldsymbol{v}}_{n + k}}} = \sum\limits_{k = 0}^1 {\Delta {t^2}{\beta _k}{{\boldsymbol{a}}_{n + k}}} \\ & \sum\limits_{k = 0}^1 {{\alpha _k}{{\boldsymbol{v}}_{n + k}}} = \sum\limits_{k = 0}^1 {\Delta t{\gamma _k}{{\boldsymbol{a}}_{n + k}}} \end{aligned}\right\} $ (33) 式中参数分别为
math-ident="5em" $ \left\{{\begin{aligned} &{\alpha _0} = - 1,\;{\rm{ }}{\alpha _1} = 1\\ &{\mu _0} = - \left( {\frac{1}{2} - {A_3} - \frac{1}{2}{A_4}} \right),\;{\rm{ }}{\mu _1} = - \left( {{A_3} + \frac{1}{2}{A_4} + \frac{1}{2}} \right)\\ &{\beta _0} = - \left( {\frac{1}{2}{A_3} + \frac{1}{6}{A_4}} \right),\;{\rm{ }}{\beta _1} = - \left( {\frac{1}{2}{A_3} + \frac{1}{3}{A_4}} \right)\\ &{\gamma _0} = 1 - A,\;{\gamma _1} = A\\ &A = {\left( {{A_1} + \frac{{{A_2}}}{2}} \right)^{ - 1}}\left( {\frac{{{A_1}}}{2} + \frac{{{A_2}}}{3}} \right) \end{aligned}}\right. $ 其次,
$ {t_{n + 1}} $ 时刻离散加速度函数为$$ \left.\begin{aligned} &{{\boldsymbol{a}}_{n + 1}} = {\dot {\boldsymbol{v}}_{n + 1}} + {{\boldsymbol{M}}^{ - 1}}{\left[ {{{\boldsymbol{f}}_\tau }} \right]_{n + 1}}\\ &{\text{ }}{\dot {\boldsymbol{v}}_{n + 1}} = {{\boldsymbol{M}}^{ - 1}}{\left( {\frac{{\partial L}}{{\partial {\boldsymbol{q}}}} + {{\boldsymbol{f}}_{{\text{ext}}}}} \right)_{n + 1}}\end{aligned}\right\} $$ (34) 为了得到更稳定的算法, 令
$$ {\left[ {{{\boldsymbol{f}}_\tau }} \right]_{n + 1}} = {B_1}\left( {{{\boldsymbol{a}}_n} - {{\dot {\boldsymbol{v}}}_n}} \right) + {B_2}\left( {{{\dot {\boldsymbol{v}}}_{n + 1}} - {{\dot {\boldsymbol{v}}}_n}} \right) $$ (35) 综上, 方程(33) ~ (35)称之为Hamel广义
$ \alpha $ 方法.注1. Hamel广义
$ \alpha $ 方法可以退化为一些经典算法, 如Newmark、HHT-$ \alpha $ 和广义$ \alpha $ 方法等. 如表1所示, 令$ {A_3} = 2\left( {3\beta - 1} \right),{A_4} = - 3\left( {4\beta - 1} \right) $ , 表格中参数$ \alpha ,{\alpha _m},{\alpha _f} $ 的含义参见文献[1-2, 4].表 1 Hamel广义$ \alpha $ 方法与经典算法比较Table 1. Comparison between Hamel generalized-$ \alpha $ method and classical algorithm$ {B_1} $ $ {B_2} $ Classical algorithm 1 0 0 Newmark 2 0 $ - \alpha $ HHT-$ \alpha $ 3 $\dfrac{ { {\alpha _m} } }{ { {\alpha _m} - 1} }$ $\dfrac{ { {\alpha _m} - {\alpha _f} } }{ {1 - {\alpha _m} } }$ generalized-$ \alpha $
method注2. 对离散形式(26)进行改进如下
$$ {\boldsymbol{q}}(t) = {{\boldsymbol{q}}_n} + \frac{{{{\boldsymbol{q}}_{n + 1}} - {{\boldsymbol{q}}_n}}}{{\Delta t}}t + \left( {\frac{1}{2}{A_3} + {A_4}} \right)\frac{{{{\boldsymbol{v}}_{n + 1}} - {{\boldsymbol{v}}_n}}}{{\Delta t}}{t^2} $$ (36) 可得到更高阶数值积分算法(33), 其中参数
${{\boldsymbol{\mu}} _0}, {{\boldsymbol{\mu}} _1},{{\boldsymbol{\gamma}} _0},{{\boldsymbol{\gamma}} _1}$ 变为math-ident="5em" $ \left\{\begin{aligned} &{{\boldsymbol{\mu}} _0} = - \left( {\frac{1}{2} - {A_3} - \frac{1}{2}{A_4}} \right){\boldsymbol{I}} + \frac{1}{{12}}{{\boldsymbol{M}}^{ - 1}}{\boldsymbol{K}}\left( {\frac{1}{2}{A_3} + {A_4}} \right)\Delta {t^2}\\ &{{\boldsymbol{\mu}} _1} = - \left( {{A_3} + \frac{1}{2}{A_4} + \frac{1}{2}} \right){\boldsymbol{I}} - \frac{1}{{12}}{{\boldsymbol{M}}^{ - 1}}{\boldsymbol{K}}\left( {\frac{1}{2}{A_3} + {A_4}} \right)\Delta {t^2}\\ &{{\boldsymbol{\gamma}} _0} = \left( {1 - A} \right){\left[ {{\boldsymbol{I}} - \frac{1}{6}{{\boldsymbol{M}}^{ - 1}}{\boldsymbol{K}}\left( {\frac{1}{2}{A_3} + {A_4}} \right)\Delta {t^2}} \right]^{ - 1}}\\ &{{\boldsymbol{\gamma}} _1} = A{\left[ {{\boldsymbol{I}} - \frac{1}{6}{{\boldsymbol{M}}^{ - 1}}{\boldsymbol{K}}\left( {\frac{1}{2}{A_3} + {A_4}} \right)\Delta {t^2}} \right]^{ - 1}} \end{aligned}\right. $ 其余参数不变, 其中
${\boldsymbol{I}}$ 为单位阵,${\boldsymbol{K}} = {{{\partial ^2}\mathcal{V}} \mathord{\left/ {\vphantom {{{\partial ^2}\mathcal{V}} ({\partial {q^i}\partial {q^j}})}} \right. } ({\partial {q^i}\partial {q^j}})}$ 为弹性矩阵. 选取特定的参数, 该数值格式会得到更高的收敛速率,具体数值算例及其收敛速率见3.3节.3. 算法性质
3.1 算法分析
考虑一个力学系统
$$ L({\boldsymbol{x}},\dot {\boldsymbol{x}}) = \frac{1}{2}{\boldsymbol{M}}{\dot {\boldsymbol{x}}^{{2}}} - \frac{1}{2}{\boldsymbol{K}}{{\boldsymbol{x}}^{{2}}} $$ (37) 相应的Euler-Lagrange方程组为
$$ \ddot {\boldsymbol{x }}+ {{\boldsymbol{\omega}} ^{{2}}}{\boldsymbol{x}} = 0,\;{\boldsymbol{\omega}} = \sqrt {{{\boldsymbol{M}}^{ {{- 1}}}}{\boldsymbol{K}}} $$ (38) Hamel广义
$ \alpha $ 方法应用于上述方程, 令$ f = 0 $ , 则离散方程可写成矩阵形式$$ {{\boldsymbol{X}}_{n + 1}} = {{\bar {\boldsymbol{A}}}}{{\boldsymbol{X}}_n},{\text{ }}n \in \left\{ {0,1,\cdots,N - 1} \right\} $$ (39) 其中
${{\boldsymbol{X}}_n} = \left[ {{{\boldsymbol{q}}_n},\Delta t{{\boldsymbol{v}}_n},\Delta {t^2}{{\boldsymbol{a}}_n}} \right]$ ,${{\bar{\boldsymbol{ A}}}}$ 为辅助矩阵. 由二阶精度条件[41]可知, 满足$$ \left( {\frac{1}{2} - A} \right)\left( {1 - {B_1}} \right) = {B_2} $$ (40) 故本文所提方法为二阶精度算法.
算法的稳定性和数值耗散取决于其辅助矩阵
${{\bar{\boldsymbol{ A}}}}$ 的特征值. 谱半径是数值耗散的一种度量, 谱半径越小, 数值耗散越大. 算法的谱半径由矩阵${{\bar{\boldsymbol{ A}}}}$ 的特征值决定, 定义为$$ \rho = \max \left( {\left| {{\lambda _1}} \right|,\left| {{\lambda _2}} \right|,\left| {{\lambda _3}} \right|} \right) $$ 其中
$ {\lambda _i} $ 为矩阵${{\bar{\boldsymbol{ A}}}}$ 的第i个特征值. 若$ \rho \leqslant 1 $ , 则对于线性问题算法是无条件稳定的.理想算法应该是严格控制高频, 但是低频区域不引起过度耗散, 在低频域的谱半径值接近于1. 随着
$ \varOmega $ 的增加, 谱半径值平滑地减小. 若辅助矩阵的主根形式为$$ {\lambda _{1,2}}\left( \varOmega \right) = a\left( \varOmega \right) + {\rm{i}}b\left( \varOmega \right) $$ 其中
$ a $ 和$ b $ 是实数, 只要$$ \mathop {\lim }\limits_{\varOmega \to \infty } b\left( \varOmega \right) = 0 $$ (41) 成立, 则高频耗散最大化. 由(41)可知
math-ident="5em" $ \left.\begin{aligned} & 12A\left( {A - 1} \right) + 8{A_5} + 3 = 0 \\ & {A_5} = {6A{A_3} + 3A{A_4} - 3{A_3} - {A_4}} \end{aligned}\right\} $ (42) 满足这个条件, Hamel广义
$ \alpha $ 方法便可实现高频耗散最大化. 此外, 当$ A = 0.5 $ 时, 可得$ {A_4} = 0 $ , 无法确定$ {A_3} $ 的值, 故而设定$ {A_3} = - 0.5 $ .令
$ {\rho _\infty } $ 表示高频极限处的谱半径, 变量$\lambda _{1,2}^\infty = \mathop {\lim }\limits_{\varOmega \to \infty } {\lambda _{1,2}}\left( \varOmega \right)$ 可表示为$$ \lambda _{1,2}^\infty = \frac{{2{A_5} + 3}}{{2(3A + {A_5} - 3)}},\;\;\lambda _3^\infty = \frac{{{B_1} + {B_2}}}{{1 + {B_2}}} $$ (43) 若
$ \lambda _{1,2}^\infty = \lambda _3^\infty $ , 则低频耗散最小, 那么可得如下关系$$ \frac{{{B_1} + {B_2}}}{{1 + {B_2}}} = \frac{{{{B_2} - 1 + {B_1}} }}{{ {{B_2} + 1 - {B_1}} }} $$ (44) 为方便描述方程(44), 借助
$ {\rho _\infty } $ , 定义$ {B_1} $ 和$ {B_2} $ $$ {B_1} = \frac{{2{\rho _\infty } - 1}}{{{\rho _\infty } - 2}},\;\;{B_2} = \frac{{{\rho _\infty } - 1}}{{2 - {\rho _\infty }}} $$ (45) 综上所述, 同时满足方程式(40)、式(42)和式(45), 则Hamel广义
$ \alpha $ 方法是无条件稳定、二阶精度的算法, 并且可实现高频耗散最大, 低频耗散最小.3.2 算法稳定性
谱分析是分析时间积分算法稳定性与精度的常用方法. 下文通过力学系统(37)来探讨算法性质. 令式(38)中
$\omega = 2\text{π}$ . 若$ 0 \leqslant \rho \leqslant 1 $ , 则Hamel广义$ \alpha $ 方法无条件稳定; 在此基础上, 若$ \rho \lt 1 $ , 该方法是数值耗散的, 当$ \rho = 1 $ 时, 该方法是非耗散的. 定义数值阻尼比$ \bar \xi $ 和周期延长率$ \left( {\bar T - T} \right)/T $ 为$$ \bar \xi = - \frac{{\ln \left( {{a^2} + {b^2}} \right)}}{{2\bar \varOmega }},{\text{ }}\frac{{\bar T - T}}{T} = \frac{{\bar \varOmega }}{\varOmega } - 1 $$ 其中
$\bar \varOmega = \arctan \left( {b/a} \right)$ . 它们用于计算时间积分算法在低频区间内的振幅和相位精度.图1展示了Hamel广义
$ \alpha $ 方法的谱半径, 发现是无条件稳定的, 并且在高频区域内可通过$ {\rho _\infty } $ 精确控制数值耗散程度.在
$ 0 \leqslant {\rho _\infty } \leqslant 1 $ 情况下, 当$\varOmega \leqslant 0.1$ 时, 谱半径接近于1, 保持了基本的低频模式, 随后快速降至$ {\rho _\infty } $ , 过滤掉不需要的高频信息. 图2和图3分别给出了该方法的数值阻尼比和周期延长率, 可以发现振幅和相位精度随着$ {\rho _\infty } $ 的增加而提高.3.3 算法收敛性
本节通过数值算例测试算法的收敛速度. 收敛速度是特定时间内随着时间步长递减, 绝对误差的下降率. 令系统方程(38)式中
$ \omega =\text{π} $ . 初始条件为$$ x(0) = 1,{\text{ }}\dot x(0) = 1 $$ 该系统的精确解为
$$ x = \cos \left( {\text{π} x} \right) + \frac{{\sin \left( {\text{π} x} \right)}}{\text{π} } $$ 通过Hamel-广义
$ \alpha $ 方法t=1时刻的绝对误差, 如图4可以发现本文方法关于位移是二阶精度, 而且随着$ {\rho _\infty } $ 增大, 算法准确度也在提升.若
$ A = 0.5, {\text{ }}{A_3} = - 1, {\text{ }}{A_4} = 1, {\text{ }}{B_1} = {B_2} = 0 $ , 通过Hamel广义$ \alpha $ 方法t = 1时刻绝对误差, 如图5, 发现注2的数值积分格式关于位移接近于四阶精度.![]()
图 4 Hamel广义$ \alpha $ 方法t = 1时刻的误差Figure 4. Errors of Hamel generalized-$ \alpha $ method at time t = 1![]()
图 5 Hamel广义$ \alpha $ 方法t = 1时刻的误差Figure 5. Errors of Hamel generalized-$ \alpha $ method at time t = 14. 数值算例
本节通过数值算例, 在
$ {\rho _\infty } = 0.8 $ 的情况下将Hamel广义$ \alpha $ 方法与广义$ \alpha $ 方法, Newmark方法($ \beta = 0.25,{\text{ }}\gamma = 0.5 $ )的结果做对比, 来验证其性质.考虑一个两自由度系统, 其Lagrange函数为
$$ L({\boldsymbol{x}},\dot {\boldsymbol{x}}) = \frac{1}{2}\left\langle {\dot {\boldsymbol{x}},{\boldsymbol{M}}\dot {\boldsymbol{x}}} \right\rangle - \frac{1}{2}\left\langle {{\boldsymbol{x}},{\boldsymbol{Kx}}} \right\rangle $$ 系数矩阵为
$$ {\boldsymbol{M}} = \left[ {\begin{array}{*{20}{c}} {{m_1}}&0 \\ 0&{{m_2}} \end{array}} \right],{\text{ }}{\boldsymbol{K}} = \left[ {\begin{array}{*{20}{c}} {{k_1} + {k_2}}&{ - {k_2}} \\ { - {k_2}}&{{k_2}} \end{array}} \right] $$ 其中
$ {m_1} = {m_2} = 1 $ 和$ {k_1} = {10^4},\;{k_2} = 1 $ .相应Euler-Lagrange方程组为
$$ {\boldsymbol{M}}\ddot {\boldsymbol{x}} + {\boldsymbol{Kx}} = {\mathbf{0}} $$ 采用如下初值
$$ {{\boldsymbol{x}}_0} = \left[ {\begin{array}{*{20}{c}} 1 \\ {10} \end{array}} \right],\;\;\;{\dot {\boldsymbol{x}}_0} = \left[ {\begin{array}{*{20}{c}} 0 \\ 0 \end{array}} \right],\;\;\;{\boldsymbol{x}} = \left[ {\begin{array}{*{20}{c}} {{x_1}} \\ {{x_2}} \end{array}} \right] $$ 图6和图7分别给出了系统1和系统2三种方法的运行结果, 通过位移、速度和加速度曲线, 可以看到Hamel广义
$ \alpha $ 方法和广义$ \alpha $ 方法在稳定性方面要比Newmark方法更具优势. 图8和图9分别给出了系统1和系统2 Hamel广义$ \alpha $ 方法和广义$ \alpha $ 方法两种方法的数值误差, 通过与经典算法对比发现, 本文提出的方法具有良好的稳定性和精度性质.![]()
图 6 系统1各种方法数值结果对比图Figure 6. Comparison of numerical results of various methods for system 1![]()
图 7 系统2各种方法数值结果对比图Figure 7. Comparison of numerical results of various methods for system 2![]()
图 8 系统1 Hamel广义$ \alpha $ 方法与广义$ \alpha $ 方法数值结果误差Figure 8. Errors of numerical results between Hamel generalized-$ \alpha $ method and generalized-$ \alpha $ method for system 1![]()
图 9 系统2 Hamel广义$ \alpha $ 方法与广义$ \alpha $ 方法数值结果误差Figure 9. Errors of numerical results between Hamel generalized-$ \alpha $ method and generalized-$ \alpha $ method for system 25. 李群空间上的Hamel-广义
$ \alpha $ 方法本节考虑李群空间中力学系统的Hamel形式,并提出相应的运动方程—–Hamel方程, 在此基础上, 构造李群空间上力学系统的时间积分算法.
5.1 Hamel方程
考虑李群空间上的力学系统, 其位形流形为李群
$G$ , 状态空间为$ TG $ ,$ {T^ * }G $ 表示相空间, 定义系统的Lagrange函数$ L:TG \to \mathbb{R} $ . 对$\forall {\boldsymbol{g}} \in G$ 令标架算子$ {\Psi _g} $ 为左不变群$ T{L_g} $ , 因此, 定义活动标架下的Lagrange函数$ l:G \times g \to \mathbb{R} $ $$ l\left( {{\boldsymbol{g}},{\boldsymbol{\xi }}} \right) = \mathcal{T}\left( {\boldsymbol{\xi}} \right) - \mathcal{V}\left({\boldsymbol{ g}} \right) $$ (46) 其中
$\mathcal{T}\left( {\boldsymbol{\xi}} \right) = \dfrac{1}{2}\left\langle {{\boldsymbol{M\xi}} ,{\boldsymbol{\xi}} } \right\rangle$ 表示系统的动能,${\boldsymbol{\xi}} = {{\boldsymbol{g}}^{ - 1}}\dot {\boldsymbol{g}}$ 表示广义速度在活动标架下的速度分量,${\boldsymbol{ M}}$ 为质量矩阵,$\mathcal{V}\left( {\boldsymbol{g }}\right)$ 表示系统的势能.给定非保守力
$$ {\boldsymbol{f}} = {{\boldsymbol{f}}_\tau } + {{\boldsymbol{f}}_{{\text{ext}}}} + {{\boldsymbol{f}}_{{\text{con}}}} $$ (47) 其中
${{\boldsymbol{f}}_\tau }$ ,${{\boldsymbol{f}}_{{\text{ext}}}}$ ,${{\boldsymbol{f}}_{{\text{con}}}}$ 分别为反馈控制力, 外力和约束力$$ {{\boldsymbol{f}}_{{\text{con}}}} = \lambda {\boldsymbol{\varPhi}} ({\boldsymbol{g}},{\boldsymbol{\xi}} ) $$ (48) $ \lambda $ 为Lagrange乘子,${\boldsymbol{\varPhi}} ({\boldsymbol{g}},{\boldsymbol{\xi}} ) = {\mathbf{0}}$ 为系统的约束方程. 由前文定理1, 可得出下面的表述是等价的.(1) 定义在
$G \times ({g} \oplus {{g}^ * })$ 上的曲线$ \left( {{\boldsymbol{g}}(t),{\boldsymbol{\xi}} (t),{\boldsymbol{\mu}} (t)} \right) $ 是作用泛函$$ s\left( {{\boldsymbol{g}},{\boldsymbol{\xi}} ,{\boldsymbol{\mu}} } \right) = \int_0^T {\left[ {l\left( {{\boldsymbol{g}},{\boldsymbol{\xi}} } \right) + \left\langle {{\boldsymbol{\mu }},\Psi _g^{ - 1}\dot {\boldsymbol{g}} - {\boldsymbol{\xi}} } \right\rangle } \right]} {\text{ d}}t $$ (49) 关于独立变分
${\boldsymbol{ \eta}} $ ,$ {\text{δ }}{\boldsymbol{\xi}} $ 和$ {\text{δ }}{\boldsymbol{\mu}} $ 的临界点, 即$$ {\text{δ }}s\left( {{\boldsymbol{g}},{\boldsymbol{\xi}} ,{\boldsymbol{\mu}} } \right) + \int_0^T {\left\langle {{\boldsymbol{f}}({\boldsymbol{g}},{\boldsymbol{\xi}} ),{\boldsymbol{\eta}} } \right\rangle } {\text{ d}}t = 0 $$ (50) 其中
${\boldsymbol{\eta }} = {T_g}{L_{{{{g}}^{ - 1}}}}{\text{δ }}{\boldsymbol{g}}$ 满足边界条件${\boldsymbol{\eta}} (0) = {\boldsymbol{\eta}} (T) = 0$ .(2) 曲线
$\left( {{\boldsymbol{g}}(t),{\boldsymbol{\xi}} (t),{\boldsymbol{\mu}} (t)} \right)$ 满足隐式Hamel方程组的弱形式$$ \begin{split} & \int_0^T {\left\langle {\Psi _g^{ - 1}\frac{{{\text{δ }}l}}{{{\text{δ }}{\boldsymbol{g}}}} - \dot {\boldsymbol{\mu}} + {{\left( {\Psi _g^{ - 1}\dot {\boldsymbol{g}},{\boldsymbol{\mu}} } \right)}^ * } + {\boldsymbol{f}},{\boldsymbol{\eta}} } \right\rangle } {\text{ d}}t+ \\ & \qquad \int_0^T {\left\langle {\frac{{{\text{δ }}l}}{{{\text{δ }}{\boldsymbol{\xi}} }} - {\boldsymbol{\mu}} ,{\text{δ }}\xi } \right\rangle } {\text{ d}}t + \int_0^T {\left\langle {\Psi _g^{ - 1}\dot {\boldsymbol{g}} - \xi ,{\text{δ }}{\boldsymbol{\mu }}} \right\rangle } {\text{ d}}t = 0 \\ \end{split} $$ (51) 其中
${\boldsymbol{\eta}} \in \Psi _g^{ - 1}({T_g}G)$ , 而且${\boldsymbol{\xi}} = {T_g}{L_{{g^{ - 1}}}}\dot {\boldsymbol{g }}$ 的变分为$$ {\text{δ }}{\boldsymbol{\xi}} = \dot {\boldsymbol{\eta }}+ \left[ {{\boldsymbol{\xi}} ,{\boldsymbol{\eta }}} \right] $$ (52) 曲线
$({\boldsymbol{g}},{\boldsymbol{\xi}} ,{\boldsymbol{\mu}} )$ 满足隐式Hamel方程组的强形式$$ \left.\begin{aligned} & \Psi _g^*\frac{{{\text{δ }}l}}{{{\text{δ }}{\boldsymbol{g}}}} - \dot {\boldsymbol{\mu}} + \left[ {{\boldsymbol{\xi}} ,{\boldsymbol{\mu}} } \right]_g^* + {\boldsymbol{f}} = {\mathbf{0}} \\ & {\boldsymbol{\mu}} = \frac{{{\text{δ }}l}}{{{\text{δ }}{\boldsymbol{\xi }}}},\;{\boldsymbol{\xi}} = \Psi _g^{ - 1}\dot {\boldsymbol{g}} \end{aligned} \right\}$$ (53) 证明参见文献[32].
由式(53)可定义加速度函数
$$ {\boldsymbol{a}}(t) = {{\boldsymbol{M}}^{ - 1}}\left\{ {\Psi _g^*\frac{{{\text{δ }}l}}{{{\text{δ }}{\boldsymbol{g}}}} + \left[ {{\boldsymbol{\xi}} (t),{\boldsymbol{\mu}} (t)} \right]_g^* + {\boldsymbol{f}}(t)} \right\} $$ (54) 设计控制力
${{\boldsymbol{f}}_\tau }$ 可获得稳定的时间积分数值算法.5.2 李群空间上的Hamel广义
$ \alpha $ 方法给定时间步长
$ \Delta t $ , 区间$\left[ {0,T} \right]$ 等分为$N$ 个子区间, 为简单起见, 令区间$\left[ {{t_n},{t_{n + 1}}} \right]$ 为$ \left[ {0,\Delta t} \right] $ , 为构建李群上数值积分算法, 定义$G \times {{g}^*}$ 上局部坐标的离散形式$$ {\boldsymbol{g}}(t) = {{\boldsymbol{g}}_n} \circ \exp \left( {\Delta {{\boldsymbol{g}}_n} \cdot \frac{t}{{\Delta t}}} \right) $$ (55) $$ {\boldsymbol{\mu}} (t) = {{\boldsymbol{\mu}} _n} + \frac{{{{\boldsymbol{\mu}} _{n + 1}} - {{\boldsymbol{\mu}} _n}}}{{\Delta t}}t $$ (56) 其中变量
$\Delta {{\boldsymbol{g}}_n}$ 表示时间区间$\left[ {{t_n},{t_{n + 1}}} \right]$ 的平均速度场, 同时令$$ {\boldsymbol{\xi}} (t) = {{\boldsymbol{\xi}} _n} + \frac{{{{\boldsymbol{\xi}} _{n + 1}} - {{\boldsymbol{\xi}} _n}}}{{\Delta t}}t $$ (57) 并且满足
${{\boldsymbol{\mu}} _0} = {\boldsymbol{M}}{{\boldsymbol{\xi}} _0}$ , 最后根据加速度函数给出相应的离散加速度math-ident="5em" $ \left.\begin{aligned} & {\boldsymbol{a}}(t) = {{\boldsymbol{a}}_n} + \frac{{{{\boldsymbol{a}}_{n + 1}} - {{\boldsymbol{a}}_n}}}{{\Delta t}}t \\ & {{\boldsymbol{a}}_n} = {{\boldsymbol{M}}^{ - 1}}\left( {{{\dot {\boldsymbol{v}}}_n} + {{\left[ {{{\boldsymbol{f}}_\tau }} \right]}_n}} \right) \\ & {{\dot {\boldsymbol{v}}}_n} = {{\boldsymbol{M}}^{ - 1}}{\left( {\Psi _g^*\frac{{{\text{δ }}l}}{{{\text{δ }}{\boldsymbol{g}}}} + {{\boldsymbol{f}}_{{\text{ext}}}} + {{\boldsymbol{f}}_{{\text{con}}}}} \right)_n} + {{\boldsymbol{M}}^{ - 1}}\left[ {{{\boldsymbol{\xi}} _n},{{\boldsymbol{\mu}} _n}} \right]_g^*{\text{ }} \end{aligned}\right\}$ (58) 式中外力, 约束力和反馈控制律离散形式为式(31).
定义辅助未知量
${{\boldsymbol{\eta}} _s}(t)$ 和$ {{\boldsymbol{\mu}} _s}(t) $ 的离散形式math-ident="5em" $ \left.\begin{aligned} & {\boldsymbol{\eta}} (t) = \left( {{A_1} + {A_2}\frac{t}{{\Delta t}}} \right)\hat {\boldsymbol{\eta}} + \left( {{A_3}\Delta t + {A_4}t} \right){{\boldsymbol{M}}^{ - 1}}{\text{δ }}\hat {\boldsymbol{\mu}} \\ & {\text{δ }}{\boldsymbol{\mu}} (t) = {\text{δ }}\hat {\boldsymbol{\mu }} \\ & {\text{δ }}{\boldsymbol{\xi}} (t) = {\text{δ }}\hat{\boldsymbol{ \xi}} \end{aligned}\right\} $ (59) 其中
$\hat {\boldsymbol{\xi }},\hat {\boldsymbol{\mu}}$ 和$\hat {\boldsymbol{\eta}}$ 为三个未知参数.将离散形式(55) ~ (59)带入方程组(51), 求解关于
$ {\text{δ }}\hat{\boldsymbol{ \xi}} ,{\text{δ }}\hat{\boldsymbol{ \mu}} $ 和$ \hat{\boldsymbol{ \eta}} $ 的驻点方程, 其中关于${\text{δ }}\hat {\boldsymbol{\mu}}$ 的驻点方程可推导出$$ \left.\begin{aligned} & \Delta {g_n} = \\ &\qquad \Delta t\left[ {\left( {{A_3} + \frac{1}{2}{A_4} + \frac{1}{2}} \right){\xi _{n + 1}} + \left( {\frac{1}{2} - {A_3} - \frac{1}{2}{A_4}} \right){\xi _n}} \right] - \\ &\qquad \Delta {t^2}\left[ {\left( {\frac{1}{2}{A_3} + \frac{1}{3}{A_4}} \right){a_{n + 1}} + \left( {\frac{1}{2}{A_3} + \frac{1}{6}{A_4}} \right){a_n}} \right] \\ & g(\Delta t) = g(0) \circ \exp \left( {\Delta {g_n}} \right) \end{aligned}\right\}$$ (60) 结合方程(55)可得出
$$ {{\boldsymbol{g}}_{n + 1}} = {{\boldsymbol{g}}_n} \circ \exp \left( {\Delta {{\boldsymbol{g}}_n}} \right) $$ (61) 关于
${\text{δ }}\hat {\boldsymbol{\xi }}$ 和$ \hat {\boldsymbol{\eta}} $ 的驻点方程可得出math-ident="5em" $ \left.\begin{aligned} & {{\boldsymbol{\xi}} _{n + 1}} = {{\boldsymbol{\xi}} _n} + \Delta t\left[ {\left( {1 - A} \right){{\boldsymbol{a}}_n} + A{{\boldsymbol{a}}_{n + 1}}} \right] \\ & A = {{\left( {\frac{{{A_1}}}{2} + \frac{{{A_2}}}{3}} \right)} \mathord{\left/ {\vphantom {{\left( {\frac{{{A_1}}}{2} + \frac{{{A_2}}}{3}} \right)} {\left( {{A_1} + \frac{{{A_2}}}{2}} \right)}}} \right. } {\left( {{A_1} + \frac{{{A_2}}}{2}} \right)}} \end{aligned}\right\}$ (62) 其次, 对于无约束系统
$ {t_{n + 1}} $ 时刻离散加速度函数为math-ident="5em" $ \left.\begin{aligned} & {{\boldsymbol{a}}_{n + 1}} = {{\dot {\boldsymbol{v}}}_{n + 1}} + {{\boldsymbol{M}}^{ - 1}}{\left[ {{{\boldsymbol{f}}_\tau }} \right]_{n + 1}} \\ & {{\dot {\boldsymbol{v}}}_{n + 1}} = {{\boldsymbol{M}}^{ - 1}}\left[ {{{\left( {\Psi _g^*\frac{{{\text{δ }}l}}{{{\text{δ }}{\boldsymbol{g}}}} + {{\boldsymbol{f}}_{{\text{ext}}}} + {{\boldsymbol{f}}_{{\text{con}}}}} \right)}_{n + 1}}} + \right. \\ &\qquad {\left[ {{{\boldsymbol{\xi}} _{n + 1}},{{\boldsymbol{\mu}} _{n + 1}}} \right]_g^*} \Biggr] \end{aligned}\right\}$ (63) 若为带约束系统, 则
$$ {\left[ {{{\boldsymbol{f}}_{{\text{con}}}}} \right]_{n + 1}} = \lambda {\boldsymbol{\varPhi}} ({{\boldsymbol{g}}_{n + 1}},{{\boldsymbol{\xi}} _{n + 1}}),{\text{ }}{\boldsymbol{\varPhi}} ({\boldsymbol{g}},{\boldsymbol{\xi}} ) = {\mathbf{0}} $$ (64) 为了得到更稳定的算法, 令
$$ {\left[ {{{\boldsymbol{f}}_\tau }} \right]_{n + 1}} = {B_1}\left( {{{\boldsymbol{a}}_n} - {{\dot {\boldsymbol{v}}}_n}} \right) + {B_2}\left( {{{\dot {\boldsymbol{v}}}_{n + 1}} - {{\dot {\boldsymbol{v}}}_n}} \right) $$ (65) 综上所述, 方程(60) ~ (65)为一般李群上的Hamel广义
$ \alpha $ 方法.注: Hamel广义
$ \alpha $ 方法可以退化为李群广义$ \alpha $ 方法等经典算法(见表2). 令${A_3} = 2\left( {3\beta - 1} \right), \;\; {A_4} = - 3\left( {4\beta - 1} \right)$ , 表格中参数$ {\alpha _m},{\alpha _f} $ 参见文献[16].表 2 Hamel广义$ \alpha $ 方法与经典算法比较Table 2. Comparison between Hamel generalized-$ \alpha $ method and classical algorithm$ {B_1} $ $ {B_2} $ Classical algorithm 1 $\dfrac{ { {\alpha _m} } }{ { {\alpha _m} - 1} }$ $\dfrac{ { {\alpha _m} - {\alpha _f} } }{ {1 - {\alpha _m} } }$ Lie group
generalized-$ \alpha $
method5.3 数值算例
本节通过空间双摆模型[18,40]的数值仿真来说明Hamel广义
$ \alpha $ 方法的性能, 模型采用的是一阶SE(3)几何精确梁单元. 最终使用Matlab2016 b实现了Hamel-广义$ \alpha $ 方法.空间双摆模型初始放置于 XOY 水平面上, 如图10所示, 仅在重力作用下运动. 两个摆臂的长度分别为0.3 m和0.5 m, 截面均为宽与高0.5 mm的矩形, 物理参数为: ρ = 2700 kg/m3, E = 2.1 × 1010 Pa, ν = 0.3, 重力加速度为g = 9.81 m/s2. 我们采用广义
$ \alpha $ 方法作为参考, 步长均为5 × 10−4 s, 仿真时长1 s, 算法参数均为$ {\rho _\infty } = 0.8 $ , 并且对两种方法的数值结果做了对比, 发现与传统广义$ \alpha $ 方法结果保持一致, 如图11和图12, 说明本方法适用于李群空间上的系统.6. 结 论
本文基于活动标架法和Hamel场变分积分子, 讨论如何发展新的数值积分算法的构造方法, 并提出了一种无条件稳定的Hamel广义
$ \alpha $ 方法. 该方法具有如下特点: 无条件稳定, 二阶精度和可以控制数值耗散; 选择特殊算法参数, 可以退化成传统的广义$ \alpha $ 方法、Newmark等; 由于活动标架法的特性, 该方法既可以用于处理一般线性空间上的问题, 也能够用于处理李群空间上的问题; 通过构造特殊的变分形式, 该方法可以演化出更高精度的数值格式. 本文通过理论分析, 给出了算法满足精度, 数值耗散等条件, 指出该算法继承了传统数值积分算法的特点, 可以通过单个参数来控制数值耗散, 可实现高频损耗并同时最小化低频损耗的目的; 甚至可以提出更高阶精度的数值格式. 算法也可以从线性空间拓展到李群空间. 通过和Newmark方法、广义$ \alpha $ 方法比较的结果表明, 本文方法具有理想的精度和稳定性, 并且可以实现更高的精度. 与李群广义$ \alpha $ 方法相比, 本文方法可以实现同等效果, 这些数值算例验证了该方法的上述特点和可行性.目前, 上述方法和算法的研究仍需进一步完善, 进而提出具无条件稳定性, 高阶精度和可以控制数值耗散以及可变时间步长的数值积分方法, 同时发展出相应的李群形式.
致谢
感谢北京理工大学宇航学院季奕博士在广义
$ \alpha $ 方法方面提供的帮助. -
![]()
图 4 Hamel广义
$ \alpha $ 方法t = 1时刻的误差Figure 4. Errors of Hamel generalized-
$ \alpha $ method at time t = 1![]()
图 5 Hamel广义
$ \alpha $ 方法t = 1时刻的误差Figure 5. Errors of Hamel generalized-
$ \alpha $ method at time t = 1![]()
图 6 系统1各种方法数值结果对比图
Figure 6. Comparison of numerical results of various methods for system 1
![]()
图 7 系统2各种方法数值结果对比图
Figure 7. Comparison of numerical results of various methods for system 2
![]()
图 8 系统1 Hamel广义
$ \alpha $ 方法与广义$ \alpha $ 方法数值结果误差Figure 8. Errors of numerical results between Hamel generalized-
$ \alpha $ method and generalized-$ \alpha $ method for system 1![]()
图 9 系统2 Hamel广义
$ \alpha $ 方法与广义$ \alpha $ 方法数值结果误差Figure 9. Errors of numerical results between Hamel generalized-
$ \alpha $ method and generalized-$ \alpha $ method for system 2表 1 Hamel广义
$ \alpha $ 方法与经典算法比较Table 1 Comparison between Hamel generalized-
$ \alpha $ method and classical algorithm$ {B_1} $ $ {B_2} $ Classical algorithm 1 0 0 Newmark 2 0 $ - \alpha $ HHT-$ \alpha $ 3 $\dfrac{ { {\alpha _m} } }{ { {\alpha _m} - 1} }$ $\dfrac{ { {\alpha _m} - {\alpha _f} } }{ {1 - {\alpha _m} } }$ generalized-$ \alpha $
method
下载: 导出CSV
表 2 Hamel广义
$ \alpha $ 方法与经典算法比较Table 2 Comparison between Hamel generalized-
$ \alpha $ method and classical algorithm$ {B_1} $ $ {B_2} $ Classical algorithm 1 $\dfrac{ { {\alpha _m} } }{ { {\alpha _m} - 1} }$ $\dfrac{ { {\alpha _m} - {\alpha _f} } }{ {1 - {\alpha _m} } }$ Lie group
generalized-$ \alpha $
method
下载: 导出CSV
-
[1] Newmark NM. A method of computation for structural dynamics. Journal of the Engineering Mechanics Division, 1959, 85(3): 67-94 doi: 10.1061/JMCEA3.0000098
[2] Chung J, Hulbert G. A time integration algorithm for structural dynamics with improved numerical dissipation: The generalized-α method. ASME Journal of Applied Mechanics, 1993, 60: 371-375
[3] Wood WL, Bossak M, Zienkiewicz OC. An alpha modification of Newmark's method. International Journal for Numerical Methods in Engineering, 1980, 15(10): 1562-1566 doi: 10.1002/nme.1620151011
[4] Hilber HM, Hughes TJR, Taylor RL. Improved numerical dissipation for time integration algorithms in structural dynamics. Earthquake Engineering and Structural Dynamics, 1977, 5(3): 283-292 doi: 10.1002/eqe.4290050306
[5] Bazilevs Y, Calo VM, Cottrell JA, et al. Variational multiscale residual-based turbulence modeling for large eddy simulation of incompressible flows. Computer Methods in Applied Mechanics and Engineering, 2007, 197(1-4): 173-201 doi: 10.1016/j.cma.2007.07.016
[6] Behnoudfar P, Calo VM, Deng Q, et al. A variationally separable splitting for the generalized-α method for parabolic equations. International Journal for Numerical Methods in Engineering, 2020, 121(5): 828-841 doi: 10.1002/nme.6246
[7] Gomez H, Hughes TJR, Nogueira X, et al. Isogeometric analysis of the isothermal Navier–Stokes–Korteweg equations. Computer Methods in Applied Mechanics and Engineering, 2010, 199(25-28): 1828-1840 doi: 10.1016/j.cma.2010.02.010
[8] Behnoudfar P, Deng Q, Calo VM. Higher-order generalized-α methods for hyperbolic problems. Computer Methods in Applied Mechanics and Engineering, 2021, 378: 113725 doi: 10.1016/j.cma.2021.113725
[9] Yen J, Petzold L, Raha S. A time integration algorithm for flexible mechanism dynamics: The DAE α-method. Computer Methods in Applied Mechanics and Engineering, 1998, 158(3-4): 341-355 doi: 10.1016/S0045-7825(97)00261-2
[10] Arnold M, Brüls O. Convergence of the generalized-α scheme for constrained mechanical systems. Multibody System Dynamics, 2007, 18(2): 185-202 doi: 10.1007/s11044-007-9084-0
[11] 张雄, 王天书. 计算动力学. 北京: 清华大学出版社, 2007 Zhang Xiong, Wang Tianshu. Computational Dynamics. Beijing: Tsinghua University Press, 2007: 236-273 (in Chinese))
[12] 田强. 基于绝对节点坐标方法的柔性多体系统动力学研究与应用. [博士论文]. 武汉: 华中科技大学, 2009 Tian Qiang. Flexible multibody dynamics research and application based on the absolute nodal coordinate method. [PhD Thesis]. Wuhan: Huazhong University of Science and Technology, 2009 (in Chinese))
[13] Behnoudfar P, Deng Q, Calo VM. High-order generalized-α method. Applications in Engineering Science, 2020, 4: 100021 doi: 10.1016/j.apples.2020.100021
[14] Ji Y, Xing Y, Wiercigroch M. An unconditionally stable time integration method with controllable dissipation for second-order nonlinear dynamics. Nonlinear Dynamics, 2021, 105(4): 3341-3358
[15] 季奕, 邢誉峰. 一种求解瞬态热传导方程的无条件稳定方法. 力学学报, 2021, 53(7): 1951-1961. (Ji Yi, , Xing Yufeng. An unconditionally stable method for transient heat conduction. Chinese Journal of Theoretical and Applied Mechanics, 2021, 53(7): 1951-1961(in Chinese)) doi: 10.1007/s11071-021-06720-9 Ji Y, Xing Y, Wiercigroch M. An unconditionally stable time integration method with controllable dissipation for second-order nonlinear dynamics[J]. Nonlinear Dynamics, 2021, 105(4): 3341-3358. (in Chinese) doi: 10.1007/s11071-021-06720-9
[16] Brüls O, Cardona A, Arnold M. Lie group generalized-α time integration of constrained flexible multibody systems. Mechanism and Machine Theory, 2012, 48: 121-137 doi: 10.1016/j.mechmachtheory.2011.07.017
[17] Arnold M, Brüls O, Cardona A. Error analysis of generalized-α Lie group time integration methods for constrained mechanical systems. Numerische Mathematik, 2015, 129(1): 149-179 doi: 10.1007/s00211-014-0633-1
[18] 刘铖, 胡海岩. 基于李群局部标架的多柔体系统动力学建模与计算. 力学学报, 2021, 53(1): 213-233 Liu Cheng, Hu Haiyan. Dynamic modeling and computation for flexible multibody systems based on the local frame of Lie group. Chinese Journal of Theoretical and Applied Mechanics, 2021, 53(1): 213-233(in Chinese))
[19] Wieloch V, Arnold M. BDF integrators for constrained mechanical systems on Lie groups. Journal of Computational and Applied Mathematics, 2021, 387: 112517 doi: 10.1016/j.cam.2019.112517
[20] Marsden JE, West M. Discrete mechanics and variational integrators. Acta Numerica, 2001, 10: 357-514 doi: 10.1017/S096249290100006X
[21] Lew A, Marsden JE, Ortiz M, et al. An overview of variational integrators. In: Finite Element Methods: 1970s and Beyond, CIMNE Barcelona, 2004: 85-146
[22] Wendlandt JM, Marsden JE. Mechanical integrators derived from a discrete variational principle. Physica D:Nonlinear Phenomena, 1997, 106(3-4): 223-246 doi: 10.1016/S0167-2789(97)00051-1
[23] Leyendecker S, Ober-Blöbaum S, Marsden JE, et al. Discrete mechanics and optimal control for constrained systems. Optimal Control Applications and Methods, 2010, 31(6): 505-528 doi: 10.1002/oca.912
[24] Kane C, Marsden JE, Ortiz M. Symplectic-energy-momentum preserving variational integrators. Journal of Mathematical Physics, 1999, 40(7): 3353-3371 doi: 10.1063/1.532892
[25] Cendra H, Marsden JE, Ratiu TS. Geometric Mechanics, Lagrangian Reduction, and Nonholonomic Systems//Mathematics Unlimited—2001 and Beyond. Berlin, Heidelberg: Springer, 2001: 221-273
[26] Hante S, Arnold M. RATTLie: A variational Lie group integration scheme for constrained mechanical systems. Journal of Computational and Applied Mathematics, 2021, 387: 112492 doi: 10.1016/j.cam.2019.112492
[27] Ober-Blöbaum S, Vermeeren M. Superconvergence of Galerkin variational integrators. IFAC-Papers on Line, 2021, 54(19): 327-333 doi: 10.1016/j.ifacol.2021.11.098
[28] Yi Z, Yue B, Deng M. Chebyshev spectral variational integrator and applications. Applied Mathematics and Mechanics, 2020, 41(5): 753-768 doi: 10.1007/s10483-020-2602-8
[29] Bloch AM, Marsden JE, Zenkov DV. Quasivelocities and symmetries in non-holonomic systems. Dynamical Systems, 2009, 24(2): 187-222 doi: 10.1080/14689360802609344
[30] Ball KR, Zenkov DV. Hamel’s formalism and variational integrators. Fields Institute Communications, 2015, 73: 477-506
[31] An Z, Wu H, Shi D. Minimum-time optimal control of robotic manipulators based on Hamel’s integrators. Meccanica, 2019, 54(15): 2521-2537 doi: 10.1007/s11012-019-01093-1
[32] Shi D, Berchenko-Kogan Y, Zenkov DV, et al. Hamel’s formalism for infinite-dimensional mechanical systems. Journal of Nonlinear Science, 2017, 27(1): 241-283 doi: 10.1007/s00332-016-9332-7
[33] An Z, Gao S, Shi D, et al. A variational integrator for the Chaplygin–Timoshenko Sleigh. Journal of Nonlinear Science, 2020, 30(4): 1381-1419 doi: 10.1007/s00332-020-09611-2
[34] Shi D, Zenkov DV, Bloch AM. Hamel’s formalism for classical field theories. Journal of Nonlinear Science, 2020, 30(4): 1307-1353 doi: 10.1007/s00332-020-09609-w
[35] 王亮, 安志朋, 史东华. 几何精确梁的 Hamel 场变分积分子[J]. 北京大学学报: 自然科学版, 2016 (4): 692-698 Wang Liang, An Zhipeng, Shi Donghua. Hamel’s field variational integrator of geometrically exact beam[J]. Acta Scientiarum Naturalium Universitatis Pekinensis, 2016, 52: 692-698(in Chinese))
[36] 高山, 史东华, 郭永新. Hamel 框架下几何精确梁的离散动量守恒律. 力学学报, 2021, 53(6): 1712-1719 doi: 10.6052/0459-1879-21-092 Gao Shan, Shi Donghua, Guo Yongxin. Discrete momentum conservation law of geometrically exact beam in Hamel’s framework. Chinese Journal of Theoretical and Applied Mechanics, 2021, 53(6): 1712-1719(in Chinese)) doi: 10.6052/0459-1879-21-092
[37] Simo JC, Tarnow N, Wong KK. Exact energy-momentum conserving algorithms and symplectic schemes for nonlinear dynamics. Computer Methods in Applied Mechanics and Engineering, 1992, 100(1): 63-116 doi: 10.1016/0045-7825(92)90115-Z
[38] Kane C, Marsden JE, Ortiz M, et al. Variational integrators and the Newmark algorithm for conservative and dissipative mechanical systems. International Journal for Numerical Methods in Engineering, 2000, 49(10): 1295-1325 doi: 10.1002/1097-0207(20001210)49:10<1295::AID-NME993>3.0.CO;2-W
[39] Bardella L, Genna F. Newmark's time integration method from the discretization of extended functionals. Journal of Applied Mechanics, 2005, 72(4): 527-537 doi: 10.1115/1.1934648
[40] 汤惠颖, 张志娟, 刘铖等. 两类基于局部标架梁单元的闭锁缓解方法. 力学学报, 2021, 53(2): 482-495 doi: 10.6052/0459-1879-20-274 Tang Huiying, Zhang Zhijuan, Liu Cheng, et al. Locking alleviation techniques of two types of beam elements based on the local frame formulation. Chinese Journal of Theoretical and Applied Mechanics, 2021, 53(2): 482-495(in Chinese)) doi: 10.6052/0459-1879-20-274
[41] Brüls O. Integrated Simulation and Reduced-order Modeling of Controlled Flexible Multibody Systems. Liège, Belgique: Université de Liège, 2005
计量
- 文章访问数: 734
- HTML全文浏览量: 349
- PDF下载量: 139
