地震工程与工程振动

计算步骤	计算项目
（1）
（2）	Δx₁ =K′^TΔu₁
（3）
（4）	Δx₃ =K′Δx₂
（5）
（6）	ΔX=Δu₁+Δu₂

计算步骤	计算项目
（1）
（2）	Δx₁ =K′^TΔu₁
（3）
（4）	Δx₃ =K′Δx₂
（5）
（6）	ΔX=Δu₁+Δu₂

方法	总时间复杂度	控制方程求解时间/s
变刚度法	1.70×10¹²	193.4
本文方法	1.20×10¹¹	37.6

方法	总时间复杂度	控制方程求解时间/s
变刚度法	1.70×10¹²	193.4
本文方法	1.20×10¹¹	37.6

基于隔离非线性法-黏弹性人工边界的结构高效抗震分析方法

PDF下载

李佳龙 ¹ , 任东兴 ¹ , 李志远 ²

地震工程与工程振动 | 2024,44(4): 208-221

收起

地震工程与工程振动 | 2024, 44(4): 208-221

基于隔离非线性法-黏弹性人工边界的结构高效抗震分析方法

全屏

李佳龙¹, 任东兴¹, 李志远²

作者信息

^1.中冶成都勘察研究总院有限公司，四川成都 610023

^2.中国水利水电科学研究院工程抗震研究中心，北京 100038

李佳龙（1991—），男，工程师，博士，主要从事结构非线性分析等研究。E-mail：wekalee@163.com

Efficient seismic analysis method for structures based on the inelasticity-separated theory and viscous-spring artificial boundary

Jialong LI¹, Dongxing REN¹, Zhiyuan LI²

Affiliations

^1.Chengdu Surveying Geotechnical Research Institute Co., Ltd., of MCC, Chengdu 610023, China

^2.Earthquake Engineering Research Center, China Institute of Water Resources and Hydropower Research, Beijing 100038, China

doi: 10.13197/j.eeed.2024.0420

文章导航

摘要

收起

黏弹性人工边界是解决土-结构动力相互作用问题的有效手段之一，但其在分析过程中通常需要进行大量的非线性动力时程计算，效率较低。建立采用黏弹性人工边界求解土-结构动力相互作用问题控制方程求解的高效计算方法，为此，引入隔离非线性法高效求解思想，构造基于隔离非线性法-黏弹性人工边界的非线性土-结构动力相互作用分析模型，推导的动力控制方程在形式上与固定边界的隔离非线性动力控制方程基本一致，仅需将黏弹性人工边界的弹簧-阻尼器元件的刚度直接集成到近场土-结构模型的初始刚度矩阵与阻尼矩阵上即可。此外，结合Woodbury公式与组合近似法提出了改进Woodbury近似法，该改进方法能同时降低控制方程求解的时间复杂度与空间复杂度，实现了黏弹性人工边界求解土-结构动力相互作用问题动力控制方程的高效求解。研究方法同时保留了黏弹性人工边界与隔离非线性法的优点，数值算例验证了所提方法的正确性与高效性。

关键词

隔离非线性法 / 黏弹性人工边界 / 土-结构动力相互作用 / 改进Woodbury近似法

Abstract

收起

Viscous-spring artificial boundary is one of the effective means to solve the problem of dynamic soil-structure interaction, but it usually requires a lot of nonlinear dynamic time-history calculation in the process of analysis, and the efficiency is low. This paper aims to establish an efficient calculation method for solving the governing equations of dynamic soil-structure interaction problem with viscoelastic artificial boundary. Therefore, the efficiently inelasticity-separated solve idea is introduced to construct the nonlinear dynamic soil-structure interaction analysis model based on the inelasticity-separated finite element method and viscous-spring artificial boundary. The derived dynamic governing equation is basically same as the governing equation of the fixed boundary that only need to directly add the spring-damper stiffness of the viscous-spring artificial boundary into the initial stiffness matrix and damping matrix of the near-field soil-structure model. In addition, an improved Woodbury approximation method is proposed by combining Woodbury formula with the combined approximation approach, which can reduce the time and space complexity in solving governing equation, and realize the efficient solution of the dynamic governing equation of dynamic soil-structure interaction problem by using viscous-spring artificial boundary. The proposed method retains the advantages of viscous-spring artificial boundary and inelasticity-separated finite element method, and the correctness and efficiency are verified by a numerical example.

Key words

inelasticity-separated finite element method / viscous-spring artificial boundary / dynamic soil-structure interaction / improved Woodbury approximation approach

引用本文

李佳龙, 任东兴, 李志远. 基于隔离非线性法-黏弹性人工边界的结构高效抗震分析方法. 地震工程与工程振动, 2024 , 44 (4) : 208 -221 . DOI: 10.13197/j.eeed.2024.0420

Jialong LI, Dongxing REN, Zhiyuan LI. Efficient seismic analysis method for structures based on the inelasticity-separated theory and viscous-spring artificial boundary[J]. Earthquake Engineering and Engineering Dynamics, 2024 , 44 (4) : 208 -221 . DOI: 10.13197/j.eeed.2024.0420

正文

收起

0　引言

收起

