欢迎访问《哈尔滨工业大学学报》编辑部网站！

一种具有可控数值阻尼的无条件稳定半显式积分算法

doi: 10.11918/202311063

长安大学建筑工程学院,西安 710061

基金项目: 国家自然科学基金(51908048,52108432,52478124) ；陕西省高校科协青年人才托举计划(20230406) ；长安大学中央高校基本科研业务费专项资金(300102283201) ；长安大学青年学者学科交叉团队建设项目(300104240923)

详细信息

作者简介

傅博(1990—),男,博士,副教授

通讯作者

陈瑾,chenjin5310@126.com

中图分类号: TU311.3

文献标识码: A

文章编号: 0367-6234(2025)02-0113-09

An unconditionally stable semi-explicit integration algorithm with controllable numerical damping

FU Bo ， ZHANG Futai ， ZHANG Qingkai ， CHEN Jin

School of Civil Engineering,Chang′an University, Xi′an 710061 , China

摘要

结构模型经过有限元方法空间离散化处理之后,可能引入虚假的高频分量,这部分高频分量会对结构的动力响应求解带来不利影响。为此,需要引入算法的数值阻尼有效地抑制这部分虚假的高频分量。使用半显式算法格式,通过匹配隐式 ρ_∞ -Bathe 算法放大矩阵的特征方程系数,提出一种具有可控数值阻尼的无条件稳定半显式积分算法,记为 NSE(New Semi-Explicit)- ρ_∞ 算法,新算法通过两个自由参数 ρ_∞ 和 γ 控制算法的数值阻尼,并且无需对结构运动方程进行加权处理。对新算法的稳定性、精度、周期延长和振幅衰减等数值特性进行分析,结果表明,新算法对于线弹性体系和非线性刚度软化体系均为无条件稳定。通过具有代表性的数值算例,将新算法与两种具有可控数值阻尼的无条件稳定显式积分算法进行对比,证明新算法能够更加有效地抑制虚假高频分量。

关键词

积分算法 / 显式 / 稳定 / 数值阻尼 / 结构动力学

Abstract

After the structural model is discretized by the finite element method, it may introduce false highfrequency components that can adversely affect the structural dynamic response of the structure. Therefore, it is necessary to introduce numerical damping into the integration algorithm to effectively suppress these false highfrequency componments. Based on a semi-explicit integration formulation, this paper proposes an unconditionally stable semi-explicit integration algorithm with controllable numerical damping by matching the characteristic equation coefficients of the amplification matrix of the implicit ρ_∞ -Bathe algorithm. The new semi-explicit (NSE)- ρ_∞ algorithm controls the numerical damping of the integration algorithm by two free coefficients ρ_∞ and γ, and does not require weighted equation of motion. The numerical characteristics of the new algorithm, such as stability, accuracy, period elongation and amplitude decay, are analyzed. It is found that the new algorithm is unconditionally stable for both linear elastic and nonlinear stiffness softening systems. Through representative numerical examples, the new algorithm is compared with two other unconditionally stable explicit integration algorithms with controllable numerical damping. The analytical results demonstrate that the new algorithm is more effective in suppressing the spurious high-frequency components.

Keywords

integration algorithm / explicit / stability / numerical damping / structural dynamics

1 新算法的建立 1.1 ρ_∞-Bathe 算法 1.2 新算法积分参数推导（SDOF 体系） 1.3 新算法积分参数推导（MDOF 体系） 2 数值特性分析 2.1 稳定性分析 2.2 精度分析 2.3 周期延长和振幅衰减 3 数值算例 4 结论

目前，在结构动力学响应的数值计算问题上，比较有效的分析方法仍然是时程分析法^[1]，即需要求解出整个动态过程中结构的位移、速度和加速度。在分析计算过程中，通常采用积分算法^[2] 求解离散化的结构运动方程。经典的积分算法如 Newmark 算法具有较好的精度和适当的计算量，被广泛用于求解结构动力学问题。然而，这些经典算法无法满足求解一些复杂结构动力学问题的需求。以实时混合试验^[3-5]为例，当数值子结构为自由度数较多的非线性结构时，如果采用常规的无条件稳定算法（如常平均加速度法），由于算法是隐式，在求解非线性问题时需要迭代，极大增加了求解时间; 如果采用常规的显式算法（如中心差分法），由于算法为条件稳定，需要极小的时间步长以满足稳定性，同样严重影响计算效率，无法满足实时混合试验的实时性要求。

为了解决上述问题，文献^[6] 最早提出了同时满足无条件稳定和显式的积分算法，称为 Chang 算法，该算法采用了常平均加速度法的速度表达式，但是位移表达式为显式表达式。由于该算法的积分参数是与结构模型固有属性（质量、阻尼和刚度）相关的变量，文献^[7]又将无条件稳定显式算法称为“基于模型的积分算法”。其他有代表性的基于模型的积分算法包括文献^[8]提出的 CR 算法、文献^[9]提出的 TL 算法。与 Chang 算法不同的是，CR 算法和 TL 算法的速度和位移表达式均为显式，文献^[7]将 Chang 算法和 CR/ TL 算法分别称为半显式和双显式算法。文献^[10]在 CR 算法的基础上，提出了一族 Gui 算法，与 γ = 1 / 2 的 Newmark 算法有相同的数值特性。为进一步提高算法的精度，文献^[11-12] 分别提出一种高阶精度双显式和半显式的基于模型的积分算法。

