Page 325 - 《振动工程学报》2025年第11期

第 11 期 李金泽，等：结构动响应的自启动单解时域积分器优化 2783

文献 [2,12] 构 造 了 自 启 动 单 解 时 域 积 分 器 框 p α 1 α 2 p−α 4 α 4
1 1 1 1 1
架。首先，利用 t n 时刻的数值解和 t n+p 时刻的二阶微 1− 1−
2 −α 6 α 6 2p 2p p p
分方程求解得到 t n+p 时刻的加速度 ¨ u n+p ： （8）
（5a） 数值放大矩阵的不变量 A i (i = 0, 1, 2, 3)在数值
M¨u n+p +C˙u n+p + Ku n+p = F n+p
˙ u n+p = ˙u n +∆t(α 3 ¨u n +α 4 ¨u n+p ) （5b） 高频极限 ( Ω := ω·∆t → ∞，其中， ω为单自由度系统
的固有频率) 下化简为：
2
u n+p = u n + p∆t˙u n +∆t (α 1 ¨u n +α 2 ¨u n+p ) （5c） 
∞
A = 1
3

式 中， p ∈ R和 α j ∈ R ( j = 1, 2, 3, 4)为 待 定 的 算 法 参 



 2pα 6 −6pα 2 + p+2α 1 +2α 2
 ∞
 A =

数；Δt 为积分步长。然后，更新 t n+1 时刻的位移、速度  2


 2pα 2


和加速度：  3pα 2 −2pα 6 +1−2α 1 −2α 2 （9）

A =
∞

 1

2
u n+1 = u n +∆t˙u n +∆t (α 5 ¨u n +α 6 ¨u n+p ) （5d）  pα 2




 2α 1 +2α 2 +2pα 6 − p−2pα 2
 ∞
 A =

˙ u n+1 = ˙u n +∆t(α 7 ¨u n +α 8 ¨u n+p ) （5e）  0
2pα 2
（5f） 需要注意的是，高频极限通常是指数值频率 Ω → ∞，
¨ u n+1 = α 9 ¨u n +α 10 ¨u n+p
式中， α j ∈ R ( j = 5, 6, ··· , 10)为待定的算法参数。式 而非物理频率的极限 ω → ∞。隐式算法实现最优可
(5) 构成了自启动单解时域积分器框架且可以利用 控数值高频耗散的条件 [2,12] 为：
扩展的 Butcher 表描述为： A = 3ρ ∞ , A = 3ρ , A = ρ 3 ∞ （10）
2
∞
∞
∞
2
1
0
∞
p α 1 α 2 α 3 α 4 将式 (9) 代入式 (10) 并求解可得：
（6）

α 5 α 6 α 7 α 8 α 9 α 10  2−ρ ∞
 p =

每个自启动单解时域积分器都对应一个具体的  1+ρ ∞





 2 3 2
Butcher 表。  2(1+ρ ∞ )(2−ρ ∞ ) α 6 +ρ −3ρ +6ρ ∞ −6

∞
∞

 α 1 = （11）
2
 2(1+ρ ∞ ) (ρ ∞ −2)

文献 [2,12] 已经对该积分算法框架进行了精度 



 1

分析。该框架实现位移、速度和加速度的一致二阶 
α 2 =

2

(1+ρ ∞ ) (2−ρ ∞ )
精度条件为：
式中， ρ ∞ ∈ [0, 1]为用户指定参数，用于控制隐式算法
α 3 = p−α 4 , α 5 = 1/2−α 6 , α 7 = 1−1/(2p),
在高频区域内的谱半径值。
α 8 = 1/(2p),α 9 = 1−1/p,α 10 = 1/p （7）
目前，隐式算法的 α 4 和 α 6 是未知的。文献 [12]
参数 α 2 和 α 4 决定了时域积分器的显式和隐式属 通过优化位移局部截断误差得到了最优的 α 4 和 α 6 参
性。当 α 2 , 0且 α 4 , 0时，框架对应了一致二阶精度 数，但 OSS21 算法存在严重的位移超调。为了降低
的隐式算法；当 α 2 = 0且 α 4 , 0时，框架对应了一致二 超调行为，本文通过最小化振幅与相位误差确定最
阶 精 度 的 速 度 隐 式 处 理 的 显 式 算 法； 当 α 2 = α 4 = 0 优的参数 α 4 ，通过降低超调响应确定最优的参数 α 6 。
时，框架对应了一致二阶精度的完全显式算法。 利用解析计算振幅与相位误差主项的一般性方

法 [19] ，本文隐式算法的振幅 δ与相位 ε误差可计算为：
2 隐 式 算 法 的 优 化  2 3
δ = δ 2 ∆t +O(∆t )


 2 3 （12）
 ε = ε 2 ∆t +O(∆t )

3
当 α 2 , 0且 α 4 , 0时，式 (6) 对应了一致二阶精度 式中，O(Δt ) 表示所有三阶及其更高阶项的总和；误
的隐式算法。此时，其 Butcher 表对应修改为： 差主项 δ 2 和 ε 2 分别为：
2
2(1+ρ ∞ )(12ρ ∞ α 4 −24α 4 −7ρ ∞ +11)×ξ −3ρ ∞ (ρ ∞ −1)(2α 4 −3)+12α 4 3
δ 2 = ξω （13）
6(1+ρ ∞ ) 2
2
2
2
4
2
8(1+ρ ∞ )(12ρ ∞ α 4 −24α 4 −7ρ ∞ +11)×ξ +(64ρ −52ρ ∞ −44)ξ +14ρ ∞ −72(ρ ∞ +1)(ρ ∞ −2)α 4 ξ −11ρ −11
ε 2 = ∞ √ ∞ ω 3
24(1+ρ ∞ ) 2 1−ξ 2
（14）
式中， ξ为单自由度系统的物理阻尼率。 求解线性结构振动时的无条件稳定。
定义误差总量 Γ = δ +ε 并依据文献 [12] 中最小 为了分析超调响应，计算在极限 ∆t → ∞时第一
2
2
2
2
个时间步的数值位移 u 为：
∞
化 Γ的技术得到最优的 α 4 ： 1
2
∞
2
ρ +2ρ ∞ −11 u ≈ c 2 ¨u 0 ∆t +c 1 ∆t +c 0 （16）
1
α 4 = ∞ （15）
8(ρ ∞ +1)(ρ ∞ −2) 式中， ¨ u 0 为单自由度系统的初始加速度； c i (i = 0, 1, 2)
需要指出的是， α 4 的上述取值也保证了隐式算法在 为算法和结构参数的表达式且 c 2 主导了位移的超调