近年来，随着各种形式复杂的大型建/构筑物大量涌现，这类结构在极端环境荷载作用下一旦发生破坏必将导致严重后果，明确这类结构的灾变机理已成为科研工作者与工程技术人员的迫切需求。地震作用下结构发生振动，同时由结构振动产生的惯性力如同新的地震反作用于场地土上从而引起新的场地震动，这一现象称之为土-结构动力相互作用或近场波动问题，该问题的求解是一个极具挑战性的课题，涉及到地震工程学、结构动力学、土动力学、材料力学以及计算机科学等众多学科的交叉^[1]。土-结构动力相互作用问题在数值分析时需考虑地震波在无限域地基中的辐射阻尼效应，常用策略是将结构和其附近土域组成近场以及外围几何无限土域组成远场，并在近场和远场交界面处设置人工边界。一般近场可采用有限元法或基于有限元法扩展得到的其他方法进行模拟，该部分可根据结构的材料属性呈现各种非线性特征，远场则通过设置合理的人工边界来模拟散射场的透射效应，并在截断边界处实现地震动的有效输入。近半个世纪以来，研究者们基于波动理论提出了多种人工边界条件，大致可分为时空耦联的全局人工边界与离散的局部人工边界两大类，全局人工边界主要包括边界元法^[2-3]和比例边界有限元法^[4]，局部人工边界主要包括无限元法^[5]、多次透射边界^[6]、完美匹配层^[7]以及黏弹性人工边界^[8-9]等。黏弹性人工边界通过在人工边界结点位置设置弹簧-阻尼器耗能元件来实现对散射波能量的吸收，该方法可以很好地处理辐射阻尼以及散射效应，物理概念清晰，程序实现简单，具有较高的计算精度与良好的稳定性，广泛应用于地震作用下土-结构动力相互作用分析^[10-15]。

目前，结构非线性分析一般采用有限元法的变刚度求解格式，一旦有单元进入非线性状态即对结构整体刚度矩阵进行更新并分解，随着结构规模的扩大以及有限元模型精细化程度的提高，其计算成本将急剧增加，计算效率偏低。为提高结构非线性分析的计算效率，LI等^[16]提出了一种结构非线性分析新方法，即隔离非线性有限元法（inelasticity-separated finite element method，IS-FEM）。该方法以有限单元法为基础，将有限元单元的材料应变分解为线弹性应变与非线性应变两部分，通过构造单元非线性应变场插值模型建立具有非线性隔离特征的单元控制方程，整体控制方程可实现结构整体刚度矩阵的弹塑性分离，无需对整体刚度矩阵进行实时更新与分解，采用Woodbury公式^[17]进行求解，可显著提高结构在发生局部材料非线性时的计算效率。目前，该方法已实现了对杆系结构^[18]、框架结构^[19]、二维平面^[20]以及三维实体结构^[21]的高效分析，其稳定性与高效性已得到验证。

鉴于此，本文将黏弹性人工边界与隔离非线性法相结合，以期提高黏弹性人工边界在求解土-结构动力相互作用问题时其控制方程求解的计算效率。考虑到黏弹性人工边界一般是以弹性围岩介质为基础构造的且在求解过程保持不变，这与隔离非线性法采用不变的初始弹性刚度矩阵进行求解的思想不谋而合，由此可建立考虑黏弹性人工边界的土-结构动力相互作用的隔离非线性运动微分方程，并基于Newmark-β法推导了隔离非线性动力控制方程。控制方程的最终形式与固定边界的土-结构模型（地震动以惯性力的形式进行加载）的控制方程基本一致，而黏弹性人工边界的弹簧-阻尼器元件的刚度可直接集成到土-结构模型的刚度矩阵与阻尼矩阵之上，地震动则通过自由场模型转变为截断边界结点上的等效结点力进行输入^[22-23]。考虑到Woodbury公式在求解非局部非线性问题时其占用的时间复杂度和空间复杂度均较大，不利于大规模问题的高效求解，为此本文将Woodbury公式与组合近似法^[24]相结合，提出了一种更为高效且使用范围更广的改进Woodbury近似法。改进方法在构造基向量时直接采用矩阵K_inf的初始形式进行计算，从而避免了计算并保存整体矩阵K_INF这一过程，降低了原Woodbury公式求解的空间复杂度；此外，为避免缩减基矩阵病态而使计算精度受到影响，提出了基向量快速正交化计算方法，该正交化过程计算量极小，可忽略不计。基于时间复杂度理论建立了改进方法的时间复杂度模型，研究表明：改进算法的时间复杂度与结构的非线性自由度数呈线性阶关系，降低了非线性问题求解的时间复杂度，对于大规模问题其计算效率相较于传统变刚度法具有明显优势。最后，以一重力坝为例验证本文所提算法的正确性与高效性。

1　隔离非线性法基本原理

收起

隔离非线性有限元法将单元的材料应变分解为线弹性应变与非线性应变两部分，通过构造单元非线性应变插值场函数并结合虚功原理建立具有非线性隔离性质的单元控制方程，实现了结构整体刚度矩阵的弹塑性分离，通过引入Woodbury公式对其求解，可显著提高结构在发生局部材料非线性时的计算效率^[17]。以图1所示的隔离非线性多边形比例边界有限元平面模型^[25]为例，在任意计算步可建立计算域的整体控制方程为

(1)

式中：K_e（维数为N×N）为结构的初始刚度矩阵，其在整个计算中保持不变，其中N为结构的整体位移自由度数；K′（维数为N×p）为结构的系数矩阵，该矩阵为稀疏矩阵，其中p为结构的非线性计算自由度数；

（维数为p×p）为包含结构材料非线性信息的分块对角稀疏矩阵，该矩阵需要在分析过程中根据结构的非线性状态变化而在计算过程中实时更新；ΔF为结构的荷载增量；

为结构中所有进入非线性状态的插值点的非线性应变增量向量。

将式（1）中的非线性应变增量向量

进行静力凝聚，可将其改写为

(2)

式中，矩阵

可看作是结构初始弹性刚度矩阵K_e的修正矩阵。基于Woodbury公式^[17]，待求位移增量ΔX可表示为

(3)

式中，

(4)

式中，矩阵

为p×p阶的满阵，在数学上可称其为初始刚度矩阵K_e的舒尔补矩阵。矩阵K_inf是仅与结构非线性区域位置相关的满阵，可在非线性分析前计算出整个结构对应矩阵K_INF并保存，非线性计算时提取相应部分即可。采用从右至左的计算顺序可将式（3）的求解过程依次转换为矩阵与向量的运算，共计6个计算步，如表1所示。

由表1可知，计算步骤（1）、（5）为初始刚度矩阵的回代，计算步骤（2）、（4）为稀疏的系数矩阵K′与向量的乘积，计算步骤（6）则为向量的加法运算，这5个计算步骤的计算量都较小，限制Woodbury公式高效计算主要在于步骤（3）的计算，因为该计算步骤涉及到一个以p×p阶矩阵为系数矩阵的方程组求解。对于局部非线性问题（即p≪N），直接采用Woodbury公式求解具有极高的效率优势。

