求解动态子结构的大质量界面虚拟激励法

孙万泉，秦浩玮，马震岳; SUN Wan-quan; QIN Hao-wei; MA Zhen-yue

使用Chrome浏览器效果最佳，继续浏览，你可能不会看到最佳的展示效果，

确定继续浏览么?

复制成功，请在其他浏览器进行阅读

求解动态子结构的大质量界面虚拟激励法 PDF

- ORCID：
孙万泉 ¹
✉
- ORCID：
秦浩玮 ¹
- ORCID：
马震岳 ²

1. 华北电力大学水利与水电工程学院, 北京 102206； 2. 大连理工大学建设工程学部水利工程学院, 辽宁大连 116023

中图分类号： TU311.3

最近更新：2023-11-09

DOI：10.16385/j.cnki.issn.1004-4523.2023.05.008

摘要

为提高大型结构随机激励的计算效率，解决不同类型组合结构的动力求解问题，提出一种动态子结构求解方法。在各子结构界面的质量矩阵上分别附加大质量单元，并直接在大质量单元上施加两次不同幅值的虚拟激励荷载，确保相邻子结构界面上的加速度响应协调；根据不同子结构界面上的内力平衡条件和线性关系建立界面力平衡方程组，求解出能够保证各子结构之间加速度响应和内力同时协调平衡的虚拟激励荷载幅值；通过求出的虚拟激励荷载分别计算各子结构的真实响应。该方法与传统的模态综合法相比，避免了模态坐标和物理坐标的反复转换，并能直接考虑子结构的高频特性；与机械导纳法相比，减小了综合方程规模，求解更为简便，并适用于外荷载激励和地震下基础加速度激励同时输入的情况。最后通过不同类型的仿真算例验证了该方法的正确性。

关键词

动态子结构; 大质量界面虚拟激励法; 随机振动; 虚拟激励

引言

随机振动中的虚拟激励法是一种既能保证计算精度，又能提高计算效率的分析方法^［

1］。该方法得到了国内外学者的广泛应用^{［参考文献 2

百度学术}2］，其中，李永华等^{［参考文献 3

百度学术}3］提出了利用大质量来直接求解绝对位移的虚拟激励法；Zhao等^{［参考文献 4

百度学术}4］从多点输入的动力学方程入手提出了一种直接求解的多点虚拟激励法；Jia等^{［参考文献 5

百度学术}5］通过在多维质量矩阵上“置大数”的方法，从理论上将一维多点激励拓展到多维多点激励。上述研究克服了传统虚拟激励法在有限元软件中求解困难的缺点，解决了传统虚拟激励法需要将绝对位移进行分解后再求解的弊端，在多点随机振动中已经有了广泛的应用。

在实际工程中，很多大型复杂结构或者不同类型的组合结构的动力计算有时会采用动态子结构方法。目前动态子结构的计算方法主要有模态综合类和机械导纳类^［

6］。其中，模态综合类可分为自由界面模态综合法、固定界面模态综合法和混合界面模态综合法三类^{［参考文献 7

百度学术}7］。上述几种模态综合方法都是建立在先分解后复原的思想上，对子结构进行高阶模态的省略，再合并建立整体模态方程求解。另外，以模态综合法为基础衍生出的新型动态子结构方法在不同的领域也有较好的应用，例如，基于波分析的动态子结构方法^{［参考文献 8

百度学术}8］和重复超单元法^{［参考文献 9

百度学术}9］等。这类方法都需要对模态坐标和物理坐标进行相互转化，较为繁琐。而机械导纳类是通过建立以子结构间的节点位移和节点力为未知数的综合方程进行联立求解。但当界面的节点数较多时，该方法的未知数较多，方程规模较大，求解困难。此外，当同时考虑结构外荷载激励和地震下基础加速度激励时，机械导纳法的理论推导和求解也比较繁琐。

为进一步提高动态子结构方法的求解效率和便捷性，本文在绝对位移直接求解的虚拟激励法^［

