力学学报 2018 , 50 (4): 863-870 https://doi.org/10.6052/0459-1879-18-111

Orginal Article

多柔体系统数值分析的模型降噪方法

齐朝晖¹^,^*^,, 曹艳¹, 王刚²^,^*^,

1大连理工大学工程力学系, 大连 116023
2大连理工大学盘锦校区海洋科学与技术学院, 辽宁盘锦 124221;

MODEL SMOOTHING METHODS IN NUMERICAL ANALYSIS OF FLEXIBLE MULTIBODY SYSTEMS

Qi Zhaohui¹^,^*^,, Cao Yan¹, Wang Gang²^,^*^,

1Department of Engineering Mechanics, Dalian University of Technology, Dalian 116023, China
2School of Ocean Science and Technology, Dalian University of Technology, Panjin 124221, Liaoning, China ;

中图分类号: O313.7

文献标识码: A

通讯作者: ^*齐朝晖, 教授, 主要研究方向: 多体系统动力学. E-mail: zhaohuiq@dlut.edu.cn ;^*王刚, 讲师, 主要研究方向: 多体系统动力学. E-mail:wanggangdut@dlut.edu.cn^*齐朝晖, 教授, 主要研究方向: 多体系统动力学. E-mail: zhaohuiq@dlut.edu.cn ;^*王刚, 讲师, 主要研究方向: 多体系统动力学. E-mail:wanggangdut@dlut.edu.cn^*齐朝晖, 教授, 主要研究方向: 多体系统动力学. E-mail: zhaohuiq@dlut.edu.cn ;^*王刚, 讲师, 主要研究方向: 多体系统动力学. E-mail:wanggangdut@dlut.edu.cn^*齐朝晖, 教授, 主要研究方向: 多体系统动力学. E-mail: zhaohuiq@dlut.edu.cn ;^*王刚, 讲师, 主要研究方向: 多体系统动力学. E-mail:wanggangdut@dlut.edu.cn

基金资助: 国家自然科学基金资助项目(91748203, 11372057).

展开

摘要

多柔体系统的动力学方程通常是一组刚性微分方程, 目前普遍采用的刚性微分方程数值解法主要通过数值阻尼滤除系统响应中的高频分量, 其求解效率难以令人满意. 为了降低多柔体系统动力学方程的刚性, 从而可采用ODE45等常规微分方程求解器进行求解, 研究了在建模过程中滤除高频振荡分量的方法. 在以当前时刻为起点的短时间内对柔性体的应力进行均匀化, 用均匀化后的应力计算柔性体的变形虚功率, 由此得到的系统动力学方程的解中不含过高频率的弹性振动, 并且可以通过调节均匀化时间区间的长度参数控制滤波的范围. 数值算例表明: 这种模型降噪方法的计算效率和精度均不低于刚性微分方程求解器, 并且在刚性微分方程求解器失效的情况下模型降噪方法仍有良好的精度和效率. 本文所提的模型降噪方法可成为求解多柔体系统动力学方程的新途径.

关键词： 多柔体系统 ; 模型光滑化 ; 刚性微分方程 ; 变形虚功率 ; 虚功率方程

Abstract

Dynamic equations of flexible multibody systems are usually a set of stiff differential equations. At present, the common numerical method for solving the stiff differential equations filters out the high frequency by using the numerical damping. The computational efficiency of this method is still unsatisfactory. In order to reduce the stiffness of dynamic equations of flexible multibody systems so greatly that the equations can be solved by regular ordinary differential equation (ODE) solvers such as MATLAB ODE45 solver, methods of filtering high frequency vibrations during the process of modeling are studied. Stresses of flexible bodies are homogenized by their mean value over a time interval from now to a short time later. The homogenized stress is then employed to replace its origin when computing the virtual deformation power. In this way, the obtained model of the flexible multibody system will not contain harmful high frequency elastic vibrations. The range of frequencies can be controlled by the length of the time interval used to homogenize stresses. As validated by the numerical examples in this paper, the precision and efficiency of the proposed method are comparable to some stiff ODE solvers. Moreover, it works well when the stiff ODE solver fails to give correct solutions in a reasonable time. Comparisons of numerical examples show that the proposed method can be a new available approach to numerical analysis of flexible multibody systems.

Keywords： flexible multibody systems ; model smoothing ; stiff ODEs ; virtual deformation power ; principle of virtual power

PDF (2374KB) 元数据多维度评价相关文章收藏文章

本文引用格式导出 EndNote Ris Bibtex

齐朝晖, 曹艳, 王刚. 多柔体系统数值分析的模型降噪方法[J]. 力学学报, 2018, 50(4): 863-870 https://doi.org/10.6052/0459-1879-18-111

Qi Zhaohui, Cao Yan, Wang Gang. MODEL SMOOTHING METHODS IN NUMERICAL ANALYSIS OF FLEXIBLE MULTIBODY SYSTEMS[J]. Acta Mechanica Sinica, 2018, 50(4): 863-870 https://doi.org/10.6052/0459-1879-18-111

多柔体系统是由多个可变形柔性部件通过铰相互连接而组成的系统. 实际工程中很多机械系统都可抽象为多柔体系统, 例如空间太阳能帆板、周边桁架式可展开天线等. 由于研究对象的复杂多样性, 多柔体系统不仅在建模理论方面不断面临新的挑战. 其动力学方程的数值求解也存在许多难点问题^[1,2,3,4].

困难之一源于铰的约束: 如果直接采用笛卡尔坐标和转动矢量描述物体在惯性坐标系中的运动, 铰的约束体现为描述参数间的代数约束方程; 如果采用物体间相对运动的参数(铰坐标)作为系统变量, 对含闭环的多体系统, 铰坐标必须满足闭环约束方程. 因而, 一般多体系统的动力学方程是微分代数混合方程组(DAEs), 它的求解本身就是计算数学的难点问题^[5,6,7,8,9]. 除此之外, 多体系统动力学对DAEs的求解方法还有更高的要求: 位置约束方程及其派生出的速度和加速度约束方程的违约量都必须被同时抑制在合理的范围内. 经国内外众多学者的长期探索, 已经提出了约束稳定化方法、广义坐标分块方法、局部参数化缩并方法、违约修正方法等一系列现实可行的方法^{[10,11,12,13]}.

即使多柔体系统中系统变量间没有约束, 求解其系统动力学方程也是十分困难的. 众所周知, 常微分方程组 $\dot{y} = f (t, y)$ 存在唯一解的条件是: 存在一个常数 $L$ 使以下不等式成立

$| | f (t, y_{2}) - f (t, y_{1}) | | L | | y_{2} - y_{1} | |$

Lipschitz常数越大, 数值求解的难度也越大^[14,15,16]. 多柔体系统中物体的高频振荡弹性运动可使系统动力学方程的Lipschitz常数高达 $1 0^{6}$ 量级以上. 因而, 在多柔体系统数值分析中有一个物理上十分不合理的现象: 当一个柔性体的弹性模量很高, 物体的运动趋近于相对简单的刚体运动时, 数值分析的难度反而大幅度提升. 分析多柔体系统时需要采用刚性微分方程求解器, 否则数值解将难以在合理的时间内获得. 刚性微分方程求解方法可划分为两类^{[17,18,19,20,21]}: 一类以标准一阶微分方程为研究对象, 其主要方法包括: BDF 方法、隐式Runge-Kutta方法、Rosenbrock方法等; 另一类以二阶微分方程为研究对象, 其主要方法包括: Newmark方法、Wilson- $θ$ 方法、HHT方法、广义 $α$ 方法等. 由于多柔体系统动力学方程属于微分代数混合方程, 还需要将这些方法与DAE解法相结合. 其主要方法有^{[22,23,24,25,26,27,28,29]}: HHT 双循环隐式积分法、BDF双循环隐式积分法、广义 $α$ 投影方法以及新型广义 $α$ 方法等. 所有这些方法都有一个共同点: 采用数值阻尼滤除了系统中的高频分量. 正是因为滤除了高频分量, 才能使积分步长远远大于最小振动周期时不引起积累误差的放大传播. 既然系统模型中的高频分量在求解时也要被滤除, 如何在建模时就滤除这些高频分量是非常值得研究的课题.

