基于全时程迭代的减震结构地震响应分析方法与ABAQUS二次开发

Seismic Response Analysis Method of Damping Structures Based on the Inter-history Iteration and Secondary Development of ABAQUS

摘要

当计算设置有黏滞阻尼器的减震结构的地震响应时，直接积分法通常需要大量计算，等效线性法又具有一定的误差，故研究既简便又高精度的分析方法具有重要意义。本文以等效线性方法为原理，以黏滞阻尼器Maxwell恢复力模型为例，提出基于全时程迭代的结构地震响应分析方法，即采用一个附加等效线性阻尼器代替这一层建筑的所有阻尼器，计算其附加等效阻尼系数，通过地震作用下结构的速度响应更新附加等效阻尼系数。然后，以更新后的附加等效阻尼系数重新计算结构的速度时程，进行迭代计算直至得到符合要求的附加等效阻尼系数；采用ABAQUS实现全时程迭代法的计算，并进行二次开发实现自动迭代计算。最后，采用直接积分法对结果进行验证。结果表明：在7层框架结构算例中，全时程迭代法具有良好的收敛性，与直接积分法的顶层位移响应曲线吻合较好，计算精度较高，误差满足工程精度的要求，可靠度高，可提高分析效率；由ABAQUS二次开发得到的自动迭代计算方法可减少人为处理数据的繁琐，降低操作误差，能简化全时程迭代方法的计算过程，使全时程迭代方法可更广泛地应用于工程实际。

Abstract

Objective Two typical approaches for assessing the seismic response of a structure equipped with viscous dampers are direct integration methods and equivalent linear methods. Although direct integration methods provide precise results, their extensive computational cost and susceptibility to numerical instability stemming from the integration step hinder their practical use in most engineering applications. The equivalent linearization method provides a more tractable alternative by linearizing the system and incorporating stiffness and damping proportional to the response level. However, the equivalent linearization method is based on a single-degree-of-freedom structure, and the seismic response of the structure is assumed to be a sine function. For this reason, practical engineering applications of the equivalent linear method can also introduce errors. Therefore, proposing a simple yet highly precise analysis method to calculate the seismic response of a structure equipped with viscous dampers is of great significance. Method This study introduces an alternative method for seismic response analysis of structures incorporating viscous dampers. Based on the equivalent linear method, the Maxwell model is used as the assumed restoring force model to account for dampers. The analysis is based on an inter-history iteration approach to simulate the response over time. In this approach, an extra equivalent linear damper replaces all dampers on each floor, and its supplementary equivalent damping coefficient is computed and updated based on the velocity response of the structure during an earthquake. Then, the velocity response of the structure is recalculated using the updated additional equivalent damping coefficient. This iterative process continues until the additional equivalent damping coefficient meets the specified criteria. The increasing adoption of finite element analysis (FEA) platforms in civil engineering enables high-precision simulation of structural dynamic behaviors through advanced numerical modeling techniques. Among the various finite element software options, ABAQUS stands out for providing several user-subroutine interfaces. These interfaces empower users to expand the functionality of the main program by writing specific function codes designed to their requirements. It becomes practical and meaningful to develop a seismic response analysis method by integrating the concept and approach of inter-history iteration and using the capabilities of ABAQUS software. This method aims to avoid complicated calculations and achieve rapid convergence, which is of great significance to the seismic design of structures equipped with viscous dampers in practical engineering. Result and Discussion This study introduces the Maxwell viscous damper model at the beginning. After the brief introduction, the inter-history iteration method, grounded in the equivalent linear method, is presented. This method is used for the seismic response analysis of structures equipped with viscous dampers. ABAQUS is utilized to model and analyze a seven-story frame structure. Simultaneously, the seven-story frame structure is also modeled in MATLAB, which calculates the structural response using the direct integration method. The outcomes obtained through the direct integration method in MATLAB are used as the standard of comparison to the modeling and analysis in ABAQUS. The error of the first five modes is small, with the error of the first and second modes being less than 0.1%, demonstrating that the two modeling approaches are very close. The comparison and error of the first five modes in the two software verify the correctness of the model in ABAQUS. Python is further utilized to perform the secondary development of ABAQUS to realize the automatic iteration of the inter-history iteration method, thus obtaining a method that efficiently calculates the additional equivalent damping coefficient. Finally, a comparative analysis is conducted between the results obtained from the automated inter-history iteration method developed in ABAQUS and those from the direct integration method. The results demonstrated that, in the case of a seven-story frame structure, the inter-history iteration method exhibits strong convergence. In addition, it aligns well with the top displacement response curve calculated using the direct integration method. The error in the results meets engineering accuracy requirements, confirming that the automatic inter-history iteration calculation method provides high calculation accuracy and reliability while enhancing analysis efficiency. The error rate between the two calculation methods remains within 5%, indicating high accuracy. Simultaneously, compared to the direct integration method, the inter-history iterative method produces a smaller calculation result. This discrepancy arises because, in the calculation formula of the additional equivalent damping coefficient, the hysteresis curve corresponding to a larger response amplitude is fuller, while the curve for a smaller response amplitude is flatter. Therefore, the calculated additional equivalent damping coefficient is larger, leading to a smaller result in the inter-history iteration method. The automatic inter-history iteration calculation method, developed through the secondary development of ABAQUS, reduces manual data processing complexity and minimizes operational errors, significantly simplifying the calculation process of the inter-history iterative method. Based on the analysis above, the inter-history iterative method demonstrates the potential for broader application in engineering practice. Conclusion The inter-history iterative method presented in this study provides several advantages: 1) Utilization of ABAQUS. The seismic response analysis of damped structures is performed using ABAQUS, eliminating the need for a complex numerical integration method. 2) Efficiency of the iterative process. Unlike the extensive computational steps required by the numerical integration method, the inter-history iterative method achieves strong convergence with only a few iterations. 3) Nonlinear problem avoidance. The model employs an equivalent linear damper, allowing for the avoidance of nonlinear challenges.

