振动工程学报

共旋框架下二维柔性多体动力学分析策略

PDF下载

季奕 ¹^,² , 伍洋 ³^,⁴ , 张慧敏 ⁵ , 孙光辉 ¹

振动工程学报 | 2025,38(10): 2288-2296

收起

振动工程学报 | 2025, 38(10): 2288-2296

共旋框架下二维柔性多体动力学分析策略

全屏

季奕¹^,², 伍洋³^,⁴, 张慧敏⁵, 孙光辉¹

作者信息

^1.哈尔滨工业大学自主智能无人系统工信部重点实验室，黑龙江　哈尔滨　150001；

^2.北京理工大学飞行器动力学与控制教育部重点实验室，北京　100081；

^3.中物院高性能数值模拟软件中心，北京　100088；

^4.北京应用物理与计算数学研究所，北京　100088；

^5.北京宇航系统工程研究所，北京　100005

通讯作者:

季奕（1995—），女，博士，副教授。E-mail： jiyi0319@outlook.com

A novel strategy for two-dimensional flexible multibody systems in the corotational frame

Yi JI¹^,², Yang WU³^,⁴, Huimin ZHANG⁵, Guanghui SUN¹

Affiliations

^1.Key Laboratory of Autonomous Intelligent Unmanned Systems, Harbin Institute of Technology, Harbin 150001, China

^2.MOE Key Laboratory of Dynamics and Control of Flight Vehicle, Beijing Institute of Technology, Beijing 100081, China

^3.CAEP Software Center for High Performance Numerical Simulation, Beijing 100088, China

^4.Institute of Applied Physics and Computational Mathematics, Beijing 100088, China

^5.Beijing Institute of Astronautical Systems Engineering, Beijing 100005, China

doi: 10.16385/j.cnki.issn.1004-4523.202310039

文章导航

摘要

收起

针对具有强非线性和大范围运动的柔性多体系统，建立了一种精确且高效的动力学分析策略。在共旋坐标系下对柔性体进行动力学建模，再利用一种优化复合方法对建立的数学模型进行求解。由于动力学模型建立在共旋坐标系下，因此现有的先进线性空间单元可以被直接使用，从而大幅降低计算成本。此外，为了准确地计算动力学响应，本文使用了一种优化复合方法，并以最小化局部截断误差为目标对算法参数进行优化设计。优化后的复合方法具有二阶精度、无条件稳定性和可控数值耗散。最后，本文对一些经典的柔性多体系统进行计算。结果表明：在相同的精度下，与目前流行的基于绝对节点坐标公式和Generalized‑α方法的动力学分析策略相比，本文策略具有显著的效率优势。

关键词

多体动力学 / 柔性多体系统 / 共旋框架 / 参数优化 / 复合方法

Abstract

收起

This paper proposes an accurate and efficient solution strategy for analyzing dynamic responses of flexible multibody systems. In the proposed strategy, flexible structures are modeled in the corotational frame, then the discrete mathematical model is solved by an optimized composite method. Due to the introduction of the corotational frame, some advanced linear elements can be directly employed, dramatically decreasing computational costs. For accurately calculating dynamic responses, an optimized three-sub-step composite method is developed wherein algorithmic parameters are optimized for minimizing local truncation errors. The optimized composite method achieves second-order accuracy, unconditional stability, and controllable stability. Some classical flexible dynamic systems are solved in this paper, and numerical results show that compared to the currently popular solution strategy based on the absolute nodal coordinate formulation and the Generalized-α method, under the same computational accuracy, our strategy has great superiorities in efficiency.

Key words

multibody dynamic / flexible multibody systems / corotational frame / parameters optimization / composite method

引用本文

季奕, 伍洋, 张慧敏, 孙光辉. 共旋框架下二维柔性多体动力学分析策略. 振动工程学报, 2025 , 38 (10) : 2288 -2296 . DOI: 10.16385/j.cnki.issn.1004-4523.202310039

Yi JI, Yang WU, Huimin ZHANG, Guanghui SUN. A novel strategy for two-dimensional flexible multibody systems in the corotational frame[J]. Journal of Vibration Engineering, 2025 , 38 (10) : 2288 -2296 . DOI: 10.16385/j.cnki.issn.1004-4523.202310039

正文

收起

具有强非线性、大变形、大转动和大量自由度的柔性多体系统在航空航天、机械制造等领域广泛存在^[1‑2]。由于此类动力学系统的复杂性，通常无法找到对应的解析解，目前基于有限单元法和时间积分方法的数值求解策略占主导地位。