本文结合多柔体系统的特点, 研究如何在建立多柔体系统动力学方程的过程中控制模型中的高频分量成分, 并通过数值实验比较模型降噪方法与典型方法在精度和效率方面的差异. 以期为大幅度改善多柔体系统动力学方程的数值性态提供一种新的途径.

1 多柔体系统数值求解的困难

多柔体系统中物体的运动可分解为物体连体基的整体运动和相对于连体基的弹性运动^[30]. 物体的弹性运动通常具有幅值小频率高的特性. 它对物体位移的贡献很小, 但却使物体加速度以可观的幅度高频振荡. 如图1所示, 由两根弹簧连接3个质量

显示原图| 下载原图ZIP| 生成PPT

图1 弹簧质量系统.块组成的弹簧质量系统, 初始静止, 在水平力 $f$ 的作用下开始运动. 其中, $m$ 是质量块的质量, $k 1$ 和 $k 2$ 是两根弹簧的刚度, $x 1$ , $x 2$ 和 $x 3$ 分别是3个质量块质心相对初始位置的位移. 在 $k 1 = k 2 = k$ 时,系统的动力学方程为

Fig.1 Spring-mass system

$m {\ddot{x}}_{1} + k (x_{1} - x_{2}) = 0$

$m {\ddot{x}}_{2} + 2 k x_{2} - k (x_{1} + x_{3}) = 0$

$m {\ddot{x}}_{3} + k x_{3} - k x_{2} = f$

如果刚度 $k$ 远远大于质量 $m$ , 整个系统可看作为一个刚体. 此时, 各物体的位移均为

$x_{0} = \frac{1}{2} a_{0} t^{2}$

其中, $a_{0} = f / (3 m)$ 为整体运动的加速度. 如果不对刚度做任何限制, 系统动力学方程的解析解为

$x_{1} = x_{0} - a_{1} (1 - \cos ω_{1} t) + a_{2} (1 - \cos ω_{2} t)$

$x_{2} = x_{0} - 2 a_{2} (1 - \cos ω_{2} t)$

$x_{3} = x_{0} + a_{1} (1 - \cos ω_{1} t) + a_{2} (1 - \cos ω_{2} t)$

其中, 弹性振动频率 $ω_{1} = \sqrt[]{k / m}$ , $ω_{2} = \sqrt[]{3} ω_{1}$ ; 弹性振动的振幅为 $a_{1} = f / (2 k)$ , $a_{2} = f / (18 k)$ . 各物体运动加速度分别为

${\ddot{x}}_{1} = a_{0} + \frac{1}{2} a_{0} (\cos ω_{2} t - 3 \cos ω_{1} t)$

${\ddot{x}}_{2} = a_{0} - a_{0} \cos ω_{2} t$

${\ddot{x}}_{3} = a_{0} + \frac{1}{2} a_{0} (3 \cos ω_{1} t + \cos ω_{2} t)$

从中可见: 弹性振动使各物体位移在整体位移 $x_{0}$ 基础上发生了扰动. 当刚度 $k$ 逐渐增大时, 扰动幅值逐渐衰减, 各物体的位移越来越趋近于整体位移 $x_{0}$ ; 然而, 各物体的加速度并没有随着刚度的增加而趋近于整体运动加速度 $a_{0}$ , 而是以与 $a_{0}$ 可比的振幅和越来越高的频率振荡(见图2), 造成系统动力学方程的数值求解反而变得越来越困难. 大量的算例表明: 多柔体系统中物体的位移和加速度都具有上述特征^[29]。

显示原图| 下载原图ZIP| 生成PPT

图2 柔性体的位移和加速度.

Fig.2 Displacement and acceleration of a flexible body

其解具有上述特征的微分方程属于刚性微分方程. 如果采用MATLAB ODE45, ODE113等常规求解器求解这类系统, 为避免累积误差放大传播, 积分步长 $Δ t$ 将受到最小振动周期 $T_{\min}$ 的限制(如ODE45要求 $Δ t 2.78 T_{\min}$ ). 多柔体系统中物体弹性振动的最高频率常常高达 $1 0^{6}$ 量级以上, 用常规求解器仿真系统1 s内的运动就需求解系统动力学方程百万次以上. 除此之外, 常规求解器大都属于变步长求解器, 它在用如此小的步长成功积分后, 会将其放大几倍作为下一步积分的尝试步长, 然后在此基础上逐步缩减直至再次找到合适的步长, 期间的大部分计算是无效的, 从而进一步大幅度地降低了求解效率^[14]. 如果采用MATLAB ODE15s, ODE23tb等刚性微分方程求解器, 可在积分步长远远大于最小振动周期的情况下仍能保证数值解的稳定性. 众所周知, 如果积分步长大于某个振动的周期, 数值解就不可能准确描述它随时间的变化^[14,15]. 一个合理的问题是: 用刚性微分方程求解器求得的解是一种什么意义下的解. 事实上, 绝大部分刚性微分方程数值方法中, 都有一个核心步骤: 用多个时刻的位移值计算 $t + Δ t$ 时刻的速度和加速度. 这相当于对系统方程进行了数值滤波, 滤除了对系统位移贡献很小但频率很高的分量. 例如, 如果用ODE15s求解图1 所示的系统(质量 $m = 10 kg$ 、刚度 $k = 1.0 \times 1 0^{5} N / m$ 、外力 $f = 20 N$ 、求解器精度参数: 相对精度 $ε_{r} = 1 0^{- 2}$ 、绝对精度 $ε_{} a = 1 0^{- 3}$ ), 则所得的各物体运动是整体运动 $x_{0}$ , 即位移均为 $t^{2} / 3$ , 加速度均为 $2 / 3 m / s^{2}$ . 刚性微分方程求解器的滤波作用大小与对解的精度要求密切相关. 如果对数值解的精度要求较高, 积分步长可能被迫缩短至某些高频振荡周期之内, 从而导致滤波作用失效. 例如, 如果将精度参数提高为 $ε_{r} = 1.0 \times 1 0^{- 6}$ , $ε_{a} = 1.0 \times 1 0^{- 10}$ , 则ODE15s 和ODE45所得结果基本相同, 但此时ODE15s 甚至比ODE45更加耗时. 2多柔体系统的模型降噪建立多柔体系统动力学方程的物理依据是虚功率方程^[29,30,31]

式中, $u$ 为物体中一点的位移; $〼$ 和 $倻$ 分别为该点处的广义应力和广义应变; $f$ 为作用在该点处的外力.

物体位移中高频弹性振动分量主要源于随时间快速变化的应力 $〼$ . 在区间 $(t, t + h)$ 内, $τ$ 时刻的应力可近似为

$\begin{equation*}\textbf{{σ}}_\tau\approx\textbf{{σ}}_t+(\tau-t)\dot{\textbf{{σ}}}_t+(\tau-t)^2\ddot{\textbf{{σ}}}_t/2\tag*{(13)}\end{equation*}$

该区间内应力平均值

为了滤除模型中的高频分量, 用平均应力$\bar{\textbf{\textit{σ}}}$替换应力$\textbf{\textit{σ}}$, 从而将变形虚功率近似为

这种方法的滤波作用可通过比较一个线性系统在修改前后的动力学方程性质说明. 当应变‒位移关系以及应力‒应变关系均为线性时, 变形虚功率可用系统刚度矩阵 $K$ 表示为

按上述方法修改后的变形虚功率

线性系统的动力学方程可以写作

$Mu + Ku = F$ (18)

将变形虚功率按式(17)近似后, 系统动力学方程成为

$(M + h^{2} K / 6) u + hKu / 2 + Ku = F$ (19)

与修改前的系统动力学方程相比, 它增加了惯性项和阻尼项, 从而可降低系统的固有频率. 参数 $h$ 越大, 保留在系统中的高频分量就越少. 这些特点可在将系统动力学方程转化为解耦的模态坐标方程后变得更加鲜明. 首先, 求解广义特征值问题

$\begin{equation*}(\textbf{{K}}-\omega^2\textbf{{M}})\textbf{{ν}}=\textbf{0}\tag*{(20)}\end{equation*}$

