Mechanical Performance Characteristics and Numerical Analysis of Consititutive Model for Circular Steel Pipes Considering Damage Accumulation
-
摘要: 由于单层网壳结构主要以圆钢管为主，因此圆钢管本构模型更适用于研究单层网壳结构的力学行为。为探究材料损伤累积对单层网壳结构动力响应的影响，本文以圆钢管考虑损伤累积的本构模型为研究对象，通过与Prandtl–Reuss材料本构模型在屈服准则、弹性刚度矩阵、弹塑性刚度矩阵、加载准则、流动法则及应力跌落等方面的对比，分析其力学行为特点，给出其数值积分格式，并基于有限元软件ANSYS开发了考虑损伤累积的用户材料子程序，采用该子程序研究单层网壳在强震作用下的动力响应。结果表明：Prandtl–Reuss材料模型高估了单层网壳结构的抗震性能，考虑材料损伤累积的结构失效时，结构动力极限荷载相比于P–R模型降低了12.24%，绝大多数杆件全截面屈服，形成大量塑性区带。一方面，损伤累积效应会降低结构的动力极限荷载，扩大损伤分布，加深损伤程度，加速结构破坏进程；另一方面，也将造成结构的失效模式由动力失稳向强度破坏转变。本文基于有限元软件ANSYS开发的考虑圆钢管损伤累积的材料子程序具有较高的计算精度，且其数值积分策略及用户材料子程序有效、可行，能够为单层网壳结构的抗震性能提供更加精确的分析与工程设计。Abstract: Since the single-layer reticulated shell structure is mainly composed of circular steel tubes, the constitutive model for circular steel tube is more suitable to study the mechanical behavior of single-layer reticulated shell structures. In this study, the constitutive model of circular steel pipe considering damage accumulation was analyzed with yield criteria, elastic stiffness matrix, elastic-plastic stiffness matrix, loading criterion, flow rule, and stress-drop, and it was compared to the constitutive model of Prandtl-Reuss material to explore the mechanical performance. By providing the format of the numerical integration, this study also explored the influence of material damage accumulation on the dynamic response of single-layer reticulated shell structure. Besides, a user-defined material subroutine considering damage accumulation was incorporated into ANSYS to investigate the dynamic response of single-layer latticed shell under earthquake action. The results show that seismic performance of single-layer reticulated shell was overestimated in Prandtl-Reuss material model. The dynamic ultimate load of structure was reduced by 12.24% compared with the Prandtl–Reuss model, and most of the members yield in the whole section, forming a large number of plastic zones, when considering the structural failure of material damage accumulation. On the other hand, the cumulative damage had the influence on reducing the dynamic limit load of the structure, expanding the damage distribution, deepening the damage degree and accelerating the process of structural failure, which also caused the failure mode of the structure to change from dynamic instability to strength failure. Moreover, the numerical integration scheme was proposed in this study showing better rationality and efficiency, and the calculation accuracy of the user-defined material subroutine is reliable to provide more accurate analysis and engineering design for the seismic performance of single-layer reticulated shell structures.
-
单层网壳结构被广泛应用于大型公共建筑中。若该类结构在地震中破坏甚至倒塌，必然会造成惨痛的人员伤亡和严重的财产损失。为此，研究人员对地震作用下网壳结构的破坏模式及失效机理展开了系列研究[1-6]。网壳结构在强震下将呈现两种失效模式，结构几何非线性所导致的动力失稳和材料塑性过度发展所导致的动力强度破坏。对于动力失稳，可根据各种稳定性判定准则取得相应的动力失稳临界荷载；而对于动力强度破坏，在合理选取特征响应量的同时，材料将因反复地震作用而产生的损伤累积效应对结构性能产生较大影响[7-8]。故理想弹塑性模型很难得到准确的强度破坏极限荷载。1958年，Kachanov[9]率先提出了连续度的概念，并用连续变量描述材料损伤的累积变化过程；在此基础上，Jean[10]在热力学框架内，提出了Lemaitre模型；Hammi等[11]在Lemaitre模型基础上进一步实现了空核的形成、增长及合并的模拟。为了更准确地预测断裂位置及断裂时刻各内变量的值，研究人员[12-16]在本构模型中考虑了温度、应力三轴比及作用角等因素的影响。然而，该类模型的提出主要基于材料试验的结果，且在结构构件层面尚未有相关数值分析与试验结果对此类模型进行验证。为准确掌握材料损伤累积对钢结构动力响应的影响，Kumar等[17]通过试验研究钢管柱的滞回性能，提出了钢材考虑循环加载的本构模型；沈祖炎[18]、董宝[19] 等对“H”型钢柱的滞回性能进行试验研究，拟合得到了钢材考虑损伤累积的本构模型；王萌等[20]研究了考虑损伤退化的钢材等效本构模型在整体框架中的应用及其用于强震作用下钢框架结构弹塑性时程分析的可行性。
由于单层网壳结构杆件主要以圆钢管为主，其受载特点和截面形式与高层建筑结构有较大区别，故以“H”型钢柱或钢管柱为试验对象而提出的本构模型并不完全适用于单层网壳结构[21]。为此，范峰等[22-23]针对空圆钢管构件进行试验研究，获得圆钢管考虑损伤累积的本构模型（以下简称圆钢管损伤本构模型）。故本文以此为研究对象，分析其力学行为特点，明确其数值积分格式，并基于有限元软件ANSYS开发了考虑损伤累积的材料子程序，采用该子程序研究单层网壳结构在强震作用下的动力响应。
1. 圆钢管损伤本构模型力学行为特点
对于延性金属材料，文献[22] 提出了考虑损伤累积材料本构模型如式（1）～（3）所示：
$$ D = \left( {1 - \beta } \right)\frac{{\varepsilon _{\text{m}}^{\text{p}}}}{{\varepsilon _{\text{u}}^{\text{p}}}} + \beta \sum\limits_{i = 1}^n {\frac{{\varepsilon _i^{\text{p}}}}{{\varepsilon _{\text{u}}^{\text{p}}}}} $$ (1) 式中：D为材料损伤因子，表征材料的损伤程度，为1个0～1之间的量，其中，0代表无损，1代表材料断裂破坏；
$ \varepsilon _{\text{m}}^{\text{p}} $ 为钢材所经历的最大塑性应变；$ \varepsilon _i^{\text{p}} $ 为钢材在第i次半循环中的塑性应变；$ \varepsilon _{\text{u}}^{\text{p}} $ 为钢材在一次拉抻时的极限塑性应变；β为待定权重系数；n为反复荷载的半循环周数。此时，对应于损伤因子D时的材料弹性模量和屈服强度为：
$$ {E_{{D}}} = \left( {1 - {\xi _1}D} \right)E $$ (2) $$ {f_{{D}}} = \left( {1 - {\xi _2}D} \right){f_{\text{y}}} $$ (3) 式（2）～（3）中，E、ED分别为无损伤和具有损伤因子D时的弹性模量，fy、fD分别为无损伤和具有损伤因子D时的屈服强度，ξ1、ξ2为待定系数。
文献[23] 针对单层网壳结构中常用的圆钢管构件开展了空间滞回试验，通过对比数值分析与试验结果，运用最小二乘法拟合了式（1）～（3）的系数β、ξ1及ξ2，得到圆钢管损伤本构模型如式（4）～（6）所示：
$$ {\qquad D=0.976\;8\frac{{\varepsilon }_{\text{m}}^{\text{p}}}{{\varepsilon }_{\text{u}}^{\text{p}}}+0.023\;2{\displaystyle \sum _{i=1}^{n}\frac{{\varepsilon }_{i}^{\text{p}}}{{\varepsilon }_{\text{u}}^{\text{p}}}}} $$ (4) $$ {E_{{D}}} = \left( {1 - 0.39D} \right)E $$ (5) $$ {f_{{D}}} = \left( {1 - 0.058D} \right){f_{\text{y}}} $$ (6) 在研究钢材的力学行为时，通常把其视为Prandtl–Reuss材料模型，即在Von Mises屈服准则和其关联的流动法则基础上推导出的理想弹塑性应力–应变关系。当采用式（4）～（6）以考虑损伤累积效应时，其本构模型与典型的Prandtl–Reuss材料模型相比，其模型具有统一性。
1.1 屈服准则
Prandtl–Reuss材料模型采用的是Von Mises屈服准则，并假定材料为理想弹塑性，故在主应力空间中屈服面为与静水压力无关的圆柱面，且仅存在一个屈服面。如果一个应力点在屈服面内时，其应力状态称之为弹性状态且只有弹性特性；另一方面，在屈服面上的应力状态称之为塑性状态，产生弹性或者弹塑性特性[24]。其屈服函数表达为：
$$ f = \sqrt {{J_2}} - k $$ (7) 式中：J2为偏应力张量第2不变量； k为材料常数，
$ k = {f_{\text{y}}}/\sqrt 3 $ ，其中，$ {f_{\rm{y}}} $ 为屈服强度。当钢管构件经历反复的地震作用时，材料中的微缺陷与微裂纹将不断发展，进而造成材料软化的现象，直接表现为其屈服强度的下降。此时，k不再为常数，而是损伤因子的函数，其函数表达式由“狗骨头”试件单轴拉伸试验与圆钢管空间滞回试验结果拟合得到[25-26]：
$$ {k_{{D}}} = {f_{{D}}}/\sqrt 3 = \left( {1 - 0.058D} \right){f_{\text{y}}}/\sqrt 3 $$ (8) 则圆钢管损伤本构模型屈服函数可表示为：
$$ f = \sqrt {{J_2}} - {k_D} = \sqrt {{J_2}} - \left( {1 - 0.058D} \right){f_{\text{y}}}/\sqrt 3 $$ (9) 由此可知，随着损伤因子的增大，屈服强度将不断减小，屈服面不断收缩，这导致多个屈服面的产生。将屈服面分为初始屈服面与后继屈服面。当对未曾屈服的应力点进行应力状态判断时，应以初始屈服面为依据；当对曾经屈服的应力点进行应力状态判断时，应以后继屈服面为依据。
1.2 弹性刚度矩阵与弹塑性刚度矩阵
对于Prandtl–Reuss材料模型，其弹性刚度矩阵为常数矩阵，用张量的形式表述为：
$$ {\boldsymbol{C}}_{ijkl}^{\text{e}} = \frac{{\nu E}}{{\left( {1 + \nu } \right)\left( {1 - 2\nu } \right)}}{\delta _{ij}}{\delta _{kl}} + \frac{E}{{2\left( {1 + \nu } \right)}}\left( {{\delta _{ik}}{\delta _{jl}} + {\delta _{ik}}{\delta _{jk}}} \right) $$ (10) 式中，
${\boldsymbol{C}}_{ijkl}^{\text{e}}$ 为4阶弹性刚度张量，E为弹性模量，ν为泊松比，δ为Kroneker符号。当考虑损伤累积效应时，弹性模量不再是常量，而随着损伤因子的增大而不断减小[27-28]。因此，在进行每个增量步时，需根据损伤程度的发展对弹性刚度矩阵进行修正。Prandtl–Reuss材料模型的弹塑性刚度矩阵用张量的形式表述为：
$$ \begin{aligned}[b] {\boldsymbol{C}}_{ijkl}^{{\text{ep}}} =& \frac{{\nu E}}{{\left( {1 + \nu } \right)\left( {1 - 2\nu } \right)}}{\delta _{ij}}{\delta _{kl}} +\\& \frac{E}{{2\left( {1 + \nu } \right)}}\left( {{\delta _{ik}}{\delta _{jl}} + {\delta _{ik}}{\delta _{jk}}} \right) - \frac{E}{{2\left( {1 + \nu } \right){k^2}}}{s_{ij}}{s_{kl}} \end{aligned} $$ (11) 式中，sij、skl分别为偏应力张量的分量。
根据式（8），修正后的圆钢管损伤本构模型的损伤弹性刚度矩阵与损伤弹塑性刚度矩阵张量表达式分别为：
$$ {\boldsymbol{C}}_{}^{\text{e}} = \frac{{\nu {E_{{D}}}}}{{\left( {1 + \nu } \right)\left( {1 - 2\nu } \right)}}{\delta _{ij}}{\delta _{kl}} + \frac{{{E_{{D}}}}}{{2\left( {1 + \nu } \right)}}\left( {{\delta _{ik}}{\delta _{jl}} + {\delta _{ik}}{\delta _{jk}}} \right) $$ (12) $$\begin{aligned}[b] {\boldsymbol{C}}_{}^{{\text{ep}}} =& \frac{{\nu {E_{{D}}}}}{{\left( {1 + \nu } \right)\left( {1 - 2\nu } \right)}}{\delta _{ij}}{\delta _{kl}} + \frac{{{E_{{D}}}}}{{2\left( {1 + \nu } \right)}}\left( {{\delta _{ik}}{\delta _{jl}} + {\delta _{ik}}{\delta _{jk}}} \right) -\\& \frac{{{E_{{D}}}}}{{2\left( {1 + \nu } \right){k_{{D}}^2}}}{s_{ij}}{s_{kl}}\\[-18pt] \end{aligned}$$ (13) 1.3 加载准则
由经典塑性理论可知，对于Prandtl–Reuss材料，当应力点处于塑性状态时，若加载，应力增量将指向屈服面切线方向；若卸载，应力增量将指向屈服面内侧。加载准则定义为：
$$ f = 0\mathop {}\limits_{}^{} ,\mathop {}\limits_{}^{} L = \frac{{\partial f}}{{\partial {\sigma _{ij}}}}{\text{d}}{\sigma _{ij}} > 0 时\text{，}加载\text{；} $$ $$ f = 0\mathop {}\limits_{}^{} ,\mathop {}\limits_{}^{} L = \frac{{\partial f}}{{\partial {\sigma _{ij}}}}{\text{d}}{\sigma _{ij}} < 0 时\text{，}卸载。 $$ 由加载函数L可简便地判断Prandtl–Reuss材料应力点的加–卸载状态。但对于圆钢管损伤本构模型，当应力点处于塑性状态时，加载将引起屈服面的收缩，应力增量将指向屈服面内侧，即加载与卸载的应力增量方向相同。这导致在应力空间中无法区分应力点的加–卸载状态。因此，对于该本构模型，需转换至应变空间[29]，用应变增量代替应力增量，才能准确判断应力点的加–卸载状态。其加载准则定义为：
$$ f = 0\mathop {}\limits_{}^{} ,\mathop {}\limits_{}^{} L = \frac{{\partial f}}{{\partial {\sigma _{ij}}}}{\boldsymbol{C}}_{ijkl}^{\text{e}}{\text{d}}{\varepsilon _{ij}} > 0 时\text{，}加载\text{；} $$ $$ f = 0\mathop {}\limits_{}^{} ,\mathop {}\limits_{}^{} L = \frac{{\partial f}}{{\partial {\sigma _{ij}}}}{\boldsymbol{C}}_{ijkl}^{\text{e}}{\text{d}}{\varepsilon _{ij}} = 0 时\text{，}中性变载\text{；} $$ $$ f = 0\mathop {}\limits_{}^{} ,\mathop {}\limits_{}^{} L = \frac{{\partial f}}{{\partial {\sigma _{ij}}}}{\boldsymbol{C}}_{ijkl}^{\text{e}}{\text{d}}{\varepsilon _{ij}} < 0 时\text{，}卸载。 $$ 1.4 流动法则
Prandtl–Reuss材料模型的流动法则要求塑性应变增量矢量
$ {\text{d}}{\boldsymbol{\varepsilon }}_{ij}^{\text{p}} $ 在塑性势能面的法线方向[30]，其计算式为：$$ {\text{d}}{\boldsymbol{\varepsilon}} _{ij}^{\text{p}} = {\text{d}}\lambda \frac{{\partial g}}{{\partial {\sigma _{ij}}}} = {s_{ij}}{\text{d}}\lambda $$ (14) 其中：
$$ {\text{d}}\lambda = \frac{1}{H}{s_{ij}}{\boldsymbol{C}}_{ijkl}^{\text{e}}{\text{d}}{{\boldsymbol{\varepsilon}} _{kl}} $$ (15) $$ H = {s_{ij}}{\boldsymbol{C}}_{ijkl}^{\text{e}}{s_{kl}} $$ (16) 式（14）～（15）中，λ为Lame常数。
损伤的引入并未对流动法则造成实质改变，故圆钢管损伤本构模型仍采用以上的流动法则。在数值计算中，由损伤弹性刚度矩阵
$ {\boldsymbol{C}}_{}^{\text{e}} $ 进一步可得等效塑性应变增量为：$$ {\text{d}}\bar {{{\boldsymbol{\varepsilon}} _{\text{p}}}} = \sqrt {\frac{2}{3}{s_{ij}}{s_{kl}}} {\text{d}}\lambda $$ (17) 1.5 应力跌落
采用Prandtl–Reuss材料模型，对处于塑性状态的应力点进行加载时，其会沿着屈服面进行塑性流动；采用圆钢管损伤本构模型，对处于塑性状态的应力点进行加载时，塑性应变增加将导致屈服面收缩，即发生类似应变软化材料的应力跌落现象[31]，且应力点从初始屈服面向后继屈服面的跌落是瞬时完成的。在此，借鉴应变软化材料对应力跌落问题的处理[30,32]：设定应力跌落过程中各应力偏量分量的原有比例保持不变，即应力跌落是从初始屈服面沿着径向向后继屈服面跌落。应力跌落并不引起球应力变化，偏应力的变化量即为应力跌落过程中应力的变化量。
设初始屈服面上的应力水平为fy，后继屈服面上的应力水平为fD，残余强度系数为α，跌落过程中的偏应力变化量与应力变化量为：
$$ {\text{d}}s_{ij}^{{\text{drop}}} = - \left( {1 - \alpha } \right){s_{ij}} $$ (18) $$ \alpha = {f_{{D}}}/{f_{\text{y}}} $$ (19) $$ {\text{d}}\sigma _{ij}^{{\text{drop}}} = {\text{d}}s_{ij}^{{\text{drop}}} $$ (20) 在应力跌落过程中，应变总量保持不变，塑性应变的增加量等于弹性应变的减少量，则此过程中的等效塑性应变增量[16]为：
$$ {\text{d}}{\overline {\boldsymbol{\varepsilon}} _{{\text{p,drop}}}} = \left( {1 - \alpha } \right){f_{\text{y}}}/{E_{{D}}} $$ (21) 由于此过程的应力跌落是瞬时完成的，故应力增量只与弹性刚度矩阵有关。弹性应变增量可由应力增量与弹性柔度矩阵求出，且塑性应变的增加量等于弹性应变的减少量，则有：
$$ {\text{d}}{\boldsymbol{\varepsilon}} _{ij}^{{\text{p,drop}}} = - {\text{d}}{\boldsymbol{\varepsilon}} _{ij}^{{\text{e,drop}}} = - D_{ijkl}^{\text{e}}{\text{d}}\sigma _{ij}^{{\text{drop}}} $$ (22) 其中：
$$ D_{ijkl}^{\text{e}} = \frac{{1 + \nu }}{{2E}}\left( { - \frac{{2\nu }}{{1 + \nu }}{\delta _{ij}}{\delta _{kl}} + {\delta _{ik}}{\delta _{jl}} + {\delta _{ik}}{\delta _{jk}}} \right) $$ (23) 由式（18）～（22）可看出，应力跌落后会产生额外的塑性应变，造成后继屈服面进一步的收缩，屈服面的收缩又会导致应力的再一次跌落。此时，采用数值方法进行迭代求解。当结果在容差范围之内时，应力跌落过程结束。此后应力点将继续在后继屈服面上进行塑性流动，该过程与Prandtl–Reuss材料模型相似。
2. 圆钢管损伤本构模型数值积分格式
本构模型数值计算的目标为：在第m–1个荷载步所有状态变量均已知的情况下，对于第m个荷载步的应变增量Δε，计算出应力增量Δσ。根据圆钢管损伤本构模型的力学行为特点，按以下步骤进行数值计算：
步骤1 根据第m–1步状态变量，计算损伤因子m–1D、损伤弹性刚度矩阵m–1Ce、损伤屈服强度m–1fD、损伤弹性刚度矩阵m–1Cep（左上标m–1表示第m–1个荷载步计算结果）。
步骤2 利用Von Mises屈服准则判断应力点的应力状态。若应力状态为塑性，转至步骤3；若应力状态为弹性，转至步骤5。
步骤3 利用加载准则判断加–卸载情况。若加载，转至步骤4；若卸载，转至步骤8。
步骤4 应力跌落与塑性流动：
1）由式（14）～（17）计算由应变增量
${\Delta {\boldsymbol{\varepsilon}} } $ 引起的塑性应变增量与等效塑性应变增量初始值及塑性应变与等效塑性应变初始值，用矩阵形式表达为：$$ {}^m{ {\delta {{\boldsymbol{\varepsilon}} ^{\text{p}}}} ^{\left( 0 \right)}} = {\text{d}}\lambda \cdot {}^{m - 1} {{{\boldsymbol{s}}}} $$ (24) $$ ^m\delta {{ {{\overline{\boldsymbol{\varepsilon}}} _{\text{p}}}} ^{\left( 0 \right)}} = \sqrt {\frac{2}{3} \cdot {}^{m - 1}{{ {\boldsymbol{s}} }^{\text{T}}} \cdot {}^{m - 1} {\boldsymbol{s}} } \cdot {\text{d}}\lambda $$ (25) $$ {\text{d}}\lambda = \frac{1}{H}\frac{{{}^{m - 1}{{ {\boldsymbol{s}} }^{\text{T}}} \cdot {}^{m - 1} {{{\boldsymbol{C}}^{\text{e}}}} \cdot {\Delta {\boldsymbol{\varepsilon}} } }}{{{}^{m - 1}{{ {\boldsymbol{s}} }^{\text{T}}} \cdot {}^{m - 1} {{{\boldsymbol{C}}^{\text{e}}}} \cdot {}^{m - 1} {\boldsymbol{s}} }} $$ (26) $$ {}^m{{{{\boldsymbol{\varepsilon}} ^{\text{p}}}} ^{\left( 0 \right)}} = {}^{m - 1} {{{\boldsymbol{\varepsilon}} ^{\text{p}}}} $$ (27) $$ ^m{ {{\overline{\boldsymbol{\varepsilon}} _{\text{p}} ^{\left( 0 \right)}}}}{ = ^{m - 1}} {{\overline{\boldsymbol{\varepsilon}} _{\text{p}}}} $$ (28) 式（24）～（28）中，
${}^{m - 1} {\boldsymbol{s}}$ 、${}^{m - 1} {{\overline {\boldsymbol{\varepsilon}}} ^{\text{p}}}$ 及$^{m - 1} {{\overline{\boldsymbol{\varepsilon}} _{\text{p}}}}$ 分别为第m–1个荷载步结束时的偏应力、塑性应变及等效塑性应变。2）建立塑性应变增量与等效塑性应变增量迭代公式。在第k次迭代中，令
$$ {}^m{ {{{\boldsymbol{\varepsilon}} ^{\text{p}}}} ^{\left( k \right)}} = {}^m{ {{{\boldsymbol{\varepsilon}} ^{\text{p}}}} ^{\left( {k - 1} \right)}} + {}^m{ {\delta {{\boldsymbol{\varepsilon }}^{\text{p}}}} ^{\left( k \right)}} $$ (29) $$ {}^m{ {{\overline{\boldsymbol{\varepsilon}} _{\text{p}} ^{\left( k \right)}}}} = {}^m{ {{\overline{\boldsymbol{\varepsilon}} _{\text{p}}}} ^{\left( {k - 1} \right)}} + {}^m\delta { {{\overline{\boldsymbol{\varepsilon}} _{\text{p}} ^{\left( k \right)}}}} $$ (30) 式中，第m个荷载步中第k次迭代的塑性应变增量
${}^m{\left\{ {\delta {{\boldsymbol{\varepsilon}} ^{\text{p}}}} \right\}^{\left( k \right)}}$ 及等效塑性应变增量${}^m\delta { {\overline{{\boldsymbol{\varepsilon}}} _{\text{p}} ^{\left( k \right)}}}$ 可由式（18）～（23）求出，右上标（k）表示第k次迭代计算结果。由式（29）、（30）以及步骤4和6，可计算第k次迭代的损伤屈服强度
${}^m{{\boldsymbol{f}}_{{D}}^{\left( k \right)}}$ 。3） 检查一致性条件：
$$ {}^m{F^{\left( k \right)}} = {}^{m - 1}\overline \sigma - \left\| {{}^m{{ {\Delta {\boldsymbol{\sigma}} _{ij}^{{\text{drop}}}} }^{\left( k \right)}}} \right\| - {}^m{{\boldsymbol{f}}_{{D}}^{\left( k \right)}} $$ (31) 其中，
$\left\| {{}^m{{ {\Delta {\boldsymbol{\sigma}} _{ij}^{{\text{drop}}}} }^{\left( k \right)}}} \right\|$ 为第m个荷载步第k次迭代的应力增量范数，由式（18）～（20）求得。若mF(k)小于等于容差，计算收敛，迭代结束；若mF(k)大于容差，则重复步骤2～4，直到收敛为止。
4） 根据迭代计算结果，更新损伤弹性刚度矩阵mCe与弹塑性刚度矩阵mCep，并用数值积分方法计算应力点在后继屈服面上的塑性流动，应力增量为：
$$ {\Delta {\boldsymbol{\sigma}} } = \int_{ {{{\boldsymbol{\varepsilon}} ^{\left( {\text{m}} \right)}}} }^{ {{{\boldsymbol{\varepsilon}} ^{\left( {\text{m}} \right)}}} + {\Delta {\boldsymbol{\varepsilon}} } } {\mathop {}\limits_{}^{} {}^m {{{\boldsymbol{C}}^{{\text{ep}}}}} } \mathop {}\limits_{}^{} {{\text{d}}{\boldsymbol{\varepsilon}} } $$ (32) 5） 塑性流动结束后，转至步骤9。
步骤5 计算试探应力：
$$ {}^m {{\boldsymbol{\sigma}} _{{\text{trial}}}^{{\text{dev}}}} = {}^m {\boldsymbol{\sigma}} + {}^{m - 1} {{{\boldsymbol{C}}^{\text{e}}}} {\Delta {\boldsymbol{\varepsilon}} } $$ (33) 步骤6 利用Von Mises屈服准则判断试探应力点是否屈服。若屈服，转至步骤7；若未屈服，转至步骤8。
步骤7 对于圆钢管损伤本构模型，借鉴J2流动理论，采用径向返回算法对应力进行更新，如图1所示。
其迭代过程为[34]：
1）计算塑性流动方向单位矢量：
$$ {\boldsymbol{ n}} = \frac{1}{{{}^m{{\overline {\boldsymbol{\sigma}} }^{\left( 0 \right)}}}}{}^m {{\boldsymbol{s}}_{{\text{trial}}}^{{\text{dev}}}} $$ (34) 式中，
${}^m {{\boldsymbol{s}}_{{\text{trial}}}^{{\text{dev}}}}$ 为试探应力${}^m {{\boldsymbol{\sigma}} _{{\text{trial}}}^{{\text{dev}}}}$ 的偏应力，${}^m{\overline {\boldsymbol{\sigma}} ^{\left( 0 \right)}}$ 为试探应力${}^m {{\boldsymbol{\sigma}} _{{\text{trial}}}^{{\text{dev}}}}$ 的等效应力。2）计算各塑性参数初始值：
$$ {}^m{\boldsymbol{\delta}} { {\overline{{\boldsymbol{\varepsilon}}} _{\text{p}} ^{\left( 0 \right)}}} = \frac{{{}^m{{\overline {\boldsymbol{\sigma}} }^{\left( 0 \right)}} - {}^m{{\boldsymbol{f}}_{{D}}^{\left( 0 \right)}}}}{{3\mu }} $$ (35) $$ {}^m\Delta { {\overline{{\boldsymbol{\varepsilon}}} _{\text{p}} ^{\left( 0 \right)}}} = {}^m{\boldsymbol{\delta}} { {\overline{{\boldsymbol{\varepsilon}}} _{\text{p}} ^{\left( 0 \right)}}} $$ (36) $$ {}^m{ {\overline{{\boldsymbol{\varepsilon}}} _{\text{p}} ^{\left( 0 \right)}}} = {}^{m{\text{-1}}} {\overline{{\boldsymbol{\varepsilon}}} _{\text{p}}} + {}^m{\boldsymbol{\delta}} { {\overline{{\boldsymbol{\varepsilon}}} _{\text{p}} ^{\left( 0 \right)}}} $$ (37) 3）计算损伤屈服强度。塑性参数初始值已知后，在第k次迭代中，由式（4）和（6）计算损伤屈服强度
${}^m{{\boldsymbol{f}}_{{D}}^{\left( k \right)}}$ 。4）建立塑性参数迭代公式。在第k次迭代中，令
$$ {}^m{\boldsymbol{\delta}} { {\overline{{\boldsymbol{\varepsilon}}} _{\text{p}} ^{\left( k \right)}}} = \frac{{{}^m{{\overline {\boldsymbol{\sigma}} }^{\left( 0 \right)}} - 3\mu \cdot {}^m\Delta {{ {\overline{{\boldsymbol{\varepsilon}}} _{\text{p}}^{\left( {k - 1} \right)}} }} \cdot {}^m{{\boldsymbol{f}}_{{D}}^{\left( k \right)}}}}{{3\mu }} $$ (38) $$ {}^m\Delta { {\overline{{\boldsymbol{\varepsilon}}} _{\text{p}} ^{\left( k \right)}}} = {}^m\Delta { {\overline{{\boldsymbol{\varepsilon}}} _{\text{p}} ^{\left( {k - 1} \right)}}} + {\boldsymbol{\delta}} { {\overline{{\boldsymbol{\varepsilon}}} _{\text{p}} ^{\left( k \right)}}} $$ (39) $$ {}^m{{\overline {{\boldsymbol{\varepsilon}}} _{\text{p}} ^{\left( k \right)}}} = {}^m{ {\overline {{\boldsymbol{\varepsilon}}} _{\text{p}} ^{\left( {k - 1} \right)}}} + {}^m{\boldsymbol{\delta}} { {\overline {{\boldsymbol{\varepsilon}}} _{\text{p}} ^{\left( k \right)}}} $$ (40) 5）检查一致性条件：
$$ {}^m{F^{\left( k \right)}} = {}^m{\overline \sigma ^{\left( 0 \right)}} - 3\mu \cdot {}^m\Delta {{\overline {{\boldsymbol{\varepsilon}}} _{\text{p}} ^{\left( {k - 1} \right)}}} \cdot {}^m{{\boldsymbol{f}}_{{D}}^{\left( k \right)}} $$ (41) 若mF(k)小于等于容差，计算收敛，迭代结束；若mF(k)大于容差，则重复步骤2 ～ 4，直到收敛为止。
6）若计算在第n次迭代收敛，则应力点在第m个荷载步结束时的应力为：
$$ {}^m {\boldsymbol{\sigma}} = {}^m{{\boldsymbol{f}}_{{D}}^{\left( n \right)}} {\boldsymbol{n}} + {}^m {p_{{\text{trial}}}^{{\text{dev}}}} $$ (42) 式中，
${}^m {p_{{\text{trial}}}^{{\text{dev}}}}$ 为试探应力的球应力。7）塑性回流结束后，转至步骤9。
步骤8 纯弹性行为，无塑性应变增量，利用弹性刚度矩阵求解，应力增量为：
$$ {\Delta {\boldsymbol{\sigma }}} = {}^{m - 1} {{{\boldsymbol{C}}^{\text{e}}}} {\Delta {\boldsymbol{\varepsilon}} } $$ (43) 步骤9 更新各状态变量，第m个荷载步结束。
3. 考虑材料损伤累积的动力响应分析
3.1 材料子程序
根据圆钢管损伤本构模型的力学行为特点与数值积分格式，基于有限元软件ANSYS，编制考虑材料损伤累积的用户材料子程序usermat。
usermat的主要任务是定义材料的应力–应变关系。具体分为两方面：由给定的应变增量
$ \Delta {\boldsymbol{\varepsilon}}$ ，需计算应力增量$ \Delta{\boldsymbol{ \sigma}} $ ，从而得到新的应力；需给出Jacobian矩阵$ {{\boldsymbol{D}}_{\text{J}}} = \partial \Delta{\boldsymbol{ \sigma}} /\partial \Delta {\boldsymbol{\varepsilon }}$ 。usermat在单元的每个高斯点上调用，且通过接口与主程序连接。增量步开始时，ANSYS主程序通过接口进入usermat，并将高斯点相关变量的初始值以输入变量的形式传递给usermat；增量步结束时，usermat将相关变量的更新值以输出变量的形式传递给ANSYS主程序。根据圆钢管损伤本构模型的力学行为特点与数值积分格式编制usermat流程如图2所示。将此材料子程序嵌入ANSYS主程序中，并结合网格结构在地震作用下的常规分析步骤，便可实现网格结构考虑材料损伤积累的地震响应分析。
图 2 usermat流程图Fig. 2 Flow chart of usermat
下载:
全尺寸图片
3.2 计算模型
以工程中常见的K8型单层球壳为例，跨度50 m，矢跨比1/4，屋面荷载2 kN/m2，以第1阶振型模态模拟结构初始几何缺陷分布模式，初始几何缺陷最大计算值取跨度的1/300，支座为周边固定铰支承。主肋杆与环杆截面φ140 mm×5 mm，斜杆截面φ121 mm×3.5 mm。采用Beam189单元对杆件进行模拟，每根杆件划分为2段，杆件截面上共划分32个积分点，为了说明截面塑性的发展程度，定义1P表示至少1个积分点进入塑性，8P表示至少8个积分点进入塑性。依次类推，32P表示全截面进入塑性。单元截面塑性发展定义如图3所示。
图 3 Beam189单元截面塑性发展定义图Fig. 3 Plastic expansion in element section of Beam189下载:
全尺寸图片
钢材初始弹性模量206×103 N/mm2，泊松比0.3，质量密度7 850 kg/m3，初始屈服强度235 N/mm2，利用开发的材料子程序模拟钢管材料的力学行为，并在计算中考虑几何非线性。分析模型如图4所示。
图 4 计算模型Fig. 4 Calculation model下载:
全尺寸图片
3.3 结果分析
选用Taft地震波三向输入，以逐步增大地面峰值加速度（PGA）的方式，通过对比考虑损伤积累结构与Prandtl–Reuss材料模型结构（以下简称P–R材料结构）在节点位移、杆件应变延性系数、进入塑性杆件比例及分布等响应指标的异同，从宏观和微观两个层面，探讨考虑材料损伤累积对网壳结构抗震性能的影响。结构各响应指标结果如图5所示。图5中，以“16P–P”为例，16P表示至少16个积分点进入塑性，P表示P–R材料结构，D表示考虑损伤累积结构。
图 5 各项响应指标结果Fig. 5 Results of responses下载:
全尺寸图片
3.3.1 节点位移
节点位移从宏观层面表征结构的变形发展与刚度变化情况，是描述强震作用下结构抗震性能最为直接的响应指标之一[35]。节点最大竖向位移–PGA曲线如图5（a）所示，节点最大竖向位移时程曲线如图5（b）所示。
由图5（a）可知，节点竖向最大位移呈现出两个阶段：PGA从100g增大至300g过程中，P–R材料结构与考虑损伤累积结构的位移响应值差距很小，两者几乎相等；PGA为400g之后，两者差距逐渐增大。这是因为PGA较小时，结构塑性发展较浅，材料损伤累积较轻，损伤累积对结构性能的影响不甚明显；PGA逐渐增大后，结构塑性深度发展，材料损伤严重积累，引起结构整体刚度下降，节点竖向最大位移响应逐渐增大。
由图5（b）可知，考虑损伤积累结构的节点竖向最大位移时程曲线与P–R材料结构相比，震荡幅值更大。这是因为考虑损伤累积后，随着地震作用的进行，材料弹性模量不断减小，强度不断软化，构件刚度随之降低，进而导致整体刚度的降低，故其位移响应更强烈，震荡更激烈。对于P–R材料结构，PGA为980g时，结构振动尚能维持收敛状态；当PGA为1 000g时，结构振动发散，且最大值已超过100 cm。再结合图5（a）可认为：P–R材料结构的动力极限荷载为980g；考虑损伤积累结构的动力极限荷载为860g，较前者下降12.24%。说明采用P–R材料模型将高估网壳结构的抗震性能，材料损伤累积的影响不可忽视。
3.3.2 杆件应变延性系数
定义杆件应变延性系数为杆件轴向总应变与弹性极限应变的比值，其可从微观层面反映杆件的塑性发展程度与延性变形能力。杆件应变延性系数–PGA曲线如图5（c）所示。
由图5（b）可知：对于P–R材料结构，当PGA为100g时，杆件最大应变延性系数小于1，说明杆件均处于弹性状态；从200g开始，杆件逐渐进入塑性，曲线斜率逐渐变小。对于考虑材料损伤累积结构，曲线呈现相同的趋势，但随着PGA的增大，与P–R材料结构的差值逐渐增大，曲线斜率变得更小；当PGA由860g小幅增长至900g时，杆件最大应变延性系数由91.90大幅增长至143.47，而P–R材料结构在失效前后的杆件最大应变延性系数分别为95.17和129.91。说明损伤累积的引入使得杆件截面的塑性发展更为充分，失效时杆件的变形也更为严重。
3.3.3 塑性杆件比例与分布
结构塑性发展程度可从塑性杆件比例及分布两方面得以体现[36]。由图5（d）可见，无论是否考虑损伤积累，各类塑性杆件比例均随PGA的增大而增加。考虑损伤积累后，当PGA较小时，各类塑性杆件比例相差不大；PGA从400g之后，两者差距逐渐显现，且1P杆件的差距较明显而32P杆件的差距较小。再结合失效时结构塑性分布（图5（e）、（f））可知：一方面，损伤因子的引入导致结构内力发生重分布，使原本的弹性杆件更容易进入塑性；另一方面，损伤累积效应削弱了杆件的刚度，使杆件的塑性程度进一步深入发展。
此外，对于考虑材料损伤累积结构，当PGA为860g时，1P、8P、16P及32P杆件的比例分别为88.47%、86.69%、83.12%及64.77%；当PGA增大为900g时，以上各值变化为97.26%、96.59%、96.20%及95.36%。说明此时的绝大多数杆件已全截面屈服，形成大量的塑性区带，结构的性能发生改变，导致位移异常增大，节点振动发散。
4. 结 论
本文以圆钢管考虑损伤累积本构模型为研究对象，分析其力学行为特点，明确其数值积分格式，并基于有限元软件ANSYS开发了考虑损伤累积的材料子程序，采用该子程序研究单层网壳结构在强震作用下的动力响应，得出如下结论：
1）本文基于圆钢管考虑损伤累积本构模型的力学行为特点提出了完整的数值积分格式，并通过计算验证了其合理性与有效性。
2）与Prandtl–Reuss材料模型结构相比，考虑材料损伤累积一方面会使结构的破坏模式由动力失稳转变为强度破坏，另一方面会降低结构的动力极限承载能力。
3）考虑材料损伤累积的结构失效时，绝大多数杆件全截面屈服，形成大量的塑性区带，导致位移迅速增大，结构破坏。
4）计算表明，本文编制的圆钢管考虑损伤积累的材料子程序计算精度较高，可用于单层网壳结构的动力响应分析。
-
图 2 usermat流程图
Fig. 2 Flow chart of usermat
下载:
全尺寸图片
图 3 Beam189单元截面塑性发展定义图
Fig. 3 Plastic expansion in element section of Beam189
下载:
全尺寸图片
图 4 计算模型
Fig. 4 Calculation model
下载:
全尺寸图片
图 5 各项响应指标结果
Fig. 5 Results of responses
下载:
全尺寸图片
-
[1] 于志伟,郑世杰,甄翠贤,等.单层柱面铝合金网壳结构强震失效机理及地震易损性[J].建筑结构学报,2020,41(增刊1):17–24. Yu Zhiwei,Zheng Shijie,Zhen Cuixian,et al.Failure mechanism and seismic vulnerability of single-layer aluminum alloy cylindrical reticulated shells[J].Journal of Building Structures,2020,41(Supp1):17–24 [2] 李玉刚,范峰,洪汉平.基于小样本记录的柱面网壳结构地震响应评估[J].工程力学,2020,37(5):228–236. doi: 10.6052/j.issn.1000-4750.2019.07.0418 Li Yugang,Fan Feng,Hong Hanping.Evaluating the seismic effects on a cylindrical lattice shell using a small number of records[J].Engineering Mechanics,2020,37(5):228–236 doi: 10.6052/j.issn.1000-4750.2019.07.0418 [3] 支旭东,张英楠,范峰,等.单层球面网壳的多重地震效应研究[J].工程力学,2018,35(9):107–113. doi: 10.6052/j.issn.1000-4750.2017.05.0370 Zhi Xudong,Zhang Yingnan,Fan Feng,et al.Responses of single-layer reticulated domes subject to multiple earthquakes[J].Engineering Mechanics,2018,35(9):107–113 doi: 10.6052/j.issn.1000-4750.2017.05.0370 [4] 段红霞,李守巨,刘迎曦.地震作用下钢结构损伤过程数值模拟[J].工程力学,2011,28(2):198–204. Duan Hongxia,Li Shouju,Liu Yingxi.Numerical simulation of damage process of steel frames under earthquake excitation[J].Engineering Mechanics,2011,28(2):198–204 [5] Yu Zhiwei,Zhi Xudong,Fan Feng,et al.Effect of substructures upon failure behavior of steel reticulated domes subjected to the severe earthquake[J].Thin-Walled Structures,2011,49(9):1160–1170. doi: 10.1016/j.tws.2011.05.002 [6] Li Q S,Chen J M.Nonlinear elastoplastic dynamic analysis of single-layer reticulated shells subjected to earthquake excitation[J].Computers & Structures,2003,81(4):177–188. doi: 10.1016/S0045-7949(02)00445-5 [7] 李永梅,郑文豪,胡琨.余震对单层网壳结构动力损伤累积效应的影响分析[J].空间结构,2016,22(4):23–28. doi: 10.13849/j.issn.1006-6578.2016.04.023 Li Yongmei,Zheng Wenhao,Hu Kun.Analysis of damage accumulation effect of single-layer latticed shells considering aftershock action[J].Spatial Structures,2016,22(4):23–28 doi: 10.13849/j.issn.1006-6578.2016.04.023 [8] 巴盼锋,张毅刚,吴金志,等.单层球面网壳动力失效全过程试验研究[J].振动与冲击,2017,36(1):31–37. doi: 10.13465/j.cnki.jvs.2017.01.005 Ba Panfeng,Zhang Yigang,Wu Jinzhi,et al.Whole process test and dynamic failure analysis for single-layer spherical lattice shells[J].Journal of Vibration and Shock,2017,36(1):31–37 doi: 10.13465/j.cnki.jvs.2017.01.005 [9] Kachanov L M.Time of the rupture process under creep conditions[J].Lzvestia AN’SSSR.Metally,1958,8(1):26–31. [10] Jean L.A continuous damage mechanics model for ductile fracture[J].Journal of Engineering Materials and Technology,1985,107(1):83–89. doi: 10.1115/1.3225775 [11] Hammi Y,Horstemeyer M F.A physically motivated anisotropic tensorial representation of damage with separate functions for void nucleation,growth,and coalescence[J].International Journal of Plasticity,2007,23(10/11):1641–1678. doi: 10.1016/j.ijplas.2007.03.010 [12] Teng X.Numerical prediction of slant fracture with continuum damage mechanics[J].Engineering Fracture Mechanics,2008,75(8):2020–2041. doi: 10.1016/j.engfracmech.2007.11.001 [13] Xue Liang,Wierzbicki T.Ductile fracture initiation and propagation modeling using damage plasticity theory[J].Engineering Fracture Mechanics,2008,75(11):3276–3293. doi: 10.1016/j.engfracmech.2007.08.012 [14] Malcher L,Andrade Pires F M,de Sá J M A C.An assessment of isotropic constitutive models for ductile fracture under high and low stress triaxiality[J].International Journal of Plasticity,2012,30/31:81–115. doi: 10.1016/j.ijplas.2011.10.005 [15] Malcher L,Mamiya E N.An improved damage evolution law based on continuum damage mechanics and its dependence on both stress triaxiality and the third invariant[J].International Journal of Plasticity,2014,56:232–261. doi: 10.1016/j.ijplas.2014.01.002 [16] Xue Liang.Damage accumulation and fracture initiation in uncracked ductile solids subject to triaxial loading[J].International Journal of Solids and Structures,2007,44(16):5163–5181. doi: 10.1016/j.ijsolstr.2006.12.026 [17] Kumar S,Usami T.An evolutionary-degrading hysteretic model for thin-walled steel structures[J].Engineering Structures,1996,18(7):504–514. doi: 10.1016/0141-0296(95)00167-0 [18] 沈祖炎,沈苏.高层钢结构考虑损伤累积及裂纹效应的抗震分析[J].同济大学学报(自然科学版),2002,30(4):393–398. doi: 10.3321/j.issn:0253-374X.2002.04.003 Shen Zuyan,Shen Su.Seismic analysis of tall steel structures with damage cumulation and fracture effects[J].Journal of Tongji University (Natural Science),2002,30(4):393–398 doi: 10.3321/j.issn:0253-374X.2002.04.003 [19] 董宝,沈祖炎,孙飞飞.考虑损伤累积影响的钢柱空间滞回过程的仿真[J].同济大学学报(自然科学版),1999,27(1):11–15. doi: 10.3321/j.issn:0367-6234.2008.12.005 Dong Bao,Shen Zuyan,Sun Feifei.Simulation of spatial hysteretic behavior of steel columns considering effects of damage accumulation[J].Journal of Tongji University(Natural Science),1999,27(1):11–15 doi: 10.3321/j.issn:0367-6234.2008.12.005 [20] 王萌,石永久,王元清.强震作用下钢材等效本构模型应用研究[J].建筑结构学报,2013,34(10):84–92. Wang Meng,Shi Yongjiu,Wang Yuanqing.Application study on equivalent constitutive model of steel subjected to strong earthquake[J].Journal of Building Structures,2013,34(10):84–92 [21] 闫翔宇,巩昊,陈志华,等.H形钢单层网壳结构形态优化研究[J].建筑钢结构进展,2021,23(12):101–108. doi: 10.13969/j.cnki.cn31-1893.2021.12.011 Yan Xiangyu,Gong Hao,Chen Zhihua,et al.Study on the shape optimization of H-shaped steel single-layer lattice shells[J].Progress in Steel Building Structures,2021,23(12):101–108 doi: 10.13969/j.cnki.cn31-1893.2021.12.011 [22] 范峰,聂桂波,支旭东.三向荷载作用下圆钢管材料本构模型研究[J].建筑结构学报,2011,32(8):59–68. Fan Feng,Nie Guibo,Zhi Xudong.Constitutive model of circular steel tubes under complicated cyclic load[J].Journal of Building Structures,2011,32(8):59–68 [23] 范峰,聂桂波,支旭东,等.圆钢管空间滞回试验及材料本构模型[J].土木工程学报,2011,44(12):18–24. Fan Feng,Nie Guibo,Zhi Xudong,et al.Spatial hysteresis experiment and constitutive model for circular steel pipes[J].China Civil Engineering Journal,2011,44(12):18–24 [24] Chen W F,Han D J.Introduction[M]//Plasticity for Structural Engineers.New York:Springer New York,1988:3–45. [25] He Sheng,Nie Yongchao,Bordas S P A,et al.Damage-plastic constitutive model of thin-walled circular steel tubes for space structures[J].Journal of Engineering Mechanics,2020,146(12):4020131. doi: 10.1061/(asce)em.1943-7889.0001844 [26] He Sheng,Wang Haosen,Bordas S P A,et al.A developed damage constitutive model for circular steel tubes of reticulated shells[J].International Journal of Structural Stability and Dynamics,2020,20(9):2050106. doi: 10.1142/s0219455420501060 [27] 李永梅,胡琨,张微敬.考虑损伤累积效应的单层球面网壳动力稳定[J].湖南大学学报(自然科学版),2014,41(6):16–21. doi: 10.1177/136943329700100105 Li Yongmei,Hu Kun,Zhang Weijing.Dynamic stability of single-layer latticed shells considering damage accumulation effect[J].Journal of Hunan University (Natural Sciences),2014,41(6):16–21 doi: 10.1177/136943329700100105 [28] 周海涛,张毅刚,吴金志.基于材料塑性损伤累积效应的梁单元数值模拟分析方法[J].空间结构,2010,16(3):13–17. Zhou Haitao,Zhang Yigang,Wu Jinzhi.Numerical simulation method considering the cumulative effect of plastic damage for beam element[J].Spatial Structures,2010,16(3):13–17 [29] 张舒原,吴运新,龚海.航空铝合金应变空间弹塑性本构模型[J].中南大学学报(自然科学版),2012,43(7):2580–2585. doi: 10.3321/j.issn:1001-9669.2008.05.020 Zhang Shuyuan,Wu Yunxin,Gong Hai.Elasto-plasticity constitutive model for aerospace aluminium alloy in strain space[J].Journal of Central South University,2012,43(7):2580–2585 doi: 10.3321/j.issn:1001-9669.2008.05.020 [30] 沈新普,岑章志,徐秉业.弹脆塑性软化本构理论的特点及其数值计算[J].清华大学学报(自然科学版),1995,35(2):22–27. Shen Xinpu,Cen Zhangzhi,Xu Bingye.The characteristics of elasto-brittle-plastic softening constitutive theory and its numerical calculation[J].Journal of Tsinghua University (Science and Technology),1995,35(2):22–27 [31] 刘杰民.二元扰动有限元法和应变软化材料的数值模拟[J].哈尔滨工业大学学报,2011,43(增刊1):213–216. Liu Jiemin.FEM of duality disturbance and numerical simulation for materials with strain softening[J].Journal of Harbin Institute of Technology,2011,43(Supp1):213–216 [32] Vignjevic R,Djordjevic N,De Vuyst T,et al.Modelling of strain softening materials based on equivalent damage force[J].Computer Methods in Applied Mechanics and Engineering,2018,335:52–68. doi: 10.1016/j.cma.2018.01.049 [33] 聂桂波.网壳结构基于损伤累积本构强震失效机理及抗震性能评估[D].哈尔滨:哈尔滨工业大学,2012. Nie Guibo.Failure mechanism and seismic performance evaluation of reticulated shell structure under strong earthquake based on damage accumulation constitutive model[D].Harbin:Harbin Institute of Technology,2012. [34] 庄茁,由小川,廖剑晖.基于ABAQUS的有限元分析和应用[M].北京:清华大学出版社,2009. [35] 马人乐,赵林,何敏娟.大型超高钢结构电视塔抗震性能分析[J].同济大学学报(自然科学版),2007,35(10):1310–1315. doi: 10.3321/j.issn:0253-374X.2007.10.004 Ma Renle,Zhao Lin,He Minjuan.Earthquake resistant behavior of super high-rising steel tower[J].Journal of Tongji University(Natural Science),2007,35(10):1310–1315 doi: 10.3321/j.issn:0253-374X.2007.10.004 [36] 汤宏伟,钟宏志.考虑杆件初弯曲的网壳弹塑性稳定性的 弱形式求积元分析[J].工程力学,2019,36(1):165–174. doi: 10.6052/j.issn.1000-4750.2017.11.0842 Tang Hongwei,Zhong Hongzhi.Elasto-plastic stability analysis of latticed sheels with initially curved members by weak-form quadrature elements[J].Engineering Mechanics,2019,36(1):165–174 doi: 10.6052/j.issn.1000-4750.2017.11.0842