对于地震作用下的结构非线性问题，采用固定边界与惯性力加载形式，隔离非线性法的动力控制方程可表示为

(5)

式中：

为结构在t时刻的等效荷载向量；

为结构的初始等效刚度矩阵，两者的具体形式为

(6)

(7)

式中：M与C分别为结构的质量矩阵与阻尼矩阵；Δt为时间步长；τ为地震动转换向量；

为t时刻的地震动加速度增量向量；Ẋ（t_k-1）与Ẍ（t_k-1）分别为结构在k-1个时间步的相对速度和相对加速度向量。采用常刚度的Rayleigh阻尼矩阵C可使等效刚度矩阵

在整个过程中保持不变。式（6）所示的动力控制方程同样可采用Woodbury公式进行求解，其求解过程与静力问题的求解过程一致。有关隔离非线性法理论更详细的介绍可以参考文献[26]。

2　考虑黏弹性人工边界的隔离非线性动力分析方法

收起

2.1　黏弹性人工边界

将隔离非线性法与黏弹性人工边界相结合，利用隔离非线性法在非线性求解方面的效率优势对土-结构动力相互作用问题的控制方程进行高效求解。以图2所示的二维平面问题为例，对黏弹性人工边界进行简要说明。数值分析时引入人工边界截去土体无限域远场部分，将近场域模型与引入的人工边界一起统称为土-结构系统。

近场域模型可采用有限元法或比例边界有限元法进行模拟，内部结点采用符号i表示，人工边界处的结点采用符号b表示，地震作用下非线性有限域的运动微分方程可表示为

(8)

式中：Mⁱⁱ、M^ib、M^bi与M^bb为近场域的质量矩阵；Cⁱⁱ、C^ib、C^bi与C^bb为近场域的阻尼矩阵；

、

与

为近场域的实时切线刚度矩阵；Xⁱ（Ẋⁱ、Ẍⁱ）与X^b（Ẋ^b、Ẍ^b）分别为结构内部结点与边界结点的位移（速度、加速度）向量；F^b为地震作用下被截去的无限土域对近场域的作用力。

人工边界处的总反应由散射场和自由场两部分组成，基于黏弹性人工边界理论^[27]可将式（8）所示的动力微分方程改写为

(9)

式中，K^b与C^b分别为黏弹性人工边界的刚度矩阵与阻尼矩阵，两者均为所有边界结点上的弹簧与阻尼器系数在法线与切线2个方向集合，

(10)

式中，K^b，j与C^b，j为边界结点j （j∈b）对应的刚度矩阵与阻尼矩阵，两者均为对角矩阵。当边界外法线方向与x轴平行时，

(11)

当边界外法线与y轴平行时，

(12)

式（11）与式（12）中所示的弹簧刚度系数

、

以及阻尼器系数

、

的计算见图2，图中参数：G为围岩材料的剪切模量，λ为拉梅常数，ρ为质量密度，R为散射源到人工边界结点的距离，c_p为纵波波速，c_s为剪切波速；α与β为修正系数，一般取α=0.8与β=1.1即可取得一个较为稳定的计算结果；A_j为边界结点j代表的等效线面积。P_g为作用于人工边界上的等效地震荷载，其形式为

(13)

可以看出，等效地震荷载P_g由三部分组成，这三部分均与自由场相关而与散射场无关，也就是说，如果对自由场进行分析获得了其人工边界处的作用力

、位移

以及速度

，那么以该自由场为近场土域的土-结构系统在各时刻的等效地震荷载也就随之确定。上述黏弹性人工边界条件是基于弹性介质条件推导而来的，已有分析表明在介质塑性区发展不是很大的情况下，采用弹性介质中构造的黏弹性边界不会带来很大误差，可满足工程精度要求^[28]。

2.2　隔离非线性动力控制方程

采用增量形式，将式（9）所示土-结构系统的运动微分方程改写为

(14)

式中：ΔXⁱ（ΔẊⁱ、ΔẌⁱ）与ΔX^b（ΔẊ^b、ΔẌ^b）分别为结构内部结点与边界结点的位移（速度、加速度）增量；ΔP_g为通过无结构规则自由场地震反应获得的等效地震荷载增量，即等效结点力增量。式（14）中等号左边第三项代表的是土-结构系统的恢复力变化项ΔF，可将其拆分为结构与近场土域的恢复力变化ΔF₁以及黏弹性人工边界引起的恢复力变化ΔF₂两项，即

(15)

结合隔离非线性法控制方程式（17）展开式的第一个表达式，结构与近场土域的恢复力变化ΔF₁又可表示为

(16)

式中：

、

与

为结构与近场土域的初始刚度矩阵；K′ⁱ与K′^b分别为内部结点i与边界结点b对应的系数向量。将式（16）代入式（15）后再代入式（14），可建立如方程（17）形式的考虑土-结构动力相互作用的隔离非线性运动微分方程为

(17)

式中，

(18)

式中：M为土-结构系统的质量矩阵；

与

分别为土-结构系统的阻尼矩阵与初始弹性刚度矩阵；

为土-结构系统在地震作用下所受的等效结点力增量；ΔX（ΔẊ、ΔẌ）为土-结构系统的位移（速度、加速度）增量向量。在式（17）中，刚度矩阵

在整个分析过程中保持不变，而土-结构系统的非线性效应则采用等号右边的虚拟荷载项

来表示。研究表明：在地震作用下基于黏弹性人工边界的土-结构动力相互作用问题，采用隔离非线性法可将非线性的运动微分方程转变为考虑非线性虚拟荷载的弹性结构运动微分方程。

采用Newmark-β平均加速度法求解式（17）为

(19)

式中：

为土-结构系统的等效弹性刚度矩阵；

为等效荷载向量，两者的具体表达分别为

(20)

(21)