一般来说，多体系统中的柔性部件可以在浮动坐标系、惯性坐标系和共旋坐标系下进行建模。浮动坐标系源于刚体动力学，随后扩展到小变形问题^[3]，目前在对柔性多体的建模上仍具有较大的局限性。对于柔性多体建模，惯性坐标系具有广泛的应用，可以通过几何精确公式（geometric exact formulation，GEF）或绝对节点坐标公式（absolute nodal coordinate formulation，ANCF）来实现。几何精确公式最初是由SIMO^[4]提出的，在过去几十年间，一些基于几何精确公式的先进梁^[5]、板/壳^[6]单元被陆续提出。但是，由于在截面大转动上描述的冗余性限制了几何精确法在柔性多体系统建模领域的发展。为了解决这些限制，绝对节点坐标法被提出。绝对节点坐标公式最早出现在SHABANA^[7]的工作中，与几何精确法不同的是，绝对节点坐标法不再将有限旋转作为节点坐标，而是使用绝对位移和全局斜率作为单元坐标。利用单元形状函数和节点坐标在全局坐标系中定义了有限元上各质点的位置和变形。总之，几何精确法和绝对节点坐标法都能够准确地分析柔性结构的大范围运动。但是，当问题局限于小应变和有限转动时，惯性坐标系不如共旋坐标系简便。

共旋坐标系^[8‑11]来源于连续介质力学把刚体运动从整体运动中分离出来，该概念首先因WEMPNER^[8]、BELYTSCHKO等^[9]提出。共旋坐标系的关键思想是利用一个局部坐标系从总运动中提取纯变形，该坐标系随每个单元平移和旋转，但不随单元变形。由于共旋坐标系中单元的纯变形较小，因此可以使用几何线性有限元公式计算共旋坐标系中的单元刚度矩阵和内力矢量，进一步通过将单元刚度矩阵和内力矢量从共旋坐标系到全局坐标系的一致变换来考虑几何非线性效应。由于共旋坐标系不依赖于几何线性单元类型，因此可以使用大量现有的具有鲁棒性和高精度的几何线性单元，并将其扩展到共旋坐标系的几何非线性分析中。

对于空间离散后系统的动力学响应分析，时间积分方法是一种非常有用的数值工具^[12‑13]，其中一些方法，比如Newmark方法和后向差分公式，已经在许多商业软件得到广泛使用。时间积分方法最初是用来求解由常微分方程（ordinary differential equation，ODE）控制的结构动力学系统。为了追求更高的精度和效率，在经典时间积分方法的基础上发展了一些更先进的时间积分方法，包括参数法^[14‑15]、保能量法^[16‑17]、配点法^[18‑19]等。由于参数法具有二阶精度、耗散可控等优点，且易于实现，因此很自然地被推广到多体系统领域。目前，Generalized‑α方法^[14，20]作为参数方法的一种，已经被广泛地用于求解多体系统。此外，为了精确地计算转动变量，学者们还构造了Lie‑Generalized‑α方法^[21]。由于参数方法的动力学平衡方程仅满足于广义的时间点，而不是时间离散点，因此他们的加速度精度只有一阶。

为了解决这一问题，学者们提出了复合方法。这类方法首次出现在BANK等^[22]的工作中，然后由BATHE等^[23]进行概念化。复合方法的思想是将不同的时间积分方法有效地组合在一个时间步长中，同时提高了精度、效率和稳定性。更重要地是，复合方法在离散的时间点上严格满足动力学平衡方程。目前已经发展出一些性能较好的复合方法，包括WEN方法^[24]、KIM方法^[25]、TTBIF^[26]等。在这些复合方法中，文献[27‑28]提出的TTBIF在精度上具有显著的优势。

在此背景下，本文为柔性多体系统开发了一种基于共旋坐标系和TTBIF的动力学响应求解策略。数值试验表明，相比当前流行的基于绝对节点坐标系和Generalized‑α方法建立的求解策略，本文策略在精度和效率上均具有明显优势。

1　共旋框架下柔性多体动力学方程

收起

柔性多体动力学建模主要研究梁、板、壳的柔性多体模型，而梁作为最简单、最常见的一种柔性多体模型在工程领域广泛存在。在此背景下，本文主要研究共旋框架下基于二维梁单元的柔性多体动力学方程。

无论采取哪种坐标系对柔性多体进行建模，对应的动力学行为都可以利用如下微分代数方程进行描述：

