建筑韧性抗震理念是国际防灾领域众多学者的广泛共识，建筑抗震韧性评价包括地震危险性分析、结构损伤分析和建筑功能损失与恢复能力计算^{[1 − 3]}，其中，基于结构非线性地震反应分析结果确定结构损伤状态是损伤分析的重要手段，而这一过程是影响建筑抗震韧性评估精度和效率的重要因素，已成为国内外学者关注的焦点。

目前，结构地震损伤分析多采用基于构件或结构宏观反应的损伤指标，如基于层间位移角^{[4 − 6]}或结构顶层最大位移^[7]等参数构造损伤指标，但对于复杂结构，此类损伤指标无法考虑地震作用下结构可能出现的平扭耦联等不规则变形，基于构件反应的损伤模型可考虑结构的局部变形，较为典型的有基于杆端转角的构件易损性模型和考虑最大变形与滞回耗能双参数线性组合的Park-Ang模型^[8]及其改进模型^{[9 − 10]}，但此类模型往往需基于大量试验数据进行标定，且缺乏清晰物理意义，计算结果具有较大不确定性^[11]，在用于韧性评估时将导致计算结果的可靠性降低。已有研究^[12]表明：直接从材料层次出发的结构损伤计算模型可以考虑材料损伤引起的整体结构力学性能劣化，具有明确的物理意义，且可以得到精细的结构损伤分布，将其应用于复杂结构的精细化抗震韧性评估具有显著的优越性，如FALEIRO等^{[13 − 14]}从混凝土的材料本构出发，采用以材料无损弹性自由能势作为权系数对损伤变量进行加权积分的损伤指标对结构进行损伤分析；STEPHENS等^[15]采用考虑材料累积塑性应变的损伤指标对结构进行损伤分析。然而，既有相关研究提出的损伤分析方法多针对特定本构模型或结构类型，适用性受限。可见，从材料层面出发，发展物理意义清晰且适用性广泛的精准化损伤分析方法，将有助于提升结构韧性评估的可靠性。

结构地震损伤分析需以非线性地震反应分析数据为基础，这通常需借助有限元法，随着现代工程结构的大型化和复杂化，对结构分析模型精细化程度的要求也越来越高，导致非线性分析所需计算耗时快速增多，效率与精度之间的矛盾十分尖锐。国内外学者利用工程结构在地震作用下的非线性变形局部化特征，发展出一系列高效数值分析方法^{[16 − 18]}，如基于重分析方法的结构局部非线性问题高效计算方法^{[19 − 22]}、子结构法^[17]、拟力法^{[23 − 25]}等。近期，LI等^{[26 − 28]}将非线性材料的应变分解为线弹性和非线性两部分，结合有限元法基本原理提出了隔离非线性分析方法，该方法将结构切线刚度矩阵转化为初始弹性刚度矩阵与由局部材料非线性引起的刚度修改项的叠加形式，进而利用Woodbury公式进行结构响应求解，能够实现地震作用下一般工程结构的高效、精细化非线性分析，然而，其分解得到的非线性应变分量不具有明确的物理意义，本质上仅为构造隔离非线性控制方程的数值处理手段，无法用来探究非线性分析过程中结构损伤演化过程。

本文基于隔离非线性法基本理论，提出了一种高效的结构精细化地震损伤分析方法，该方法首先将非线性材料的应变分解为线弹性、损伤和塑性三部分，并基于连续损伤力学和有效应力空间中的塑性力学理论确定损伤和塑性应变的演化法则，随后构造单元损伤和塑性应变的分布场模型，引入损伤和塑性自由度表征分离出的损伤和塑性应变，并结合虚功原理建立隔离损伤-塑性控制方程，进而引入Woodbury公式进行快速求解。本文方法可使非线性分析的主要计算量集中于小规模损伤和塑性矩阵的更新和分解，避免了大规模整体切线刚度矩阵耗时的迭代更新，引入的损伤与塑性应变物理意义清晰，且可在每个增量步的控制方程求解过程中实时获取，从而可基于此在分析过程中对结构任意部位的损伤状态进行实时量化计算，实现对结构在地震作用下的精细化损伤分布及其演化过程进行精准描述与追踪。

1   材料线弹性、损伤和塑性应变分解

连续损伤力学^{[29 − 30]}通过定义损伤内变量来描述由材料内部微缺陷导致的刚度、强度退化。在等温绝热条件下，基于不可逆热力学第二定律可推导含损伤变量的应力-应变关系，对于单轴材料本构，当采用标量损伤变量d描述材料损伤时，应力-应变关系式为：

式中：${E_{\mathrm{e}}}$为材料初始弹性模量；${E_{\mathrm{d}}}(d) = (1 - d){E_{\mathrm{e}}}$为材料弹性损伤模量：${\varepsilon ^{\rm{e}}}$为材料弹性应变；d为损伤变量，其描述材料弹性模量的退化。

连续损伤力学理论认为材料产生损伤后，损伤材料塑性变形的计算应采用基于有效应力空间的塑性力学方法，以图1所示单轴材料应力-应变关系进行说明，图中$\overline \sigma$-$\varepsilon$曲线为损伤材料的有效应力-应变曲线，损伤材料的有效应力$\overline \sigma$服从经典弹塑性力学关系，材料塑性变形的计算采用有效应力$\overline \sigma$和表征材料塑性特性的内变量${q^{\rm{p}}}$，而有效应力基于“等效应变假设”^{[29 − 32]}计算，该假设认为应力$\sigma$作用在损伤材料上引起的应变与有效应力$\overline \sigma$作用在相应无损伤材料上引起的应变等价，如图1所示，由此可得材料应力与有效应力的关系为：

