考虑土-结相互作用的大型结构高效地震分析方法

胡静静，余丁浩，李 钢，王 睿，张 晗，苏 璞

(大连理工大学海岸和近海工程国家重点实验室，辽宁，大连 116024)

随着经济社会的迅速发展和科学技术的不断突破，近年来我国逐渐涌现大量的大型复杂工程结构，如大型城市综合体、地下核电站、超高层建筑物等。这些结构具有如下特点[1]：① 体量大且构成复杂；② 土与结构接触部位刚度突变明显、相互作用效应显著且不容忽略。这些大型工程结构在地震下的安全问题关系着国计民生[2]，因此对此类结构进行地震反应分析具有重大意义，特别是开展考虑土-大型结构相互作用的高效非线性分析方法研究显得尤为重要。

目前，考虑土-结动力相互作用的数值分析方法主要有直接法和子结构法[3]。子结构法是先把地基和上部结构分别离散分析，再根据接触界面上的连续性条件考虑两个子结构间的相互作用，该方法的优点是灵活性大、计算量小，国内外相关学者已基于此方法进行了大量研究[4-5]。然而，该方法前期主要限于线性或者等效线性问题求解，为将其扩展到非线性地震反应分析方面，WOLF[6]总结了在时域中进行土-结相互作用非线性分析的步骤，提出线性和非线性体系均适用的频域-时域混合分析程序，但对于非线性体系来说该方法计算结果仅是一种有限近似；傅敏红等[7]基于子结构法提出了数值解和解析解耦合的土-结相互作用求解方法，可将其拓展到非线性计算领域，但该方法求解过程中基础和地基土必须分别假定为刚性和弹性，因而其适用性有限。由于子结构法是在叠加原理的基础上建立的，仅可用于某些特定条件下的非线性地震反应分析，这使得此类方法的应用受到较大限制[8]。直接法也称整体有限元法，该方法将结构与较大范围的近场地基土域看作整体进行建模分析，可以更完整、真实的考虑结构在地震作用下的动力响应，能够很好的模拟土与结构界面的分离、滑移等强动力非线性接触行为。整体有限元法[9]已成为处理土-结相互作用问题的最有效方法，该方法最关键的环节是如何处理土与地下结构之间接触面的动力学行为。地震作用下，土与周围结构的变形不一致，土-结接触面上发生的相对位移、张开、重新闭合等强接触非线性行为[10]直接影响接触面附近土体与结构的动力响应。目前有限元法在处理土与结构接触面的动力接触行为时主要有两种方法[11]：一是接触力学法，即通过直接定义不同介质之间接触表面对的力学传递特性，建立接触面力传递的力学模型和接触方程，通过接触算法求解接触方程，代表性求解算法有Lagrange 乘子法[12-13]、Penalty法[14-15]等；另一种是接触面单元法，即通过在两种接触的介质之间建立接触面单元，同时赋予接触面单元特殊的应力-应变关系来模拟接触面的力学行为，代表性方法有Goodman 单元[16]和薄层单元[17]等。随着有限元技术的不断突破，接触面单元法凭借其概念简单、能够与普通单元一样进行程序化实现等优势已逐渐发展为土-结接触面数值分析的主流方法[18]。

随着工程结构规模的扩大，采用整体有限元法求解考虑土-结相互作用的大型结构非线性地震反应时，由于需同时对主体结构、大范围土域和接触面进行建模，其单元数与结点自由度数通常非常巨大，这就对计算效率提出了更高要求。目前，结构非线性分析主要采用变刚度求解方法，该方法需要在迭代求解过程中对结构整体刚度矩阵进行实时更新和分解，对于大型结构而言，大规模刚度矩阵的时变特性使得非线性分析计算效率较低。为提高非线性分析效率，众多学者分别从不同角度提出了多种加速计算方法，如方程组高效求解方法、快速迭代方法、并行计算方法等，这些方法或者通过提高大规模线性方程组求解速度来改善计算效率，或者通过减少刚度矩阵更新次数来提高效率，或者通过降低方程组规模来提高计算效率。地震作用下结构通常并非全部进入非线性状态，而是仅在某些局部区域出现非线性，而结构的其他大部分区域则始终保持为弹性状态，众多学者基于工程结构的局部非线性特征提出了多种高效分析方法，主要有多尺度方法[19-21]、子结构方法[22-24]、重分析方法[25-26]、拟力法[27-29]等，然而，这些方法适用性通常受限，难以直接用来进行考虑土-结构相互作用的大型结构地震反应分析。近年来，通过将拟力法中的变形分解思想进行拓展，李钢等提出了一种求解局部材料非线性问题的新型高效通用有限元算法，即隔离非线性有限元法[30-33]，该方法通过对材料应变进行分解并在单元内构造非线性应变场函数，能够在非线性迭代求解过程中保持结构弹性刚度矩阵不变，避免了传统方法由于结构整体刚度矩阵时变[34-35]导致的计算效率低下，在保证精细化模拟的同时大幅提高计算效率，目前已开发出基于该方法的多种单元模型，如纤维梁单元[36]、壳单元[37]、实体单元[38]等，但还未涉及土-结构相互作用问题的高效非线性分析。虽然土-结构接触区域在地震作用下可能出现较大范围的强接触非线性变形，但此区域通常仅占整体计算模型的小部分，因此，考虑土-结构相互作用的非线性分析问题具有显著的局部非线性特征，结合适当的局部非线性算法有助于显著提升此类问题的分析效率。