利用其特征向量组成的矩阵 $T = [퀻$ 对位移进行模态坐标变换: $u = Tw$ . 由于矩阵 $T$ 满足等式: $T_{} TKT = diag (ω_{1}^{2}, ω_{2}^{2}, \cdot \cdot \cdot, ω_{n}^{2})$ 以及 $T_{} TMT = E$ , 因而, 修改前后系统的模态坐标分别满足方程

${\ddot{w}}_{i} + ω_{i}^{2} w_{i} = f_{i}, i = 1,2, \cdot \cdot \cdot, n$

$(1 + \frac{h^{2}}{6} ω_{i}^{2}) {\ddot{\overset{̅}{w}}}_{i} + ω_{i}^{2} (\frac{h}{2} {\dot{\overset{̅}{w}}}_{i} + {\overset{̅}{w}}_{i}) = f_{i}, i = 1,2, \cdot \cdot \cdot, n$

其中, 模态力 $f_{i}$ 为向量 $T^{T} F$ 的分量. 修改前后模态坐标变化频率与振幅的比值分别为

$λ_{ω} ≜ (1 + h^{2} ω_{i}^{2} / 6)^{- 1} (1 + 5 h^{2} ω_{i}^{2} / 48)^{\frac{1}{2}}$

$λ_{a} ≜ (1 + h^{2} ω_{i}^{2} / 6)^{- 1}$

这两个比值随系统频率 $ω_{i}$ 和模型降噪参数 $h$ 的变化如图3所示. 从图3可以看到: (1)无论参数 $h$ 取何值, 修改后的模型不影响系统的刚体运动(频率为零的模态); (2)修改后的模型降低了高频分量的振幅和频率, 低频分量衰减幅度小, 高频分量衰减幅度大; (3)随着 $h$ 值的增大高频分量逐渐减少, 反之 $h$ 值逐渐减小, 模型降噪对系统的扰动也随之降低.

显示原图| 下载原图ZIP| 生成PPT

图3 频率比与振幅比.

Fig.3 Amplitudes and frequencies ratios

模型降噪相当于对原系统进行了一次模态减缩, 但又不必具体求解系统的模态, 只需调整降噪参数 $h$ 的大小就可调整保留在系统中的高频分量.

3 数值算例

利用模型降噪方法滤除系统高频分量后, 系统动力学方程可用任何现有微分方程求解器求解. 如果先降噪再用ODE45求解并取 $h = 0.1$ , 则将这种方法简写为MS $_$ ODE45 $_$ 0.1; 如果用ODE15s 求解并取 $h = 0.01$ , 则将其缩写为MS $_$ ODE15s $_$ 0.01.

例1图1所示的系统, 其系统参数如表1所示.

表1 系统参数

Table 1 System parameters

$m$ /kg	$k_{1}$ /(N $\cdot$ m $^{- 1})$	$k_{2}$ /(N $\cdot$ m $^{- 1})$	$f$ /N
	$1 0^{7}$	$1 0^{8}$	$20$

新窗口打开

用ODE45, ODE15s和MS $_$ ODE45 $_$ 0.1分别对其在0 $~$ 5 s内的运动进行数值求解. 当精度控制参数分别为 $ε_{} r = 1.0 \times 1 0^{- 4}$ 和 $ε_{} a = 1.0 \times 1 0^{- 6}$ 时, 各方法所得结果与系统位移解析解

$x_{1} = x_{2} = x_{3} = t^{2} / 3$

解析解和数值解之间的绝对误差以及求解所耗时间如表2所示. 从中可见: 3种方法所得结果均具有很高的精度, 但计算效率相差很大. 其中,ODE15s最差, 说明刚性微分方程求解器不适合采用高精度控制参数. 模型降噪方法MS $_$ ODE45 $_$ 0.1在这种情况下计算效率明显高于其他两种方法.

表2 精度与效率比较( $ε r = 1.0 × 1 0 - 4$ , $ε a = 1.0 × 1 0 - 6$ )

Table 2 Comparisons of precision and efficiency ( $ε r = 1.0 × 1 0 - 4$ , $ε a = 1.0 × 1 0 - 6$ )

Solvers	Maximum absolute error	Time/s
ODE45	$1.5 \times 1 0^{- 5}$	376.9
ODE15s	$2.3 \times 1 0^{- 5}$	884.4
MS $_$ ODE45 $_$ 0.1	$1.6 \times 1 0^{- 5}$	0.27

新窗口打开

当精度控制参数分别取 $ε_{r} = 1.0 \times 1 0^{- 2}$ 和 $ε_{a} = 1.0 \times 1 0^{- 3}$ 时, 解析解和数值解之间的绝对误差以及求解所耗时间如表3所示.

表3 精度与效率比较( $ε r = 1.0 × 1 0 - 2$ , $ε a = 1.0 × 1 0 - 3$ )

Table 3 Comparisons of precision and efficiency ( $ε r = 1.0 × 1 0 - 2$ , $ε a = 1.0 × 1 0 - 3$ )

Solvers	Maximum absolute error	Time/s
ODE45	$8.8 \times 1 0^{11}$	376.7
ODE15s	$1.6 \times 1 0^{- 4}$	0.24
MS $_$ ODE45 $_$ 0.1	$2.2 \times 1 0^{- 3}$	0.24

新窗口打开

在这种情况下, 非刚性微分方程求解器ODE45因积累误差放大传播导致数值解发散. ODE15s和本文所提方法MS $_$ ODE45 $_$ 0.1都在较短的时间内得到了具有合理精度的数值解.

例2 为检验本文所提方法处理刚柔耦合问题的能力, 将图1所示系统中两个刚度间的差异增大, 令 $k_{1} = 2.0 \times 1 0^{3}$ N/m和 $k_{2} = 1.0 \times 1 0^{8}$ N/m, 同时保持质量 $m$ 和外力 $f$ 不变. 在这种情况下, 系统位移的解析解为

$x_{1} = t^{2} / 3 - 2 [1 - \cos (10 \sqrt[]{3} t)] / 900$

$x_{2} = x_{3} = t^{2} / 3 + [1 - \cos (10 \sqrt[]{3} t)] / 900$

图4为采用MS $_$ ODE45 $_$ 0.001方法所求得的数值解与解析解之间的比较. 从中可见: 模型降噪法给出的质点位移 $x_{1}$ 和 $x_{2}$ 与解析解式(26)和式(27)是吻合的.

显示原图| 下载原图ZIP| 生成PPT

图4 位移数值解与解析解比较.

Fig.4 Comparisons of numerical displacements and their analytical solution

解析解和数值解之间的绝对误差以及求解所耗时间如表4所示. 从中可以看到: 模型降噪法在 $h$ 取不同值时都得到了具有足够精度的数值解; 降噪后采用ODE15s求解时效率更高, 但精度一般低于采用ODE45时的精度; 参数 $h$ 值越小, 解的精度越高; $h$ 值越大, 耗时越短.

表4为取精度控制参数 $ε_{} r = 1.0 \times 1 0^{- 3}$ 和 $ε_{a} = 1.0 \times 1 0^{- 6}$ 时, 系统在^[0,5] s内的运动数值解的精度和效率.

表4 精度与效率比较(刚柔耦合系统)

Table 4 Comparisons of precision and efficiency (rigid-flexible coupling system)

Solvers	Maximum absolute error	Time/s
MS $_$ ODE45 $_$ 0.1	$2.34 \times 1 0^{- 3}$	0.31
MS $_$ ODE45 $_$ 0.01	$2.16 \times 1 0^{- 3}$	1.06
MS $_$ ODE45 $_$ 0.001	$6.83 \times 1 0^{- 4}$	5.80
MS $_$ ODE15s $_$ 0.001	$2.35 \times 1 0^{- 2}$	0.57

新窗口打开

例3 本算例将给出本文所提方法在柔性多体系统动力学中的应用. 基于文献^[32]所提的应变插值大变形平面柔性梁的建模方法, 将梁的轴向应变和曲率修改为

$\begin{array}{l} \overset{̅}{ε} = ε_{t} + \frac{1}{2} h {\dot{ε}}_{t} + \frac{1}{6} h^{2} {\ddot{ε}}_{t} \\ \overset{̅}{κ} = κ_{t} + \frac{1}{2} h {\dot{κ}}_{t} + \frac{1}{6} h^{2} {\ddot{κ}}_{t} \end{array}\}$

因此, 修正后的节点内力为

$f_{int} = K_{e} 瀼$