3，5］的启发下，提出一种求解动态子结构的大质量界面虚拟激励法。通过对各个子结构界面上添加大质量点，并在大质量点上施加虚拟激励荷载，在满足子结构间加速度相等和内力平衡的条件下建立界面平衡方程组，求解能够保持所有子结构界面协调平衡的虚拟激振力幅值，并最终求解出各个子结构的响应。该方法与传统的模态综合法相比，避免了模态坐标和物理坐标的反复转换，并能直接考虑子结构的高频特性；与机械导纳法相比，综合方程里减少了系统响应的未知数，从而减小了求解方程的规模，使求解更为简便。

1 绝对位移直接求解的虚拟激励法

对于一个有 $N$ 个支座和 $n$ 个自由度的离散结构，其多点地震激励运动方程可写成如下的分块矩阵形式^［

1］：

[\begin{matrix} M_{s} & M_{s b} \\ M_{b s} & M_{b} \end{matrix}] \{\begin{matrix} {\ddot{X}}_{s} \\ {\ddot{X}}_{b} \end{matrix}\} + [\begin{matrix} C_{s} & C_{s b} \\ C_{b s} & C_{b} \end{matrix}] \{\begin{matrix} {\dot{X}}_{s} \\ {\dot{X}}_{b} \end{matrix}\} + [\begin{matrix} K_{s} & K_{s b} \\ K_{b s} & K_{b} \end{matrix}] \{\begin{matrix} X_{s} \\ X_{b} \end{matrix}\} = \{\begin{matrix} 0 \\ P_{b} \end{matrix}\}

（1）

式中 m维列向量 $X_{b}$ ， ${\dot{X}}_{b}$ 和 ${\ddot{X}}_{b}$ 分别代表N个支座的地面强迫位移、速度和加速度；n维列向量 $X_{s}$ ， ${\dot{X}}_{s}$ 和 ${\ddot{X}}_{S}$ 分别代表结构系统所有非支座节点位移、速度和加速度；m维列向量 $P_{b}$ 代表地面作用于N个支座的力； $M$ ， $C$ 和 $K$ 分别为结构的质量、阻尼和刚度矩阵；下标“s”和“b”分别对应结构的非支座自由度和支座自由度。

将上式按第二项展开有：

\begin{matrix} M_{b s} \cdot {\ddot{X}}_{s} + M_{b} \cdot {\ddot{X}}_{b} + C_{b s} \cdot {\dot{X}}_{s} + C_{b} \cdot {\dot{X}}_{b} + \\ K_{b s} \cdot X_{s} + K_{b} \cdot X_{b} = P_{b} \end{matrix}

（2）

因连接地面处的支座地震力 $P_{b}$ 等于支座质量 $M_{b}$ 与地面加速度 ${\ddot{U}}_{b}$ 的乘积，即 $P_{b} = M_{b} \cdot {\ddot{U}}_{b}$ ，此时有：

\begin{matrix} M_{b s} \cdot {\ddot{X}}_{s} + M_{b} \cdot {\ddot{X}}_{b} + C_{b s} \cdot {\dot{X}}_{s} + C_{b} \cdot {\dot{X}}_{b} + \\ K_{b s} \cdot X_{s} + K_{b} \cdot X_{b} = M_{b} \cdot {\ddot{U}}_{b} \end{matrix}

（3）

在式（3）左右两边同乘 ${M_{b}}^{- 1}$ 可得：

{M_{b}}^{- 1} (M_{b s} \cdot {\ddot{X}}_{s} + M_{b} \cdot {\ddot{X}}_{b} + C_{b s} \cdot {\dot{X}}_{s} + C_{b} \cdot {\dot{X}}_{b} + K_{b s} \cdot X_{s} + K_{b} \cdot X_{b}) = {M_{b}}^{- 1} \cdot M_{b} \cdot {\ddot{U}}_{b}

（4）