图  1  非线性材料(弹塑性损伤本构关系)增量应变分解

Figure  1.  Incremental strain decomposition of nonlinear stress-strain relationship: the elastic plastic damage constitutive law

下载: 全尺寸图片 幻灯片

式中，$\overline \sigma = {E_{\rm{e}}}{\varepsilon ^{\rm{e}}}$为材料的有效应力。

LI等^{[26 − 27]}通过引入拟力法^{[23 − 25]}中的位移分解原理，将非线性材料的应变分解为线弹性和非线性两部分，提出了隔离非线性有限元法。受其概念启发，本文根据材料初始弹性刚度和弹性损伤刚度(即图1中材料卸载刚度${E_{{\mathrm{d}},k}}$)将非线性材料的应变分解为弹性、损伤和塑性三部分。

图1以材料单轴弹塑性损伤本构关系为例给出本文应变分解的基本原理，假设在第k个增量步内，材料从B点状态到达C点状态，此时，应力为${\sigma _{{C}}} = {\sigma _{{B}}} + \Delta \sigma$，塑性内变量为$q_{{C}}^{\rm{p}} = q_{{B}}^{\rm{p}} + \Delta q$，损伤内变量为${d_{{C}}} = {d_{{B}}} + \Delta d$，其中$\Delta$代表增量，下标B、C代表B点或C点状态。应变增量$\Delta \varepsilon$是可恢复应变增量$\Delta {\varepsilon ^{\rm{e}}}$和不可恢复塑性应变增量$\Delta {\varepsilon ^{\rm{p}}}$之和：

式中，不可恢复塑性应变增量$\Delta {\varepsilon ^{\rm{p}}}$为：

式中，${\varepsilon ^{\rm{p}}}$为不可恢复塑性应变。可恢复应变增量$\Delta {\varepsilon ^{\rm{e}}}$为应变增量$\Delta \varepsilon$与塑性应变增量$\Delta {\varepsilon ^{\rm{p}}}$的差值，本文将可恢复应变增量$\Delta {\varepsilon ^{\rm{e}}}$分解为线弹性应变增量$\Delta \varepsilon '$和损伤应变增量$\Delta {\varepsilon ^{\rm{d}}}$之和：

式中，线弹性应变增量$\Delta \varepsilon '$是假设材料无损时(即弹性模量不发生退化)的可恢复应变增量。损伤应变增量$\Delta {\varepsilon ^{\rm{d}}}$是由材料弹性模量损伤引起的可恢复应变增量，其与塑性应变增量$\Delta {\varepsilon ^{\rm{p}}}$之和称为材料应变增量的非线性部分，$F(\cdot)$代表函数。

考虑更为一般的多轴材料本构模型，将式(3)、式(4)和式(5)中的应变增量$\Delta \varepsilon$、可恢复应变增量$\Delta {\varepsilon ^{\rm{e}}}$、弹性应变增量$\Delta \varepsilon '$、损伤应变增量$\Delta {\varepsilon ^{\rm{d}}}$、塑性应变增量$\Delta {\varepsilon ^{\rm{p}}}$和应力增量$\Delta \sigma$替换为相应的向量或张量形式$\Delta {\boldsymbol{\varepsilon}}$、$\Delta {{\boldsymbol{\varepsilon}} ^{\rm{e}}}$、$\Delta {\boldsymbol{\varepsilon}} '$、$\Delta {{\boldsymbol{\varepsilon}} ^{\rm{d}}}$、$\Delta {{\boldsymbol{\varepsilon}} ^{\rm{p}}}$和$\Delta \sigma$。此时，基于应变分解的材料应力计算表达式为：

式中，${{\boldsymbol{D}}_{\rm{e}}}$为材料初始弹性刚度矩阵，且式(6)表明，无论何种弹塑性损伤材料，应变分解后任意时刻的材料应力可表示为恒常的初始弹性刚度矩阵与线弹性应变增量向量的乘积。

本文在式(6)中新引入变量$\Delta {{\boldsymbol{\varepsilon}} ^{\rm{d}}}$和$\Delta {{\boldsymbol{\varepsilon}} ^{\rm{p}}}$，需定义其与材料应力之间的关系式，首先讨论损伤应变，将含损伤变量的材料本构关系式(2)拓展至多轴，并给出其增量形式：

式中：${\boldsymbol{ R}}$反映了损伤变量演化引起的材料切线刚度退化；${\bf\textit{ω}}$为与有效应力增量$\Delta \overline{\boldsymbol{ \sigma}}$相关的量^[29]。

考虑式(5)与线弹性应变增量$\Delta {\boldsymbol{\varepsilon}} ' = {\boldsymbol{D}}_{\rm{e}}^{ - 1}\Delta {\boldsymbol{\sigma}}$，由式(7)得材料应力增量与损伤应变增量的关系式为：

进一步由“弹性应变等效假设”^[33]可知，对于产生塑性应变的材料，材料有效应力增量和应变增量之间的关系为：

式中：${{\boldsymbol{D}}_{\rm{ep}}}$为材料有效弹塑性切线刚度矩阵，即有效应力-应变曲线的切线刚度。基于式(2)和式(7)可得材料应力增量与塑性应变增量的关系式为：