其中, $瀼$ 是单元应变参数列阵, 包括单元轴向应变和两端的节点处曲率; $K_{e}$ 是单元刚度矩阵, 其具体形式为

$\begin{array}{l} K_{e} = [\begin{array}{l} EAL & 0 & 0 \\ 0 & EIL / 3 & EIL / 6 \\ 0 & EIL / 6 & EIL / 3 \end{array}] \end{array}$

其中, $EA$ 和 $EI$ 分别为梁的抗拉和抗弯刚度; $EI$ 是梁单元长度.

对于图5所示的双摆机构, 材料密度为 $ρ = 7801 kg / m^{3}$ , 弹性模量为 $E = 2.1 \times 1 0^{11}$ Pa. 两根摆长均为 $L = 1.8 m$ , 截面积均为 $A = 2.5 \times 1 0^{- 3} m^{2}$ . 机构在 $g = - 9.81 m / s^{2}$ 重力加速度作用下, 从水平静止位置落下.

显示原图| 下载原图ZIP| 生成PPT

图5 自由下落的柔性双摆机构.

Fig.5 Free-falling flexible double pendulum mechanism

分别采用MATLAB的刚性方程求解器radau5 $(h = 0)$ 和ODE45求解器 $(h = 0.001; h = 0.01)$ 分别对该问题进行数值积分, 取精度控制参数 $ε_{r} = 1.0 \times 1 0^{- 3}$ 和 $ε_{a} = 1.0 \times 1 0^{- 6}$ . 得到双摆机构上末端点在竖直方向和水平方向的位移, 分别如图6和图7所示.

显示原图| 下载原图ZIP| 生成PPT

图6 末端点的轴向变形.

Fig.6 Axial deformation of the end point

显示原图| 下载原图ZIP| 生成PPT

图7 末端点的横向挠度.

Fig.7 Lateral deflection of the end point

由图6和图7可以看出, 当柔性双摆的材料刚度较大时 $(E = 2.1 \times 1 0^{11} Pa)$ , 柔性梁的变形很小. 采用本文所提的模型降噪方法可以使用非刚性微分方程求解器ODE45求解该问题, 并且结果与刚性微分方程求解器radau5的求解结果很接近. 此外, 采用模型降噪方法能够更好地滤除系统中的高频, 得到更为平滑的变形曲线.

表5对比了上述不同方法的计算时间, 数据显示: 如果不使用模型降噪方法, 普通的微分方程求解器ODE45无法在有限的时间内获得系统的数值解; 使用模型降噪技术后, 对于刚度较大的柔性多体系统, 本文所提方法的计算效率更高.

表5 效率比较(柔性双摆机构)

Table 5 Efficiency comparison (flexible double pendulum mechanism)

Solvers	Time/s
radau5	22 680
MS $_$ ODE45 $_$ 0.01	14
MS $_$ ODE45 $_$ 0.001	116
MS $_$ ODE45	- -

新窗口打开

4 结论

在多柔体系统的虚功率方程中用 $(t, t + h)$ 区间内的平均应力替代 $t$ 时刻应力, 所建立的系统动力学方程可有效滤除解中的高频振荡分量, 从而大幅度降低了系统动力学方程的刚性. 降噪后的系统动力学方程可采用如ODE45等常规求解器求解也可采用如ODE15s, radau5等刚性微分方程求解器求解. 参数 $h$ 可控制保留在系统中的高频分量, $h$ 值越大滤波效果越明显. 在滤除了高频分量的前提下, $h$ 值越小数值解的精度越高.

The authors have declared that no competing interests exist.

参考文献

原文顺序

文献年度倒序

文中引用次数倒序

被引期刊影响因子