本文提出一种考虑土-结相互作用的大型结构高效地震反应分析方法，该方法以隔离非线性有限元法为理论基础，采用无厚度三维Goodman 单元模拟土-结接触行为，通过将单元相对位移分解为初始弹性和非线性两部分，并引入额外非线性自由度建立单元非线性相对位移场，建立了能够模拟土-结构相互作用的隔离非线性接触面单元控制方程，随后对主体结构、土体和土-结接触模型的控制方程进行集成建立了隔离非线性的整体结构动力控制方程，最后引入Woodbury 公式求解，在非线性分析时仅对代表结构局部非线性行为的小规模Schur 补矩阵进行不断更新和分解，大大提高了结构非线性分析效率，数值算例结果验证了本文方法的有效性和高效性。

1 隔离非线性的土-结相互作用计算模型

为实现土-结构接触行为的高效模拟，本节基于变形分解思想建立隔离非线性的接触面单元模型，考虑到Goodman 单元简单易于实现，并且与界面本构模型结合可较好地模拟土-结构接触区域的滑移、张开和压密等接触非线性行为，目前在工程中得到广泛应用[39-43]，本文以该单元为基础对土-结构的接触非线性行为进行描述。

1.1 基本假设

八结点Goodman 接触面单元是由两个四结点平面组成，如图1 所示面Ⅰ和面Ⅱ，单元厚度e=0，图中ξηζ 代表单元等参坐标系，两接触面之间任意一接触点对(点1 和点2，以下简称点对)的受力变形关系可认为通过一组剪切方向和法向弹簧进行模拟，通过对这些弹簧定义适当的应力-相对位移本构关系即可实现接触面间水平相对错动、法向张开和闭合等行为的计算分析。

图1 Goodman 单元工作原理Fig.1 Principle of the Goodman element

如图2(a)所示，空间八结点无厚度Goodman接触面单元的两个接触面对应结点的坐标完全相同，图中xyz代表整体坐标系，即：

图2 三维Goodman 单元Fig.2 Three-dimensional Goodman element

式中：xi、yi、zi是接触面Ⅰ (结构表面)的单元结点i(i=1, 2, 3, 4)在整体坐标系内的坐标值；xi+4、yi+4、zi+4是接触面Ⅱ (土体表面)的单元结点i+4(i=1, 2, 3, 4)在整体坐标系内的坐标值。平面问题的Goodman 单元本质上是一维单元[44]，因而本文中三维Goodman 单元本质上可以看作某种特殊的二维单元，基于此可将空间八结点接触面单元进行等参处理，如图2(b)所示，其中ξηζ 代表等参坐标系，等参坐标系下单元形函数可表示为：

式中，ξi和ηi分别为结点i(i=1, 2, 3, 4)对应的等参坐标值，即ξi、ηi=±1。

为确定整体坐标系xyz与等参坐标系ξηζ 之间的对应关系，将(ξ,η,ζ)空间内形状规则的单元映射为(x,y,z)空间的形状扭曲的单元，需引入如下等参变换方程：

八结点接触面单元每个结点具有3 个自由度，在迭代过程中的任一线性化增量求解步中，单元的结点位移增量(整体坐标系下)可表示为：