式中：Δt为时间步长；Ẋ（t_k-₁）与MẌ（t_k-₁）分别为土-结构系统在第k-1个时间步的速度和加速度向量。若采用常阻尼矩阵且时间步长Δt在整个分析过程中保持不变，则等效弹性刚度矩阵

可保持不变。将式（19）中虚拟荷载项

移到等式左边，可将其改写为

(22)

式中：ΔX与

均为未知，可见式（22）需再建立一个补充方程方可实现求解。由于隔离非线性法控制方程（1）的第二个表达式表示的是单元内的平衡关系，与边界形式无关，联立该表达式可建立考虑黏弹性人工边界的隔离非线性法的动力分析控制方程为

(23)

可以看出，上式（23）与考虑惯性力加载的隔离非线性动力控制方程在形式上基本一致，仅是在等效初始刚度

中附加了黏弹性人工边界的刚度矩阵K^b与阻尼矩阵C^b，由于构造的边界刚度K^b与阻尼矩阵C^b的对角特性，因此构造的等效初始刚度矩阵

同样具有稀疏带宽特性。同样的，可采用Woodbury公式对式（23）进行求解，迭代求解过程与静力问题的迭代求解过程基本一致，在此不做过多赘述。

2.3　分析流程

图3给出了采用隔离非线性法求解土-结构动力相互作用问题的整体分析流程，主要过程可分为人工边界构造、边界等效结点力计算以及迭代求解三部分。人工边界根据弹性围岩的材料参数进行构造，这一过程并不需要对黏弹性边界进行修改，其只是借用该边界模型来处理土-结构相互作用问题，因此人工边界的准确性与稳定性完全取决于已有黏弹性边界模型的参数选择，等效结点力则采用波场分离法通过无自重自由场模型进行计算。通过式（8）~式（23）的推导可以看出，基于隔离非线性方法的求解格式对基于黏弹性人工边界的土-结构相互作用控制方程没有额外修正，因此黏弹性人工边界的应用范围广（能够方便地反映介质非均质、几何不规则、软土夹杂、捕虏体等地基复杂因素）、计算程序简单、数值稳定性高、时域直接求解满足非线性分析要求等优点，基于隔离非线性法的计算模型也都具备，仅需要近场模拟时在地基模型中反映各种复杂因素即可。本文虽然是以二维单一地基模型为例构造的隔离非线性动力求解格式，但理论上对于三维以及分层或者任意分层地基模型均适用。此外，对于任意入射角输入的地震动，相关研究已经获得了这类情况的等效结点力计算方法^[22]，可直接应用。迭代求解过程与静力问题的迭代求解过程基本一致，采用Woodbury公式求解具有与完全Newton-Raphson迭代法同样的二阶收敛速度，但考虑到Woodbury公式使用的限制，迭代求解时建议直接采用本文第三章提出的改进Woodbury近似法进行求解，其准确性与高效性并不会因为增加黏弹性人工边界而改变。

3　改进Woodbury高效求解方法

收起

通过上述分析，可以看出不论是静力问题，亦或固定边界或考虑黏弹性人工边界的动力问题，整个计算域的整体控制方程基本一致，且控制方程均可采用Woodbury公式进行求解，唯一的差别是采用初始刚度矩阵还是等效初始刚度矩阵进行计算。然而，由于隔离非线性法在单元模型构造时在单元中引入了较多的非线性应变插值点，这使得当结构的非线性区域较大时采用Woodbury公式求解的时间复杂度与空间复杂度占用较高，效率优势并不明显。这是因为Woodbury公式在求解时需提前计算并保存矩阵K_inf，而矩阵K_inf对应的整体矩阵K_INF是一个满阵，这对于大型结构而言计算并保存整体矩阵K_INF所需要的存储空间是无法接受的；此外，Woodbury公式在求解时其时间复杂度与结构的非线性自由度数p呈立方阶关系，即O（p³）^[29]，当非线性自由度数较高时（即非局部非线性问题），其计算效率反而不如传统变刚度法。

鉴于此，本文将Woodbury公式与组合近似法^[24]相结合，提出一种高效的近似求解方法，称之为改进Woodbury近似法（improved Woodbury approximation approach，IWAA）。本文所提出的改进Woodbury近似法首先通过近似求解策略获得矩阵

的基向量，后又通过Gram-Schmidt法获得正交化基向量，从而保证计算精度。下面以式（1）所示的静力控制方程为例对该近似求解算法为例进行说明，对于动力问题仅需将初始刚度矩阵K_e更换为初始等效刚度

即可，而求解时的迭代过程均一致。

3.1　近似求解

组合近似法^[24]（combined approximation approach，CA法）作为一种目前最为流行的结构重分析方法，该方法将二项式级数展开和多项式拟合相结合，利用局部近似效率上的优势和全局近似精度上的优势估算结构响应。下面采用组合近似法对Woodbury公式的计算步骤（3）进行求解（表1），该计算步骤可改写为

(24)

基于组合近似法的求解策略，待求未知量Δx₂可转换为若干个基向量线性组合的形式，即

(25)

式中：下标s为选择的基向量个数；R_B =[r₁，r₂，…，r_s]是由s个基向量组成的基向量矩阵，其阶数为p×s；y={y₁，y₂，…，y_s}^T为各基向量对应系数组成的系数向量列阵。一般地，基向量个数s将远远小于舒尔补矩阵的维数，即s≪p。将式（25）代入式（24）并在两端同时左乘矩阵

得

(26)

式中：K_R为s×s阶的缩减基矩阵；F_R为s×1阶的缩减荷载向量，两者分别为

(27)

(28)

也就是说，对于Woodbury公式计算步骤（3）的求解，如果获得了基向量矩阵R_B以及系数向量y，则待求未知量Δx₂即可获得。基于组合近似法的求解流程，式（25）所示的基向量的构造过程可表述为

(29)