[1]	田强, 刘铖, 李培等. 多柔体系统动力学研究进展与挑战 . 动力学与控制学报, 2017, 15(5): 385-405 [本文引用: 1] (Tian Qiang, Liu Cheng, Li Pei, et al. Advances and challenges in dynamics of flexible multibody systems .Journal of Dynamics and Control, 2017, 15(5): 385-405 (in Chinese)) [本文引用: 1]
[2]	王琪, 庄方方, 郭易圆等. 非光滑多体系统动力学数值算法的研究进展 . 力学进展, 2013, 43(1): 101-111 [本文引用: 1] (Wang Qi, Zhuang Fangfang, Guo Yiyuan, et al. Advances in the research on numerical methods for non-smooth dynamics of multibody systems .Advances in Mechanics, 2013, 43(1): 101-111 (in Chinese)) [本文引用: 1]
[3]	刘才山, 陈滨, 彭瀚等. 多体系统多点碰撞接触问题的数值求解方法 . 动力学与控制学报, 2003, 1(1): 59-65 [本文引用: 1] (Liu Caishan, Chen Bin, Peng Han, et al. Numerical resolution of multibody systems with multiple contact/impact points .Journal of Dynamics and Control, 2003, 1(1): 59-65 (in Chinese)) [本文引用: 1]
[4]	姚文莉, 岳嵘. 有争议的碰撞恢复系数研究进展 . 振动与冲击, 2015, 34(19): 43-48 [本文引用: 1] (Yao Wenli, Yue Rong. Advance in controversial restitution coefficient study for impact problems .Journal of Vibration and Shock, 2015, 34(19): 43-48 (in Chinese)) [本文引用: 1]
[5]	Stechlinski PG, Barton PI. Dependence of solutions of nonsmooth differential-algebraic equations on parameters .Journal of Differential Equations, 2017, 262(3): 2254-2285 [本文引用: 2]
[6]	Dhamacharoen A. Efficient numerical methods for solving differential algebraic equations .Journal of Applied Mathematics ＆ Physics, 2016, 4(1): 1-9 [本文引用: 1]
[7]	Petzold L. Differential/algebraic equations are not ODE’S .Siam Journal on Scientific ＆ Statistical Computing, 2012, 3(3): 367-384 [本文引用: 1]
[8]	Shampine LF, Reichelt MW, Kierzenka JA. Solving index-I DAEs in MATLAB and simulink .Siam Review, 1999, 41(3): 538-552 [本文引用: 1]
[9]	Gear CW, Petzold LR. ODE methods for the solution of differential /algebraic systems .Siam Journal on Numerical Analysis, 1984, 21(4): 716-728 [本文引用: 1]
[10]	Haddouni M, Acary V, Garreau S, et al. Comparison of several formulations and integration methods for the resolution of DAEs formulations in event-driven simulation of nonsmooth frictionless multibody dynamics .Multibody System Dynamics, 2017, 41(3): 201-231 [本文引用: 1]
[11]	Marques F, Souto A P, Flores P. On the constraints violation in forward dynamics of multibody systems .Multibody System Dynamics, 2016, 39(4): 1-35 [本文引用: 1]
[12]	Schweizer B, Lu D, Li P. Co-simulation method for solver coupling with algebraic constraints incorporating relaxation techniques .Multibody System Dynamics, 2016, 36(1): 1-36 [本文引用: 1]
[13]	潘振宽, 赵维加, 洪嘉振等. 多体系统动力学微分/代数方程组数值方法 . 力学进展, 1996, 26(1): 83-96 [本文引用: 1] (Pan Zhenkuan, Zhao Weijia, Hong Jiazhen, et al. On numerical algorithms for differential/algebraic equations of motion of multibody system .Advances in Mechanics, 1996, 26(1): 26-40 (in Chinese)) [本文引用: 1]
[14]	袁新鼎. 刚性常微分方程初值问题的数值解法. 北京: 科学出版社, 1987 [本文引用: 3] (Yuan Xinding.Numerical Methods for Solving Initial Value Problems of Stiff Ordinary Differential Equations. Beijing: Science Press, 1987 (in Chinese)) [本文引用: 3]
[15]	程正兴. 数值逼近与常微分方程数值解. 西安: 西安交通大学出版社, 2000 [本文引用: 2] (Cheng Zhengxing.Numerical Approximation and Numerical Solution of Ordinary Differential Equation. Xi’an: Xi’an Jiaotong University Press, 2000 (in Chinese)) [本文引用: 2]
[16]	Hairer E, Lubich C. Numerical Analysis of Ordinary Differential Equations. Berlin, Heidelberg: Springer 2015 [本文引用: 1]
[17]	Ibrahim Z, Nasir NM, Othman K, et al. Adaptive order of block backward differentiation formulas for stiff ODEs .Numerical Algebra, 2017, 7(1): 95-106 [本文引用: 1]
[18]	El-Zahar ER, Habib HM, Rashidi MM, et al. A comparison of explicit semi-analytical numerical integration methods for solving stiff ODE systems .American Journal of Applied Sciences, 2015, 12(5): 304-320 [本文引用: 1]
[19]	Ariel G, Engquist B, Tsai R. A multiscale method for highly oscillatory ordinary differential equations with resonance .Mathematics of Computation, 2009, 78(266): 929-956 [本文引用: 1]
[20]	Tokman M. Efficient integration of large stiff systems of ODEs with exponential propagation iterative (EPI) methods .Journal of Computational Physics, 2006, 213(2): 748-776 [本文引用: 1]
[21]	Meijaard JP. Application of Runge-Kutta-Rosenbrock methods to the analysis of flexible multibody systems .Multibody System Dynamics, 2003, 10(3): 263-288 [本文引用: 1]
[22]	Shabana AA, Hussein BA. A two-loop sparse matrix numerical integration procedure for the solution of differential/algebraic equations: Application to multibody systems .Journal of Sound ＆ Vibration, 2009, 327(3): 557-563 [本文引用: 1]
[23]	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): 341-355 [本文引用: 1]
[24]	Arnold M, Brüls O. Convergence of the generalized-α, scheme for constrained mechanical systems .Multibody System Dynamics, 2007, 18(2): 185-202 [本文引用: 1]
[25]	丁洁玉, 潘振宽. 多体系统动力学微分‒代数方程广义<inline-formula><mml:math xmlns:mml="http://www.w3.org/1998/Math/MathML" id="Mml184-0459-1879-50-4-863"><mml:mo>-</mml:mo><mml:mi>α</mml:mi></mml:math></inline-formula>投影法 . 工程力学, 2013, 30(4): 380-384 [本文引用: 1] (Ding Jieyu, Pan Zhenkuan. generalized<inline-formula><mml:math xmlns:mml="http://www.w3.org/1998/Math/MathML" id="Mml185-0459-1879-50-4-863"><mml:mo>-</mml:mo><mml:mi>α</mml:mi></mml:math></inline-formula> projection method for differential-algebraic equations of multibody dynamics .Engineering Mechanics, 2013, 30(4): 380-384 (in Chinese)) [本文引用: 1]
[26]	郭晛, 章定国, 陈思佳. Hilber-Hughes-Taylor<inline-formula><mml:math xmlns:mml="http://www.w3.org/1998/Math/MathML" id="Mml186-0459-1879-50-4-863"><mml:mo>-</mml:mo><mml:mi>α</mml:mi></mml:math></inline-formula> 法在接触约束多体系统动力学中的应用 . 物理学报, 2017, 66(16): 144-154 [本文引用: 1] (Guo Xian, Zhang Dingguo, Chen Sijia. Application of Hilber-Hughes-Taylor<inline-formula><mml:math xmlns:mml="http://www.w3.org/1998/Math/MathML" id="Mml187-0459-1879-50-4-863"><mml:mo>-</mml:mo><mml:mi>α</mml:mi></mml:math></inline-formula> method to dynamics of flexible multibody system with contact and constraint .Acta Physica Sinica, 2017, 66(16): 144-154 (in Chinese)) [本文引用: 1]
[27]	姚廷强, 迟毅林, 黄亚宇. 柔性多体系统动力学新型广义<inline-formula><mml:math xmlns:mml="http://www.w3.org/1998/Math/MathML" id="Mml188-0459-1879-50-4-863"><mml:mo>-</mml:mo><mml:mi>α</mml:mi></mml:math></inline-formula>数值分析方法 . 机械工程学报, 2009, 45(10): 53-60 [本文引用: 1] (Yao Tingqiang, Chi Yilin, Hu Yayu. New generalized-α algorithms for multibody dynamics .Journal of Mechanical Engineering, 2009, 45(10): 53-60 (in Chinese)) [本文引用: 1]
[28]	郭晛, 章定国. 柔性多体系统动力学典型数值积分方法的比较研究 . 南京理工大学学报(自然科学版), 2016, 40(6): 726-733 [本文引用: 1] (Guo Xian, Zhang Dingguo. Comparative study of typical numerical integration methods of flexible multi-body systems dynamics .Journal of Nanjing University of Science and Technology, 2016, 40(6): 726-733 (in Chinese)) [本文引用: 1]
[29]	Banerjee A. Flexible Multibody Dynamics: Efficient Formulations and Applications. John Wiley ＆ Sons , Inc., 2016 [本文引用: 3]
[30]	洪嘉振. 计算多体系统动力学. 北京: 高等教育出版社, 1999 [本文引用: 2] (Hong Jiazhen.Computational Dynamics of Multibody Systems. Beijing: Higher Education Press, 1999 (in Chinese)) [本文引用: 2]
[31]	齐朝晖. 多体系统动力学. 北京: 科学出版社, 2008 [本文引用: 1] (Qi Zhaohui.Dynamics of Multibody Systems. Beijing: Science Press, 2008 (in Chinese)) [本文引用: 1]
[32]	张志刚, 齐朝晖, 吴志刚. 一种基于应变插值大变形梁单元的刚‒柔耦合动力学分析 . 振动工程学报, 2015, 28(3): 337-344 [本文引用: 1] (Zhang Zhigang, Qi Zhaohui, Wu Zhigang. Rigid-flexible dynamics analysis of a large deformation beam element based on interpolation of strains .Journal of Vibration Engineering, 2015, 28(3): 337-344 (in Chinese)) [本文引用: 1]

多柔体系统动力学研究进展与挑战

2017

... 多柔体系统是由多个可变形柔性部件通过铰相互连接而组成的系统. 实际工程中很多机械系统都可抽象为多柔体系统, 例如空间太阳能帆板、周边桁架式可展开天线等. 由于研究对象的复杂多样性, 多柔体系统不仅在建模理论方面不断面临新的挑战. 其动力学方程的数值求解也存在许多难点问题^[1,2,3,4]. ...

多柔体系统动力学研究进展与挑战

2017

非光滑多体系统动力学数值算法的研究进展

2013

非光滑多体系统动力学数值算法的研究进展

2013

多体系统多点碰撞接触问题的数值求解方法

2003

多体系统多点碰撞接触问题的数值求解方法

2003

有争议的碰撞恢复系数研究进展

2015

有争议的碰撞恢复系数研究进展

2015

Dependence of solutions of nonsmooth differential-algebraic equations on parameters

2017

... 困难之一源于铰的约束: 如果直接采用笛卡尔坐标和转动矢量描述物体在惯性坐标系中的运动, 铰的约束体现为描述参数间的代数约束方程; 如果采用物体间相对运动的参数(铰坐标)作为系统变量, 对含闭环的多体系统, 铰坐标必须满足闭环约束方程. 因而, 一般多体系统的动力学方程是微分代数混合方程组(DAEs), 它的求解本身就是计算数学的难点问题^[5,6,7,8,9]. 除此之外, 多体系统动力学对DAEs的求解方法还有更高的要求: 位置约束方程及其派生出的速度和加速度约束方程的违约量都必须被同时抑制在合理的范围内. 经国内外众多学者的长期探索, 已经提出了约束稳定化方法、广义坐标分块方法、局部参数化缩并方法、违约修正方法等一系列现实可行的方法^{[10,11,12,13]}. ...

... 表4为取精度控制参数

ε_{} r = 1.0 \times 1 0^{- 3}

和

ε_{a} = 1.0 \times 1 0^{- 6}

时, 系统在^[0,5] s内的运动数值解的精度和效率. ...

Efficient numerical methods for solving differential algebraic equations

2016

Differential/algebraic equations are not ODE’S

2012

Solving index-I DAEs in MATLAB and simulink

1999

ODE methods for the solution of differential /algebraic systems

1984

Comparison of several formulations and integration methods for the resolution of DAEs formulations in event-driven simulation of nonsmooth frictionless multibody dynamics

2017

On the constraints violation in forward dynamics of multibody systems

2016

Co-simulation method for solver coupling with algebraic constraints incorporating relaxation techniques

2016

多体系统动力学微分/代数方程组数值方法

1996

多体系统动力学微分/代数方程组数值方法

1996

1987