单元两片接触面之间任意点对的相对位移增量为Δu=[ΔuΔvΔw]T，其可由单元变形矩阵B和结点位移增量Δq表示：

式中，I为3×3 的单位矩阵。

在接触面单元中，需考虑沿其剪切方向和法向的相对位移和应力，取接触面单元上任一点沿接触面的法线方向为ye轴，如图2(b)所示，建立局部坐标系xeyeze，xeze即为过该点的切平面[45]。在局部坐标系下，单元相对位移增量表示为Δue=[ΔueΔveΔwe]T，其中Δue和Δve对应局部xe向和ze向的剪切相对位移分量，用来描述土-结接触面之间的相对错动变形；Δwe对应局部ye向的法向相对位移分量，用来描述土-结接触面之间的张开、闭合等变形。局部坐标系下的Δue可由如下表达式计算：

式中，L为坐标转换矩阵。令LB=B'，则有：

在局部坐标系xeyeze下，设ksη和ksξ为当前计算步下单元中任意点对在xe向和ze向(剪切方向)的切线刚度系数，kζ为该点对ye向(法向)的切线刚度系数，则局部坐标系下该点的增量应力-相对位移关系可表示为：

记作：

式中：Δσi=[ΔτηΔτξΔσζ]T为应力增量向量，其中，Δτη和Δτξ分别为局部xe向和ze向的剪应力分量，Δσζ为局部ye向的正应力分量；Dt为该点对在当前计算步下的瞬时切线刚度矩阵。

1.2 相对位移分解

在任一线性化增量迭代计算步中，接触面单元中任意点对的相对位移增量向量Δue可根据该点的初始弹性刚度分解为初始弹性和非线性两部分：

相对位移增量向量中包含2 个剪切相对位移分量Δue、Δve和法向相对位移分量Δwe(即Δue=[ΔueΔveΔwe])。单元两个剪切方向的行为并无耦合，但法向行为会对剪切方向行为产生影响。图3 和图4分别为单元法向和剪切方向的典型应力-相对位移关系图，当两接触面闭合受压时(法向相对位移we＜0)，单元法向行为可通过一较大的刚度kp进行模拟(如图3 所示)，此时单元存在剪切摩擦行为，相应的剪切应力-相对位移关系如图4 所示；当两接触面打开时(法向相对位移we＞0)，单元法向正应力和剪切应力均为0 值。

图3 单元法向应力-相对位移关系Fig.3 Normal stress-relative displacement of Goodman element

图4 单元剪切应力-相对位移关系Fig.4 Shear stress-relative displacement of Goodman element

图5 剪切相对位移分解Fig.5 Decomposition of shear relative displacement

由于法线方向上的应力-相对位移关系在原点附近存在拐点，导致单元内任意点对在不同加载方向上初始刚度不唯一，不同初始接触状态下单元变形分解的具体实现方式存在差异，相应的，对于接触面非线性状态的判定标准也不相同。为此，本文建立了不同初始接触状态下的单元位移分解方法，具体论述如下：

1) 当单元中某点对的法向初始状态为受压，且在当前计算步下其法向依然受压，但两个剪切方向的瞬时切线刚度偏离其初始刚度，则认为该点对在剪切方向进入非线性状态，即在剪切方向产生非零的非线性相对位移，但法向依然处于初始弹性状态，即法向的非线性相对位移分量为0，此时单元非线性相对位移增量表示为：

应当指出的是，本文方法中初始刚度主要用来作为变形分解的依据，从而使得在计算求解过程中整体刚度保持不变。当考虑卸载刚度退化或强度下降时，本文方法依然适用，此时，每个增量步中的变形增量均可以初始弹性刚度为基准分解为两部分，分解出的非线性变形将为非零值，对应的非线性自由度将依然处于激活状态，并参与控制方程和相应求解公式的构造与计算。

2) 当单元中某个点对的法向初始状态为受压，且在当前计算步下其法向变为受拉，则此时该点对三个方向的切线刚度均偏离其初始刚度(均为0)，因而均会判定为进入非线性状态，相应非线性相对位移增量的三个分量均为非零值。图6给出了此时法向相对位移增量的分解示意图，假设当前增量计算步下其应力状态从a点变化到b点，则对应的非线性相对位移增量计算表达式为：

图6 法向相对位移分解Fig.6 Decomposition of normal relative displacement