式中，${({\boldsymbol{D}}_{\rm{ep}}^{ - 1} - {\boldsymbol{D}}_{\rm{e}}^{ - 1})^{ - 1}}$为有效塑性刚度。

2   单元塑性、损伤应变场

图2所示为一个一般结构体V的有限单元离散示意图，其中X、Y、Z代表笛卡尔坐标系，u、v、w代表沿三个坐标轴的位移分量，单元域$\varOmega$内的位移场可通过插值方法近似表示，任意线性化增量计算步内的单元增量位移场可表示为：

图  2  一般结构体V的有限元离散及四结点四边形单元$\varOmega$

Figure  2.  Finite-element mesh for an arbitrary volume V and four-node quadrilateral element $\varOmega$

下载: 全尺寸图片 幻灯片

式中，$\Delta {\boldsymbol{u}} = {[ {\Delta u}\;\;\;{\Delta v}\;\;\;{\Delta w} ]^{\rm{T}}}$代表单元内任一点的位移增量；单元的结点位移增量为$\Delta\boldsymbol{a}= [\Delta u_1\; \; \; \Delta v_1\; \; \; \Delta w_1\; \; \; \cdots\; \; \; \Delta u_i\; \; \; \Delta v_i\; \; \; \Delta w_i\; \; \; \cdots]^{\rm{T}}$；N为单元的位移插值函数矩阵。进一步将式(11)代入几何方程，可得单元的增量应变场为：

式中：L代表几何方程中的微分算子^[34]；B为应变矩阵。

式(6)中新引入的变量$\Delta {{\boldsymbol{\varepsilon}} ^{\rm{d}}}$和$\Delta {{\boldsymbol{\varepsilon}} ^{\rm{p}}}$在描述单元域$\varOmega$内应力分布时不可或缺，需建立单元的塑性、损伤应变场模型。本文参考单元位移场和应变场的模拟策略，在单元$\varOmega$内预先设置若干损伤、塑性应变插值点，使用插值方法建立单元的损伤、塑性应变场函数。为保证后文所提控制方程刚度矩阵的对称性，减小非线性分析过程中的存储量和计算量，本文方法令单元内预设的损伤和塑性应变插值点与高斯积分点位置重合。据此可建立单元内任一点损伤应变增量$\Delta {{\boldsymbol{\varepsilon}} ^{\rm{d}}}$和塑性应变增量$\Delta {{\boldsymbol{\varepsilon}} ^{\rm{p}}}$统一计算式：

其中：

式中，当in = d时代表损伤应变增量计算表达式，此时$\Delta {{\boldsymbol{\varepsilon}} ^{\rm{d}}}$代表单元内任一点的损伤应变增量，$\Delta {\boldsymbol{\varepsilon}} _{\rm{chp}}^{\rm{d}}$包含了单元中各个损伤应变插值点处的材料损伤应变增量，${h_{\rm{d}}}$为单元$\varOmega$内预设的损伤应变插值点个数，${{\boldsymbol{C}}^{\rm{d}}}$代表损伤应变插值函数矩阵，其中${\boldsymbol{C}}_{{i}}^{\rm{d}}$为第i个插值函数矩阵。同理，当in = p时，式(13)中的变量代表与塑性相关的向量和矩阵。

在某个线性化增量计算步中，若结构局部区域${V_{\rm{d}}}$仅发生材料损伤，局部区域${V_{\rm{p}}}$内材料仅出现塑性行为，局部区域${V_{\rm{dp}}}$内材料同时出现塑性和损伤行为，如图2所示，则${V_{\rm{d}}}$内仅损伤应变插值点的应变增量为非零值，${V_{\rm{p}}}$内仅塑性应变插值点的应变增量为非零值，${V_{\rm{dp}}}$内损伤、塑性应变插值点的应变增量均为非零值。对于图2所示结构模型中的单元$\varOmega$，区域${\varOmega _{\rm{d}}}$内损伤应变插值点的应变增量为非零值，${\varOmega _{\rm{p}}}$内塑性应变插值点的应变增量为非零值，${\varOmega _{\rm{dp}}}$内塑性、损伤应变插值点的应变增量均为非零值，而区域${\varOmega _{\rm{e}}}$处于线弹性状态。若单元内某个塑性(或损伤)应变插值点的应变增量为非零值，认为其在当前增量步下处于激活状态，反之认为其处于未激活状态，对于单元$\varOmega$，假设区域${\varOmega _{\rm{d}}}$和${\varOmega _{\rm{dp}}}$中共有l_d个损伤应变插值点处于激活状态(${l_{\rm{d}}} {\leqslant} {h_{\rm{d}}}$)，其中区域${\varOmega _{\rm{dp}}}$中有l_dp个损伤应变插值点，单元域中其余损伤应变插值点处于未激活状态(即这些点的损伤应变增量为零)；同样地，假设区域${\varOmega _{\rm{p}}}$和${\varOmega _{\rm{dp}}}$中共有l_p个塑性应变插值点处于激活状态(${l_{\rm{p}}} {\leqslant} {h_{\rm{p}}}$)，其中区域${\varOmega _{\rm{dp}}}$中有l_pd个塑性应变插值点，单元域中其余塑性应变插值点处于未激活状态(即这些点的塑性应变增量为零)，据此可对式(13)进行消元处理并将其改写为如下形式：

其中：