... Lipschitz常数越大, 数值求解的难度也越大^[14,15,16]. 多柔体系统中物体的高频振荡弹性运动可使系统动力学方程的Lipschitz常数高达

1 0^{6}

量级以上. 因而, 在多柔体系统数值分析中有一个物理上十分不合理的现象: 当一个柔性体的弹性模量很高, 物体的运动趋近于相对简单的刚体运动时, 数值分析的难度反而大幅度提升. 分析多柔体系统时需要采用刚性微分方程求解器, 否则数值解将难以在合理的时间内获得. 刚性微分方程求解方法可划分为两类^{[17,18,19,20,21]}: 一类以标准一阶微分方程为研究对象, 其主要方法包括: BDF 方法、隐式Runge-Kutta方法、Rosenbrock方法等; 另一类以二阶微分方程为研究对象, 其主要方法包括: Newmark方法、Wilson-

θ

方法、HHT方法、广义

α

方法等. 由于多柔体系统动力学方程属于微分代数混合方程, 还需要将这些方法与DAE解法相结合. 其主要方法有^{[22,23,24,25,26,27,28,29]}: HHT 双循环隐式积分法、BDF双循环隐式积分法、广义

α

投影方法以及新型广义

α

方法等. 所有这些方法都有一个共同点: 采用数值阻尼滤除了系统中的高频分量. 正是因为滤除了高频分量, 才能使积分步长远远大于最小振动周期时不引起积累误差的放大传播. 既然系统模型中的高频分量在求解时也要被滤除, 如何在建模时就滤除这些高频分量是非常值得研究的课题. ...

... 其解具有上述特征的微分方程属于刚性微分方程. 如果采用MATLAB ODE45, ODE113等常规求解器求解这类系统, 为避免累积误差放大传播, 积分步长

Δ t

将受到最小振动周期

T_{\min}

的限制(如ODE45要求

Δ t  2.78 T_{\min}

). 多柔体系统中物体弹性振动的最高频率常常高达

1 0^{6}

量级以上, 用常规求解器仿真系统1 s内的运动就需求解系统动力学方程百万次以上. 除此之外, 常规求解器大都属于变步长求解器, 它在用如此小的步长成功积分后, 会将其放大几倍作为下一步积分的尝试步长, 然后在此基础上逐步缩减直至再次找到合适的步长, 期间的大部分计算是无效的, 从而进一步大幅度地降低了求解效率^[14]. 如果采用MATLAB ODE15s, ODE23tb等刚性微分方程求解器, 可在积分步长远远大于最小振动周期的情况下仍能保证数值解的稳定性. 众所周知, 如果积分步长大于某个振动的周期, 数值解就不可能准确描述它随时间的变化^[14,15]. 一个合理的问题是: 用刚性微分方程求解器求得的解是一种什么意义下的解. 事实上, 绝大部分刚性微分方程数值方法中, 都有一个核心步骤: 用多个时刻的位移值计算

t + Δ t

时刻的速度和加速度. 这相当于对系统方程进行了数值滤波, 滤除了对系统位移贡献很小但频率很高的分量. 例如, 如果用ODE15s求解图1 所示的系统(质量

m = 10 kg

、刚度

k = 1.0 \times 1 0^{5} N / m

、外力

f = 20 N

、求解器精度参数: 相对精度

ε_{r} = 1 0^{- 2}

、绝对精度

ε_{} a = 1 0^{- 3}

), 则所得的各物体运动是整体运动

x_{0}

, 即位移均为

t^{2} / 3

, 加速度均为

2 / 3 m / s^{2}

. 刚性微分方程求解器的滤波作用大小与对解的精度要求密切相关. 如果对数值解的精度要求较高, 积分步长可能被迫缩短至某些高频振荡周期之内, 从而导致滤波作用失效. 例如, 如果将精度参数提高为

ε_{r} = 1.0 \times 1 0^{- 6}

ε_{a} = 1.0 \times 1 0^{- 10}

, 则ODE15s 和ODE45所得结果基本相同, 但此时ODE15s 甚至比ODE45更加耗时. 2多柔体系统的模型降噪建立多柔体系统动力学方程的物理依据是虚功率方程^[29,30,31] ...

... [14,15]. 一个合理的问题是: 用刚性微分方程求解器求得的解是一种什么意义下的解. 事实上, 绝大部分刚性微分方程数值方法中, 都有一个核心步骤: 用多个时刻的位移值计算

t + Δ t

时刻的速度和加速度. 这相当于对系统方程进行了数值滤波, 滤除了对系统位移贡献很小但频率很高的分量. 例如, 如果用ODE15s求解图1 所示的系统(质量

m = 10 kg

、刚度

k = 1.0 \times 1 0^{5} N / m

、外力

f = 20 N

、求解器精度参数: 相对精度

ε_{r} = 1 0^{- 2}

、绝对精度

ε_{} a = 1 0^{- 3}

), 则所得的各物体运动是整体运动

x_{0}

, 即位移均为

t^{2} / 3

, 加速度均为

2 / 3 m / s^{2}

ε_{r} = 1.0 \times 1 0^{- 6}

ε_{a} = 1.0 \times 1 0^{- 10}

1987

... Lipschitz常数越大, 数值求解的难度也越大^[14,15,16]. 多柔体系统中物体的高频振荡弹性运动可使系统动力学方程的Lipschitz常数高达

1 0^{6}

θ

方法、HHT方法、广义

α

α

投影方法以及新型广义

α

... 其解具有上述特征的微分方程属于刚性微分方程. 如果采用MATLAB ODE45, ODE113等常规求解器求解这类系统, 为避免累积误差放大传播, 积分步长

Δ t

将受到最小振动周期

T_{\min}

的限制(如ODE45要求

Δ t  2.78 T_{\min}

). 多柔体系统中物体弹性振动的最高频率常常高达

1 0^{6}

t + Δ t

时刻的速度和加速度. 这相当于对系统方程进行了数值滤波, 滤除了对系统位移贡献很小但频率很高的分量. 例如, 如果用ODE15s求解图1 所示的系统(质量

m = 10 kg

、刚度

k = 1.0 \times 1 0^{5} N / m

、外力

f = 20 N

、求解器精度参数: 相对精度

ε_{r} = 1 0^{- 2}

、绝对精度

ε_{} a = 1 0^{- 3}

), 则所得的各物体运动是整体运动

x_{0}

, 即位移均为

t^{2} / 3

, 加速度均为

2 / 3 m / s^{2}

ε_{r} = 1.0 \times 1 0^{- 6}

ε_{a} = 1.0 \times 1 0^{- 10}

t + Δ t

时刻的速度和加速度. 这相当于对系统方程进行了数值滤波, 滤除了对系统位移贡献很小但频率很高的分量. 例如, 如果用ODE15s求解图1 所示的系统(质量

m = 10 kg

、刚度

k = 1.0 \times 1 0^{5} N / m

、外力

f = 20 N

、求解器精度参数: 相对精度

ε_{r} = 1 0^{- 2}

、绝对精度

ε_{} a = 1 0^{- 3}

), 则所得的各物体运动是整体运动

x_{0}

, 即位移均为

t^{2} / 3

, 加速度均为

2 / 3 m / s^{2}

ε_{r} = 1.0 \times 1 0^{- 6}

ε_{a} = 1.0 \times 1 0^{- 10}

2000

... Lipschitz常数越大, 数值求解的难度也越大^[14,15,16]. 多柔体系统中物体的高频振荡弹性运动可使系统动力学方程的Lipschitz常数高达

1 0^{6}

θ

方法、HHT方法、广义

α

α

投影方法以及新型广义

α

... 其解具有上述特征的微分方程属于刚性微分方程. 如果采用MATLAB ODE45, ODE113等常规求解器求解这类系统, 为避免累积误差放大传播, 积分步长

Δ t

将受到最小振动周期

T_{\min}

的限制(如ODE45要求

Δ t  2.78 T_{\min}

). 多柔体系统中物体弹性振动的最高频率常常高达

1 0^{6}

t + Δ t

时刻的速度和加速度. 这相当于对系统方程进行了数值滤波, 滤除了对系统位移贡献很小但频率很高的分量. 例如, 如果用ODE15s求解图1 所示的系统(质量

m = 10 kg

、刚度

k = 1.0 \times 1 0^{5} N / m

、外力

f = 20 N

、求解器精度参数: 相对精度

ε_{r} = 1 0^{- 2}

、绝对精度

ε_{} a = 1 0^{- 3}

), 则所得的各物体运动是整体运动