3) 当单元中某点对的法向初始状态为受拉，且在当前计算步下其法向依然受拉，则此时该点对的三个方向均处为零刚度状态，与初始状态相同，因而均认为处于初始弹性状态，相应非线性相对位移增量向量中的三个分量均为0，即=[0 0 0]。

4) 当单元中某点对的法向初始状态为受拉，且在当前计算步下其法向变为受压，则该点对三个方向的瞬时切线刚度均不再是初始状态时的零刚度，因而均认为进入非线性状态，相应各分量的非线性相对位移计算表达式的形式与式(13)和式(14)相同，为防止初始零刚度引起的计算奇异问题，本文在计算时使用一个接近于0 的小数近似代替各方向初始时刻的刚度值。

本文中对于非线性相对位移的定义是相对于初始状态而言的，代表了单元任意点对当前状态与其初始状态的数值偏离程度，其作用主要是为了能够始终使用恒定的初始刚度对单元任意时刻的应力状态进行描述，以便于下文推导出能够使结构弹性刚度矩阵在分析过程中始终保持不变的隔离非线性控制方程。

1.3 单元控制方程

任意状态下单元的应力增量场σin可以用单元弹性相对位移增量场与弹性矩阵De的乘积来计算，结合式(10)可将其进一步表达为：

相较于传统有限元方法，隔离非线性方法需在单元中额外设置若干个非线性变形插值点，以便于对单元非线性变形场进行模拟，本文以单元高斯点作为该单元的非线性相对位移插值点，并在每个接触单元中设置4 个高斯积分点，如图7所示，则单元的非线性相对位移场可通过如下插值形式给出：

图7 单元非线性变形插值点Fig.7 Inelastic deformation interpolation point of the element

式中，I3为3×3 的单位矩阵。

将式(7)和式(16)代入式(15)可将任意状态下单元的应力增量场Δσin进一步表示为：

在计算分析过程中，对于弹性单元，其相应的应力增量可表示为单元初始刚度矩阵与单元相对位移的乘积，而对于进入非线性状态的单元，其任意点对在任意线性化增量迭代步中的应力增量Δσin又可通过材料瞬时切线本构矩阵Dt与总相对位移增量Δue乘积的形式求得(如式(9)所示)，结合式(9)和式(10)可得出如下表达式：

将式(11)与式(16)代入式(19)，整理后可将单元应力增量场的计算表达式重写为如下形式：

引入虚功原理建立单元控制方程，其中本文接触面单元的虚功方程为：

根据矩阵运算相关规则，将式(23)进行整理，即可建立接触面单元的控制方程：

其中：

式中：ke为单元的初始弹性刚度矩阵，在计算过程中始终不变，其具体数值与单元初始状态有关；k′和为与单元非线性有关的系数矩阵，可通过高斯积分方法进行数值求解。由于本文令非线性变形插值点与高斯积分点重合，可以证明插值函数矩阵C满足Kronecker delta 性质[46]，因此将单元的第j个高斯积分点坐标代入非线性变形插值函数Ci后可得：

式中：Dt,i代表当前计算步下单元非线性插值点处的材料切线本构矩阵；Hi为第i个插值点处的积分权系数；为第i个插值点处的变形矩阵；|J|i为插值点处的雅可比矩阵行列式[45]，其计算表达式为：

式中：

在实际计算过程中，可能存在某些插值点的全部或个别非线性变形分量被判定处于初始弹性状态的情况，如1.2 节相对位移分解情况1)和情况3)所述，此时这些插值点中相应非线性变形分量为0，在单元控制方程推导时应将这些变形分量直接剔除，在实际程序实现时，只需要将单元控制方程式(27)中相应非线性自由度及系数矩阵k′和中对应行、列直接删除即可，这一过程可由程序根据单元状态自动完成，当整个单元都处于线弹性状态时，单元控制方程将变为keΔq=ΔFe。

2 考虑土-结相互作用结构地震反应分析

本文采用整体式模拟方法建立考虑土-结构相互作用的结构分析模型，以图8 所示框架剪力墙高层结构为例进行说明，计算模型中包含主体结构、土体和土-结构接触区域三部分，其中，主体结构的梁柱构件采用隔离非线性纤维梁单元[47]进行模拟；主体结构中的剪力墙采用隔离非线性壳单元[37]进行模拟；土体部分采用隔离非线性实体单元模型[48]进行模拟，该模型将材料应变分解为线弹性应变与非线性应变两部分，并通过插值形式构造非线性应变场，可实现土体部分局部非线性行为的隔离表达和高效分析；土-结构接触区域采用本文前述建立的隔离非线性接触单元进行模拟。通过将不同单元的控制方程进行集成，可得如下形式的整体计算模型控制方程：