式中：当in = d时，${\boldsymbol{C}}_{\rm{a}}^{\rm{d}}$代表单元$\varOmega$中处于激活状态的损伤应变插值点处的插值函数矩阵，$\Delta {\boldsymbol{\varepsilon}} _{\rm{achp}}^{\rm{d}}$为相应损伤应变插值点处的应变增量，该向量中的每个元素代表一个损伤自由度；当in = p时，${\boldsymbol{C}}_{\rm{a}}^{\rm{p}}$代表单元$\varOmega$中处于激活状态的塑性应变插值点处的插值函数矩阵，$\Delta {\boldsymbol{\varepsilon}} _{\rm{achp}}^{\rm{p}}$为相应塑性应变插值点处的应变增量，该向量中的每个元素代表一个塑性自由度。为建立如式(14)所示的应变场函数，需在计算过程中实时确定单元$\varOmega$中各损伤应变和塑性应变插值点的状态，为此，本文方法在每次迭代时对各损伤应变插值点进行损伤状态判断，对各塑性应变插值点进行塑性状态判断。

将式(12)、式(14)代入式(6)，则单元内任一点的应力增量可通过下式进行计算：

式(15)表明，单元应力增量的分布与结点位移增量$\Delta {\boldsymbol{a}}$、损伤应变增量向量$\Delta {\boldsymbol{\varepsilon}} _{\rm{achp}}^{\rm{d}}$和塑性应变增量向量$\Delta {\boldsymbol{\varepsilon}} _{\rm{achp}}^{\rm{p}}$有关。

3   单元隔离损伤-塑性控制方程

为构建单元$\varOmega$的隔离损伤-塑性控制方程，需首先建立单元的虚功方程，即：

式中：$\delta {\boldsymbol{\varepsilon}}$和$\delta {\boldsymbol{a}}$分别代表单元虚应变和虚结点位移向量；$\Delta {\boldsymbol{f}}$为结点荷载增量；$\Delta {\boldsymbol{\sigma}}$为与荷载$\Delta f$相平衡的单元应力增量。将式(16)中的$\delta {\boldsymbol{\varepsilon}}$替换为($\delta {\boldsymbol{\varepsilon}} - \delta {{\boldsymbol{\varepsilon}} ^{\rm{d}}} - \delta {{\boldsymbol{\varepsilon}} ^{\rm{p}}} + \delta {{\boldsymbol{\varepsilon}} ^{\rm{d}}} + \delta {{\boldsymbol{\varepsilon}} ^{\rm{p}}}$)，其中，$\delta {{\boldsymbol{\varepsilon}} ^{\rm{d}}}$和$\delta {{\boldsymbol{\varepsilon}} ^{\rm{p}}}$为单元虚塑性、损伤应变向量，并将单元应变场式(12)、单元损伤和塑性应变场式(14)和应力增量计算式(15)代入式(16)中，进一步，在单元结点位移增量未知量$\Delta {\boldsymbol{a}}$的基础上引入损伤应变增量$\Delta {\boldsymbol{\varepsilon}} _{\rm{achp}}^{\rm{d}}$和塑性应变增量$\Delta {\boldsymbol{\varepsilon}} _{\rm{achp}}^{\rm{p}}$作为未知量，可建立如下方程：

其中：

式中：${{\boldsymbol{k}}_{\rm{e}}}$为单元$\varOmega$的初始弹性刚度矩阵，其在计算过程中保持不变；当in = d时，$\Delta {{\boldsymbol{f}}_{\rm{d}}}$是与单元${\varOmega _{\rm{d}}}$和${\varOmega _{\rm{dp}}}$区域中的材料损伤状态相关的应力增量，${{\boldsymbol{k}}'_{\rm{d}}}$和${{\boldsymbol{k}}_{\rm{d}}}$与单元$\varOmega$的初始弹性属性和${\varOmega _{\rm{d}}}$与${\varOmega _{\rm{dp}}}$区域(处于激活状态)中的l_d个损伤应变插值点有关；当in = p时，$\Delta {{\boldsymbol{f}}_{\rm{p}}}$是与单元${\varOmega _{\rm{p}}}$和${\varOmega _{\rm{dp}}}$区域中的材料塑性状态相关的应力增量；${{\boldsymbol{k}}'_{\rm{p}}}$和${{\boldsymbol{k}}_{\rm{p}}}$与单元$\varOmega$的初始弹性属性和${\varOmega _{\rm{p}}}$与${\varOmega _{\rm{dp}}}$区域(处于激活状态)中的l_p个塑性应变插值点有关；${{\boldsymbol{k}}'_{\rm{dp}}}$和${{\boldsymbol{k}}'_{\rm{pd}}}$与单元$\varOmega$的初始弹性属性和${\varOmega _{\rm{dp}}}$区域(处于激活状态)中的l_dp个损伤应变插值点和l_pd个塑性应变插值点有关。

根据本文定义的材料应力与塑性应变的关系式(10)和材料应力与损伤应变的关系式(8)，在单元域$\varOmega$中考虑损伤与塑性应变场式(14)，可得单元任意点处的材料应力增量为：

根据插值基函数${{\boldsymbol{C}}_i }$的正交性，可采用式(19)和式(20)计算各插值点处材料应力增量$\Delta {{\boldsymbol{f}}_{\rm{in}}}$。将其代入式(17)后单元控制方程可重新写为：

式中：

其中：