x_{0}

, 即位移均为

t^{2} / 3

, 加速度均为

2 / 3 m / s^{2}

ε_{r} = 1.0 \times 1 0^{- 6}

ε_{a} = 1.0 \times 1 0^{- 10}

2000

... Lipschitz常数越大, 数值求解的难度也越大^[14,15,16]. 多柔体系统中物体的高频振荡弹性运动可使系统动力学方程的Lipschitz常数高达

1 0^{6}

θ

方法、HHT方法、广义

α

α

投影方法以及新型广义

α

... 其解具有上述特征的微分方程属于刚性微分方程. 如果采用MATLAB ODE45, ODE113等常规求解器求解这类系统, 为避免累积误差放大传播, 积分步长

Δ t

将受到最小振动周期

T_{\min}

的限制(如ODE45要求

Δ t  2.78 T_{\min}

). 多柔体系统中物体弹性振动的最高频率常常高达

1 0^{6}

t + Δ t

时刻的速度和加速度. 这相当于对系统方程进行了数值滤波, 滤除了对系统位移贡献很小但频率很高的分量. 例如, 如果用ODE15s求解图1 所示的系统(质量

m = 10 kg

、刚度

k = 1.0 \times 1 0^{5} N / m

、外力

f = 20 N

、求解器精度参数: 相对精度

ε_{r} = 1 0^{- 2}

、绝对精度

ε_{} a = 1 0^{- 3}

), 则所得的各物体运动是整体运动

x_{0}

, 即位移均为

t^{2} / 3

, 加速度均为

2 / 3 m / s^{2}

ε_{r} = 1.0 \times 1 0^{- 6}

ε_{a} = 1.0 \times 1 0^{- 10}

Numerical Analysis of Ordinary Differential Equations. Berlin,

2015

... Lipschitz常数越大, 数值求解的难度也越大^[14,15,16]. 多柔体系统中物体的高频振荡弹性运动可使系统动力学方程的Lipschitz常数高达

1 0^{6}

θ

方法、HHT方法、广义

α

α

投影方法以及新型广义

α

Adaptive order of block backward differentiation formulas for stiff ODEs

2017

... Lipschitz常数越大, 数值求解的难度也越大^[14,15,16]. 多柔体系统中物体的高频振荡弹性运动可使系统动力学方程的Lipschitz常数高达

1 0^{6}

θ

方法、HHT方法、广义

α

α

投影方法以及新型广义

α

A comparison of explicit semi-analytical numerical integration methods for solving stiff ODE systems

2015

... Lipschitz常数越大, 数值求解的难度也越大^[14,15,16]. 多柔体系统中物体的高频振荡弹性运动可使系统动力学方程的Lipschitz常数高达

1 0^{6}

θ

方法、HHT方法、广义

α

α

投影方法以及新型广义

α

A multiscale method for highly oscillatory ordinary differential equations with resonance

2009

... Lipschitz常数越大, 数值求解的难度也越大^[14,15,16]. 多柔体系统中物体的高频振荡弹性运动可使系统动力学方程的Lipschitz常数高达

1 0^{6}

θ

方法、HHT方法、广义

α

α

投影方法以及新型广义

α

Efficient integration of large stiff systems of ODEs with exponential propagation iterative (EPI) methods

2006

... Lipschitz常数越大, 数值求解的难度也越大^[14,15,16]. 多柔体系统中物体的高频振荡弹性运动可使系统动力学方程的Lipschitz常数高达

1 0^{6}

θ

方法、HHT方法、广义

α

α

投影方法以及新型广义

α

Application of Runge-Kutta-Rosenbrock methods to the analysis of flexible multibody systems

2003

... Lipschitz常数越大, 数值求解的难度也越大^[14,15,16]. 多柔体系统中物体的高频振荡弹性运动可使系统动力学方程的Lipschitz常数高达

1 0^{6}

θ

方法、HHT方法、广义

α

α

投影方法以及新型广义

α

A two-loop sparse matrix numerical integration procedure for the solution of differential/algebraic equations: Application to multibody systems

2009

... Lipschitz常数越大, 数值求解的难度也越大^[14,15,16]. 多柔体系统中物体的高频振荡弹性运动可使系统动力学方程的Lipschitz常数高达

1 0^{6}

θ

方法、HHT方法、广义

α

α

投影方法以及新型广义

α

A time integration algorithm for flexible mechanism dynamics: the DAE α-method

1998

... Lipschitz常数越大, 数值求解的难度也越大^[14,15,16]. 多柔体系统中物体的高频振荡弹性运动可使系统动力学方程的Lipschitz常数高达

1 0^{6}

θ

方法、HHT方法、广义

α

α

投影方法以及新型广义

α

Convergence of the generalized-α, scheme for constrained mechanical systems

2007

... Lipschitz常数越大, 数值求解的难度也越大^[14,15,16]. 多柔体系统中物体的高频振荡弹性运动可使系统动力学方程的Lipschitz常数高达

1 0^{6}

θ

方法、HHT方法、广义

α

α

投影方法以及新型广义

α

多体系统动力学微分?代数方程广义<inline-formula><mml:math xmlns:mml="http://www.w3.org/1998/Math/MathML" id="Mml184-0459-1879-50-4-863"><mml:mo>-</mml:mo><mml:mi>α</mml:mi></mml:math></inline-formula>投影法

2013

... Lipschitz常数越大, 数值求解的难度也越大^[14,15,16]. 多柔体系统中物体的高频振荡弹性运动可使系统动力学方程的Lipschitz常数高达

1 0^{6}

θ

方法、HHT方法、广义

α

α

投影方法以及新型广义

α

2013

... Lipschitz常数越大, 数值求解的难度也越大^[14,15,16]. 多柔体系统中物体的高频振荡弹性运动可使系统动力学方程的Lipschitz常数高达

1 0^{6}

θ

方法、HHT方法、广义

α

α

投影方法以及新型广义

α

Hilber-Hughes-Taylor<inline-formula><mml:math xmlns:mml="http://www.w3.org/1998/Math/MathML" id="Mml186-0459-1879-50-4-863"><mml:mo>-</mml:mo><mml:mi>α</mml:mi></mml:math></inline-formula> 法在接触约束多体系统动力学中的应用

2017

... Lipschitz常数越大, 数值求解的难度也越大^[14,15,16]. 多柔体系统中物体的高频振荡弹性运动可使系统动力学方程的Lipschitz常数高达

1 0^{6}

θ

方法、HHT方法、广义

α

α

投影方法以及新型广义

α

2017

... Lipschitz常数越大, 数值求解的难度也越大^[14,15,16]. 多柔体系统中物体的高频振荡弹性运动可使系统动力学方程的Lipschitz常数高达

1 0^{6}

θ

方法、HHT方法、广义

α

α

投影方法以及新型广义

α

柔性多体系统动力学新型广义<inline-formula><mml:math xmlns:mml="http://www.w3.org/1998/Math/MathML" id="Mml188-0459-1879-50-4-863"><mml:mo>-</mml:mo><mml:mi>α</mml:mi></mml:math></inline-formula>数值分析方法

2009

... Lipschitz常数越大, 数值求解的难度也越大^[14,15,16]. 多柔体系统中物体的高频振荡弹性运动可使系统动力学方程的Lipschitz常数高达

1 0^{6}

θ

方法、HHT方法、广义

α

α

投影方法以及新型广义

α

2009

... Lipschitz常数越大, 数值求解的难度也越大^[14,15,16]. 多柔体系统中物体的高频振荡弹性运动可使系统动力学方程的Lipschitz常数高达

1 0^{6}

θ

方法、HHT方法、广义

α

α

投影方法以及新型广义

α

柔性多体系统动力学典型数值积分方法的比较研究

2016

... Lipschitz常数越大, 数值求解的难度也越大^[14,15,16]. 多柔体系统中物体的高频振荡弹性运动可使系统动力学方程的Lipschitz常数高达

1 0^{6}

θ

方法、HHT方法、广义

α

α

投影方法以及新型广义

α

柔性多体系统动力学典型数值积分方法的比较研究