图8 隔离非线性单元模型Fig.8 Element model of inelasticity-separated method

式中：Ke为结构整体初始弹性刚度矩阵，在计算过程中保持恒定不变；K'和由各单元控制方程中的系数矩阵k'和集成而来；ΔQ为结构的结点位移增量；ΔF为结构的荷载增量；集合了整个结构中所有非线性单元中的非线性插值点的非零非线性变形。K'与材料的非线性程度无关，仅与产生非零增量非线性变形的区域位置有关，矩阵不仅与产生非零非线性变形的区域位置有关，还与材料的非线性本构模型相关。结构控制方程表现为整体弹性刚度矩阵Ke和代表局部非线性变形的刚度项K'和相互分离的特征，设结构整体自由度数为n、非线性自由度数为m，则矩阵Ke、K'和的规模分别为n×n、n×m、m×m，由于控制方程中非线性相关系数矩阵在集成时仅考虑了局部区域内进入非线性的相关单元，m代表了结构非线性区域的规模，其值通常远小于n(即m＜＜n)，这表明控制方程式(31)中非线性自由度和相关系数矩阵的规模通常极小。

对于考虑土体-结构相互作用的整体结构动力非线性分析，可列出增量形式的运动微分方程为：

式中：M为结构质量矩阵；C为结构阻尼矩阵；ΔQ¨(t) 、 ΔQ˙(t)分别为结构各时刻的相对加速度以及速度； ΔPg(t)为地震作用力；ΔF(t)为结构的恢复力向量。将式(31)代入式(32)，并结合Newmarkβ 法求解运动微分方程，可基于隔离非线性法建立考虑土体-结构相互作用的整体结构动力分析控制方程：

其中：

算法实现流程如图9 所示。可以看出，在每次迭代时的单元状态确定过程中，需对各单元的非线性状态进行判定，这可通过将单元每个非线性插值点处的瞬时切线刚度与其初始弹性刚度对比实现，对于某个单元，若全部插值点处的瞬时切线刚度与初始弹性刚度相等，就表明该单元在当前增量计算步下处于初始弹性状态，否则就代表该单元在当前计算步下处于非线性状态。

图9 算法流程图Fig.9 Algorithm flow chart

此外，在每次迭代时，根据确定的非线性单元位置和数量，可直接确定非线性自由度数目，从而确定求解方程中非线性相关系数矩阵的规模及其中的非零元位置和数量，本文方法采用Fortran 编程，利用动态内存分配机制实时开辟和销毁相关系数矩阵所需存储空间。可以看出，相较于传统方法，本文方法虽然计算效率更高，但需增加额外的存储开销，然而，由于计算求解过程中涉及到的非线性相关系数矩阵均具有高度的稀疏性[34]，且本文方法采用仅存储非零元的行压缩格式(compressed sparse row，CSR)进行矩阵存储，因此额外增加的存储开销并不显著。

3 数值验证

为验证本文方法准确性和高效性，建立如图10所示的接触摩擦算例模型，上部块体作用有给定竖向荷载P1和单调递增的水平荷载P2，考虑两个计算工况，其中工况1 中P1荷载为50 kN，工况2中P1荷载为100 kN。上部块体与下部块体均采用实体单元模拟，弹性模量分别为2.1×1010Pa 和5×107Pa，泊松比均为0.3。上、下块体接触面采用本文建立的接触面单元模拟，接触本构采用Clough-Duncan 双曲线模型[50]，具体本构模型参数取值详见表1。

表1 接触面本构参数Table 1 Constitutive parameters of interface

图10 接触摩擦块模型 /mFig.10 Contact friction block model

采用本文方法对该模型进行分析，两个工况下计算所得接触面平均剪切应力-相对位移曲线如图11 所示，图中与理论解进行了对比，其中理论解为基于Clough-Duncan 模型计算出的接触面极限剪切应力(即剪切相对位移u→∞时的剪切应力值)，表2 给出数值解(即本文方法加载至最大位移处对应的平均剪应力)与理论解的相对误差。由图11 和表2 可见，本文方法计算结果与理论解吻合良好。