{F I (q, q ˙, q ¨) + Φ q T (q) λ + F E (q) = Q (t) Φ (q, t) = 0

(1)

式中，q为广义坐标向量，包括平动和转动变量；F_I为惯性力；F_E为弹性力；

Φ

（q，t）为包含运动和几何约束的向量；

Φ q T

（q，t）为向量

Φ

（q，t）关于q求导的矩阵；λ为拉格朗日乘子；Q（t）为外部激励项。采用不同坐标系，式（1）中F_I和F_E的构造方式是不同的，下面介绍共旋框架下基于一般二维梁单元柔性多体系统动力学方程的建立。

共旋坐标系的核心是利用一个局部坐标系从总运动中提取纯变形。如图1所示，为了描述一个梁模型在平面的运动，首先对每个单元建立一个随单元进行平动和转动的局部坐标系x_LO_Ly_L，并规定：坐标原点O_L位于梁单元左节点1；x_L沿着梁单元两节点连线，由节点1指向节点2。β和β₀为当前构型和初始构型下x_L轴与x_G轴的夹角；α=β-β₀。

几何线性的梁单元在全局坐标系x_GO_Gy_G下的广义坐标向量为：

q = [u 1 w 1 θ 1 u 2 w 2 θ 2] T

(2)

式中，u_i、w_i和θ_i（i=1，2）分别为节点i沿x_G方向的变形量、沿y_G方向的变形量和截面转角。对应其在局部坐标系x_LO_Ly_L下的广义坐标向量为：

q L = [u ¯ θ ¯ 1 θ ¯ 2] T

(3)

式中，

u ¯

=l_c-l₀，

θ ¯

₁=

θ

₁-β，

θ ¯

₂=

θ

₂-β，其中l_c和l₀分别表示当前和初始时刻的梁长。

共旋坐标法建模的核心就是建立式（2）和（3）中的投影关系。对式（3）进行变分可得：

δ q L = B δ q

(4)

式中，B为二者之间的投影矩阵，对应的元素为：

B = [b 1 b 2 b 3] = [- c - s 0 c s 0 - s / l c c / l c 1 s / l c - c / l c 0 - s / l c c / l c 0 s / l c - c / l c 1]

(5)

其中：

{c = c o s β = (x 2 + u 2) - (x 1 + u 1) l c s = s i n β = (y 2 + w 2) - (y 1 + w 1) l c

(6)

式中，x_i（i=1，2）分别为节点i在完成刚体转动后在全局坐标系下的坐标。

梁单元的弹性势能为：

V = δ q T F E = δ q L T F E, L = δ q T B T F E, L

(7)

式中，F_E和F_E，L分别为全局和局部坐标系下的弹性力。

因为式（7）对任意δq均成立，因此可以得到F_E和F_E，L之间的关系为：

F E = B T F E, L

(8)

其中：

F E, L = [N M 1 M 2] T

(9)

式中，N、M₁和M₂分别表示轴力、节点1和2的弯矩。

梁单元的动能为：

K = 12 [∫ l 0 A ρ (u ˙ G 2 + w ˙ G 2) d l + ∫ l 0 I ρ θ ˙ G 2 d l]

(10)

式中，

A ρ

和

I ρ

分别表示单位长度质量和转动惯量；l为梁长。

进一步地将

u ˙

_G、

w ˙

_G和

θ ˙

_G关于x的表达式代入式（10），积分后可得矩阵形式的动能表达式，即

K = 12 q ˙ T T T M L T q ˙

(11)

其中，旋转矩阵T的表达式如下：

T = [c s 0000 - s c 0000 001000 000 c s 0 000 - s c 0 000001]

(12)

局部坐标系下质量矩阵

M L

的表达式如下：

M L = M L, 1 + M L, 2

(13)

其中：

M L, 1 = I ρ 30 l 0 [000000 036 3 l 0 0 - 36 3 l 0 0 3 l 0 4 l 02 0 - 3 l 0 - l 02 000000 0 - 36 - 3 l 0 036 - 3 l 0 0 3 l 0 - l 02 0 - 3 l 0 4 l 02]

(14)

M L, 2 = A ρ l 0 420 [140 m 1 070 - m 1 0 m 1 156 22 l 0 m 2 54 - 13 l 0 0 22 l 0 4 l 02 0 13 l 0 - 3 l 02 70 m 2 0140 - m 2 0 - m 1 54 13 l 0 m 2 156 - 22 l 0 0 - 13 l 0 - 3 l 02 0 22 l 0 4 l 02]

(15)

m 1 = 21 θ ¯ 1 - 14 θ ¯ 2, m 2 = 14 θ ¯ 1 - 21 θ ¯ 2

(16)

由此可得梁单元的惯性力为：

F I = d d t ∂ K ∂ q ˙ - ∂ K ∂ q

(17)

本文利用Newton迭代方法求解非线性微分代数方程（1）。与常微分方程的求解类似，式（1）中第一行表达式里惯性力F_I和弹性力F_E关于状态变量q、

q ˙

和

q ¨

的求导可定义为切线质量矩阵、切线阻尼矩阵和切线刚度矩阵，即

M = ∂ F I ∂ q ¨

(18)

C = ∂ F I ∂ q ˙

(19)

K = ∂ F E ∂ q + ∂ F I ∂ q

(20)

此外，式（1）中约束方程

Φ

（q，t）关于状态变量

q

、

q ˙

和

q ¨

的求导可见后文。

2　基于复合方法的时域求解公式

收起

针对由常微分方程控制的结构动力学系统，与现有的时间积分方法相比，TTBIF在位移、速度和加速度的计算方面均展现出优势。在这种情况下，本工作计划将TTBIF拓展到动力学行为更加复杂的由微分代数方程控制的柔性多体系统。本节对TTBIF的时间步进格式和算法参数进行描述。

在TTBIF中，一个时间步[t，t+γΔt]被划分成3个部分，分别为[t，t+γΔt]、[t+γΔt，t+2γΔt]和[t+2γΔt，t+Δt]，如图2所示。

前2个分步为了最大程度保留重要的低频信息，均采用梯形法则（trapezoidal rule，TR），对应的位移、速度和加速度之间的假设关系如下：

{q t + γ Δ t = q t + 12 γ Δ t (q ˙ t + q ˙ t + γ Δ t) q ˙ t + γ Δ t = q ˙ t + 12 γ Δ t (q ¨ t + q ¨ t + γ Δ t)

(21)

{q t + 2 γ Δ t = q t + γ Δ t + 12 γ Δ t (q ˙ t + γ Δ t + q ˙ t + 2 γ Δ t) q ˙ t + 2 γ Δ t = q ˙ t + γ Δ t + 12 γ Δ t (q ¨ t + γ Δ t + q ¨ t + 2 γ Δ t)

(22)

为了能够过滤掉由于空间离散等引入的虚假高频信息，最后一个分步使用四点欧拉后向插值格式，对应的时间步进方程为：

{q t + Δ t = q t + Δ t (θ 0 q ˙ t + θ 1 q ˙ t + γ Δ t + θ 2 q ˙ t + 2 γ Δ t + θ 3 q ˙ t + Δ t) q ˙ t + Δ t = q ˙ t + Δ t (θ 0 q ¨ t + θ 1 q ¨ t + γ Δ t + θ 2 q ¨ t + 2 γ Δ t + θ 3 q ¨ t + Δ t)

(23)

从式（21）~（23）中可以看出，TTBIF的数值性能完全由5个算法参数控制，分别为：γ、θ₀、θ₁、θ₂和θ₃。

首先利用一些基本数值性能建立部分算法参数之间的关系，包括：

①二阶精度，即局部截断误差^[12]σ=A（ωΔt）- A_exact（ωΔt）=O（Δt²），其中A（ωΔt）表示TTBIF的数值传递矩阵，解析传递矩阵的表达式为：A_exact（ωΔt）=exp（-ωΔt）；

②可控高频耗散：即谱半径ρ（ωΔt→∞）=ρ_∞，ρ_∞为人为引入的算法参数，取值区间为[0，1]，其中谱半径ρ^[12]的定义为TTBIF特征根模的最大值，即ρ=max|λ_i|≤1（i=1，2，3）；

③无条件稳定性：即谱半径0≤ρ（0≤ωΔt≤∞）≤1。

利用上述三个条件，可以得到如下参数关系，即

θ 3 = 4 γ θ 0 - 3 γ + 1 ρ ∞ γ - 3 γ + 2

(24)

θ 2 = 2 γ (θ 0 + θ 3 - 1) - 2 θ 3 + 1 2 γ

(25)

θ 1 = 4 γ (1 - θ 0 - θ 3) + 2 θ 3 - 1 2 γ

(26)

θ 0 = 4 c 2 + c 1 2 (ρ ∞ + 1) c 3 4 c 3

(27)

其中：

c 1 = - 2 + 5 γ - 3 γ 2 - ρ ∞ γ + ρ ∞ γ 2

(28)

c 2 = (2 + 2 γ - 11 γ 2 + 3 γ 3) + 2 ρ ∞ (1 - 3 γ + 3 γ 2 + γ 3) + γ 2 ρ ∞ 2 (1 - γ)

(29)

c 3 = 8 (2 - 4 γ + γ 2 + ρ ∞ γ 2)

(30)

进一步，以最小化局部截断误差为优化目标确定最后一个算法参数γ。根据二阶精度的要求，TTBIF的局部截断误差表达式变为σ=s₀（γ，ρ_∞）。O（Δt³）+s₁（γ，ρ_∞）O（Δt⁴）+⋯+s_m（γ，ρ_∞）·O（Δt^m⁺³），其中s_i（i=0，1，2，3，⋯）为关于γ和ρ_∞的函数。为了最大化计算精度，这里利用∂s₀（γ，ρ_∞）/∂γ=0找到s₀的最小值，从而建立γ和ρ_∞的关系，即

γ = - 3241297 ρ ∞ 5 + 195493 ρ ∞ 4 - 15519434 ρ ∞ 3 + 50232711 ρ ∞ 2 - 102227201 ρ ∞ + 5771599

(31)

下面介绍共旋框架下基于TTBIF在一个时间步t→t+γΔt→t+2γΔt→t+Δt内求解式（1）的流程。

步骤1：假设x_t₊_γ_Δ_t= x_t，进一步利用式（21）预测τ=t+γΔt时刻的速度和加速度。

步骤2：计算τ=t+γΔt时刻动力学方程（1）的残余向量，即

G = [G 1 G 2] T

(32)

{G 1 = F I (q τ, q ˙ τ, q ¨ τ) + Φ q T (q τ) λ + F E (q τ) - Q (τ) G 2 = 4 γ 2 Δ t 2 Φ (q τ, τ)

(33)

第二个等式中，4/（γ²Δt²）为额外引入的系数。系数4/（γ²Δt²）的引入既不破坏原有动力学方程的平衡条件，同时可以抵消掉雅可比矩阵求解中由于∂q/∂

q ¨

=（γ²Δt²）/4引入的系数，从而保证雅可比矩阵的所有元素均为时间步长Δt的零次幂。

步骤3：计算τ=t+γΔt时刻修正量，即

(34)

其中：

J = [∂ G 1 ∂ q ¨ ∂ G 1 ∂ λ ∂ G 2 ∂ q ¨ ∂ G 2 ∂ λ]

(35)

∂ G 1 ∂ q ¨ = [M + γ 2 Δ t 2 4 (K + ∂ Φ q T ∂ q λ) + γ Δ t 2 C] | τ, ∂ G 1 ∂ λ = Φ q T | τ, ∂ G 2 ∂ q ¨ = Φ q | τ, ∂ G 2 ∂ λ = 0

(36)

步骤4：对x_t₊_γ_Δ_t进行修正，即

x t + γ Δ t = x t + γ Δ t + Δ x

(37)

将式（37）代入步骤2，若残量不满足精度要求，则对步骤2~4进行循环，否则进入下一个分步。

步骤5：假设x_t₊₂_γ_Δ_t= x_t₊_γ_Δ_t，进一步利用式（22）预测τ=t+2γΔt时刻的速度和加速度。

步骤6：计算τ=t+2γΔt时刻动力学方程（1）的残余向量，对应的表达式同式（32）和（33）。

步骤7：计算τ=t+2γΔt时刻修正量，对应的表达式同式（34）~（36）。

步骤8：对x_t₊₂_γ_Δ_t进行修正，即

x t + 2 γ Δ t = x t + 2 γ Δ t + Δ x

(38)

将式（38）代入步骤6，若残量不满足精度要求，则对步骤6~8进行循环，否则进入下一个分步。

步骤9：假设x_t_+Δ_t=x_t₊₂_γ_Δ_t，进一步利用式（23）预测τ=t+Δt时刻的速度和加速度。

步骤10：计算τ=t+Δt时刻动力学方程（1）的残余向量，其中G₂的表达式变为：

G 2 = 1 θ 32 Δ t 2 Φ (q τ, τ)

(39)

步骤11：计算τ=t+Δt时刻修正量，其中

∂ G 1

∂ q ¨

的表达式变为：

∂ G 1 ∂ q ¨ = [M + θ 32 Δ t 2 (K + ∂ Φ q T ∂ q λ) + θ 3 Δ t C] | τ

(40)

步骤12：对x_t_+Δ_t进行修正，即

x t + Δ t = x t + Δ t + Δ x

(41)

将式（41）代入步骤10，若残量不满足精度要求，则对步骤10~12进行循环，否则进入下一个时间步。

从式（33）和（39）中可以看到，系数4/（γ²∆t²）和1/（

θ 32

∆t²）的引入可以确保雅可比矩阵的所有元素都是O（Δt⁰），从而避免了由于时间步长取值过小导致的雅可比矩阵奇异性问题。此外，为了方便读者执行本文提出的策略，图3提供了Co‑TTBIF的计算流程。

3　数值算例

收起

本节利用2个算例来验证本文策略（记为Co‑TTBIF）在柔性多体系统动力学仿真中的优势。作为对比，基于绝对节点坐标单元和共旋坐标单元与Generalized‑α方法的组合求解策略也被考虑在本节，分别记为ANCF‑Gα和Co‑Gα。考虑到TTBIF在一个时间步长内包含3个分步，而Generalized‑α方法是单分步的。因此，在所有算例中，本文策略的时间步长取为ANCF‑Gα和Co‑Gα的3倍。此外，为了过滤柔性体空间离散引入的虚假振荡，所有策略采用ρ_∞=0。

3.1　自由下落单摆模型

首先考虑受重力作用下自由下落的柔性单摆模型^[29]，如图4所示。物理参数假设为：摆长L=1.2 m、横截面积A=a×b=0.015×0.01=0.00015 m²、杨氏模量E=10 MPa、泊松比ν=0.3、密度ρ_o=5540 kg/m³。在本算例中，考虑3种求解策略，分别是：（1）4个考虑剪切变形的二维ANCF梁单元（12个自由度）与Generalized‑α方法；（2）4个考虑剪切变形的二维共旋梁单元（6个自由度）与Generalized‑α方法；（3）4个考虑剪切变形的二维共旋梁单元与TTBIF。假设∆t（ANCF‑Gα）=10^-4 s、∆t（Co‑Gα）=10^-4 s和∆t（Co‑TTBIF）=3×10^-4 s。参考解由取更小∆t的ANCF‑Gα提供。

对于该保守系统，在当前初始条件下，系统的总能量理论上恒为0。图5~7提供了该模型的能量和位移结果。从图5中可以看到，ANCF‑Gα、Co‑Gα和Co‑TTBIF都较好地保守系统的能量。进一步，从图6和7中可以看到：对于具有大转动的柔性多体系统，ANCF‑Gα、Co‑Gα和Co‑TTBIF对动力学响应都给出了准确的预测。

为了展示Co‑TTBIF在柔性多体系统动力学仿真中的优势，图8对ANCF‑Gα、Co‑Gα和Co‑TTBIF的Newton迭代次数和CPU时间进行了比较。从中发现：（1）相同单元数下，共旋单元用更少的自由度数实现了与ANCF单元相同的精度；（2）因为TTBIF的动力学平衡方程严格满足时间离散点，而Generalized‑α方法的动力学平衡方程只在广义的时间点满足，因此TTBIF更容易收敛，Newton迭代次数更少，计算量更小；（3）相同精度下，Co‑TTBIF具有显著的效率优势，相比当前流行的ANCF‑Gα，大约可以节约95%的计算量。

3.2　受迫激励下的柔性悬臂梁模型

第二个算例考虑一个集中力作用下的柔性悬臂梁模型^[30]，如图9所示，物理参数为：梁长L₀=10 m、横截面积A=a×b=0.25×0.5=0.125 m²、杨氏模量E=210 GPa、密度ρ=7850 kg/m³。在悬臂梁的自由端受到正弦激励P=P₀sin（ωt），其中P₀=10 MN、ω=50 rad/s。该悬臂梁模型利用20个共旋梁单元进行空间离散。Co‑Gα和Co‑TTBIF的时间步长分别取为10^-4 s和3×10^-4 s，参考解由取更小∆t的Co‑Gα提供。

图10绘制了Co‑Gα和Co‑TTBIF在自由端的位移随时间的演化规律，Ref表示参考解。图中给出了两种求解策略的平均绝对误差（average absolute errors，AAE）。可以发现：这两种求解策略都能较好地预测动力学响应，但Co‑TTBIF的精度更高一些。图11比较了Co‑Gα和Co‑TTBIF在自由端的速度随时间的演化规律，可以发现由于柔性体的空间离散导致速度中存在数值振荡，特别是在转角的计算中。随着计算时间增加，相比Co‑TTBIF，Co‑Gα的数值振荡更加剧烈。

进一步，为了讨论人工阻尼量对柔性多体系统计算结果的影响，图12对比了Co‑TTBIF和Co‑Gα在不同ρ_∞取值下自由端速度的结果。数值结果显示：随着耗散程度的逐渐降低，Co‑TTBIF和Co‑Gα的数值震荡越来越剧烈。因此，在柔性连续体的动力学仿真中，一定程度的数值耗散量是必备的，其可以有效地过滤掉不需要的虚假信息。

为了进一步比较Co‑Gα和Co‑TTBIF的效率，图13给出了两种求解策略的Newton迭代次数和CPU时间。结果表明：在Co‑TTBIF比Co‑Gα略有精度优势时，其计算量仅为Co‑Gα的1/8。Co‑TTBIF比Co‑Gα针对柔性多体采用了完全相同的建模方式，但由于TTBIF的动力学平衡方程严格满足时间离散点，导致其Newton迭代次数更少，计算量显著低于Generalized‑α方法。

综上所述，本节通过对两个具有大旋转和大变形的柔性多体系统的动力学分析发现：与当前柔性多体系统动力学分析领域广泛使用的求解策略ANCF‑Gα相比，本文提出的Co‑TTBIF在求解同类动力学系统时具有可观的精度和效率优势。

4　结论

收起

本文针对柔性多体系统提出了一种共旋框架下基于TTBIF的动力学分析策略，记为Co‑TTBIF。在本文策略中，柔性体的动力学模型被建立在共旋坐标系下，因此任意先进的几何线性单元均可以被使用，进而降低动力学方程的非线性程度。随后，利用优化的复合方法TTBIF对建立的动力学模型进行时域分析，在有效地过滤由于柔性体空间离散引入的数值振荡的同时可以高精度地计算重要的低频响应。与已有工作相比，本文的创新之处在于利用共旋理论在保证动力学模型精度的前提下显著降低模型自由度，同时利用复合型方法求解动力学模型，在保证计算精度和人工可控阻尼的前提下有效降低迭代次数，从而实现柔性多体系统动力学响应的精确和快速仿真。

利用两个经典的柔性多体动力学模型分析了本文策略Co‑TTBIF在处理这类具有强非线性和大范围运动的动力学系统的表现。计算结果显示：与当前流行的绝对坐标框架下基于Generalized‑α方法的动力学求解策略（ANCF‑Gα）相比，在精度相近的前提下，本文策略具有显著的效率优势。

基金

收起

国家自然科学基金资助项目(12202058; 12302044)
中国博士后科学基金资助项目(2023T160051; 2022M710386)
北京市科协青年人才托举工程项目(BYESS2023344)

参考文献

收起

文献

收起

参考文献引证文献

排序方式：

[1]

杜超凡，周晓婷，章定国，等.无网格法在中心刚体-旋转柔性梁系统动力学分析中的应用[J].振动工程学报，2022,35(1):178-187.

Chaofan

，ZHOU

Xiaoting

，ZHANG

Dingguo

，et al. Meshless methods applied to the dynamic analysis for hub-beam system[J]. Journal of Vibration Engineering，2022,35(1): 178-187.

[2]

孙加亮，孙佳昊，金栋平，等.充气自旋结构的动力学特性研究[J].振动工程学报，2023,36(3):662-670.

SUN

Jialiang

，SUN

Jiahao

，JIN

Dongping

，et al. Study on dynamic characteristics of inflatable spin structure[J]. Journal of Vibration Engineering，2023,36(3): 662-670.

[3]

DE VEUBEKE

B F

. The dynamics of flexible bodies[J]. International Journal of Engineering Science，1976,14(10): 895-913.

[4]

SIMO

J C

. A finite strain beam formulation. The three-dimensional dynamic problem. Part Ⅰ[J]. Computer Methods in Applied Mechanics and Engineering，1985,49(1): 55-70.

[5]

MOURID

，MAMOURI

，KOULI

，et al. An efficient implicit time stepping scheme for nonlinear dynamics analysis of 3D geometrically exact beam using exponential mapping[J]. Mechanics of Advanced Materials and Structures，2022,29(28): 7189-7203.

[6]

ZHANG

，LIU

，TANG

H Y

. Geometrically exact shell with drilling rotations formulated on the special Euclidean group SE（3）[J]. International Journal for Numerical Methods in Engineering，2021,122(18): 4886-4921.

[7]

SHABANA

A A

. Flexible multibody dynamics： review of past and recent developments[J]. Multibody System Dynamics，1997,1(2): 189-222.

[8]

WEMPNER

. Finite elements，finite rotations and small strains of flexible shells[J]. International Journal of Solids and Structures，1969,5(2): 117-153.

[9]

BELYTSCHKO

，SCHWER

，KLEIN

M J

. Large displacement，transient analysis of space frames[J]. International Journal for Numerical Methods in Engineering，1977,11(1): 65-84.

[10]

WANG

，QI

Z H

，XU

J S

. A high-precision co-rotational formulation of 3D beam elements for dynamic analysis of flexible multibody systems[J]. Computer Methods in Applied Mechanics and Engineering，2020,360: 112701.

[11]

SHI

J B

，LIU

Z Y

，HONG

J Z

. Multibody dynamic analysis using a rotation-free shell element with corotational frame[J]. Acta Mechanica Sinica，2018,34(4): 769-780.

[12]

邢誉峰，张慧敏，季奕.动力学常微分方程的时间积分方法[M].北京：科学出版社，2022.

[13]

HUGHES

T J R

. The Finite Element Method： Linear Static and Dynamic Finite Element Analysis[M]. New Jersey： Prentice Hall Press，1987.

[14]

CHUNG

，HULBERT

G M

. A time integration algorithm for structural dynamics with improved numerical dissipation： the generalized-α method[J]. Journal of Applied Mechanics，1993,60(2): 371-375.

[15]

ZHANG

H M

，XING

Y F

. A three-parameter single-step time integration method for structural dynamic analysis[J]. Acta Mechanica Sinica，2019,35(1): 112-128.

[16]

ZHANG

H M

，XING

Y F

，JI

. An energy-conserving and decaying time integration method for general nonlinear dynamics[J]. International Journal for Numerical Methods in Engineering，2020,121(5): 925-944.

[17]

LIU

C Y

，LI

J Y

，YANG

Z Q

，et al. Two high-order energy-preserving and symmetric Gauss collocation integrators for solving the hyperbolic Hamiltonian systems [J]. Mathematics and Computers in Simulation，2023,205: 19-32.

[18]

KIM

，REDDY

J N

. A new family of higher-order time integration algorithms for the analysis of structural dynamics[J]. Journal of Applied Mechanics，2017,84(7): 071008.

[19]

FUNG

T C

. Unconditionally stable collocation algorithms for second order initial value problems[J]. Journal of Sound and Vibration，2001,247(2): 343-365.

[20]

邵慧萍，蔡承文. 结构动力学方程数值积分的三参数算法[J]. 应用力学学报，1988,5(4): 76-81.

SHAO

Huiping

，CAI

Chengwen

. A three parameters algorithm for numerical integration of structural dynamic equations[J]. Chinese Journal of Applied Mechanics，1988,5(4): 76-81.

[21]

BRÜLS

，CARDONA

，ARNOLD

. Lie group generalized-α time integration of constrained flexible multibody systems[J]. Mechanism and Machine Theory，2012,48: 121-137.

[22]

BANK

R E

，COUGHRAN

W M

，FICHTNER

，et al. Transient simulation of silicon devices and circuits[J]. IEEE Transactions on Computer-Aided Design of Integrated Circuits and Systems，1985,4(4): 436-451.

[23]

BATHE

K J

，BAIG

M M I

. On a composite implicit time integration procedure for nonlinear dynamics[J]. Computers & Structures，2005,83(31-32): 2513-2524.

[24]

WEN

W B

，LEI

，WANG

，et al. An efficient hybrid implicit time integration method with high accuracy for linear and nonlinear dynamics[J]. European Journal of Mechanics-A，Solids，2023,97: 104811.

[25]

KIM

. An improved implicit method with dissipation control capability： the simple generalized composite time integration algorithm[J]. Applied Mathematical Modelling，2020,81: 910-930.

[26]

，XING

Y F

. An optimized three-sub-step composite time integration method with controllable numerical dissipation[J]. Computers & Structures，2020,231: 106210.

[27]

LIU

T H

，HUANG

F L

，WEN

W B

，et al. Further insights of a composite implicit time integration scheme and its performance on linear seismic response analysis[J]. Engineering Structures，2021,241: 112490.

[28]

，XING

Y F

. A three-sub-step composite method for the analysis of rigid body rotations with Euler parameters[J]. Nonlinear Dynamics，2023,111(15): 14309-14333.

[29]

BERZERI

，SHABANA

A A

. Development of simple models for the elastic forces in the absolute nodal co-ordinate formulation[J]. Journal of Sound and Vibration，2000,235(4): 539-565.

[30]

T N

，BATTINI

J M

，HJIAJ

. Efficient formulation for dynamics of corotational 2D beams[J]. Computational Mechanics，2011,48(2): 153-161.

2025年第38卷第10期

PDF下载

引用本文

BibTeX

文章信息

doi: 10.16385/j.cnki.issn.1004-4523.202310039

接收时间：2023-10-18
首发时间：2026-02-04

补充材料

相关文章

文章信息

作者

出版历史

收稿日期：2023-10-18
修回日期：2024-01-19

基金

国家自然科学基金资助项目(12202058; 12302044)

中国博士后科学基金资助项目(2023T160051; 2022M710386)

北京市科协青年人才托举工程项目(BYESS2023344)

作者信息

^1.哈尔滨工业大学自主智能无人系统工信部重点实验室，黑龙江　哈尔滨　150001；

^2.北京理工大学飞行器动力学与控制教育部重点实验室，北京　100081；

^3.中物院高性能数值模拟软件中心，北京　100088；

^4.北京应用物理与计算数学研究所，北京　100088；

^5.北京宇航系统工程研究所，北京　100005

通讯作者:

季奕（1995—），女，博士，副教授。E-mail： jiyi0319@outlook.com

参考文献

分享链接

https://castjournals.cast.org.cn/joweb/zdgcxb/CN/10.16385/j.cnki.issn.1004-4523.202310039

分享至

全文二维码

扫描看全文

引用本文

BibTeX

本文的引用情况

2种不同金属材料的力学参数

科 Family	属数 Number of genus	种数 Number of species	占总种数比例 Percentage of total species (%)	属 Genus	种数 Number of species	占总种数比例 Percentage of total species (%)
鹅膏菌科Amanitaceae	2	11	5.26	鹅膏菌属 Amanita	10	4.78
小菇科 Mycenaceae	2	12	5.74	丝盖伞属 Inocybe	5	2.39
多孔菌科 Polyporaceae	8	14	6.70	蜡蘑属 Laccaria	5	2.39
红菇科 Russulaceae	3	23	11.00	小皮伞属 Marasmius	6	2.87
				小菇属 Mycena	11	5.26
				光柄菇属 Pluteus	5	2.39
				红菇属 Russula	17	8.13
				栓菌属 Trametes	5	2.39

关闭全屏

BibTeX
EndNote
RefWorks
TxT