对于式（29），在构造r₂至r_s+₁的基向量过程中，不再将矩阵K_inf看作是一个整体，而是保留其原有形式进行分步计算，因此无须提前计算并保存整体矩阵K_INF，可大大降低求解的前处理过程并降低对计算机存储性能的需求，即算法的空间复杂度得到降低。为保证基向量之间线性无关，需满足式（30）所示关系：

(30)

一般地，设定E_er=0.95即可满足计算精度要求。由于式（29）构造基向量过程中存在过渡变量K″_pr1=Δx₁以及

，每次基向量计算时将这些过渡变量存储起来，从而在计算缩减基矩阵K_R时无须重复计算。K_R的最终形式可表示为

(31)

其中的各元素分别为

(32)

由于舒尔补矩阵

为对称矩阵，可以证明缩减基矩阵K_R同样为对称矩阵，即元素

。根据缩减基矩阵R_B构造缩减荷载向量F_R，则待求的系数向量y可表示为

，将该系数向量y与基向量矩阵R_B代入式（25），则待求未知量Δx₂即可表示为

(33)

采用上述方法，可将Woodbury公式的计算步骤（3）p×p阶的线性方程组求解转换为一个s×s阶的线性方程组求解，可明显减小问题的计算规模，降低计算量。

3.2　正交优化

为避免式（31）所示的缩减基矩阵病态而使得式（33）所示缩减方程的求解精度受到影响，可利用Gram-Schmidt法对3.1节所述中获得的基向量r_i（i=1，2，…，s）进行正交化处理^[24]，该计算过程可表述为

(34)

(35)

(36)

式中，

为舒尔补矩阵；

与q_i分别为第i个非标准化与标准化基向量。采用上述表达虽可直接进行基向量正交化计算，但过程复杂且需进行大型舒尔补矩阵与向量的重复计算，可借用4.1节已有的计算结果进行计算。

采用归纳总结法可获得标准化基向量q_i的计算过程。第i个非标准化基向量

与标准化基向量q_i可分别表示为

(37)

式中：R_i =[r₁ r₂…r_i ]为由前i个基向量组成的基向量矩阵；L_i ={l_i，₁，l_i，₂，…，l_i，i}^T与H_i ={h_i，₁，h_i，₂，…，h_i，i}^T分别为第i个非标准化基向量

以及标准化基向量q_i对应的系数向量，l_i，j与h_i，j则分别为

与q_i中原基向量r_j对应的系数，其中j≤i≤s。L_i与H_i的通项公式则分别为

(38)

(39)

式中：变量

为缩减基矩阵K_R中第i行中前1到j个元素组成的行向量；

为缩减基矩阵K_R中前i行、前i列元素组成的子矩阵。对于系数向量L_i，始终有l_i，j=1，上述通项公式对所有基向量均使用，因此不用再将第一个标准化基向量看作是一个特例进行计算。可以看出，L_i与H_i这2个系数向量的计算是层层递进且相互耦合的，但又仅与缩减基矩阵K_R相关而与原基向量无关。也就是说，如果获得了缩减基矩阵K_R，那么通过简单运算即可获得全部非标准化基向量L_i与标准化基向量的系数向量H_i。

当基向量正交化后，式（24）所示的待求未知量Δx₂可重新表示为

(40)

式中：Q_B=[q₁q₂…q_s]为由s个标准化基向量组成的p×s阶基向量矩阵；u={u₁，u₂，…，u_s}^T为各标准化基向量对应系数组成的系数向量。将式（40）代入式（24），并同时在等式两端左乘矩阵

并考虑到标准化基向量的正交性，可知标准化基向量的系数向量

，那么未知变量Δx₂可表示为

(41)

综上所述，正交化后的基向量q_i可直接用以构造以舒尔补矩阵为系数矩阵的方程组的解Δx₂，若要通过增加基向量的个数来提高求解精度时，只需通过式（29）构造新的基向量r_k，然后再通过式（32）计算新增基向量r_k对应缩减基矩阵K_R中的元素并通过式（38）与式（39）构造系数向量H_k，从而可通过式（37）获得新增基向量对应的正交化向量q_k。基于上述操作，可直接在式（41）计算结果的基础上进行修改，无须重新求解，同时保证了计算精度和求解效率。

3.3　时间复杂度

时间复杂度是一种能客观评价某种算法计算效率的工具，可以定量评估算法在计算机中执行所需要的理论时间，不受计算机性能、编程技巧等外在因素的影响。采用传统变刚度法求解结构非线性问题的控制方程（采用LDL^T分解法求解），其一次迭代求解的时间复杂度T_L约为^[30]

(42)

式中，m为刚度矩阵带宽，若对结构自由度编号进行优化，则刚度矩阵的半带宽一般可取m_k = N¹^/2。

经统计，改进Woodbury近似法一次迭代求解主要有s+2次常刚度矩阵回代，s+1次稀疏矩阵K′（以及矩阵K′^T、矩阵K″_p的逆）与向量的乘积，以及一个s×p维的矩阵与一个p×s维的矩阵乘积。对改进方法的时间复杂度进行统计，其在一个迭代步骤的时间复杂度约为

(43)

式中，参数α与β分别为结构控制方程中系数矩阵K′和非线性状态矩阵K″_p中每列非零元素个数的最大值，且参数α、β取值越大其对应的时间复杂度也就高。可以看出，本文所提出改进方法的时间复杂度与结构的非线性自由度数呈线性阶关系。已有研究表明，通常选取较少基向量也可获得较为满意的计算结果，对于小规模问题且矩阵修改量适中的情况，选取2~3个基向量即可，对于中等规模问题且矩阵修改量较大时需要5~6个基向量，而对于大规模修改问题需要9~10个基向量。