当 ${M_{b}}^{- 1} \to 0$ 时， ${M_{b}}^{- 1} (M_{b s} \cdot {\ddot{X}}_{s} + C_{b s} \cdot {\dot{X}}_{s} + C_{b} \cdot {\dot{X}}_{b} + K_{b s} \cdot X_{s} + K_{b} \cdot X_{b})$ 也为0，则可以得到：

{\ddot{X}}_{b} = {\ddot{U}}_{b}

（5）

由此可得，只要支撑处的质量矩阵足够大，即在结构总体质量矩阵中支撑质量矩阵处置入一个足够大的数，就能保证支撑处的地面加速度与支撑地震响应加速度相等^［

10］，从而使得式（1）的求解容易很多。

将式（5）代入式（1）第一项后展开可得：

\begin{matrix} M_{s} \cdot {\ddot{X}}_{s} + M_{s b} \cdot {\ddot{U}}_{b} + C_{s} \cdot {\dot{X}}_{s} + C_{s b} \cdot {\dot{X}}_{b} + \\ K_{s} \cdot X_{s} + K_{s b} \cdot X_{b} = 0 \end{matrix}

（6）

式（6）中未知量 ${\ddot{U}}_{b}$ ， ${\dot{X}}_{b}$ 和 $X_{b}$ 可以由下式求得^［

1］：

S (i ω) = P^{*} P^{T}

（7）

P = e^{- i ω t} S_{X} Q

（8）

{\ddot{U}}_{b} = P e^{i ω t}

（9）

{\dot{X}}_{b} = {\ddot{U}}_{b} / (ω i)

（10）

X_{b} = - {\ddot{U}}_{b} / ω^{2}

（11）

式中 $S (i ω)$ 为加速度功率谱密度； $ω$ 为角频率；上标“T”和“*”分别为矩阵的转置和伴随； $e$ 为自然对数的底数； $i$ 为虚数单位； $S_{X}$ 为忽略局部场地效应下的地面节点位移； $t$ 为加速度功率谱分解后得到的时间分量； $P$ 为n $\times$ r（r $\leq n$ ）实阵；Q为与行波效应相关的n $\times$ r实阵。

式（6）整理后可得：

\begin{matrix} M_{s} \cdot {\ddot{X}}_{s} + C_{s} \cdot {\dot{X}}_{s} + K_{s} \cdot X_{s} = - (M_{s b} \cdot {\ddot{U}}_{b} + \\ C_{s b} \cdot {\dot{X}}_{b} + K_{s b} \cdot X_{b}) \end{matrix}

（12）

将式（9）~（11）代入式（12）整理可得：

\begin{matrix} M_{s} \cdot {\ddot{X}}_{s} + C_{s} \cdot {\dot{X}}_{s} + K_{s} \cdot X_{s} = - (M_{s b} + \\ C_{s b} / (ω i) - K_{s b} / ω^{2}) \cdot P e^{i ω t} \end{matrix}

（13）

根据虚拟激励法理论^［

1］假定阻尼力不是与绝对速度而是与相对速度成正比，即将公式（12）中的

{\dot{X}}_{s}

表示为动态相对速度，并用0矩阵代替

{\dot{U}}_{b}

。此时相当于将公式（13）中等效荷载项中的阻尼项忽略，因此可以得到虚拟激励动力方程：

M_{s} \cdot {\ddot{X}}_{s} + C_{s} \cdot {\dot{X}}_{s} + K_{s} \cdot X_{s} = - (M_{s b} - K_{s b} / ω^{2}) \cdot P e^{i ω t}

（14）

式（14）可以看作在每一个频率点 $ω$ 处的简谐运动方程。

由此可得，在结构边界处添加大质量M后，对边界施加 $M \cdot P e^{i ω t}$ 的加速度激励，即可得到边界处加速度为 $P e^{i ω t}$ 的结构内部动力方程。

2 大质量界面虚拟激励法的提出

以最简单的两个子结构为例进行该方法的说明，具体实现步骤如下：