由于在推导式(14)时引入了插值点状态判断和自由度消元操作，式(21)中仅包含发生增量损伤应变和增量塑性应变(处于激活状态)的局部材料信息，其方程规模与单元$\varOmega$中的损伤自由度数目l_d和塑性自由度数目l_p有关；若单元$\varOmega$仅发生局部增量损伤应变，其方程规模仅与损伤自由度数目l_d有关，经自由度消元后式(14)中的矩阵${\boldsymbol{C}}_{\rm{a}}^{\rm{p}}$和向量$\Delta {\boldsymbol{\varepsilon}} _{\rm{achp}}^{\rm{p}}$为空，式(21)退化为弹性损伤单元的刚度方程，同样地，若仅发生局部增量塑性应变，其方程规模仅与塑性自由度数目l_p有关，式(14)中的矩阵${\boldsymbol{C}}_{\rm{a}}^{\rm{d}}$和向量$\Delta {\boldsymbol{\varepsilon}} _{\rm{achp}}^{\rm{d}}$为空，式(21)退化为弹塑性单元的刚度方程，此时式(21)左端的方程均为$2 \times 2$的块矩阵；若单元$\varOmega$未发生增量损伤应变和增量塑性应变，式(14)中的矩阵${\boldsymbol{C}}_{\rm{a}}^{\rm{in}}$和向量$\Delta {\boldsymbol{\varepsilon}} _{\rm{achp}}^{\rm{in}}$为空，式(21)退化为弹性单元的刚度方程(即${{\boldsymbol{k}}_{\rm{e}}}\Delta {\boldsymbol{a }}= \Delta {\boldsymbol{f}}$)；若分析过程中单元$\varOmega$域中的所有损伤应变插值点和塑性应变插值点始终处于激活状态，则${{\boldsymbol{k}}'_{\rm{in}}}$、${{\boldsymbol{k}}_{\rm{in}}}$、${{\boldsymbol{k}}'_{\rm{dp}}}$和${{\boldsymbol{k}}'_{\rm{pd}}}$为常数矩阵。

4   结构控制方程求解及损伤状态分析

4.1   结构控制方程

将单元控制方程进行集成，可得到结构的隔离损伤-塑性控制方程：

式中：$\Delta {\boldsymbol{X}}$和$\Delta {\boldsymbol{F}}$分别为结构的结点位移增量和结点荷载增量；$\Delta {\boldsymbol E}_{\rm{achp}}^{\rm{d}}$和$\Delta {\boldsymbol E}_{\rm{achp}}^{\rm{p}}$为结构的损伤应变增量和塑性应变增量，其中分别包含了结构中处于激活状态的m_d个损伤自由度和m_p个塑性自由度；${{\boldsymbol{K}}_{\rm{e}}}$为结构的初始弹性刚度矩阵，其阶数等于结构位移自由度数n，其在计算过程中保持不变；${{\boldsymbol{K}}'_{\rm{d}}}$、${{\boldsymbol{K}}'_{\rm{p}}}$、${{\boldsymbol{K}}'_{\rm{dp}}}$和${{\boldsymbol{K}}'_{\rm{pd}}}$的规模分别为n$\times$m_d阶、n$\times$m_p阶、m_d$\times$m_p阶和m_p$\times$m_d阶，其与材料损伤、塑性本构关系无关，仅与发生增量损伤、塑性应变的区域位置相关，${{\boldsymbol{K}}''_{\rm{d}}}$和${{\boldsymbol{K}}''_{\rm{p}}}$矩阵的规模分别为m_d$\times$m_d阶和m_p$\times$m_p阶，其与发生增量损伤、塑性应变的区域位置以及区域内材料的损伤、塑性状态均有关，计算过程中该矩阵的确定依赖于材料塑性、损伤本构模型，矩阵维度m_d和m_p在计算过程中随荷载条件和插值点的损伤、塑性状态的变化而变化，可用来实时反映结构中发生增量损伤、塑性应变的区域规模。式(25)中${\Delta }{\boldsymbol{X}}$、${\Delta }{\boldsymbol E}_{\rm{achp}}^{\rm{d}}$和${\Delta }{\boldsymbol E}_{\rm{achp}}^{\rm{p}}$为未知量，相应的自由度数目与方程数目相等，方程可解。

式(25)提出的隔离损伤-塑性控制方程可直接用于结构的静力非线性分析，而针对结构的运动状态，本文结合运动微分方程，建立隔离损伤-塑性的动力分析控制方程，实现高效的结构非线性地震反应分析。结构在地震作用下的增量运动微分方程为：

式中：${\boldsymbol{M}}$为质量矩阵；${\boldsymbol{C }}$为阻尼矩阵；$\Delta\ddot{\boldsymbol{X}}\mathit{_g }$为地震动加速度增量向量；$\Delta \ddot {\boldsymbol{X}}$和$\Delta \dot {\boldsymbol{X}}$分别为结构加速度增量向量和速度增量向量；$\Delta {\boldsymbol{F}}$为结构整体恢复力的变化。根据隔离非线性的动力分析理论^[27]，首先将式(26)转化为考虑损伤和塑性虚拟荷载的弹性结构运动微分方程，接着引入Newmark-$\beta$法^[35]对弹性结构运动微分方程进行数值求解，可得与式(25)形式相同的隔离损伤与塑性的动力分析控制方程，控制方程与式(25)相比仅是将弹性刚度${{\boldsymbol{K}}_{\rm{e}}}$和荷载项$\Delta {\boldsymbol{F}}$替换为包含动力分析相关参量的有效刚度和有效荷载^[27]。

从式(25)的推导过程可以看出，通过对非线性材料的应变进行分解并增设额外的损伤、塑性自由度，可在结构控制方程中将整体弹性刚度矩阵(即${{\boldsymbol{K}}_{\rm{e}}}$)与代表损伤变形的刚度项(${{\boldsymbol{K}}'_{\rm{d}}}$和${{\boldsymbol{K}}''_{\rm{d}}}$)和代表塑性变形的刚度项(${{\boldsymbol{K}}'_{\rm{p}}}$和${{\boldsymbol{K}}''_{\rm{p}}}$)相分离，且引入的损伤应变$\Delta {\boldsymbol{\varepsilon}} _{\rm{achp}}^{\rm{d}}$和塑性应变$\Delta {\boldsymbol{\varepsilon}} _{\rm{achp}}^{\rm{p}}$分别描述损伤引起的材料弹性刚度退化和塑性引起的材料不可恢复变形，相比隔离非线性控制方程^[26]在数值意义上将材料应变分解为线弹性和非线性具有更加明确的物理意义。

4.2   控制方程求解

为高效求解本文控制方程，将结构控制方程进行凝聚后可写为：

式中：

上式表明：当(${m_{\rm{d}}} + {m_{\rm{p}}}$)远小于n时，任意时刻结构的切线刚度矩阵可表示为初始弹性刚度矩阵的低秩摄动形式。

从式(27)的推导过程可知，本文所提方法将非线性材料的应变分解为弹性、损伤和塑性三部分并增设额外的损伤和塑性自由度，最终可在控制方程中将代表局部损伤变形的刚度矩阵${{\boldsymbol{K}}''_{\text{d}}}$和代表局部塑性变形的刚度矩阵${{\boldsymbol{K}}''_{\text{p}}}$从整体刚度矩阵中隔离表达出来，鉴于此，本文将所提方法命名为“隔离损伤-塑性方法”。

采用完全Newton-Raphson迭代法进行结构非线性分析，引入Woodbury公式^[26](式(28))在每个迭代步中对式(27)进行高效求解：

式中：

由式(28)和式(29)可以看出，结构整体初始弹性刚度矩阵${{\boldsymbol{K}}_{\rm{e}}}$在分析过程中始终保持不变，其仅需在分析前进行一次分解即可。由于地震作用下工程结构的材料损伤和塑性行为往往集中于局部区域，具有典型的局部化特征，因此，本文所提控制方程在求解过程中损伤和塑性自由度数目(${m_{\rm{d}}} + {m_{\rm{p}}}$)远小于结构位移自由度数目n (即(m_d + m_p) << n)，此时Woodbury公式在求解结构位移响应时的主要计算开销仅集中于小规模的${{\boldsymbol{K}}_{\rm{Sch}}}$矩阵((${m_{\rm{d}}} + {m_{\rm{p}}}$)$\times$(${m_{\rm{d}}} + {m_{\rm{p}}}$)阶)的更新与分解，避免了传统变刚度法中整体切线刚度矩阵(n$\times$n阶)的反复合成与分解运算，实现了计算效率的提升。

4.3   结构损伤分析

在结构动力非线性分析的每一增量步中，本文控制方程求解式(28)的求解原则为从右至左依次进行，使得求解过程始终是矩阵和向量之间的运算，如图3所示。

图  3  Woodbury公式求解步骤

Figure  3.  Computational process of Woodbury formula

下载: 全尺寸图片 幻灯片

由图3可知，结构的损伤应变增量$\Delta {\boldsymbol E}_{\rm{achp}}^{\rm{d}}$和塑性应变增量$\Delta {\boldsymbol E}_{\rm{achp}}^{\rm{p}}$为控制方程求解过程中的中间变量(图3中红色字体标注的求解步)，可在计算过程中直接获得，其包含了各单元中处于激活状态的相应插值点处的损伤应变增量$\Delta {\boldsymbol{\varepsilon}} _{\rm{achp}}^{\rm{d}}$与塑性应变增量$\Delta {\boldsymbol{\varepsilon}} _{\rm{achp}}^{\rm{p}}$。将上述增量损伤应变和增量塑性应变进行累加，并结合式(14)定义的材料损伤应变场和塑性应变场模型，即可确定地震作用下结构中任意部位在任意时刻t的损伤应变和塑性应变(使用${\boldsymbol{\varepsilon}} _{\rm{achp}}^{\rm{d}}(t)$和${\boldsymbol{\varepsilon}} _{\rm{achp}}^{\rm{p}}(t)$表示)。进一步通过构造基于损伤及塑性应变的损伤指标，即可对任意时刻结构任意部位损伤状态量化分析，从而实现对结构在地震作用下的精细化损伤分布及其演化过程进行精准描述与追踪。这一过程可表述为如下数学形式：

式中，${\rm DA{M}_{{\text{pt}}}}$为损伤指数。应当指出的是，本文方法并不针对特定的结构或模型，因此上式的具体表达形式并不唯一，可根据具体结构形式或材料特性进行构造，如对于弹塑性材料，通常采用塑性应变来对材料损伤进行量化，此时材料损伤指数可具体表示为：

式中：${\varepsilon ^{{\mathrm{p,eq}}}}$为材料塑性应变的等效单轴形式；${\varepsilon ^{{\mathrm{pu,eq}}}}$为材料极限等效塑性应变^[36]。若以地震反应末期材料的塑性应变为依据计算上式(31)，则计算结果代表了震后的不可恢复变形，可用来在韧性评价中表征结构不同部位的可修复性，若以累计塑性应变为依据计算上式(31)，则计算结果可用来表征结构中任意部位的损伤演化过程和程度。针对一般的各向同性均质连续损伤材料，标量损伤变量足够描述多轴应力状态下的卸载刚度退化效应，通常可根据材料多轴应力状态下的损伤应变与标量损伤变量的关系计算材料的损伤指数。

基于式(30)获取任意时刻结构的精细化损伤分布状态后，可以此为基础进一步对构件的各构成单元的材料损伤指数通过加权积分方式计算各个构件的实时损伤指数^[37]：

式中：$\psi$为权重系数，例如，依据文献[37]该系数可取为材料的弹性应变能密度；V为构件域；${\rm DA{M}_{\rm{comp}}}$为构件的损伤指数。

对于结构整体损伤指数的计算，可通过加权求和的方式对构件、楼层和结构的损伤指数进行逐级计算^[38]。首先对各楼层中的结构构件按照构件类别、布置方向等属性划分为不同的组，然后将每一组各构件的损伤指数进行加权求和(如以构件耗能作为权系数^[38])得到各组的损伤指数，再对同一楼层中各组的损伤指数进行加权求和得到各楼层的损伤指数，最后对各层损伤指数进行加权求和得到整体结构的损伤指数^{[39 − 40]}。

5   数值算例

为了验证本文所提方法的正确性与高效性，对文献[41]中4层约束砌体结构进行地震反应分析。结构层高为3 m，总高度为12 m，纵墙上窗洞尺寸为1.5 m$\times$1.8 m，门洞尺寸为1.0 m$\times$1.8 m。该结构在纵横墙交接处设置钢筋混凝土构造柱，且所有纵横墙在楼层标高处设置钢筋混凝土圈梁，构造柱和圈梁采用C20混凝土，截面尺寸为250 mm$\times$250 mm，纵筋采用4根HRB335直径16 mm的钢筋。砌体墙厚度为250 mm，材料为烧结普通砖，强度等级为MU10，混合砂浆强度等级为M10。结构楼面恒载为5.49 kN/m²，活载为2.0 kN/m²，屋面恒载为6.371 kN/m²。设防烈度为8度(0.2 g)，Ⅱ类场地，由《建筑抗震设计规范》^[4]得场地覆盖层厚度${\geqslant}$5 m的等效剪切波速$500 {\geqslant} {{V_{\rm se}}} > 250$ m/s，据此选取1978年9月16日伊朗东部TABAS地区发生地震中Ferdows台站记录的TABAS\FER-L1方向的地震波分量(所在场地的剪切波速V_se = 302.64 m/s)，加速度时程如图4所示，并以PGA = 400 cm/s²进行调幅。

图  4  TABAS地震波加速度记录

Figure  4.  Acceleration record of TABAS earthquake

下载: 全尺寸图片 幻灯片

图5为该结构的数值分析模型，模型假设在结构首层底端设置固定约束，圈梁和构造柱采用纤维梁单元^[27]进行模拟，共划分1894个纤维梁单元，圈梁和构造柱中的钢筋通过在纤维梁柱模型中添加钢筋纤维来模拟。本文方法在分析纤维梁单元时，首先对梁柱截面变形进行分解，分解过程引入的截面损伤变形由产生卸载刚度退化的混凝土材料提供，而引入的截面塑性变形由钢筋和混凝土材料共同提供。砌体墙采用离散宏单元^[41]进行模拟，该单元由剪切子单元和界面子单元构成，通过耦合两种子单元的变形模拟砌体墙的弯曲和剪切非线性行为，单元详细特征可参考文献[41]，整体结构共划分4198个剪切单元，剪切单元之间采用9770个无厚界面单元相互连接，模型共有39934个位移自由度。本文假设砌体与圈梁和构造柱的拉结钢筋具有足够的强度，采用离散宏单元中界面子单元的一侧与纤维梁单元共结点的方式模拟，界面子单元的另一侧与离散宏单元中的剪切子单元连接，以实现界面子单元模拟砌体墙的弯曲变形和轴向变形。整体结构模型在楼板的每一区格设置两根斜向刚性杆传递水平荷载和协调楼面结点位移。为简化模拟，本算例模型中仅考虑离散宏单元中剪切子单元斜向弹簧的非线性变形和界面子单元的轴向和弯曲非线性变形，对界面子单元的剪切和扭转变形假设为弹性。结构第一振型是一阶Y向平动，一阶振型的周期为T₁=0.14 s，第二振型是一阶X向平动，二阶振型的周期是T₂=0.1 s，本文分析结构在Y向(弱轴)水平地震作用下的动力响应。

图  5  约束砌体结构数值分析模型

Figure  5.  The discrete macro-element model of the masonry structure

下载: 全尺寸图片 幻灯片

结构顶层外纵墙与外横墙交接处(图5示A点)的位移时程曲线如图6所示，对比隔离损伤-塑性方法与隔离非线性方法^{[41 − 42]}的计算结果，两者吻合良好，验证了本文所提分析方法的正确性。由各楼层位移时程曲线发现该结构底层层间位移角最大，峰值达到0.4%，依据约束砌体结构0.4%倒塌限值^[43]判定该结构已达到倒塌的破坏状态。

图  6  顶层位移时程曲线

Figure  6.  History responses at top floor

下载: 全尺寸图片 幻灯片

对于本算例，采用材料的卸载刚度退化系数作为损伤指数，该系数与最大损伤应变(即${\max_{\rm{h}}}{\varepsilon ^{\rm{d}}}$)存在一一对应关系，可直接依据计算出的最大损伤应变获取。图7给出了本文砌体结构1号外横墙关键部位A、B、C和D剪切单元的弹塑性损伤本构模型^[41]的损伤变量d与最大损伤应变${\max_{\rm{h}}}{\varepsilon ^{\rm{d}}}$的关系曲线，图8给出了根据本文控制方程求解过程中计算出的相应1号外横墙关键部位A、B、C和D剪切单元的损伤变形$\Delta {\boldsymbol E}_{\rm{achp}}^{\rm{d}}$与塑性变形$\Delta {\boldsymbol E}_{\rm{achp}}^{\rm{p}}$的时程曲线，1号外横墙关键部位A、B、C和D剪切单元的位置如图9(a)所示，其中A单元、C单元和D单元分别为结构1号外横墙一层、二层和三层的刚度突变位置，B单元为结构1号外横墙的一层中部开门洞填充墙的门洞上部墙肢。结合图7所示损伤变量曲线和图8所示损伤变形时程曲线，图9和图10分别对结构1号外横墙墙肢(在图5中标注)与2号内横墙墙肢(在图5中标注)的损伤指数分布时程演化过程进行了精细化的描述与追踪，基于此可实现对结构损伤状态的精准刻画。

图  7  损伤变量与损伤变形对应关系

Figure  7.  Relation of damage deformation and damage variable

下载: 全尺寸图片 幻灯片

图  8  结构关键部位非线性应变时程曲线

Figure  8.  Inelastic deformation history responses of key elements

下载: 全尺寸图片 幻灯片

图  9  1号外横墙损伤云图

Figure  9.  Damage nephogram of No.1 outer transversal wall

下载: 全尺寸图片 幻灯片

图  10  2号内横墙损伤云图

Figure  10.  Damage nephogram of No.2 inner transversal wall

下载: 全尺寸图片 幻灯片

由图9和图10可知：结构各层剪切单元发生不同程度的剪切破坏，下部1、2层损伤明显比上部3、4层严重，这是由于在水平地震作用下，结构发生了典型的剪切破坏模式^[44]，结构底部一层所受地震剪力最大，二层次之，破坏主要集中在结构一、二层，而三层和四层所受地震剪力小，仅出现轻微损伤。图11得到了1号外横墙关键部位(图9(a)给出了关键部位的位置)的损伤指数时程曲线，由图可知损伤指数逐渐增大，结构性能水平不断降低。

图  11  关键部位损伤指数时程曲线

Figure  11.  Damage evolutions of key elements

下载: 全尺寸图片 幻灯片

图12给出了结构损伤自由度数和塑性自由度数以及非线性自由度数(损伤自由度与塑性自由度的总和)随增量步的变化曲线。非线性分析过程中激活的最大非线性自由度数为10343，占结构位移自由度数的25.9%，而每个增量步的平均非线性自由度为3597，仅占结构位移自由度数的9%，说明结构具有显著的局部非线性特征，采用式(28)进行求解时，由于需进行实时更新和分解的K_Sch矩阵阶数等于非线性自由度数目，局部非线性引发的较少的非线性自由度可保证本文方法效率的大幅度提升。

图  12  非线性、损伤和塑性自由度时程曲线

Figure  12.  Time-history curves of inelastic, damage and plastic degrees of freedom

下载: 全尺寸图片 幻灯片

为了进一步论证本文方法的高效性，将传统非线性分析方法(矩阵分解采用LDL^T法)^[45]和本文方法的时间复杂度^{[46 − 48]}进行统计对比，如图13所示，本文隔离损伤-塑性方法每一迭代步的时间复杂度与产生非线性自由度数目有关，本算例平均时间复杂度为0.37$\times$10¹⁰，在本文方法计算过程的每一迭代步中采用传统非线性分析方法求解时，传统非线性分析方法的时间复杂度与非线性数目无关，计算得传统非线性分析方法的平均时间复杂度为1.39$\times$10¹⁰，是本文方法的3.76倍，证明了本文方法相对传统非线性分析方法的高效性。

图  13  时间复杂度时程曲线

Figure  13.  Time-history curves of time complexity

下载: 全尺寸图片 幻灯片

6   结论

本文方法将非线性材料的应变分解为线弹性、损伤和塑性三部分，引入损伤和塑性自由度表征分离出的损伤和塑性应变，结合有限元法理论建立了隔离损伤-塑性方法，得到如下结论：

(1) 通过非线性材料的应变分解，使控制方程中结构弹性刚度矩阵保持不变，结构中损伤和塑性信息集中于控制方程中的损伤和塑性矩阵块中，实现了结构线弹性、损伤、塑性状态的分离。

(2) 隔离损伤-塑性方法对结构进行非线性分析时仅对小规模损伤和塑性矩阵进行迭代更新，避免了传统变刚度法中对整体切向刚度矩阵的迭代更新，实现结构非线性问题的高效求解，且损伤和塑性矩阵对应的损伤和塑性自由度的位置可衡量结构中发生实时损伤与塑性变形的区域规模。

(3) 隔离损伤-塑性方法引入的损伤与塑性应变具有明确的物理意义，且可在每一增量步的控制方程求解过程中实时获取，从而可基于结构损伤与塑性应变场建立结构中任意部位的损伤量化指标，实现对结构在地震作用下的精细化损伤分布及其演化过程进行精准描述与追踪。