表2 本文方法数值解与理论解对比Table 2 Comparison between numerical solution and theoretical solution

图11 剪切应力-相对位移曲线Fig.11 Shear stress-relative displacement curve

为了验证本文方法的高效性，利用有限元计算软件ABAQUS 对该算例进行分析，模型规模相同条件下两个工况计算时间分别为344 s 和342 s，而本文方法的计算时间分别为35 s 和47 s，仅为ABAQUS 的10.2%和13.7%。

4 数值算例

4.1 算例概况

采用本文方法对某大型城市综合体结构进行考虑土体-结构相互作用的地震反应分析，结构类型为框架核心筒体系，结构示意图如图12 所示。该结构地上31 层、地下2 层，裙房地上部分共5 层，地上结构高度99.45 m，地下埋深-10.5 m，平面尺寸约为126 m×67.2 m，建筑总面积为115 506.72 m2。截取的地基土域平面尺寸为378 m×201.6 m，深度取60 m。对结构进行双向地震动加载，为使分析过程中模型结构出现较为显著的非线性变形，以便于对本文方法进行验证，x方向和z方向的地震波峰值加速度分别为2.2 m/s2和1.87 m/s2，两个方向的地震加速度时程如图13 所示。

图12 整体土-结有限元模型Fig.12 soil-structure finite element model

图13 地震波加速度记录Fig.13 Acceleration record of the earthquake

该综合体的整体土-结计算分析模型如图12 所示，其中主体结构的梁、柱构件采用隔离非线性纤维梁单元模拟，其单元数为80 658；剪力墙结构采用隔离非线性壳单元模拟，其单元数为76 391；土体采用隔离非线性实体单元模型模拟，其单元数为138 160；土-结相互作用采用本文隔离非线性接触面单元模拟，其单元数为1116。整体模型的单元总数为296 325，总自由度数为615 513。模型中梁、柱构件的混凝土采用修改后的修正Kent-Park 本构模型[34]进行模拟，剪力墙混凝土材料采用修正斜压场模型，混凝土材料参数见表3；钢筋采用Dodd 和Cooke[51]提出的本构模型；地基土域的材料采用Drucker-Prager 模型[52-53]，钢筋及土体相关材料参数见表4；土域边界采用滚轴边界进行处理[54]；土体与结构接触面的本构关系选用Clough-Duncan 双曲线模型[50]，其中界面摩擦角取值为31°，其他参数取值均同表1。本文算例模型对于桩基础的模拟采用“复合式分析”方法，即将桩弥散于土体中，把群桩基础视为桩、土合一的复合材料[55]，以降低建模和分析的复杂度。此外，本文算例中，主体结构的地下部分外侧均为框架-剪力墙，其中梁柱构件采用隔离非线性纤维梁单元模拟，梁柱与剪力墙通过共结点方法整体建模，接触单元仅连接墙单元和土体单元。

表3 混凝土材料参数Table 3 Material parameters of concrete

表4 钢筋及土体材料参数Table 4 Material parameters of steel and soil

4.2 结果分析

对建立的考虑土-结相互作用的整体模型分别采用传统非线性分析方法和本文方法进行地震反应分析，其中传统分析方法通过在每次迭代时对整体刚度进行更新分解计算结构位移响应。

任选土体与结构接触面上四个参考点对，分别为图14 中所示A、B、C和D四个位置，这四个点对的法向相对位移时程响应如图15 所示，图中法向相对位移是指结构与土体之间的相对开度，接触面张开为正，可以看出，四个点对的初始状态均为受压接触。由图15 可以看出，在前5 s内，这四个参考点处接触点对的法向相对位移均在0 附近，说明结构与土体在该时间段内是接触状态；在5 s 之后，接触点对的相对位移数值呈现“忽高忽低”的波动现象，且其最大值开始逐渐增大，说明这些接触点对附近的土-结构接触面表现出开、合交替出现的行为，并在10 s 左右时法向相对位移达到峰值，此外，由图15(d)可以看出，点D处接触点对的法向相对位移在9.5 s 之后始终大于0，说明该点附近的接触面在9.5 s 之后始终处于脱开状态。

图14 接触面参考点对Fig.14 Reference point pair of interface

图15 法向相对位移反应Fig.15 Normal relative displacement response