将结构的非线性自由度p与结点自由度数N的比值定义为非线性自由度占比γ，即γ=p/N，该参数可用以衡量结构进入非线性状态程度。为更直观地对比IWAA相较于传统变刚度法的计算效率提升情况，将IWAA在一个迭代步的时间复杂度T_IWAA与传统变刚度法求解的时间复杂度T_L的比值定义为时间复杂度比率η，即η=T_IWAA/T_L。若η<1表示其计算效率高于传统变刚度法，η=1表示两者效率相当，η>1则表示该算法效率不及变刚度法。以多边形比例边界单元模型为例，图4给出了在不同非线性自由度占比γ（γ=0.01、0.10、0.50、1.00、5.00、10.00）以及不同基向量s（s=3、6、10）情况下，其时间复杂度比率η随结点自由度数目N的变化情况。可以看出，在任意基向量情况下，当非线性自由度占比γ<0.1时，该IWAA对于任意规模结构均具有高效性；当非线性自由度占比γ逐渐增加时，其具有高效性的临界结点自由度数N虽然也在增加，当N>3 000时（图4（c）所示：γ=10，s=10），IWAA相较于传统变刚度法均具有高效性。当结点自由数目N=10⁵时（取γ=10，s=10），IWAA的时间复杂度不到传统变刚度法的1/10；而当N=10⁸（取γ=10，s=10）时，其时间复杂度则仅为变刚度法的1/1000。也就是说，当结点自由度数目N越大，IWAA相较于传统变刚度法的计算效率提升也就越高，具有明显的效率优势。

4　数值验证

收起

选取一建于弹性基岩上的二维混凝土重力坝为例进行地震反应分析，分析时基岩选择为300 m×100 m的范围，计算模型整体几何尺寸如图5（a）所示。本例采用多边形比例边界单元^[25]对大坝与地基进行网格划分，其中重力坝共计13666个多边形单元，弹性地基共计1 448个单元，分别如图5（b）与（c）所示。重力坝材料参数：弹性模量E_e=30.0 GPa，泊松比为v=0.2，密度ρ=2630 kg/m³，抗压强度为24.1 MPa，抗拉强度为2.9 MPa；基岩考虑为弹性，其材料参数：弹性模量E_e=20.0 GPa，泊松比v=0.2，密度ρ=2 500 kg/m³。使用双向El Centro波对结构进行非线性地震反应分析，采用Drucker-Prager模型^[31-32]模拟大坝的非线性效应，为了获得大坝更为明显的非线性反应，将地震动强度调整为原始强度的一倍，其中地震持续时间为30 s，分析时步长选择为0.02 s，共计1500个荷载步。

4.1　边界验证

首先计算用于该问题分析的等效结点力并验证所构造的黏弹性人工边界的正确性。以图5（c）所示的300 m×100 m无自重地基作为自由场模型，基于波场分离法将调整后的双向地震动加速度记录转换为等效结点力施加于具有黏弹性人工边界的自由场模型上并进行弹性动力有限元分析。图6为自由场模型顶部中点在水平与竖向2个方向的位移时程响应，其中解析解采用传递矩阵法^[33]进行计算，可以看出解析解与数值解完全吻合，表明构造的黏弹性人工边界正确，通过自由场模型获得的等效结点力可直接用于本问题分析。

4.2　计算结果分析

动力分析时，首先施加一个代表大坝与地基自重的静力荷载，将结构自重引起的单元应力作为初应力代入计算，由自重引起的结构位移以及单元应变则不考虑，然后再在此基础上进行动力非线性分析，从而保证地应力平衡。分别采用变刚度法（即使直接采用Newmark-β法进行求解）与本文所提算法对该土-结构动力相互作用模型进行计算。上述分析已说明基于自由场模型获得的边界结点处的等效结点力正确，那么将该等效结点力进行变刚度求解可认为是该算例的精确解。

图7~图9给出了该重力坝顶部的绝对位移、速度以及加速度时程曲线，可以看出本文算法与变刚度法的计算结果基本完全吻合，相对误差极小，完全满足工程精度要求。此外，由图7可知，考虑上部混凝土大坝后绝对位移变化曲线与自由场曲线大致趋势一致，自由场模型位移曲线相对光滑，而图7锯齿现象明显，这是由于大坝自身相较于地面本身振动而引起的。由于大坝较高且地震动强度增加了一倍，因为鞭梢效应，顶部的最大加速度反应较强。

通过分析可知，大坝发生非线性行为主要集中在前5 s，图10为大坝前5 s非线性发展过程，可以明显看出，当地震动强度较大时大坝颈部明显出现塑性贯穿区，在实际结构中可能已出现贯穿裂缝。本算例采用DP模型来模拟混凝土的非线性行为，虽然该模型无法考虑混凝土开裂后的非线性行为，但也可在一定程度上模拟混凝土的抗拉压行为，具有一定的参考意义。

采用本文方法与变刚度法求解的累计时间复杂度对比，如表2所示。由表可知，本文方法相较于传统变刚度法降低超过93%。基于MATLAB平台，变刚度法求解控制方程的累计计算时间约为193.4 s，而采用本文方法求解整体控制方程的求解时间约为37.6 s，减少了约80%，不论是从算法理论还是实际求解时间方面均显示出本文方法的高效性。

5　结论

收起

本文将隔离非线性法与黏弹性人工边界相结合，提出了黏弹性人工边界在求解土-结构动力相互作用问题时控制方程求解的高效计算方法。基于弹性围岩介质构造的并联式弹簧-阻尼器元件在分析过程中可保持弹簧刚度与阻尼器系数不变，这与隔离非线性法采用不变的初始刚度矩阵进行分析不谋而合，为此建立了考虑黏弹性人工边界的隔离非线性土-结构系统的运动微分方程，并通过Newmark-β法推导了动力分析控制方程。此外，考虑到已有隔离非线性法在求解控制方程时占用的时间复杂度与空间复杂度均较大，为此本文将Woodbury公式与组合近似法相结合，提出了改进Woodbury近似法，拓展了隔离非线性法的应用范围。得出如下结论：

1）考虑黏弹性人工边界的隔离非线性动力控制方程与固定边界加惯性力输入的隔离非线性动力控制方程在表现形式上基本一致，不同点仅在于等效初始刚度矩阵附加了黏弹性人工边界的刚度矩阵项与阻尼矩阵项，而等效刚度矩阵可在整个计算过程中保持不变且仍具有稀疏带宽特性，结构系数矩阵以及非线性状态矩阵则与黏弹性人工边界无关。