Graphical abstract

关键词

黏滞阻尼器 / 全时程迭代 / 附加等效阻尼系数 / ABAQUS / 二次开发

Key words

seismic response / linearization / damping / ABAQUS / secondary development

本刊网刊

在土木工程结构振动控制领域中，黏滞阻尼器是一种高效的耗能减震装置并被广泛应用于实际工程项目^[1]。通过合理的减震设计，在结构中设置黏滞阻尼器可耗散地震传递的能量，从而减小主体结构在多遇地震下的地震响应，使主体结构能保持线性的工作状态，仅局部布置的黏滞阻尼器表现出非线性状态，形成了典型的局部非线性问题^[2‒3]。

黏滞阻尼器作为一种广泛应用的消能减震技术，具有构造简单、便于施工和减震机理明确等优点^[4]，其恢复力计算模型在实际工程中常被简化，相关学者已提出了一系列有效的模型，如线性、非线性、Maxwell、Kevin和Wiechert模型^[5‒8]，这些模型基于不同的假设，采用不同的公式描述黏滞阻尼器的力‒位移行为。

目前，针对装有黏滞阻尼器的结构地震响应研究的方法主要分为两种：直接积分法和等效线性法。直接积分法需将黏滞阻尼器恢复力模型代入结构的运动方程，采用如Newmark^[9]、Rosenbrock^[10‒11]等计算法获得结构的地震响应，但其计算量大，稳定性受到积分步长的影响，在普通工程中不适用。等效线性方法将非线性的阻尼器转换为等效的线性阻尼器，通过等效的阻尼系数，利用时程分析方法即可求解结构的地震响应^[12]。但等效线性方法基于单自由度结构且需要假定结构的地震响应为正弦函数，而实际上地震运动是随机的，因此等效线性方法在实际工程的应用中也存在着一定的误差^[13]。贾传果等^[2]以等效线性方法为基本原理，提出全时程迭代思想，推导出多自由度结构中黏滞阻尼器的附加等效阻尼系数的计算公式，可减少等效线性方法的误差。

有限元软件在土木工程的应用日渐丰富，可帮助使用者更好地分析结构响应。其中，ABAQUS提供了若干用户子程序接口，允许使用者通过子程序编写满足具体功能的代码，扩展主程序的功能^[14‒15]。罗维刚等^[16]采用ABAQUS自身提供的二次开发接口，研发了针对钢筋混凝土框架结构前处理的建模插件，可有效减少重复性建模工作。高喆等^[17]使用MATLAB和Python语言针对性地输出ABAQUS模型中特定节点或杆件的应力应变，并运用到实际相应节点或杆件的强度分析中，证明了后处理二次开发技术的可靠性，提高了用户的工作效率。

以全时程迭代的思想与方法为基础，结合ABAQUS软件的功能寻找一种可避免大量复杂计算、尽可能快速收敛的地震响应分析方法，对于实际工程中黏滞阻尼器减震结构的抗震设计具有重要意义。

本文首先介绍选用的Maxwell黏滞阻尼器模型，然后基于全时程迭代方法，给出黏滞阻尼器减震结构的地震响应分析方法；针对7层黏滞阻尼器减震框架结构中的一榀框架，使用ABAQUS建模分析并结合MATLAB验证此分析方法的可靠性^[18]；再使用Python对ABAQUS进行二次开发实现该计算方法的自动迭代，得到可高效计算附加等效阻尼系数的方法。

1 Maxwell黏滞阻尼器恢复力模型