第一步，将结构划分为子结构1和2，其中结构外荷载位于子结构1的内部节点上，记为 $F_{1}^{m}$ 。此时对于子结构1和2分别有：

[\begin{matrix} M_{1}^{m m} & M_{1}^{m s} \\ M_{1}^{s m} & M_{1}^{s s} \end{matrix}] \{\begin{matrix} {\ddot{X}}_{m 1} \\ {\ddot{X}}_{s 1} \end{matrix}\} + [\begin{matrix} C_{1}^{m m} & C_{1}^{m s} \\ C_{1}^{s m} & C_{1}^{s s} \end{matrix}] \{\begin{matrix} {\dot{X}}_{m 1} \\ {\dot{X}}_{s 1} \end{matrix}\} +

[\begin{matrix} K_{1}^{m m} & K_{1}^{m s} \\ K_{1}^{s m} & K_{1}^{s s} \end{matrix}] \{\begin{matrix} X_{m 1} \\ X_{s 1} \end{matrix}\} = \{\begin{matrix} F_{1}^{m} \\ F_{1}^{s} \end{matrix}\}

（15）

[\begin{matrix} M_{2}^{m m} & M_{2}^{m s} \\ M_{2}^{s m} & M_{2}^{s s} \end{matrix}] \{\begin{matrix} {\ddot{X}}_{m 2} \\ {\ddot{X}}_{s 2} \end{matrix}\} + [\begin{matrix} C_{2}^{m m} & C_{2}^{m s} \\ C_{2}^{s m} & C_{2}^{s s} \end{matrix}] \{\begin{matrix} {\dot{X}}_{m 2} \\ {\dot{X}}_{s 2} \end{matrix}\} +

[\begin{matrix} K_{2}^{m m} & K_{2}^{m s} \\ K_{2}^{s m} & K_{2}^{s s} \end{matrix}] \{\begin{matrix} X_{m 2} \\ X_{s 2} \end{matrix}\} = \{\begin{matrix} 0 \\ F_{2}^{s} \end{matrix}\}

（16）

式中下标“1”和“2”为子结构的编号；角标“m”和“s”分别为子结构的内部自由度和界面自由度； $F_{1}^{s}$ 和 $F_{2}^{s}$ 分别为子结构1和2的交界面上的界面力。

第二步，分别在子结构1和2的交界面节点处添加大质量单元M₀（结构整体质量的10⁵~10⁸倍）。并对外荷载 $F_{1}^{m}$ 做Fourier变换，得到外荷载在不同频率 $ω$ 下的Fourier谱 $F_{1}^{m} (ω)$ 。

第三步，分别在子结构1和2界面的大质量单元上施加与其实际运动方向相同，且大小、频率一致的单位简谐虚拟激励力 $F_{1}^{s} (ω)$ 和 $F_{2}^{s} (ω)$ （ $F_{1}^{s} (ω) = F_{2}^{s} (ω)$ ）。并与外荷载 $F_{1}^{m} (ω)$ 一起在频域内对两个子结构进行求解。此时，由于大质量单元的存在，在相同的界面虚拟激励下，两个子结构界面的加速度响应是相同的，即 ${\ddot{X}}_{s 1} = {\ddot{X}}_{s 2} = F_{1}^{s} (ω) / M_{0}$ ，原理如上节所述。此时两个子结构界面处的内力方向相反，绝对值并不相同，即 $N_{s 1} \neq N_{s 2}$ 。