2）提出的改进Woodbury近似法用于求解隔离非线性控制方程，可降低Woodbury公式求解的时间复杂度与空间复杂度，同时具备良好的计算精度与收敛性，使得隔离非线性法不再局限于局部非线性问题的求解，对于非局部非线性问题甚至是全局非线性问题均适用。

3）本文工作虽然是以二维平面问题为例进行分析，但是对于三维问题的求解理论上同样适用；此外，本文虽然是以单一地基模型为例进行分析，理论上该求解思想对于多层或任意成层地基模型以及任意入射角度的地震波输入均适用。

基金

收起

国家自然科学基金项目(52038002; 52008075)
国家自然科学青年基金项目(52209168; 52209157)
四川省自然科学基金项目(2024NSFSC0156)

参考文献

收起

文献

收起

参考文献引证文献

排序方式：

[1]

宋崇民，渠艳龄，刘磊，等. 土-结构动力相互作用远场问题数值分析方法综述[J]. 水力发电学报，2019，38(9):1-17.

SONG

Chongmin

, QU

Yanling

, LIU

Lei

, et al. A review of numerical methods for far-field modeling in dynamic soil-structure interaction systems[J]. Journal of Hydroelectric Engineering, 2019, 38(9): 1-17. (in Chinese)

[2]

BANERJEE

P K

, BUTTERFIELD

. Boundary element methods in engineering science[M]. London: McGraw-Hill, 1981.

[3]

LACHAT

J C

, WATSON

J O

. Effective numerical treatment of boundary integral equations: A formulation for three-dimensional elastostatics[J]. International Journal for Numerical Methods in Engineering, 1976, 10(5): 991-1005.

[4]

SONG

C M

. The scaled boundary finite element method: Introduction to theory and implement[M]. New Jersey: Wiley, 2018.

[5]

UNGLESS

R F

. An infinite finite element[D]. Vancouver: University of British Columbia, 1973.

[6]

廖振鹏，黄孔亮，杨柏坡，等. 暂态波透射边界[J]. 中国科学A辑，1984，14(6):556-564.

LIAO

Zhenpeng

, HUANG

Kongliang

, YANG

Baipo

, et al. Transient wave transmission boundary[J]. Science in China, Ser A, 1984, 14(6): 556-564.(in Chinese)

[7]

BERENGER

J P

. A perfectly matched layer for the absorption of electromagnetic waves[J]. Journal of Computational Physics, 1994, 114(2): 185-200.

[8]

刘晶波，王振宇，杜修力，等. 波动问题中的三维时域粘弹性人工边界[J]. 工程力学，2005，22(6):46-51.

LIU

Jingbo

, WANG

Zhenyu

, DU

Xiuli

, et al. Three-dimensional visco-elastic artificial boundaries in time domain for wave motion problems[J]. Engineering Mechanics, 2005, 22(6): 46-51. (in Chinese)

[9]

赵密，王利涛，高志懂，等. 土-结构相互作用分析中几类地震输入方式的比较[J]. 震灾防御技术，2017，12(3):589-598.

ZHAO

, WANG

Litao

, GAO

Zhidong

, et al. Comparison of earthquake input methods in soil-structure interaction analysis[J]. Technology for Earthquake Disaster Prevention, 2017, 12(3): 589-598. (in Chinese)

[10]

赵密，苏成坤，王丕光，等. 地震作用下水-结构-土动力相互作用分析[J]. 工程力学，2022，39(3):51-63.

ZHAO

, SU

Chengkun

, WANG

Piguang

, et al. Analasis of water-structure-soil dynamic interaction under earthquakes[J]. Engineering Mechanics, 2022, 39(3): 51-63. (in Chinese)

[11]

赵密，欧阳文龙，黄景琦，等.P波作用下跨断层隧道轴线地震响应分析[J]. 岩土力学，2019，40(9):3645-3655.

ZHAO

, OUYANG

Wenlong

, HUANG

Jingqi

, et al. Analysis of axis dynamic response of rock tunnels through fault fracture zone under P waves of earthquake[J]. Rock and Soil Mechanics, 2019, 40(9): 3645-3655. (in Chinese)

[12]

许成顺，戴金，豆鹏飞，等. 液化场地-群桩-上部结构动力相互作用简化分析方法[J]. 地震工程与工程振动，2020，40(6):25-35.

Chengshun

, DAI

Jin

, DOU

Pengfei

, et al. Simplified analysis method for dynamic interaction of liquefiedsite-group pile-superstructure[J]. Earthquake Engineering and Engineering Dynamics, 2020, 40(6): 25-35. (in Chinese)

[13]

文学章，戴进峰. 瑞利波作用下桩筏基础的动力响应[J]. 地震工程与工程振动，2019，39(2):177-182.

WEN

Xuezhang

, DAI

Jinfeng

. Kinematic response of piled raft foundations to Rayleigh waves[J]. Earthquake Engineering and Engineering Dynamics, 2019, 39(2): 177-182. (in Chinese)

[14]

梁建文，杨晶丽，韩冰，等. 土体非线性对土-结构动力相互作用的影响———以Millikan图书馆为例[J]. 地震工程与工程振动，2022，42(3):1-10.

LIANG

Jianwen

, YANG

Jingli

, HAN

Bing

, et al. Effects of the soil nonlinearity on dynamic soil-structure interaction: A case study of Millikan library building[J]. Earthquake Engineering and Engineering Dynamics, 2022, 42(3): 1-10. (in Chinese)

[15]

陈少林，朱学江，赵宇昕，等. 考虑土骨架非线性的饱和土-结构相互作用分析[J]. 地震工程与工程振动，2019，39(1):114-127.

CHEN

Shaolin

, ZHU

Xuejiang

, ZHAO

Yuxin

, et al. Analysis of saturated soil-structure interaction considering soil skeleton nonlinearity[J]. Earthquake Engineering and Engineering Dynamics, 2019, 39(1): 114-127. (in Chinese)