当黏滞阻尼器具有刚度特点且其力学性能与频率相关时，Maxwell模型可更为精确地进行阻尼器恢复力计算^[19]，其模型如图1所示。图1中，u\left(t\right)为阻尼器两端的相对位移，C_{\mathrm{0}}为零频率时的线性阻尼系数，K为“无限大”频率区域内的刚度系数，F_d(t)为恢复力。

Maxwell模型中恢复力为：

F_{\mathrm{d}}=C_{\mathrm{d}}{\mathit{\nu }}^{\mathit{a}}（1）

\left\{\begin{array}{l}{u}_{\mathrm{1}}\left(t\right)+u_{\mathrm{2}}\left(t\right)=u\left(t\right),\\ C_{\mathrm{0}}{\dot{u}}_{\mathrm{1}}\left(t\right)=K{u}_{\mathrm{2}}\left(t\right)=F_{\mathrm{d}}\left(t\right)\end{array}\right.（2）

式（1）、（2）中，\mathit{v}为速度向量，\mathit{a}为加速度向量，C_{\mathrm{d}}为阻尼系数，u_{\mathrm{1}}\left(t\right)、u_{\mathrm{2}}\left(t\right)分别为阻尼单元、弹簧单元，{\dot{u}}_{\mathrm{1}}\left(t\right)为u_{\mathrm{1}}\left(t\right)\mathrm{的导}\mathrm{数}，t为时间。

在进行非线性时程分析时，黏滞阻尼器的力学行为通常采用Maxwell模型进行数值模拟^[20]。当前主流有限元分析软件（包括本文采用的ABAQUS软件）均采用Maxwell模型，其优势在于能够准确地表征此类装置，且兼具刚度特性和频率依赖性的特殊力学响应。因此，在下文第2.1节基于全时程迭代的黏滞阻尼器减震结构地震响应分析中推导附加等效阻尼系数时均采用非线性模型中的阻尼力，其计算公式为F_{\mathrm{d}}=C_{\mathrm{d}}{\mathit{\nu }}^{\mathit{\alpha }}。

2 基于ABAQUS全时程迭代过程二次开发

2.1 黏滞阻尼器等效线性化

本文围绕非线性黏滞阻尼器对结构地震响应的影响展开，不考虑因结构非线性滞回特性所产生的附加等效阻尼比。假设减震结构的恢复力是线性的，则黏滞阻尼器减震结构的运动方程为：

\begin{array}{c}\mathit{M}\mathit{a}\left(t\right)+\mathit{C}\mathit{v}\left(t\right)+{\mathit{F}}_{\mathrm{d}}\left(\mathit{v}\left(t\right)\right)+\mathit{K}\mathit{d}\left(t\right)=\mathit{P}\left(t\right)\end{array}（3）

式中，\mathit{M}为减震结构的质量矩阵，\mathit{C}为减震结构的固有阻尼矩阵，{\mathit{F}}_{\mathrm{d}}为黏滞阻尼器的非线性阻尼力向量，\mathit{K}为刚度矩阵，\mathit{d}为位移向量，\mathit{P}为外部激励向量。

黏滞阻尼器减震结构模型如图2所示。图2中，n、m分别为平面框架的总楼层数、总跨数，{\theta }_{ij}为阻尼器D_ij 的轴线方向与水平方向的夹角，i、j分别为平面框架的楼层数、跨数。为方便求解黏滞阻尼器的非线性阻尼力向量，以图2的平面框架结构为例，只考虑各节点的水平移动，忽略竖向移动和转动，即只考虑水平的自由度。

图3为黏滞阻尼器减震结构的等效模型。

图3中，S_i 为等效刚度。根据等效线性化方法，结构每层均用一个水平的线性阻尼器D_i进行等效替换，令C_i为阻尼器的等效阻尼系数，K_i 为等效刚度系数。

将式（3）进行等效线性化，可得线性运动方程为：

\mathit{M}\mathit{a}\left(t\right)+C_{\mathrm{0}}\mathit{v}\left(t\right)+{\mathit{C}}_{\mathrm{a}}\mathit{v}\left(t\right)+\mathit{K}\mathit{d}\left(t\right)=\mathit{P}\left(t\right)（4）

式中，{\mathit{C}}_{\mathrm{a}}为附加等效阻尼矩阵，

{\mathit{C}}_{\mathrm{a}}=\left[\begin{array}{ccccc}{C}_{\mathrm{1}}+C_{\mathrm{2}}& -C_{\mathrm{1}}& \cdots & \mathrm{0}& \mathrm{0}\\ -C_{\mathrm{1}}& C_{\mathrm{2}}+C_{\mathrm{3}}& \cdots & \mathrm{0}& \mathrm{0}\\ ⋮& ⋮& \ddots & ⋮& ⋮\\ \mathrm{0}& \mathrm{0}& \cdots & C_{n-\mathrm{1}}+C_n& -C_n\\ \mathrm{0}& \mathrm{0}& \cdots & -C_n& C_n\end{array}\right]。

比较式（3）和（4），可得两种计算方法中黏滞阻尼器阻尼力的误差为：

\mathit{d}={\mathit{F}}_{\mathrm{d}}\left(\mathit{v}\right)-{\mathit{C}}_{\mathrm{a}}\mathit{v}（5）

要实现等效线性化，需根据误差最小的原则确定等效阻尼矩阵，为此需求解每层楼的最佳等效阻尼系数。由最小二乘法的思想可知，当误差平方的数学期望E最小时，可求得最佳的等效阻尼系数为：

E\left({\mathit{d}}^{\mathrm{T}}\mathit{d}\right)=E\left({\mathit{F}}_{\mathrm{d}}^{\mathrm{T}}{\mathit{F}}_{\mathrm{d}}+{\mathit{v}}^{\mathrm{T}}{\mathit{C}}_{\mathrm{a}}^{\mathrm{T}}{\mathit{C}}_{\mathrm{a}}\mathit{v}-\mathrm{2}{\mathit{F}}_{\mathrm{d}}^{\mathrm{T}}{\mathit{C}}_{\mathrm{a}}\mathit{v}\right)\doteq \mathrm{m}\mathrm{i}\mathrm{n}（6）

参考文献[2]经过推导可得，附加等效阻尼系数C_i 为：

C_i={\displaystyle \sum _{k=\mathrm{1}}^r}{\displaystyle \sum _{j=\mathrm{1}}^m}{C}_{ij}{\left(v_i^k\mathrm{s}\mathrm{i}\mathrm{n}{\theta }_{ij}\right)}^{\mathrm{1}+{\alpha }_{ij}}/{\displaystyle \sum _{k=\mathrm{1}}^r}{\left(v_i^k\right)}^{\mathrm{2}}（7）

式中，v_i^k为层间相对速度，C_{ij}为第i层第j跨非线性阻尼器的阻尼系数，\mathrm{s}\mathrm{i}\mathrm{n}{\theta }_{ij}为非线性阻尼器的安装角度，{\alpha }_{ij}为第i层第j跨非线性阻尼器的阻尼指数，r为积分过程总步数，k为积分步数。

一般黏滞阻尼器的阻尼指数取值范围为0.3~1.0，阻尼指数越小，结构的等效阻尼系数越大。由式（7）可知，当其他条件相同的情况下，减震结构的速度响应越小，对应的附加等效阻尼系数越大。根据已有资料可知：阻尼器的阻尼指数与阻尼系数对结构频率的影响很小，黏滞阻尼器对结构频率的影响很小；求得等效阻尼器与等效之前的非线性阻尼器对结构地震响应的影响相似。这些均为全时程迭代计算黏滞阻尼器减震结构地震响应的前提。

2.2 引入全时程迭代的等效方法

根据附加等效阻尼系数计算公式可知迭代计算过程：首先，假设已知每层楼的相对速度v_i，通过式（7）求得每层的附加等效阻尼系数C_i，再求得附加等效阻尼矩阵{\mathit{C}}_{\mathrm{a}}；然后，对等效线性化后的公式进行数值积分得到每层楼的相对速度v_i；最后，通过式（7）求得新的附加等效阻尼系数C_i，进入迭代计算。

虽然使用数值积分方法计算的结果较为准确，但如果自行计算或使用MATLAB编程计算仍会存在计算过程复杂、计算量大的问题。ABAQUS中计算结构地震响应的方式即为数值积分法，本文采用ABAQUS计算，再将最后计算得到的附加等效阻尼系数与一般的数值积分方法所得结果进行对比，验证本文方法的可行性。全时程迭代计算流程可简述为：

第1步：按式（4）建立黏滞阻尼器减震结构运动方程，并假设每层楼的初始附加等效阻尼系数均为0；

第2步：使用ABAQUS建立模型，输入地震波进行作业分析，输出减震结构的地震响应；

第3步：根据ABAQUS输出的结构速度响应，代入式（7）中计算新的附加等效阻尼系数；

第4步：将第3步新计算的附加等效阻尼系数输入ABAQUS模型，重复第2、3步，依次迭代，判断当连续两次迭代中减震结构地震响应最大值的差值绝对值小于允许值时，即可停止迭代得到最终的附加等效阻尼系数。

采用全时程迭代方法的优点为：1）采用ABAQUS进行减震结构的地震响应分析，避免使用数值积分方法的繁琐过程；2）比起数值积分方法的大量计算过程，本方法仅需迭代几次即可收敛；3）模型采用等效线性阻尼器，可避免非线性问题。

2.3 基于ABAQUS的二次开发

ABAQUS处理分析作业经历3个阶段，其中，前处理和后处理是实现二次开发的关键阶段^[21‒22]。由于实际中工程问题差异较大且偏重点不同，通用软件无法满足用户在专业方面的需求。为解决此种情况，ABAQUS提供了若干用户子程序接口，允许使用者通过子程序编写能满足具体功能的代码扩展主程序的功能，以此实现二次开发以帮助使用者提升工作效率^[23‒25]。目前，利用编程语言基于通用软件平台进行二次开发，是各领域研究的一个重要发展方向^[26‒27]。ABAQUS通过Python语言提供库函数，使用Python编写的函数便可在ABAQUS中直接进行材料特性定义、截面定义、建立模型、指派截面方向、划分网格、创建作业、提交作业、后处理分析作业等，可避免在ABA- QUS交互界面中手动操作，实现自动化处理，节省使用者大量时间，提高分析效率^[28‒30]。

2.3.1 实现过程

基于全时程迭代的方法进行计算，需要对输出的减震结构地震响应大量数据进行提取和运算，如果人工提取和计算，再将计算出的附加等效阻尼系数结果重新输回ABAQUS的模型，这个过程较为繁琐且无意义；迭代运算是基于循环运算，所以大量计算不可避免。因此，针对完成分析的作业后处理进行二次开发非常必要。后处理主要任务包括：提取分析作业结果中所需节点的时程响应（速度、位移等），对提取的数据进行清洗和处理，计算每层楼的附加等效阻尼系数，改变ABAQUS中每层楼的阻尼器系数，创建新作业并提交分析，自动实现迭代计算。当前后两次结构响应差值的绝对值小于允许值且迭代计算次数在设定的允许次数以内时，迭代结束并输出各层楼的最大层间相对速度和附加等效阻尼系数。

结合ABAQUS二次开发流程，再考虑在第2.2节中引入全时程迭代后得到的计算流程，可整理得到此次二次开发总体操作流程图，如图4所示。

2.3.2 主要代码

1）ABAQUS功能调用代码

图4中，需在第2步“在ABAQUS中建立结构模型”开始就实现输入代码后自动调用和计算。之后每一步骤涉及到ABAQUS中的操作都需实现自动调用相应功能。实现ABAQUS后处理的代码如附录A所示。附录表A1展示了部分实现调用ABAQUS或CAE中打开模型、打开输出数据库、打开模型阻尼器管理器、在阻尼器信息栏中输入阻尼系数、创建新作业并提交、打开新作业的odb文件，以及绘制所需节点的时程曲线并输出时程数据等功能的程序脚本，为接下来顺利进行阻尼系数计算公式的编写做好准备工作。

2）功能函数的编写

在以上代码的基础上进行等效阻尼系数的计算，需分析的节点单元较多且每个节点的时程响应数据较多，所以在计算附加等效阻尼系数时根据数据特点和运算流程编写相应的函数，以便提高计算效率，使整个代码更加简洁，增加可读性。以7层框架结构为例，解释主要函数功能。

a. getV函数

通过读取ABAQUS生成的数据文件，获取减震结构7层相对速度，具体实现方法如下：

首先，使用pandas库中的read_fwf函数读取ABA-QUS生成的存有减震结构7层顶部速度及地面运动速度的8个*. rpt文件。然后，使用pandas库中的drop系列函数清洗数据，去掉nan、重复出现的时间点及对应的速度。接着，重新生成索引，并转换其数据类型，从string转换为float。最后，求出最大速度。

b. ft_iterate函数

先定义初始阻尼c₀=10^‒10－10，通过getV得到一次迭代过程中减震结构7层顶部的相对速度时程和最大速度。然后，计算出附加等效阻尼系数。最后，根据阻尼器的角度，计算出下一次需要输入结构的附加等效阻尼系数。

c. creatName函数

通过输入自己自定义的部分uniqueName和迭代次数c_nt，生成一系列的变量名和函数名，如：ABAQUS作业名、生成的odb文件名、8个存储速度的rpt文件名和绘制速度的绘图变量名称。

d. changeC函数

通过输入要改变的模型名detailModelName、阻尼器所在的部分的名称changePartName、阻尼器名称changeCName、阻尼系数具体值c1_n、改变黏滞阻尼器的阻尼系数。

e. outPut函数

通过输入一系列通过creatName生成的文件名与变量名：odbName、vv0f、vv1f、vv2f、vv3f、vv4f、vv5f、vv6f、vv7f、xyPlotName，实现打开odb文件，绘制需要的点、V曲线，保存点的横纵坐标数据到rpt文件。

f. openModel函数

通过输入模型名称modelName，实现模型的打开。

g. mainProcess函数

此函数整合了其他所有函数，是本代码的核心函数。

首先，赋值最大速度的误差为1（即100%），最大迭代次数为10，调用openModel函数打开模型；然后，通过判断最大速度差百分比（即前后两次迭代速度差与其两者平均值的比值）是否小于设置的允许速度差百分比，迭代次数是否达到最大迭代次数，决定是否进行循环内的操作。循环内，先调用creatName函数生成所需要的变量名、文件名；再调用changeC改变阻尼器的阻尼系数；接着调用submitWork提交作业并进行计算；计算完成后调用outPut获取需要的数据并保存结果；再调用ft_iterate函数，在这个过程中，调用getV函数获取减震结构各层相对速度，计算非线性恢复力和新的附加等效阻尼系数，迭代次数+1，计算最大的速度差百分比，并打印相关参数的值。当前述循环结束后，判断迭代次数是否达到最大迭代次数，并输出有关信息。

附录附录表A2汇总了部分编写的函数及其说明。

根据编写好的函数，用mainProcess函数调用计算所需的所有函数。图5为函数运算流程图。

3 案例分析

为验证本文提出的ABAQUS二次开发的可靠性，选取叶正强等^[18]研究的结构实例作计算模型。此模型为黏滞阻尼器减震的钢筋混凝土框架结构，在验证计算中选取其中一榀框架进行相应的地震响应分析。平面框架及阻尼器布置如图6所示。图6中，结构每层层高为3.3 m，梁的截面尺寸为300 mm×600 mm，柱的截面尺寸为500 mm×500 mm。梁、柱材料的等效弹性模量取3.0×10¹⁰ Pa，等效质量密度为3 500 kg/m³。采用瑞雷阻尼计算结构的固有阻尼，附加阻尼器的阻尼系数均为134 kN·s/m，阻尼指数为0.6。地震波分别采用El‒Centro波和Taft波（以下分别简称E波、T波），峰值加速度P_GA取0.20g和0.62g^[31]。

3.1 模态对比

按照总体操作流程图4，首先，根据上述模型信息在ABAQUS中建立图6的7层平面框架结构，用以采用全时程迭代法计算结构地震响应。由于复杂的阻尼器模型通常不能直接在ABAQUS中定义，一般需要基于Fortran的二次开发过程实现，为简化本文方法的验证过程，选择在MATLAB中进行此7层框架结构的建模，以采用直接积分法计算结构响应。在验证迭代方法前先检查模型建立的正确性，展示两种软件中建模后同种模型的前5阶模态对比和误差，如表1所示。

由表1可知，前5阶模态的误差较小，尤其是第1、2阶的模态误差在0.1%以内，证明两者的建模非常接近，可开始后续迭代计算。

3.2 收敛性证明

以地震波为E波、P_GA取0.20g为例，每次迭代计算过程中最大层间相对速度如表2所示。

为了便于进一步分析，绘制最大层间相对速度随迭代次数的变化规律，如图7所示。

每次迭代计算过程中每层楼的附加等效阻尼系数如表3所示。附加等效阻尼系数随迭代次数的变化规律如图8所示。

由表3、图8可知：从第1次迭代到第2次迭代，每层楼的最大层间相对速度和附加等效阻尼系数变化较大；到第3次迭代时变化减小；至第4、5次迭代时，这些参数基本保持不变，表明全时程迭代法的收敛速度较快。

3.3 模型地震响应分析

全时程迭代法与直接积分法的结构顶层最大位移对比如表4所示。

表4中，两种计算方法的误差率控制在5%以内，精度较高。P_GA取0.62g时的误差均大于P_GA取0.20g时。同时，全时程迭代法得到的计算结果比直接积分法的偏小，其原因是在附加等效阻尼系数的计算公式中，较大的反应量对应的滞回曲线较为饱满而较小的反应量对应的滞回曲线较为扁平，因此附加等效阻尼系数计算结果偏大，最终导致全时程迭代法的计算结果偏小。

平面框架结构顶层位移响应如图9所示。

图9中，当P_GA分别为0.20g和0.62g时，在E波和T波两条地震波作用下，分别用直接积分法和全时程迭代法计算出此框架黏滞阻尼器减震结构顶层的位移响应。

由图9可知，本文提出的全时程迭代方法与直接积分方法所计算的7层框架结构顶层位移响应曲线重合度较高。

4 结 论

以等效线性化为基础推导出的全时程迭代方法，结合ABAQUS建模与后处理，提出基于全时程迭代的减震结构地震响应分析和ABAQUS二次开发以便更好地进行全时程迭代计算。以7层框架单榀结构为例进行分析，验证了其收敛性和精度。主要结论为：

1）全时程迭代方法是基于等效线性方法和多自由度结构推导的，具有充分的理论基础和物理意义。

2）全时程迭代方法计算减震结构地震响应分析收敛速度快、精度高，满足工程要求。

3）全时程迭代方法考虑了等效阻尼系数与结构相应的耦合作用。

4）ABAQUS二次开发的自动迭代计算程序可简便地实现全时程迭代法，促进了其在实际工程中的推广和应用。

附录见本刊网络版，扫描标题旁的二维码可阅读网络全文。

参考文献

[1] Liang Shahe, Zhang Zhiqiang, Li Aiqun.Seismic reduction principle and response analysis of variable damping viscous damper system[J].International Journal of Structural Integrity,2024,15(1):31‒60. doi:10.1108/ijsi-08-2023-0085

[2] Jia Chuanguo, Gou Yingqi, Li Jianguang,et al.Seismic resp-onse analysis method for structures with viscous dampers based on inter-history iteration[J].Journal of Vibration and Shock,2020,39(6):89‒97.

[3] 贾传果,苟英旗,李建广,等.基于全时程迭代的黏滞阻尼器减震结构地震响应分析方法研究[J].振动与冲击,2020,39(6):89‒97.

[4] Jia Chuanguo, Li Yutao, Su Hongchen,et al.The inter-history iteration method for seismic response analysis of seismically isolated structures and its secondary development ba-sed on ABAQUS[J].The Structural Design of Tall and Sp-ecial Buildings,2024,33(6):e2078. doi:10.1002/tal.2078

[5] Siami Kaleybar R, Tehrani P.Effects of using different arrangements and types of viscous dampers on seismic performance of intermediate steel moment frames in comparison with different passive dampers[J].Structures,2021,33:3382‒3396. doi:10.1016/j.istruc.2021.06.079

[6] Wang Zilong, Ren Wenjie.Research on optimization design for structure with viscous dampers based on genetic algorithm[J].World Earthquake Engineering,2021,37(2):123‒131. doi:10.3969/j.issn.1007-6069.2021.02.015

[7] 王子龙,任文杰.基于遗传算法的黏滞阻尼器减震结构优化设计研究[J].世界地震工程,2021,37(2):123‒131. doi:10.3969/j.issn.1007-6069.2021.02.015

[8] Chen Gaoyuan.Application research on viscous fluid dampers of a medical building in high intensity area[D].Zhenzhou:North China University of Water Resources and Electric Power,2023.

[9] 陈高源.液体黏滞阻尼器在高烈度区医疗建筑中的应用研究[D].郑州:华北水利水电大学,2023.

[10] Morro A.Modelling of viscoelastic materials and creep behaviour[J].Meccanica,2017,52(13):3015‒3021. doi:10.1007/s11012-016-0585-x

[11] Bochnia J, Blasiak S.Stress relaxation and creep of a poly-mer-aluminum composite produced through selective laser sintering[J].Polymers,2020,12(4):830. doi:10.3390/polym12040830

[12] Tyler R G.Preliminary tests on all energy absorbing element for braced structures under earthquake loading[J].Bulletin of the New Zealand National Society for Earthquake Engineering,1983,16(3):201‒212. doi:10.5459/bnzsee.16.3.201-212

[13] Tyler R G.Further notes on a steel energy-absorbing element for braced frameworks[J].Bulletin of the New Zeal and National Society for Earthquake Engineering,1985,18(3):270‒279. doi:10.5459/bnzsee.18.3.270-279

[14] Jia Chuanguo, Li Yutao, Su Hongchen,et al.Numerical simulation and experimental investigation for the dynamic rocking response of two stacked rigid blocks[J].Journal of Sound and Vibration,2024,569:118082. doi:10.1016/j.jsv.2023.118082

[15] Jia Chuanguo, Pan Jiafu, Yu Chenyang,et al.Research on the seismic response analysis method for the structures with metallic dampers based on inter-history iteration[J].Progress in Steel Building Structures,2021,23(3):64‒73.

[16] 贾传果,潘家富,于晨洋,等.基于全时程迭代的金属阻尼器减震结构地震响应分析方法研究[J].建筑钢结构进展,2021,23(3):64‒73.

[17] American Society of Civil Engineers. Minimum design loa- ds for buildings and other structures:A SCE7—10 [S].Was-hington D C:American Society of Civil Engineers,2010.

[18] Wu Xingyi, Vanapalli S K.Three-dimensional modeling of the mechanical behavior of a single pile in unsaturated expansive soils during infiltration[J].Computers and Geotechnics,2022,145:104696. doi:10.1016/j.compgeo.2022.104696

[19] Liang Jiali, Liang Jianwen.A user-defined element for dynamic analysis of saturated porous media in ABAQUS[J].Computers and Geotechnics,2020,126:103693. doi:10.1016/j.compgeo.2020.103693

[20] Luo Weigang, Tao Yuru, Song Yirong.The RC frame structure failure mode based on ABAQUS secondary developme-nt[J].Journal of Civil Engineering and Management,2016,33(5):1‒5.

[21] 罗维刚,陶昱儒,宋一容.基于ABAQUS二次开发的RC框架结构屈服机制[J].土木工程与管理学报,2016,33(5):1‒5.

[22] Gao Zhe, Yu Chaoshuai, Liu Zhaobin,et al.Application of Abaqus post-processing technique in diesel engine finite element analysis based on Python‒Matlab[J].Tractor & Fa-rm Transporter,2017,44(4):46‒49.

[23] 高喆,禹朝帅,刘钊宾,等.基于Python‒Matlab的Abaqus后处理技术在柴油机有限元分析中的应用[J].拖拉机与农用运输车,2017,44(4):46‒49.

[24] Ye Zhengqiang, Li Aiqun, Cheng Wenrang,et al.Study on vibration energy dissipation design of structures with fluid viscous dampers[J].Journal of Building Structures,2001,22(4):61‒66.

[25] 叶正强,李爱群,程文瀼,等.采用黏滞流体阻尼器的工程结构减振设计研究[J].建筑结构学报,2001,22(4):61‒66.

[26] Lu Lyan Ywan, Lin Ging Long,Shih Ming Hsiang.An expe-rimental study on a generalized Maxwell model for nonli-near viscoelastic dampers used in seismic isolation[J].Engineering Structures,2012,34:111‒123. doi:10.1016/j.engstruct.2011.09.012

[27] Du Yian.Study on performance of viscous damper applied to frame-supported tower[D].Tianjin:Tianjin University,2018.

[28] 杜怡安.应用于框架支撑式塔器的黏滞阻尼器性能研究[D].天津:天津大学,2018.

[29] Winowlin Jappes J T, Khan M A, Vignesh V,et al.Design and simulation of skate scooter made of composite material[J].Materials Today:Proceedings,2022,52:255‒259. doi:10.1016/j.matpr.2021.07.463

[30] Ganesh R, Dude D P, Kuna M,et al.ProCrackPlast:A finite element tool to simulate 3D fatigue crack growth under large plastic deformations[J].International Journal of Fracture,2023,243(1):65‒90. doi:10.1007/s10704-023-00732-9

[31] Ding Li, Yu Chao, Zhang Xu,et al.Incremental strain gradient plasticity model and torsion simulation of copper micro-wires[J].International Journal of Mechanical Sciences,2023,239:107891. doi:10.1016/j.ijmecsci.2022.107891

[32] He Miao, Li Qi, Li Xiaying.A new simulator for hydromechanical coupling analysis of injection-induced fault activ-ation[J].Geomechanics and Geophysics for Geo-energy and Geo-resources,2022,8(2):43. doi:10.1007/s40948-022-00353-x

[33] Wang Huakun, Gao Jing, Liu Tao,et al.Axial buckling behavior of H-piles considering mechanical-electrochemical interaction induced damage[J].Marine Structures,2022,83:103157. doi:10.1016/j.marstruc.2022.103157

[34] Gao Jingxiang, Cao Yongzhi, Wang Kaijie,et al.Polycrystalline simulation and experimental study of spatiotemporal anisotropy aluminum alloy irradiated by nanosecond laser[J].Optics & Laser Technology,2021,142:107244. doi:10.1016/j.optlastec.2021.107244

[35] Zhao Ying, Wang Yicheng, Hao Jibo,et al.Study on mechanical properties of cellular structures with negative Poisson's ratio based on the development of Abaqus plug-in tool[J].Composite Structures,2023,322:117348. doi:10.1016/j.compstruct.2023.117348

[36] Shen Yang, Moore R H, Deo A.Visualizing Abaqus output database in ParaView:A universal converter in Python and C++[J].SoftwareX,2023,22:101331. doi:10.1016/j.softx.2023.101331

[37] Gholami M, Afrasiab H, Baghestani A M,et al.Hygrothermal degradation of elastic properties of fiber reinforced co-mposites:A micro-scale finite element analysis[J].Composite Structures,2021,266:113819. doi:10.1016/j.compstruct.2021.113819

[38] Yuan Weihai, Wang Haocheng, Zhang Wei,et al.Particle finite element method implementation for large deformation analysis using Abaqus[J].Acta Geotechnica,2021,16(8):2449‒2462. doi:10.1007/s11440-020-01124-2

[39] 中华人民共和国国住房和城乡建设部. 建筑抗震设计规范:GB50011—2010 [S].北京:中国建筑工业出版社,2010.

基金资助

国家自然科学基金项目(51408080)