图16 给出这四个点对的剪切相对位移时程曲线。综合各点对的法向和剪切相对位移时程曲线，可知接触面发生了法向的脱开和重新闭合、剪切方向的相对摩擦和错位滑移等非线性行为，这也验证了本文方法模拟接触面变形的有效性。此外，图15 和图16 表明本文方法与传统方法的分析结果基本一致，表5 给出本文方法和传统方法的峰值点相对位移数据及其相对误差值，可以看出，本文方法具有较高的计算精度。

表5 相对位移相对误差Table 5 Relative error of relative displacement

图16 剪切相对位移反应Fig.16 Shear relative displacement response

图17 和图18 分别给出了结构在x和z方向的结构顶点位移时程曲线和层间位移角曲线，为分析土-结构相互作用对结构动力响应的影响，建立不考虑土-结相互作用的计算模型并进行分析(图中模型a 为图12(a)所示主体结构模型，模型b 为图12(c)所示土-结相互作用整体模型)，从图17 和图18 可以看出，采用本文方法和传统方法对模型b 分析的结果基本一致。此外，由图17(a)和图17(b)可以看出，在第10 s～12 s 内两种模型分析得到的结构顶点位移达到峰值，且在x方向模型a 的峰值位移小于模型b；由图18(a)可以看出，模型a 除第7 层～第13 层的x方向层间位移角大于模型b 外，其他楼层层间位移角均小于模型b；由图18(b)可以看出，模型a 的层间位移角均小于模型b。上述结果表明：不考虑土-结相互作用会低估结构的位移反应。

图17 顶点位移时程曲线Fig.17 Time-history curves of top displacement

图18 层间位移角时程曲线Fig.18 Time-history curves of inter-story displacement angle

图19 给出了本文方法在整体土-结相互作用模型计算过程中的非线性自由度变化过程。从图中可以看出，结构在前200 分析步内已进入了非线性，但非线性自由度数目很小。整个分析过程中产生的非线性自由度最大值为135 206，占总自由度数(615 513)的22%，而每个增量步的平均非线性自由度为59 643，仅占总自由度数(615 513)的9.7%。这说明，本文方法在进行该结构地震反应分析过程中需要实时更新和分解的Schur 补矩阵阶数远小于结构整体刚度矩阵阶数，Woodbury 公式的使用可以显著提升结构分析效率。

图19 非线性自由度时程曲线Fig.19 Time-history curves of inelastic degrees of freedom

本文方法的高效性可通过与传统非线性分析方法的计算耗时进行对比来验证。其中，传统分析方法需在每个迭代步进行一次整体刚度矩阵的更新和分析，两种方法计算过程均采用相同平台，计算采用的硬件设备为小型服务器(CPU 为2 个Intel Xeon Gold 6254，内存为384 GB)。采用本文方法对该算例进行分析共耗时约6311 s；而采用传统非线性分析方法则需耗费约42 072 s，为本文方法所需计算时间的6.7 倍。

5 结论

为实现大型结构考虑土体-结构相互作用的高效地震反应分析，本文建立了基于隔离非线性的土-结构接触行为非线性计算模型：首先根据初始接触状态建立了接触面单元的相对位移分解方法；随后引入额外非线性自由度并结合Goodman 接触单元格式对接触面上的非线性变形进行模拟，能够保证单元刚度矩阵在分析过程中恒定不变；进一步采用Woodbury 公式进行求解，有效提升了此类问题的分析效率。本文得到结论如下：

(1) 隔离非线性法要求根据初始弹性刚度进行变形分解，以便于保持结构刚度的恒常性及实现局部非线性行为的隔离表达，然而，单元法向应力-相对位移关系在零点附近存在拐点，导致单元内任意点对在不同加载方向上的初始刚度不唯一，为在分析过程中始终使用恒定的初始刚度对单元任意时刻的应力状态进行描述，本文分别对不同初始状态下的单元变形分解的非线性相对位移进行了定义。数值算例分析结果表明，本文建立的三维接触面单元模型可对土与结构接触面处的非线性行为进行有效模拟，且不考虑土-结相互作用会低估上部结构的位移响应。

(2) 计算结果显示，分析过程中出现的非线性自由度较少，采用Woodbury 公式对结构控制方程进行求解时能够避免传统非线性分析方法在迭代求解过程中所需的大规模整体刚度矩阵的实时更新和分解，仅需对一个代表局部非线性行为的小规模Schur 补矩阵进行迭代更新，在保证精度的同时显著降低求解时间，对大型结构分析问题具有较高的效率优势。