[16]

, YU

D H

. Efficient inelasticity-separated finite-element method for material nonlinearity analysis[J]. Journal of Engineering Mechanics, 2018, 144(4): 0001426.

[17]

WOODBURY

M A

. Inverting modified matrices[R]. Princeton: Princeton University, 1950.

[18]

靳永强. 基于隔离非线性理论的杆系结构弹塑性大位移高效分析方法[D]. 大连：大连理工大学，2021.

JIN

Yongqiang

. Efficient analysis method for skeletal structures with elasto-plastic large displacement based on the inelasticity-separated theory[D]. Dalian: Dalian University of Technology, 2021. (in Chinese)

[19]

余丁浩，李钢，李宏男. 采用自适应Woodbury公式的框架结构高效动力分析与易损性计算方法[J]. 振动工程学报，2023，36(1):32-43.

Dinghao

, LI

Gang

, LI

Hongnan

. Efficient structural dynamic analysis and fragility solution methods for frame using adaptive Woodbury method[J]. Journal of Vibration Engineering, 2023, 36(1): 32-43. (in Chinese)

[20]

李佳龙，李钢，于龙. 隔离非线性平面应变单元模型及其在Drucker-Prager模型中的应用[J]. 岩土力学，2020，41(5):1492-1501, 1509.

Jialong

, LI

Gang

, YU

Long

. Inelasticity-separated plane-strain element model and its application to Drucker-Prager model[J]. Rock and Soil Mechanics, 2020, 41(5): 1492-1501, 1509. (in Chinese)

[21]

李佳龙，李钢，李宏男. 基于隔离非线性的实体单元模型与计算效率分析[J]. 工程力学，2019，36(9):40-49, 59.

Jialong

, LI

Gang

, LI

Hongnan

. The inelasticity-separated solid element model and computational efficiency analysis[J]. Engineering Mechanics, 2019, 36(9): 40-49, 59. (in Chinese)

[22]

杜修力. 工程波动理论与方法[M]. 北京：科学出版社，2009.

Xiuli

. Theories and methods of wave motion for engineering[M]. Beijing: Science Press, 2009. (in Chinese)

[23]

HUANG

J Q

, ZHAO

, DU

X L

. Non-linear seismic responses of tunnels within normal fault ground under obliquely incident P waves[J]. Tunnelling and Underground Space Technology, 2017, 61: 26-39.

[24]

KIRSCH

. Reanalysis of structures: A unified approach for linear, nonlinear, static, and dynamic systems[M]. Dordrecht, the Netherlands: Springer, 2008.

[25]

李佳龙，李钢，余丁浩. 多边形比例边界有限元非线性高效分析方法[J]. 工程力学，2020，37(9):8-17.

Jialong

, LI

Gang

, YU

Dinghao

. Polygon scaled boundary finite element nonlinear efficient analysis method[J]. Engineering Mechanics, 2020, 37(9): 8-17. (in Chinese)

[26]

余丁浩. 隔离非线性有限元法及在RC结构抗震分析中的应用[D]. 大连：大连理工大学，2019.

Dinghao

. Inelasticity-separated finite element method and its application in seismic response analysis of RC structurs[D]. Dalian: Dalian University of Technology, 2019. (in Chinese)

[27]

杜修力，赵密，王进廷. 近场波动模拟的人工应力边界条件[J]. 力学学报，2006，38(1):49-56.

Xiuli

, ZHAO

, WANG

Jinting

. A stress artificial boundary in FEA for near-field wave problem[J]. Chinese Journal of Theoretical and Applied Mechanics, 2006, 38(1): 49-56. (in Chinese)

[28]

吴艺，房营光. 弹塑性地基中黏性与黏弹性人工边界条件有效性的验证[J]. 岩石力学与工程学报，2006，25(增刊2):3468-3473.

, FANG

Yingguang

. Validity of viscous and viscoplastic artificial boundary conditions for elastoplastic foundation[J]. Chinese Journal of Rock Mechanics and Engineering, 2006, 25(S2): 3468-3473. (in Chinese)

[29]

贾硕，李钢，李宏男. 基于Woodbury非线性方法的迭代算法对比分析[J]. 工程力学，2019，36(8):16-29, 58.

JIA

Shuo

, LI

Gang

, LI

Hongnan

. Comparative analysis of the iterative algorithms for nonlinear method based on the Woodbury formula[J]. Engineering Mechanics, 2019, 36(8): 16-29, 58. (in Chinese)

[30]

GOLUB

G H

, VAN LOAN

C F

. Matrix computations[M]. 4th ed. Beijing: Posts & Telecom Press, 2014.

[31]

JIANG

J F

, WU

Y F

. Identification of material parameters for Drucker-Prager plasticity model for FRP confined circular concrete columns[J]. International Journal of Solids and Structures, 2012, 49(3/4): 445-456.

[32]

LEE

. Large deflections of cantilever beams of non-linear elastic material under a combined loading[J]. International Journal of Non-Linear Mechanics, 2002, 37(3): 439-443.

[33]

KENNETT

B L N

. Seismic wave propagation in stratified media[M]. Cambridge: Cambridge University Press, 1983.

2024年第44卷第4期

PDF下载

引用本文

BibTeX

文章信息

doi: 10.13197/j.eeed.2024.0420

接收时间：2023-02-07
首发时间：2026-03-30

补充材料

相关文章

文章信息

作者

出版历史

收稿日期：2023-02-07
修回日期：2023-05-11

基金

国家自然科学基金项目(52038002; 52008075)

国家自然科学青年基金项目(52209168; 52209157)

四川省自然科学基金项目(2024NSFSC0156)

作者信息

^1.中冶成都勘察研究总院有限公司，四川成都 610023

^2.中国水利水电科学研究院工程抗震研究中心，北京 100038

参考文献

分享链接

https://castjournals.cast.org.cn/joweb/dzgcygczd/CN/10.13197/j.eeed.2024.0420

分享至

全文二维码

扫描看全文

引用本文

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