需要注意的是，上述基于模型的积分算法均不具有数值阻尼，算法的数值阻尼与结构的阻尼功能类似，可以减小结构响应。数值阻尼通常与 Ω = ωΔt（ω 和 Δt 分别为圆频率和时间积分步长）正相关，当算法用于高频时，ω 与 Ω 较大，导致数值阻尼较大，进而算法可以抑制或消除虚假高阶/ 高频响应。其抑制效果与数值阻尼的大小直接相关，数值阻尼也会影响低阶响应的精度，所以，与无数值阻尼的算法相比，有数值阻尼的算法优势在于可以抑制或消除虚假高频分量的影响，缺点在于低阶响应的精度相对较低。因此，设计有数值阻尼算法时，希望低频具有尽量小的数值阻尼，而高频具有尽量大的数值阻尼，即可控的数值阻尼。文献^[13] 采用 CR 算法的表达式，通过匹配广义 α（G-α）算法^[14] 的数值特性，提出一族具有可控数值阻尼的 KR-α 双显式算法。类似地，文献^[15]提出了一族具有可控数值阻尼的 Chang-α 半显式算法。 KR-α 算法和 Chang-α 算法均是引入了额外的加权结构运动方程，应用起来相对复杂，并且数值阻尼仅为单参数控制，不够灵活。文献^[16]提出一种与 Newmark 算法数值特性匹配的 GCR 双显式算法，与 Newmark 算法一样，GCR 算法在 γ >1 / 2 且 β≥（γ + 0.5）² / 4 具有数值阻尼，不过此时算法精度降为 1 阶。文献^[17] 在 TL 算法的基础上，提出一种具有可控数值阻尼的 TL-φ 算法，该算法同样存在精度较低的问题。此外，该算法对于无阻尼结构体系并不具有数值阻尼。文献^[18]基于隐式 Bathe 算法^[19]，提出一种具有数值阻尼的无条件稳定显式算法，但是该算法的数值阻尼不可控。

综上，本文采用半显式 Chang 算法的表达式，基于具有可控数值阻尼的隐式 ρ_∞-Bathe 算法^[20]，利用极点配置理论，提出一种具有可控数值阻尼的无条件稳定半显式的 NSE-ρ_∞ 算法，与 KR-α 算法和 Chang-α 算法相比，新算法无需引入额外的加权结构运动方程，使用更加方便，具有更加灵活且优异的可控数值特性。

1 新算法的建立

1.1 ρ_∞-Bathe 算法

隐式 ρ_∞-Bathe 算法具有二阶精度，无条件稳定，且具有可控数值阻尼。需使用两个子步计算 t + Δt 时刻的未知位移、速度和加速度。

在第 1 个子步中，使用梯形规则计算 t + γΔt 时刻的结构响应，格式如下:

m {\ddot{x}}_{i + γ} + c {\dot{x}}_{i + γ} + k x_{i + γ} = f_{i + γ}

(1)

x_{i + γ} = x_{i} + \frac{γ Δ t}{2} ({\dot{x}}_{i} + {\dot{x}}_{i + γ})

(2)

{\dot{x}}_{i + γ} = {\dot{x}}_{i} + \frac{γ Δ t}{2} ({\ddot{x}}_{i} + {\ddot{x}}_{i + γ})

(3)

在第 2 个子步中，使用 3 点后向欧拉公式计算 t + Δt时刻的结构响应，格式如下:

m {\ddot{x}}_{i + 1} + c {\dot{x}}_{i + 1} + k x_{i + 1} = f_{i + 1}

(4)

x_{i + 1} = x_{i} + Δ t (q_{0} {\dot{x}}_{i} + q_{1} {\dot{x}}_{i + γ} + q_{2} {\dot{x}}_{i + 1})

(5)

{\dot{x}}_{i + 1} = {\dot{x}}_{i} + Δ t (s_{0} {\ddot{x}}_{i} + s_{1} {\ddot{x}}_{i + γ} + s_{2} {\ddot{x}}_{i + 1})

(6)

式中:m，c，k 分别为质量、阻尼系数和刚度; c = 2mξω，ξ 为阻尼比，ω 为圆频率;

\ddot{x}

，

\dot{x}

，x 分别为加速度、速度和位移; f 为外界激励; q₀，q₁，q₂，s₀，s₁，s₂ 和 γ 为控制算法性能的待定参数。为了使 ρ_∞-Bathe 算法具有二阶精度，这些待定参数需要满足下列关系:

\{\begin{matrix} q_{0} = (γ - 1) q_{1} + \frac{1}{2} \\ q_{2} = - γ q_{1} + \frac{1}{2} \\ s_{0} = (γ - 1) s_{1} + \frac{1}{2} \\ s_{2} = - γ s_{1} + \frac{1}{2} \end{matrix}

(7)

为了在二阶精度的基础上满足无条件稳定的特性，待定参数还需要满足下列关系:

s_{1} = q_{1}

(8)

为了直接规定高频范围内的数值阻尼，建立了 q₁ 与高频谱半径 ρ_∞之间的关系:

q_{1} = \frac{ρ_{\infty} + 1}{2 γ (ρ_{\infty} - 1) + 4}, ρ_{\infty} \in [0,1]

(9)

至此，ρ_∞-Bathe 算法实际上是关于两个自由参数 γ 和 ρ_∞的两子步无条件稳定隐式算法。此外，γ 通常选用范围是（0，2），且应避免取 1 和 2 /（1-ρ_∞）。

1.2 新算法积分参数推导（SDOF 体系）

新算法（NSE-ρ_∞算法）采用 Chang 算法的表达式，对于单自由度（ SDOF）体系，其速度、位移差分方程分别如下:

{\dot{x}}_{i + 1} = {\dot{x}}_{i} + \frac{1}{2} Δ t {\ddot{x}}_{i} + \frac{1}{2} Δ t {\ddot{x}}_{i + 1}

(10)

x_{i + 1} = x_{i} + α_{1} Δ t {\dot{x}}_{i} + α_{2} Δ t^{2} {\ddot{x}}_{i}

(11)

将式（4）、（10）和（11）写成递推矩阵的形式:

[\begin{matrix} x_{i + 1} \\ Δ t {\dot{x}}_{i + 1} \\ Δ t^{2} {\ddot{x}}_{i + 1} \end{matrix}] = A [\begin{matrix} x_{i} \\ Δ t {\dot{x}}_{i} \\ Δ t^{2} {\ddot{x}}_{i} \end{matrix}] + [\begin{matrix} 0 \\ 0 \\ Δ t^{2} / m \end{matrix}] f_{i + 1}

(12)

式中

A

为放大矩阵，表达式如下:

A = [\begin{matrix} 1 & α_{1} & α_{2} \\ - \frac{Ω^{2}}{2 ξ Ω + 2} & - \frac{α_{1} Ω^{2} - 2}{2 ξ Ω + 2} & - \frac{α_{2} Ω^{2} - 1}{2 ξ Ω + 2} \\ - \frac{Ω^{2}}{ξ Ω + 1} & - \frac{α_{1} Ω^{2} + 2 ξ Ω}{ξ Ω + 1} & - \frac{α_{2} Ω^{2} + ξ Ω}{ξ Ω + 1} \end{matrix}]

(13)

对于 3 × 3 的矩阵，其特征方程可以表示为

λ^{3} - 2 A_{1} λ^{2} + A_{2} λ - A_{3} = 0

(14)

式中:λ 为放大矩阵

A

的特征值，根据 |

A

-λI | = 0 可以求得 NSE-ρ_∞ 算法的特征方程系数为

\{\begin{matrix} A_{1} = (4 - α_{1} Ω^{2} - 2 α_{2} Ω^{2}) / 4 (ξ Ω + 1) \\ A_{2} = (α_{1} Ω^{2} - 2 α_{2} Ω^{2} - 2 ξ Ω + 2) / 2 (ξ Ω + 1) \\ A_{3} = 0 \end{matrix}

(15)

令 NSE-ρ_∞算法与 ρ_∞-Bathe 算法的特征方程系数完全相等，即

A_{i}^{NSE - ρ_{\infty}} = A_{i}^{ρ_{\infty} - Bathe}, i = 1,2, 3

(16)

通过上式可以求解得到 NSE-ρ_∞ 算法在 SDOF 体系下的积分参数。具体表达形式如下:

α_{1} = \frac{c_{7} ξ Ω^{3} + (c_{8} ξ^{2} + c_{9}) Ω^{2} + c_{10} ξ Ω + c_{11}}{c_{1} Ω^{4} + c_{2} ξ Ω^{3} + (c_{3} ξ^{2} + c_{4}) Ω^{2} + c_{5} ξ Ω + c_{6}}

(17)

α_{2} = \frac{c_{12} ξ Ω^{3} + (c_{13} ξ^{2} + c_{14}) Ω^{2} + c_{15} ξ Ω + c_{16}}{c_{1} Ω^{4} + c_{2} ξ Ω^{3} + (c_{3} ξ^{2} + c_{4}) Ω^{2} + c_{5} ξ Ω + c_{6}}

(18)

式中常系数 c_i（ i = 1~16）取值如下:c₁ =（2q₁γ-1）² γ²，c₂ = 4γ [（2q₁γ-1）γ-1]（2q₁γ-1），c₃ =-16（2q₁γ-1）γ，c₄ =（16

q_{1}^{2}

+ 4）γ²-16q₁γ + 4，c₅ = 16-16（2q₁-1）γ，c₆ = 16，c₇ = 4（2q₁-1）² γ ²，c₈ =-16（ 2q₁-1）γ，c₉ = 4（ 2q₁-1）² γ²，c₁₀ = 16-16（2q₁-1）γ，c₁₁ = 16，c₁₂ =-8

q_{1}^{2}

（ γ-1）²γ²，c₁₃ =-16（2q₁-1）q₁（γ-1）γ，c₁₄ = 2（2q₁-1）（4q₁γ-2q₁-1）γ²，c₁₅ = 64q₁γ（γ-1）ξ² + 8q₁γ（-6γ + 4）+ 8γ，c₁₆ = 8。

1.3 新算法积分参数推导（MDOF 体系）

对于线弹性多自由度（MDOF）体系，NSE-ρ_∞ 算法的速度、位移差分方程和运动方程如下:

{\dot{X}}_{i + 1} = {\dot{X}}_{i} + \frac{1}{2} Δ t {\ddot{X}}_{i} + \frac{1}{2} Δ t {\ddot{X}}_{i + 1}

(19)

X_{i + 1} = X_{i} + α_{1} Δ t {\dot{X}}_{i} + α_{2} Δ t^{2} {\ddot{X}}_{i}

(20)

M {\ddot{X}}_{i + 1} + C {\dot{X}}_{i + 1} + K X_{i + 1} = F_{i + 1}

(21)

令 Y 为模态坐标系的位移向量，与 X 的关系为

X = Φ Y ， Φ = [ϕ_{1} ， ϕ_{2} ， \dots ， ϕ_{n}]

为模态矩阵;

ϕ_{1}

（ i = 1，…，n）为第 i 阶振型向量，利用振型的正交性，则式（19）~（21）可以改写至模态坐标系中:

{\dot{Y}}_{i + 1} = {\dot{Y}}_{i} + \frac{1}{2} Δ t {\ddot{Y}}_{i} + \frac{1}{2} Δ t {\ddot{Y}}_{i + 1}

(22)

Y_{i + 1} = Y_{i} + α_{1}^{*} Δ t {\dot{Y}}_{i} + α_{2}^{*} Δ t^{2} {\ddot{Y}}_{i}

(23)

M^{*} {\ddot{Y}}_{i + 1} + C^{*} {\dot{Y}}_{i + 1} + K^{*} Y_{i + 1} = Φ^{T} F_{i + 1}

(24)

式中

α_{1}^{*} = Φ^{- 1} α_{1} Φ 和 α_{2}^{*} = Φ^{- 1} α_{2} Φ

为对角型积分参数矩阵。假定阻尼矩阵符合经典阻尼假设，

M^{*} = Φ^{T} M Φ ， C^{*} = Φ^{T} C Φ ， K^{*} = Φ^{T} K Φ

分别为模态质量、阻尼和刚度矩阵。

第 i 阶模态的积分参数可以由 SDOF 体系的积分参数（式（25）、（26））表示，具体如下:

α_{1 i}^{*} = \frac{\frac{c_{7}}{2} Δ t^{3} C_{i}^{*} K_{i}^{*} + \frac{c_{8}}{4} Δ t^{2} C_{i}^{* 2} + c_{9} Δ t^{2} M_{i}^{*} K_{i}^{*} + \frac{c_{10}}{2} Δ t M_{i}^{*} C_{i}^{*} + c_{11} M_{i}^{* 2}}{c_{1} Δ t^{4} K_{i}^{* 2} + \frac{c_{2}}{2} Δ t^{3} C_{i}^{*} K_{i}^{*} + \frac{c_{3}}{4} Δ t^{2} C_{i}^{* 2} + c_{4} Δ t^{2} M_{i}^{*} K_{i}^{*} + \frac{c_{5}}{2} Δ t M_{i}^{*} C_{i}^{*} + c_{6} M_{i}^{* 2}}

(25)

α_{2 i}^{*} = \frac{\frac{c_{12}}{2} Δ t^{3} C_{i}^{*} K_{i}^{*} + \frac{c_{13}}{4} Δ t^{2} C_{i}^{* 2} + c_{14} Δ t^{2} M_{i}^{*} K_{i}^{*} + \frac{c_{15}}{2} Δ t M_{i}^{*} C_{i}^{*} + c_{16} M_{i}^{* 2}}{c_{1} Δ t^{4} K_{i}^{* 2} + \frac{c_{2}}{2} Δ t^{3} C_{i}^{*} K_{i}^{*} + \frac{c_{3}}{4} Δ t^{2} C_{i}^{* 2} + c_{4} Δ t^{2} M_{i}^{*} K_{i}^{*} + \frac{c_{5}}{2} Δ t M_{i}^{*} C_{i}^{*} + c_{6} M_{i}^{* 2}}

(26)

式中:

M_{i}^{*} ， K_{i}^{*} = M_{i}^{*} ω_{i}^{2} 和 C_{i}^{*} = 2 M_{i}^{*} ω_{i} ξ_{i}

分别为第 i 阶模态的质量、刚度和阻尼矩阵，ω_i 和 ξ_i分别为第 i阶模态的圆频率和阻尼比。从而，模态坐标系下的积分参数矩阵

α_{1}^{*} 和 α_{2}^{*}

表达式如下:

α_{1}^{*} = α^{- 1} [\frac{c_{7}}{2} Δ t^{3} C^{*} K^{*} + \frac{c_{8}}{4} Δ t^{2} C^{* 2} + c_{9} Δ t^{2} M^{*} K^{*} + \frac{c_{10}}{2} Δ t M^{*} C^{*} + c_{11} M^{* 2}]

(27)

α_{2}^{*} = α^{- 1} [\frac{c_{12}}{2} Δ t^{3} C^{*} K^{*} + \frac{c_{13}}{4} Δ t^{2} C^{* 2} + c_{14} Δ t^{2} M^{*} K^{*} + \frac{c_{15}}{2} Δ t M^{*} C^{*} + c_{16} M^{* 2}]

(28)

式中

α = c_{1} Δ t^{4} K^{* 2} + \frac{c_{2}}{2} Δ t^{3} C^{*} K^{*} + \frac{c_{3}}{4} Δ t^{2} C^{* 2} + c_{4} Δ t^{2} M^{*} K^{*} + \frac{c_{5}}{2} Δ t M^{*} C^{*} + c_{6} M^{* 2}

。

最后将

α_{1}^{*} 和 α_{2}^{*}

左乘

Φ

，右乘

Φ

^-1得到 NSEρ_∞ 算法用于 MDOF 体系的积分参数矩阵:

α_{1} = Φ α^{- 1} [\frac{c_{7}}{2} Δ t^{3} C^{*} K^{*} + \frac{c_{8}}{4} Δ t^{2} C^{* 2} + c_{9} Δ t^{2} M^{*} K^{*} + \frac{c_{10}}{2} Δ t M^{*} C^{*} + c_{11} M^{* 2}] Φ^{- 1}

(29)

α_{2} = Φ α^{- 1} [\frac{c_{12}}{2} Δ t^{3} C^{*} K^{*} + \frac{c_{13}}{4} Δ t^{2} C^{* 2} + c_{14} Δ t^{2} M^{*} K^{*} + \frac{c_{15}}{2} Δ t M^{*} C^{*} + c_{16} M^{* 2}] Φ^{- 1}

(30)

2 数值特性分析

2.1 稳定性分析

对于线弹性体系，如果放大矩阵的谱半径

ρ = m a x \{|λ_{1}| ， |λ_{2}| ， |λ_{3}|\} ⩽ 1

，则积分算法是稳定的。图1为无阻尼 SDOF 体系时 NSE-ρ_∞ 算法的谱半径，可知无论参数 ρ_∞ 和 γ 取何值，NSE-ρ_∞算法的谱半径 ρ≤1，即算法在线弹性体系中始终无条件稳定。

接着探究 NSE-ρ_∞ 算法在非线性情况下的数值稳定性。利用根轨迹法对算法用于非线性 SDOF 体系的稳定性进行评估。当切线刚度 k_t 从 0 增加到 + ∞ 时，根轨迹可以反映不同结构状态的稳定性，包括线弹性、非线性刚度软化和非线性刚度硬化。在 m = 1，ω = 2π，ξ = 0. 02 和 Δt = 0. 02 的条件下，γ 取 0. 05 时，得到 NSE-ρ_∞算法用于非线性 SDOF 体系的根轨迹图，如图2所示。可以看出，根轨迹均从开环的极点开始，到开环的零点或者无穷远处结束。根轨迹与单位圆交点z =-1 为算法的稳定极限。

式（31）中为 NSE-ρ_∞ 算法用于非线性 SDOF 体系的离散闭环传递函数，令其分母的表达式为 0，则为特征方程。将 z =-1 代入式（31）的特征方程中，可以得到 k_t/ k 的稳定极限，如式（32）所示。

G_{c l}^{N L} (z) = \frac{X (z)}{F (z)} = \frac{G_{1} (z)}{1 + k_{t} G_{1} (z)}

(31)

图1无阻尼 SDOF 体系时 NSE-ρ_∞算法的谱半径

Fig.1Spectral radius of NSE-ρ_∞ algorithm applied to undamped SDOF system

\frac{k_{t}}{k} ⩽ \frac{4 \{- c_{1} Ω^{4} + ξ (c_{7} - c_{2}) Ω^{3} + [(c_{8} - c_{3}) ξ^{2} - c_{4}] Ω^{2} + ξ (c_{9} - c_{5}) Ω - c_{6}\}}{Ω^{2} [(c_{7} - 2 c_{10}) Ω^{2} + ξ (c_{8} - 2 c_{11}) Ω + c_{9} - 2 c_{12}]} = {(\frac{k_{t}}{k})}_{l i m}

(32)

图2NSE-ρ_∞算法用于非线性 SDOF 体系时的根轨迹（γ = 0. 05）

Fig.2Root loci of NSE-ρ_∞ algorithm applied to nonlinear SDOF system (γ = 0. 05)

在 γ 分别取 0. 05、1.99 且 ξ = 0 时求得 NSE-ρ_∞算法的稳定极限，如图3所示，由图3（a）可以看出，当 γ = 0. 05∈（0，1）且 ρ_∞= 1 时，NSE-ρ_∞ 的稳定极限与（k_t/ k）lim = 1 相切，即 NSE-ρ_∞算法存在一个不稳定点为 Ω≈9.176，但是其余 ρ_∞ 取值的 NSE-ρ_∞ 算法对于非线性刚度软化体系和线弹性体系均为无条件稳定。由图3（b）可以看出，当 γ∉（0，1）时，无论 ρ_∞ 取何值，NSE-ρ_∞算法的稳定极限与（k_t / k）lim = 1均不存在相切或者相交点，因此，ρ_∞ 取任何值时，NSE-ρ_∞ 算法对于非线性刚度软化体系和线弹性体系均为无条件稳定。此外，对于非线性刚度硬化体系，相比 KR-α 和 Chang-α 算法，NSE-ρ_∞算法具有更大的无条件稳定区间，即

\frac{k_{t}}{k} \in [0, \frac{- 4}{ρ_{\infty}^{2} + 2 ρ_{\infty} - 3}], ρ_{\infty} \in [0,1)

(33)

图3NSE-ρ_∞ 算法的稳定性极限

Fig.3Stability limits of NSE-ρ_∞algorithm

2.2 精度分析

为了验证 NSE-ρ_∞算法的数值精度，分析有阻尼线弹性 SDOF 体系在受迫振动时的收敛速率。圆频率 ω = 2π，初位移和初速度均为 1，计算持续时间 t = 1 s，阻尼比 ξ = 0. 05，外界荷载 f（t）= sin（πt）。

图4为当 γ = 1.5 时，NSE-ρ_∞算法的收敛速率。可以看出，新算法对于位移、速度和加速度均能保持二阶数值精度，ρ_∞ 的取值影响计算误差的大小，基本的趋势为 ρ_∞ 越大，计算误差越小。

图4NSE-ρ_∞ 算法的收敛速率（γ = 1.5）

Fig.4Convergence rates of NSE-ρ_∞ algorithm (γ = 1.5)

2.3 周期延长和振幅衰减

除收敛速率外，周期延长（period elongation，E_P）和振幅衰减（ amplitude decay，D_A）是两项可以衡量积分算法准确性的指标。具体计算公式如下:

E_{P} = \frac{\bar{T} - T_{0}}{T_{0}} = \frac{Ω}{\bar{Ω}} - 1

(34)

D_{A} = 1 - e x p (- 2 π \bar{ξ} \frac{Ω}{\bar{Ω}})

(35)

式中:

\bar{T} = 2 π / \bar{ω} ， \bar{ω} = \bar{Ω} / Δ t

，T₀为无阻尼 SDOF 结构体系的自振周期。数值阻尼比

\bar{ξ}

和表观频率

\bar{Ω}

的表达式如下:

\bar{ξ} = \frac{- l n (σ^{2} + ε^{2})}{2 \bar{Ω}}

(36)

\bar{Ω} = \frac{{t a n}^{- 1} (ε / σ)}{\sqrt{1 - {\bar{ξ}}^{2}}}

(37)

式中:σ 和 ε 分别为算法放大矩阵 A 特征值的实部与虚部。

由于算法涉及到两个自由参数 γ 和 ρ_∞，采用控制变量的方法分别对 γ 和 ρ_∞对 E_P和 D_A的影响进行分析。

图5给出了当 ρ_∞ = 0.5 且 γ∈（1，2）时，NSEρ_∞ 算法的 E_P和 D_A。可以看出，当 ρ_∞ 取定值时，E_P和 D_A随着 γ 的增大均一直增大。

图6为当 γ = 1.5 且 ρ_∞ ∈[0，1]时，NSE-ρ_∞ 算法的 E_P和 D_A。可以看出，当 γ 取定值时，E_P和 D_A随着 ρ_∞ 的增大均一直减小。

图5NSE-ρ_∞算法的周期延长和振幅衰减（ρ_∞ = 0.5，γ∈（1，2））

Fig.5Period elongation (E_P) and amplitude decay (D_A) of NSE-ρ_∞ algorithm (ρ_∞ = 0.5, γ∈ (1, 2) )

图6NSE-ρ_∞ 算法的周期延长和振幅衰减（γ = 1.5，ρ_∞∈[0，1]）

Fig.6Period elongation (E_P) and amplitude decay (D_A) of NSE-ρ_∞ algorithm (γ = 1.5, ρ_∞ ∈[0, 1])

通过上述分析可以得到，NSE-ρ_∞ 算法在 γ = 1.99（非常接近 γ 选值范围的上限）且 ρ_∞ = 0 时，抑制虚假高频分量的能力非常强，将在下一节通过数值算例进一步验证。

图7对比了 NSE-ρ_∞算法与 GCR、KR-α、Chang-α、ρ_∞-Bathe 等算法的 E_P 和 D_A，可以看出: KR-α 和 Chang-α 具有相同的 E_P和 D_A，因为他们的数值特性均与 G-α 算法匹配; NSE-ρ_∞ 和 ρ_∞-Bathe 的 E_P和 D_A 相同; GCR 的 E_P最大，而 NSE-ρ_∞ 和 ρ_∞-Bathe 的 E_P 最小，具有最优的精度。

图7几种积分算法的周期延长和振幅衰减对比

Fig.7Comparisons of period elongation (E_P) and amplitude decay (D_A) for different integration algorithms

3 数值算例

采用文献^[21] 承受轴向荷载的一维自由固支杆件的轴向波动问题作为数值算例，该模型被广泛用来验证各种积分算法抑制虚假高频分量的能力^[22-24]，其示意如图8所示。

图8轴向荷载激励的自由固支杆件

Fig.8Free-fixed bar subjected to axial load excitation

该模型参数均设置为无量纲数，轴向荷载 F（t）= 10⁴，材料的杨氏模量 E = 3 × 10⁷，材料密度 ρ = 7.3 × 10^-4，横截面面积 A = 1，杆件长度 L = 200。有限元数值模型使用 N = 200 个等长 l = L / N 的两节点线性单元集成。该模型的整体集中质量矩阵 M 和刚度矩阵 K 分别计算如下:

M = \frac{ρ A l}{2} {[\begin{matrix} 2 & 0 \\ 2 \\ 2 \\ ⋱ \\ 2 \\ 0 & 1 \end{matrix}]}_{N \times N}

(38)

K = \frac{E A}{l} {[\begin{matrix} 2 - 1 0 \\ - 1 2 - 1 \\ - 1 2 ⋱ \\ ⋱ ⋱ - 1 \\ - 1 2 - 1 \\ 0 - 1 1 \end{matrix}]}_{N \times N}

(39)

计算该固支杆件中心处的位移和速度时程响应，其位移精确解为

u (x, t) = \frac{F (t)}{E A} x - \frac{8 F (t) L}{π^{2} E A} \sum_{s = 1}^{\infty} [\frac{(- 1)^{s - 1}}{(2 s - 1)^{2}} . s i n \frac{(2 s - 1) π x}{2 L} c o s \frac{(2 s - 1) π c t}{2 L}]

(40)

c = \sqrt{\frac{E}{ρ l}}

(41)

式中 c 表示波在杆件中的传播速度。

使用时间步长 Δt = 9.88 × 10^-7 的 KR-α 算法（ρ_∞ =0）、Chang-α 算法（ρ_∞=0）、G-α算法（ ρ_∞ = 0）、 NSE-ρ_∞算法（γ = 1.99，ρ_∞ = 0）、ρ_∞-Bathe 算法（γ = 1.99，ρ_∞= 0）作为数值解，各算法均选择对应的最大数值阻尼参数，以便于进一步比较各算法对虚假的高频分量抑制能力的强弱。此外，选用 NSE-ρ_∞ 算法（γ = 1.99，ρ_∞ = 1）作为对比，此时，该算法的数值阻尼最小。

图9和图10分别给出了固支杆件中心处的位移和速度时程响应。可以看出，无论对于位移的求解或是速度的求解，当选用最大数值阻尼参数时，KR-α、Chang-α 和 G-α 算法结果接近，而 NSE-ρ_∞和 ρ_∞-Bathe 算法结果接近，相比 KR-α、Chang-α 和 G-α 算法对虚假的高频分量有更强的抑制能力。而对于具有最小数值阻尼的 NSE-ρ_∞ 算法（γ = 1.99，ρ_∞ = 1），其抑制虚假高频分量的能力明显弱于其他算法，特别是对于速度的求解（图10）。

图9杆件中心的位移-时程响应

Fig.9Displacement-time history response of bar center

图10杆件中心的速度-时程响应

Fig.10Velocity-time history response of bar center

4 结论

本文采用半显式的算法格式，提出了一种具有可控数值阻尼的基于模型的积分算法，即 NSE-ρ_∞ 算法。该算法与两子步无条件稳定的 ρ_∞-Bathe 隐式算法具有相同的特征根。 NSE-ρ_∞ 算法不需要引入额外的加权运动方程使算法具有可控的数值阻尼，因此，算法格式更精简。此外，新算法由两个自由参数 γ 和 ρ_∞ 控制算法的数值阻尼，相比 KR-α 算法和 Chang-α 算法仅使用一个自由参数 ρ_∞控制算法的数值阻尼，新算法无疑具有更加灵活的数值阻尼可控性。通过具有代表性的数值算例说明了 NSE-ρ_∞ 算法能够更为有效地抑制虚假高频分量。根据研究结果得到如下结论:

1）NSE-ρ_∞ 算法兼具显式、无条件稳定和可控数值阻尼的优点。以非常简单的算法格式，给算法引入了更加灵活可控的数值阻尼。

2）NSE-ρ_∞ 算法的数值阻尼由两个自由参数 γ 和 ρ_∞控制，其中，γ 的取值范围为（0，2），但应避免取 1 和 2 /（1-ρ_∞）; ρ_∞ 的取值则更为简单，仅需满足[0，1]即可。

3）NSE-ρ_∞ 算法对于线弹性体系和非线性刚度软化体系均为无条件稳定，对于非线性刚度硬化体系，相比 KR-α 和 Chang-α 算法，NSE-ρ_∞算法具有更大的无条件稳定区间。

4）无论是对于位移的求解还是速度的求解，NSE-ρ_∞ 算法对虚假的高频分量具有非常强的抑制能力。

图1无阻尼 SDOF 体系时 NSE-ρ_∞算法的谱半径

Fig.1Spectral radius of NSE-ρ_∞ algorithm applied to undamped SDOF system

下载: 全尺寸图片

图2NSE-ρ_∞算法用于非线性 SDOF 体系时的根轨迹（γ = 0. 05）

Fig.2Root loci of NSE-ρ_∞ algorithm applied to nonlinear SDOF system (γ = 0. 05)

下载: 全尺寸图片

图3NSE-ρ_∞ 算法的稳定性极限

Fig.3Stability limits of NSE-ρ_∞algorithm

下载: 全尺寸图片

图4NSE-ρ_∞ 算法的收敛速率（γ = 1.5）

Fig.4Convergence rates of NSE-ρ_∞ algorithm (γ = 1.5)

下载: 全尺寸图片

图5NSE-ρ_∞算法的周期延长和振幅衰减（ρ_∞ = 0.5，γ∈（1，2））

Fig.5Period elongation (E_P) and amplitude decay (D_A) of NSE-ρ_∞ algorithm (ρ_∞ = 0.5, γ∈ (1, 2) )

下载: 全尺寸图片

图6NSE-ρ_∞ 算法的周期延长和振幅衰减（γ = 1.5，ρ_∞∈[0，1]）

Fig.6Period elongation (E_P) and amplitude decay (D_A) of NSE-ρ_∞ algorithm (γ = 1.5, ρ_∞ ∈[0, 1])

下载: 全尺寸图片

图7几种积分算法的周期延长和振幅衰减对比

Fig.7Comparisons of period elongation (E_P) and amplitude decay (D_A) for different integration algorithms

下载: 全尺寸图片

图8轴向荷载激励的自由固支杆件

Fig.8Free-fixed bar subjected to axial load excitation

下载: 全尺寸图片

图9杆件中心的位移-时程响应

Fig.9Displacement-time history response of bar center

下载: 全尺寸图片

图10杆件中心的速度-时程响应

Fig.10Velocity-time history response of bar center

下载: 全尺寸图片

图1无阻尼 SDOF 体系时 NSE-ρ_∞算法的谱半径

Fig.1Spectral radius of NSE-ρ_∞ algorithm applied to undamped SDOF system

图2NSE-ρ_∞算法用于非线性 SDOF 体系时的根轨迹（γ = 0. 05）

Fig.2Root loci of NSE-ρ_∞ algorithm applied to nonlinear SDOF system (γ = 0. 05)

图3NSE-ρ_∞ 算法的稳定性极限

Fig.3Stability limits of NSE-ρ_∞algorithm

图4NSE-ρ_∞ 算法的收敛速率（γ = 1.5）

Fig.4Convergence rates of NSE-ρ_∞ algorithm (γ = 1.5)

图5NSE-ρ_∞算法的周期延长和振幅衰减（ρ_∞ = 0.5，γ∈（1，2））

Fig.5Period elongation (E_P) and amplitude decay (D_A) of NSE-ρ_∞ algorithm (ρ_∞ = 0.5, γ∈ (1, 2) )

图6NSE-ρ_∞ 算法的周期延长和振幅衰减（γ = 1.5，ρ_∞∈[0，1]）

Fig.6Period elongation (E_P) and amplitude decay (D_A) of NSE-ρ_∞ algorithm (γ = 1.5, ρ_∞ ∈[0, 1])

图7几种积分算法的周期延长和振幅衰减对比

Fig.7Comparisons of period elongation (E_P) and amplitude decay (D_A) for different integration algorithms

图8轴向荷载激励的自由固支杆件

Fig.8Free-fixed bar subjected to axial load excitation

图9杆件中心的位移-时程响应

Fig.9Displacement-time history response of bar center

图10杆件中心的速度-时程响应

Fig.10Velocity-time history response of bar center

图1无阻尼 SDOF 体系时 NSE-ρ_∞算法的谱半径

Fig.1Spectral radius of NSE-ρ_∞ algorithm applied to undamped SDOF system

图2NSE-ρ_∞算法用于非线性 SDOF 体系时的根轨迹（γ = 0. 05）

Fig.2Root loci of NSE-ρ_∞ algorithm applied to nonlinear SDOF system (γ = 0. 05)

图3NSE-ρ_∞ 算法的稳定性极限

Fig.3Stability limits of NSE-ρ_∞algorithm

图4NSE-ρ_∞ 算法的收敛速率（γ = 1.5）

Fig.4Convergence rates of NSE-ρ_∞ algorithm (γ = 1.5)

图5NSE-ρ_∞算法的周期延长和振幅衰减（ρ_∞ = 0.5，γ∈（1，2））

Fig.5Period elongation (E_P) and amplitude decay (D_A) of NSE-ρ_∞ algorithm (ρ_∞ = 0.5, γ∈ (1, 2) )

图6NSE-ρ_∞ 算法的周期延长和振幅衰减（γ = 1.5，ρ_∞∈[0，1]）

Fig.6Period elongation (E_P) and amplitude decay (D_A) of NSE-ρ_∞ algorithm (γ = 1.5, ρ_∞ ∈[0, 1])

图7几种积分算法的周期延长和振幅衰减对比

Fig.7Comparisons of period elongation (E_P) and amplitude decay (D_A) for different integration algorithms

图8轴向荷载激励的自由固支杆件

Fig.8Free-fixed bar subjected to axial load excitation

图9杆件中心的位移-时程响应

Fig.9Displacement-time history response of bar center

图10杆件中心的速度-时程响应

Fig.10Velocity-time history response of bar center

CHOPRA A K. Dynamics of structures:theory and applications to earthquake engineering[M].4th ed. Upper Saddle River: Prentice Hall,2012

邢誉峰, 季奕, 张慧敏, 等. 时间积分方法的研究进展与挑战[J]. 北京航空航天大学学报,2022,48(9):1692.XING Yufeng, JI Yi, ZHANG Huimin,et al. Advances and challenges in time integration methods[J]. Journal of Beijing University of Aeronautics and Astronautics,2022,48(9):1692. DOI:10.13700/j.bh.1001-5965.2022.0288

吴斌, 王贞, 许国山, 等. 工程结构混合试验技术研究与应用进展[J]. 工程力学,2022,39(1):1. WU Bin, WANG Zhen, XU Guoshan,et al. Research and application progress in hybrid testing of engineering structures[J]. Engineering Mechanics,2022,39(1):1. DOI:10.6052/j.issn.1000-4750.2021.10. ST06

TIAN Yingpeng, SHAO Xiaoyu, ZHOU Huimeng,et al. Advances in real-time hybrid testing technology for shaking table substructure testing[J]. Frontiers in Built Environment,2020,6:123. DOI:10.3389/fbuil.2020.00123

DONG X, TANG Z, DU X. State of the art and development trends in numerical simulation for real-time hybrid simulation[J]. Earthquake Engineering and Resilience,2022,1:245. DOI:10.1002/EER2.25

CHANG S Y. Explicit pseudodynamic algorithm with unconditional stability[J]. Journal of Engineering Mechanics,2002,128(9):935. DOI:10.1061/(ASCE)0733-9399(2002)128:9(935)

KOLAY C, RICLES J M. Assessment of explicit and semi-explicit classes of model-based algorithms for direct integration in structural dynamics[J]. International Journal for Numerical Methods in Engineering,2016,107(1):49. DOI:10.1002/nme.5153

CHEN C, RICLES J M. Development of direct integration algorithms for structural dynamics using discrete controltheory[J]. Journal of Engineering Mechanics,2008,134(8):676. DOI:10.1061/(ASCE)0733-9399(2008)134:8(676)

TANG Yu, LOU Menglin. New unconditionally stable explicit integration algorithm for real-time hybrid testing[J]. Journal of Engineering Mechanics,2017,143(7):04017029. DOI:10.1061/(ASCE)EM.1943-7889.0001235

GUI Yao, WANG Jinting, JIN Feng,et al. Development of a family of explicit algorithms for structural dynamics with unconditional stability[J]. Nonlinear Dynamics,2014,77(4):1157. DOI:10.1007/s11071-014-1368-3

FU Bo, ZHANG Futai. A dual-explicit model-based integration algorithm with higher-order accuracy for structural dynamics[J]. Applied Mathematical Modelling,2022,110:513. DOI:10.1016/J. APM.2022.06.005

傅博, 张付泰, 陈瑾. 一种新型高精度半显式基于模型的积分算法[J]. 同济大学学报(自然科学版),2023,51(5):738. FU Bo, ZHANG Futai, CHEN Jin. A new semi-explicit modelbased integration algorithm with high accuracy[J]. Journal of Tongji University(Natural Science Edition),2023,51(5):738. DOI:10.11908/j.issn.0253-374x.21547

KOLAY C, RICLES J M. Development of a family of unconditionally stable explicit direct integration algorithms with controllable numerical energy dissipation[J]. Earthquake Engineering and Structural Dynamics,2014,43(9):1361. DOI:10.1002/eqe.2401

CHUNG J, HULBERT G M. A time integration algorithm for structural dynamics with improved numerical dissipation:the generalized-alpha method[J]. Journal of Applied Mechanics,1993,60(2):371. DOI:10.1115/1.2900803

CHANG S Y. A dual family of dissipative structure-dependent integration methods for structural nonlinear dynamics[J]. Nonlinear Dynamics,2019,98(1):703. DOI:10.1007/s11071-019-05223-y

FU Bo, FENG Decheng, JIANG Huanjun. A new family of explicit model-based integration algorithms for structural dynamic analysis[J]. International Journal of Structural Stability and Dynamics,2019,19(6):1950053. DOI:10.1142/S0219455419500536

TANG Yu, REN Dawei, QIN Hui,et al. New family of explicit structure-dependent integration algorithms with controllable numerical dispersion[J]. Journal of Engineering Mechanics,2021,147(3):04021001. DOI:10.1061/(ASCE)EM.1943-7889.0001901

孟凡涛, 赵建锋. 基于Bathe隐式算法的结构动力学显式算法[J]. 振动与冲击,2019,38(6):226. MENG Fantao, ZHAO Jianfeng. Explicit algorithm for structural dynamics based on Bathe implicit algorithm[J]. Vibration and Shock,2019,38(6):226. DOI:10.13465/j.cnki.jvs.2019.06.034

BATHE K J, BAIG M M I. On a composite implicit time integration procedure for nonlinear dynamics[J]. Computers and Structures,2005,83:2513. DOI:10.1016/j.compstruc.2005.08.001

NOH G, BATHE K. The Bathe time integration method with controllable spectral radius:the ρ_∞-Bathe method[J]. Computers and Structures,2019,212:299. DOI:10.1016/j.compstruc.2018.11.001

GERADIN M. Mechanical vibrations[M].3rd ed. Hoboken: Wiley,2015

WEN Weibin, TAO Yong, DUAN Shengyu,et al. A comparative study of three composite implicit schemes on structural dynamic and wave propagation analysis[J]. Computers & Structures,2017,190:126. DOI:10.1016/j.compstruc.2017.05.006

MALAKIYEH M M, SHOJAEE S, BATHE K J. The Bathe time integration method revisited for prescribing desired numerical dissipation[J]. Computers & Structures,2019,212:289. DOI:10.1016/j.compstruc.2018.10.008

JI Y, XING Y. An optimized three-sub-step composite time integration method with controllable numerical dissipation[J]. Computers and Structures,2020,231:106210. DOI:10.1016/j.compstruc.2020.106210

LI Jinze, LI Hua, LIAN Yiwei,et al. A suite of second-order composite sub-step explicit algorithms with controllable numerical dissipation and maximal stability bounds[J]. Applied Mathematical Modelling,2023,114:601. DOI:10.1016/J. APM.2022.10.012

1 新算法的建立

2 数值特性分析

3 数值算例

4 结论

友情链接LINKS