2016

... Lipschitz常数越大, 数值求解的难度也越大^[14,15,16]. 多柔体系统中物体的高频振荡弹性运动可使系统动力学方程的Lipschitz常数高达

1 0^{6}

θ

方法、HHT方法、广义

α

α

投影方法以及新型广义

α

Flexible Multibody Dynamics: Efficient Formulations and Applications. John Wiley ＆ Sons

2016

... Lipschitz常数越大, 数值求解的难度也越大^[14,15,16]. 多柔体系统中物体的高频振荡弹性运动可使系统动力学方程的Lipschitz常数高达

1 0^{6}

θ

方法、HHT方法、广义

α

α

投影方法以及新型广义

α

... 从中可见: 弹性振动使各物体位移在整体位移

x_{0}

基础上发生了扰动. 当刚度

k

逐渐增大时, 扰动幅值逐渐衰减, 各物体的位移越来越趋近于整体位移

x_{0}

; 然而, 各物体的加速度并没有随着刚度的增加而趋近于整体运动加速度

a_{0}

, 而是以与

a_{0}

可比的振幅和越来越高的频率振荡(见图2), 造成系统动力学方程的数值求解反而变得越来越困难. 大量的算例表明: 多柔体系统中物体的位移和加速度都具有上述特征^[29]. ...

... 其解具有上述特征的微分方程属于刚性微分方程. 如果采用MATLAB ODE45, ODE113等常规求解器求解这类系统, 为避免累积误差放大传播, 积分步长

Δ t

将受到最小振动周期

T_{\min}

的限制(如ODE45要求

Δ t  2.78 T_{\min}

). 多柔体系统中物体弹性振动的最高频率常常高达

1 0^{6}

t + Δ t

时刻的速度和加速度. 这相当于对系统方程进行了数值滤波, 滤除了对系统位移贡献很小但频率很高的分量. 例如, 如果用ODE15s求解图1 所示的系统(质量

m = 10 kg

、刚度

k = 1.0 \times 1 0^{5} N / m

、外力

f = 20 N

、求解器精度参数: 相对精度

ε_{r} = 1 0^{- 2}

、绝对精度

ε_{} a = 1 0^{- 3}

), 则所得的各物体运动是整体运动

x_{0}

, 即位移均为

t^{2} / 3

, 加速度均为

2 / 3 m / s^{2}

ε_{r} = 1.0 \times 1 0^{- 6}

ε_{a} = 1.0 \times 1 0^{- 10}

1999

... 多柔体系统中物体的运动可分解为物体连体基的整体运动和相对于连体基的弹性运动^[30]. 物体的弹性运动通常具有幅值小频率高的特性. 它对物体位移的贡献很小, 但却使物体加速度以可观的幅度高频振荡. 如图1所示, 由两根弹簧连接3个质量 ...

... 其解具有上述特征的微分方程属于刚性微分方程. 如果采用MATLAB ODE45, ODE113等常规求解器求解这类系统, 为避免累积误差放大传播, 积分步长

Δ t

将受到最小振动周期

T_{\min}

的限制(如ODE45要求

Δ t  2.78 T_{\min}

). 多柔体系统中物体弹性振动的最高频率常常高达

1 0^{6}

t + Δ t

时刻的速度和加速度. 这相当于对系统方程进行了数值滤波, 滤除了对系统位移贡献很小但频率很高的分量. 例如, 如果用ODE15s求解图1 所示的系统(质量

m = 10 kg

、刚度

k = 1.0 \times 1 0^{5} N / m

、外力

f = 20 N

、求解器精度参数: 相对精度

ε_{r} = 1 0^{- 2}

、绝对精度

ε_{} a = 1 0^{- 3}

), 则所得的各物体运动是整体运动

x_{0}

, 即位移均为

t^{2} / 3

, 加速度均为

2 / 3 m / s^{2}

ε_{r} = 1.0 \times 1 0^{- 6}

ε_{a} = 1.0 \times 1 0^{- 10}

1999

... 其解具有上述特征的微分方程属于刚性微分方程. 如果采用MATLAB ODE45, ODE113等常规求解器求解这类系统, 为避免累积误差放大传播, 积分步长

Δ t

将受到最小振动周期

T_{\min}

的限制(如ODE45要求

Δ t  2.78 T_{\min}

). 多柔体系统中物体弹性振动的最高频率常常高达

1 0^{6}

t + Δ t

时刻的速度和加速度. 这相当于对系统方程进行了数值滤波, 滤除了对系统位移贡献很小但频率很高的分量. 例如, 如果用ODE15s求解图1 所示的系统(质量

m = 10 kg

、刚度

k = 1.0 \times 1 0^{5} N / m

、外力

f = 20 N

、求解器精度参数: 相对精度

ε_{r} = 1 0^{- 2}

、绝对精度

ε_{} a = 1 0^{- 3}

), 则所得的各物体运动是整体运动

x_{0}

, 即位移均为

t^{2} / 3

, 加速度均为

2 / 3 m / s^{2}

ε_{r} = 1.0 \times 1 0^{- 6}

ε_{a} = 1.0 \times 1 0^{- 10}

2008

... 其解具有上述特征的微分方程属于刚性微分方程. 如果采用MATLAB ODE45, ODE113等常规求解器求解这类系统, 为避免累积误差放大传播, 积分步长

Δ t

将受到最小振动周期

T_{\min}

的限制(如ODE45要求

Δ t  2.78 T_{\min}

). 多柔体系统中物体弹性振动的最高频率常常高达

1 0^{6}

t + Δ t

时刻的速度和加速度. 这相当于对系统方程进行了数值滤波, 滤除了对系统位移贡献很小但频率很高的分量. 例如, 如果用ODE15s求解图1 所示的系统(质量

m = 10 kg

、刚度

k = 1.0 \times 1 0^{5} N / m

、外力

f = 20 N

、求解器精度参数: 相对精度

ε_{r} = 1 0^{- 2}

、绝对精度

ε_{} a = 1 0^{- 3}

), 则所得的各物体运动是整体运动

x_{0}

, 即位移均为

t^{2} / 3

, 加速度均为

2 / 3 m / s^{2}

ε_{r} = 1.0 \times 1 0^{- 6}

ε_{a} = 1.0 \times 1 0^{- 10}

2008

... 其解具有上述特征的微分方程属于刚性微分方程. 如果采用MATLAB ODE45, ODE113等常规求解器求解这类系统, 为避免累积误差放大传播, 积分步长

Δ t

将受到最小振动周期

T_{\min}

的限制(如ODE45要求

Δ t  2.78 T_{\min}

). 多柔体系统中物体弹性振动的最高频率常常高达

1 0^{6}

t + Δ t

时刻的速度和加速度. 这相当于对系统方程进行了数值滤波, 滤除了对系统位移贡献很小但频率很高的分量. 例如, 如果用ODE15s求解图1 所示的系统(质量

m = 10 kg

、刚度

k = 1.0 \times 1 0^{5} N / m

、外力

f = 20 N

、求解器精度参数: 相对精度

ε_{r} = 1 0^{- 2}

、绝对精度

ε_{} a = 1 0^{- 3}

), 则所得的各物体运动是整体运动

x_{0}

, 即位移均为

t^{2} / 3

, 加速度均为

2 / 3 m / s^{2}

ε_{r} = 1.0 \times 1 0^{- 6}

ε_{a} = 1.0 \times 1 0^{- 10}

一种基于应变插值大变形梁单元的刚?柔耦合动力学分析

2015

... 例3 本算例将给出本文所提方法在柔性多体系统动力学中的应用. 基于文献^[32]所提的应变插值大变形平面柔性梁的建模方法, 将梁的轴向应变和曲率修改为 ...

一种基于应变插值大变形梁单元的刚?柔耦合动力学分析

2015

〈

〉

多柔体系统数值分析的模型降噪方法

MODEL SMOOTHING METHODS IN NUMERICAL ANALYSIS OF FLEXIBLE MULTIBODY SYSTEMS

1 多柔体系统数值求解的困难

3 数值算例

4 结 论

参考文献 文献选项 原文顺序 文献年度倒序 文中引用次数倒序 被引期刊影响因子

4 结论

参考文献

原文顺序

文献年度倒序

文中引用次数倒序

被引期刊影响因子