第四步，与上一步操作相同，分别在两个子结构的界面大质量单元上重新施加另一个不同幅值的简谐虚拟激振力 $F_{1}^{s'} (ω)$ 和 $F_{2}^{s'} (ω)$ （ $F_{1}^{s^{'}} (ω) = F_{2}^{s'} (ω)$ ）。同样地，求得该虚拟力下的界面加速度和内力绝对值结果为： ${\ddot{X}}_{s 1}^{'} = {\ddot{X}}_{s 2}^{'} = F_{1}^{s'} (ω) / M_{0}$ ， $N_{s 1}^{'} \neq N_{s 2}^{'}$ 。

第五步，为保证两个子结构之间的界面加速度和界面内力同时协调平衡，也就是要求在哪个虚拟激励幅值（ $F^{*} (ω)$ ）下，可以使得两个子结构的界面加速度都等于 ${\ddot{X}}^{*}$ ，界面内力绝对值也都等于 $N^{*}$ （但方向相反）。因此，可以根据线弹性关系联立方程如下：

\{\begin{array}{l} \frac{F^{*} (ω) - F_{1}^{s} (ω)}{F_{1}^{s^{'}} (ω) - F_{1}^{s} (ω)} = \frac{N^{*} - N_{s 1}}{N_{s 1}^{'} - N_{s 1}} \\ \frac{F^{*} (ω) - F_{2}^{s} (ω)}{F_{2}^{s^{'}} (ω) - F_{2}^{s} (ω)} = \frac{N^{*} + N_{s 2}}{N_{s 2}^{'} - N_{s 2}} \\ F_{1}^{s} (ω) = F_{2}^{s} (ω) \\ F_{1}^{s^{'}} (ω) = F_{2}^{s'} (ω) \end{array}

（17）

以上方程类似两条直线求交点，比较容易理解。整理方程（17）后可求得 $F^{*} (ω)$ 如下式所示：

F^{*} (ω) = \frac{(N_{s 1} + N_{s 2}) * [F_{1}^{s} (ω) - F_{1}^{s^{'}} (ω)]}{(N_{s 1}^{'} - N_{s 1}) - (N_{s 2}^{'} - N_{s 2})} + F_{1}^{s} (ω)

（18）

可见，当系统只有一个交界面一个节点时，该方法也只有一个待求未知数 $F^{*} (ω)$ 。

另外，当子结构界面上有n个节点时，式（17）和（18）中各变量应写为矩阵形式。界面上所有节点应全部施加大质量单元和虚拟激励，考虑各节点内力的相互关联，然后联立求解方程。例如： $F^{*} (ω) = [\begin{matrix} F_{1}^{*} \\ ⋮ \\ F_{n}^{*} \end{matrix}]$ ，其中 $F_{i}^{*}$ 表示该界面上第i个节点上的待求激励； $F_{1}^{s} (ω) = [\begin{matrix} F_{1 ‑ 1}^{s} \\ ⋮ \\ F_{1 ‑ n}^{s} \end{matrix}]$ ，其中 $F_{1 ‑ i}^{s}$ 表示在第一个子结构界面的第i个节点上施加的虚拟激励； $N_{s 1} = [\begin{matrix} N_{s 1 ‑ 1}^{1} & \dots & N_{s 1 ‑ n}^{1} \\ ⋮ & ⋮ \\ N_{s 1 ‑ n}^{n} & \dots & N_{s 1 ‑ n}^{n} \end{matrix}] \times [\begin{matrix} 1 \\ ⋮ \\ 1 \end{matrix}]$ = $[\begin{matrix} N_{s 1}^{1} \\ ⋮ \\ N_{s 1}^{n} \end{matrix}]$ ，其中， $N_{s 1}^{j}$ 表示子结构1上第j个节点处的内力， $N_{s 1 ‑ i}^{j}$ 表示子结构1在第i个节点处施加虚拟激励后在第j个节点处产生的内力。

第六步，在求出能够保证两个子结构的界面加速度和界面内力都协调平衡的虚拟激励 $F^{*} (ω)$ 后，就可以根据公式（15）和（16）求出各自子结构的真实内部响应。

第七步，采用以上方法对全频段进行遍历求解，即可得到结构在随机激励下的全部频域响应，根据频域解再应用Fourier逆变换，就可得到时域解。

同理，如果系统包含n个子结构交界面，利用以上方法可以联立方程组求解n个未知数 $F_{n}^{*} (ω)$ 。由于篇幅所限，不再赘述。

3 工程算例

首先以相对简单的多自由度集中质量模型进行分析，验证这一求解方法。然后再以一个相对复杂的地震和外荷载同时输入的桥梁模型，以及地震作用下的大型复杂水电站厂房结构模型进行进一步计算验证。其中，从工程角度，子结构的划分应遵循以下原则^［

11］：按照实际复杂结构几何形状和装配关系划分；尽量割断较少的联系，即尽量用较少的修改就能取得化整为零的最大效果等等。

3.1 集中质量模型

为了验证本文提出的界面虚拟激励法求解动态子结构的准确性和计算精度，先以六个自由度的集中质量杆件系统在水平简谐激励作用下的响应计算为例进行分析。计算模型和子结构划分如图1所示，水平简谐激励 $F = s i n (10 π \cdot t)$ kN，作用于第一层质量单元。各单元间距为2 m，各单元质量为M=10 kg。

图1 集中质量模型

Fig.1 Lumped-mass model

利用动态子结构模型的计算结果与整体结构计算结果进行相互对比验证，记结果误差值 $P = \frac{| S 2 - S 1 |}{S 2} \times 100 %$ ，其中，本文的子结构求解方法的幅值解为S1，整体模型求解的幅值结果为S2。计算结果对比如表1所示。

表1 利用界面虚拟激励法的子结构求解与整体结构求解结果对比

Tab.1 Comparison of the substructure calculating results and the global structure calculating results using the interface virtual excitation method

单元号（从下至上）	振幅S2/m	振幅S1/m	误差/%
1	0.21988	0.21995	0.03
2	0.08784	0.08787	0.03
3	0.01591	0.01608	1.06
4	-0.18593	-0.18497	0.52
5	-0.21534	-0.21422	0.52
6	-0.20626	-0.20653	0.13

由表1结果可知，利用本文方法的计算结果与整体直接计算结果几乎相等，误差最大为1.06%，因此可以有效证明本文方法的准确性和精度。

3.2 地震和外荷载同时作用下的桥梁结构分析

以某混凝土桥梁为例，全桥长为360 m，跨径布置为110 m+160 m+100 m，宽为13 m，厚度为3.5 m，两墩分别为高65和105 m、宽3.5 m的实心矩形墩。该桥三维模型和子结构划分如图2所示，模型采用六面体实体网格剖分，共计376个单元和805个节点。子结构界面上的节点全部施加大质量单元和虚拟激励求解。墩底处完全约束，桥台处约束竖向位移。桥墩处受到垂直于桥向的地震加速度荷载，同时，桥面中心受到竖直向下，均方根值为10 kN的白噪声随机激励的外荷载F。

图2 三维有限元模型

Fig. 2 Three-dimensional finite element model

地震加速度模型选用胡聿贤和周锡元的修正Kanai‑Tajimi加速度随机模型：

S (ω) = \frac{1 + 4 ξ_{g}^{2} (ω / ω_{g})^{2}}{[1 - (ω / ω_{g})^{2}]^{2} + 4 ξ_{g}^{2} (ω / ω_{g})^{2}} \cdot \frac{ω^{6}}{ω^{6} + ω_{c}^{6}} \cdot S_{0}

（19）

式中 S₀为谱强度因子； $ω_{g}$ 和 $ξ_{g}$ 分别为场地特征频率和阻尼比； $ω_{c}$ 为控制低频含量参数。根据文献［

12］进行参数取值，

ω_{g} = 13.96 r a d / s

，

ξ_{g} = 0.8

，

ω_{c} = 0.6 π

，

S_{0} = 8.6697 c m^{2} / s^{3}

。

根据结构的自振频率，取频率积分区间为 $ω \in [0.5,60] r a d / s$ ，步长为0.5 $r a d / s$ 。地震虚拟力荷载为地震加速度与大质量的乘积，其中地震加速度是由地震谱转换而来的，如公式（7）~（9）所示。为更清晰地验证本文方法在不同频段下都具有足够的准确性，取本文动态子结构方法的计算结果与整体结构直接求解的结果进行比较。选取子结构1上的节点A，子结构1和2连接界面上的节点B，以及子结构3上的节点C进行频域结果比较，并单独列出高频区的结果进行对比。在地震加速度与桥面外荷载共同激励下的计算结果如图3~5所示。

图3 节点A的竖向位移幅值

Fig.3 Vertical displacement amplitudes of node A

图4 节点B的竖向位移幅值

Fig.4 Vertical displacement amplitudes of node B

图5 节点C的竖向位移幅值

Fig.5 Vertical displacement amplitudes of node C

从图3~5中可以看出，求解动态子结构的界面虚拟激励法在全频段内，无论是子结构界面还是子结构内部，在求解随机激励问题时与整体结构直接求解结果都非常吻合，该方法的子结构计算结果与整体求解的误差在0.6%左右。极少数点的误差虽能达到6%，但是从计算结果来看已经达到10^-6量级，对整体响应结果的影响很小。

动力子结构方法能够大幅度地缩减动力分析的规模，实现大规模复杂模型的分块并行求解，避免了整体高阶动力学方程的求解，提高了计算效率。计算效率的定量分析与具体动力模型的规模和子结构的划分有关。本算例中，整体模型计算时间为9.3 s，运用子结构计算时间为4 s，计算效率提高了57%。随着模型规模和计算步数的增加，该子结构法的计算优势也更加突出。

3.3 地震作用下的大型复杂结构分析

水电站厂房是由大体积不规则混凝土结构和上部框架板肋结构组成的大型复杂结构体系。以某水电站厂房为例，根据其结构形式特点，划分为上部和下部两个子结构，如图6所示。厂房总长为65 m、宽为35 m、高为48.5 m，模型共计83050个单元和53108个节点。边界条件为基础底面与尾水管上游侧完全约束，采用与3.2节相同的加速度模型在基础底面施加与厂房横河向呈45°角的水平地震加速度荷载。

图6 三维有限元模型

Fig.6 Three-dimensional finite element model

根据结构的自振频率，取频率区间为 $ω \in [0,35] H z$ ，取步长为1 Hz。为验证本文方法在大型复杂结构下也具有足够的先进性与准确性，取本文动态子结构方法的计算结果与整体结构直接求解的结果进行比较。选取子结构1上的节点a和b，子结构2上的节点c和d进行横向与纵向的频域结果比较，如图7~10所示。

图7 节点a在0~35 Hz频段的位移幅值

Fig.7 Displacement amplitude of node a in 0~35 Hz frequency band

图8 节点b在0~35 Hz频段的位移幅值

Fig.8 Displacement amplitude of node b in 0~35 Hzfrequency band

图9 节点c在0~35 Hz频段的位移幅值

Fig.9 Displacement amplitude of node c in 0~35 Hzfrequency band

图10 节点d在0~35 Hz频段的位移幅值

Fig.10 Displacement amplitude of node d in 0~35 Hz frequency band

从图7~10中可以看出，求解动态子结构的大质量界面虚拟激励法在整体模型较大、较复杂时也具有足够的精度。其中，各个点位的求解误差均小于5%。该大型复杂结构在整个频率区间的整体模型计算时间为560 s，分为两个子结构的情况下计算时间为455 s，计算效率提高了19%。

4 结论

在直接求解的虚拟激励法的原理上提出了一种求解动态子结构的大质量界面虚拟激励法，通过理论推导与计算分析，可得到以下结论：

借助大质量点可以保证相邻子结构的界面节点在相同的虚拟激励下拥有相同的加速度。在保证界面加速度相同，界面内力数值相等，且方向相反的情况下，可以推导出保证各个子结构界面运动协调的虚拟激励荷载。从而求解各个子结构的运动方程，实现大型结构的简化计算。

本文提出的方法在界面单节点与界面多节点以及一维与多维的情况下都具有较高的准确性，并且大幅减小了模型整体计算求解的规模，实现了各个子结构模型单独求解，提高了模型的计算效率。相较于传统动态子结构方法，在处理混合荷载作用下的结构时也能保证足够的精度和简便。

参考文献

林家浩，张亚辉. 随机振动的虚拟激励法［M］. 北京：科学出版社， 2004： 126-143. [百度学术]

Lin Jiahao， Zhang Yahui. Pseudo Excitation Method for Random Vibration［M］. Beijing： Science Press， 2004： 126-143. [百度学术]

Zhou Bo， Guo Wei， Han Xiaoshuang， et al. Random seismic response analysis of jacket structure with Timoshenko's beam theory［J］. Ships and Offshore Structures， 2016， 11（4）： 438-444. [百度学术]

李永华，李思明. 绝对位移直接求解的虚拟激励法［J］. 振动与冲击， 2009， 28（10）： 185-190. [百度学术]

LI Yonghua， LI Siming. Pseudo excitation method based on solving absolute displacement［J］. Journal of Vibration and Shock， 2009， 28（10）： 185-190. [百度学术]

Zhao B， Wang Y Q， Chen Z H， et al. Research on the random seismic response analysis for multi- and large-span structures to multi-support excitations［J］. Earthquake Engineering and Engineering Vibration， 2015， 14（3）： 527-538. [百度学术]

Jia H Y， Zhang D Y， Zheng S X， et al. Local site effects on a high-pier railway bridge under tridirectional spatial excitations： nonstationary stochastic analysis［J］. Soil Dynamics and Earthquake Engineering， 2013， 52： 55-69. [百度学术]

Su Jintao， Lou Jianping， Jiang Xiaolu. Review of the theory of hybrid modeling［J］. IOP Conference Series： Earth and Environmental Science， 2021， 769（4）： 042077. [百度学术]

邹明松，吴有生，孙建刚，等. 动态子结构高效集成的主动力和界面力处理方法［J］. 振动工程学报， 2019， 32（3）： 439-443. [百度学术]

ZOU Mingsong， WU Yousheng， SUN Jiangang， et al. A method to deal with external and interface forces for efficient integration of dynamic substructure analysis［J］. Journal of Vibration Engineering， 2019， 32（3）： 439-443. [百度学术]

孙宜强，陶江平，王波. 基于波分析的动态子结构法及其在船舶管系振动计算中的应用［J］. 船舶工程， 2021， 43（增刊1）： 299-302. [百度学术]

SUN Yiqiang， TAO Jiangping， WANG Bo. Dynamic substructure method based on wave approach and its application in vibration calculation of ship pipe system［J］. Ship Engineering， 2021， 43（Sup1）： 299-302. [百度学术]

姜东，徐宇，曹鑫，等. 行星轮系动态特性分析的重复超单元方法［J］. 机械设计与制造， 2020（7）： 67-70. [百度学术]

Jiang Dong， Xu Yu， Cao Xin， et al. Repetitive super-element method of dynamic analysis of planetary gear system［J］. Machinery Design & Manufacture， 2020（7）： 67-70. [百度学术]

郑史雄，贾宏宇，张金，等. 直接求解多维多点非平稳激励的虚拟激励法［J］. 工程力学， 2014， 31（1）： 85-90. [百度学术]

ZHENG Shixiong， JIA Hongyu， ZHANG Jin， et al. Pseudo excitation method of direct solving multi-dimensional and multi-support non-stationary excitation［J］. Engineering Mechanics， 2014， 31（1）： 85-90. [百度学术]

王永岩. 动态子结构方法理论及应用［M］. 北京：科学出版社， 1999： 13-14. [百度学术]

Wang Yongyan. Theory and Application of Dynamic Substructure Method［M］. Beijing： Science Press， 1999： 13-14. [百度学术]

DING Yang， LIN Wei， LI Zhongxian. Non-stationary random seismic analysis of long-span spatial structures under multi-dimensional and multi-support earthquake excitations［J］. Engineering Mechanics， 2007， 24（3）： 97-103. [百度学术]