常见的工程材料(如混凝土、钢材)为依赖应变历史材料(SHDM)，此类材料的正应变增量过程(应变加载)中的切线模量，不同于负应变增量过程(应变卸载)中的切线模量，即材料点应变在加载和卸载(下文简称“应变加卸载”)时具有不同的应力变化规律。当不考虑其他非线性因素时，SHDM结构(SHDM制成的结构)的有限元静力问题是一个材料非线性问题。此静力问题完成求解时，会获得满足静力平衡的荷载-位移曲线(以下简称平衡路径)。平衡路径具有硬化段，但当材料兼具应变软化性质时(如混凝土)，它可能还具有极值点、软化段。此静力问题进行求解时：解法一般采用增量迭代法，即总体上采用增量法求解，增量步采用迭代法求解[1]；迭代时，为保证收敛解为真实静力解，用以更新单元内高斯积分点(以下简称高斯点)应力的应变须仅经历过加载或卸载。
隐式静力算法[2](ISA)、动力松弛法[3](DRM)是两种求解静力问题增量步的迭代算法。对迭代中的应力更新步骤，必须采用增量步内的结点位移累积增量计算高斯点应力增量[1]，否则会引入虚假应变历史，从而使得满足静力平衡的迭代收敛解为虚假静力解[4]。此外，静力问题也可采用显式拟静力算法(EQSA)近似求解[5]。
采用ISA迭代求静力增量步时，需在预判高斯点应变加卸载状态后进行总体刚度矩阵K的组装和求逆，但迭代是否收敛依赖于材料特征和K特征。平衡路径硬化段内，K正定；极值点处，K奇异；软化段内，K负定。荷载控制类ISA在ABAQUS有限元软件中为一般非线性问题的推荐算法[6]，但它不能求解K非正定的静力问题[7]。位移控制类ISA、整体弧长类ISA不适合求解应变软化结构的静力问题[2－8]，前者还难以合适选定控制位移分量[9]，后者还难以跨越材料非线性引起的极值点[10]。局部弧长类ISA目前算例主要考虑受拉应变软化情形[11－13]。因此，对于复杂静力分析问题，文献[2]建议先将其转化为动力问题，再用EQSA给出静力近似解。
EQSA本质上是动力问题的显式时间积分算法；当通过低加载速率使动力系统的惯性效应满足拟静力判据时，显式时间积分算法可视为EQSA[2]。采用EQSA近似求解静力问题时，会在真实时间域向量式地更新结点位移a。因此，EQSA的优点为：无需预判高斯点应变加卸载状态和构造K[6]。但是，EQSA的局限性为：① 计算成本具有随结点质量的减小呈指数比例增大的规律[8]；② 合理确定粘滞阻尼系数时的经验依赖性；③ a依时间序列产生，由此给出的近似静力解中存在虚假应变历史[14]。
DRM是一种静力分析算法，其首先构造一个包含粘滞阻尼[3]或动态阻尼[15]的虚拟动力系统，然后通过此系统的动能阻尼耗散计算静力解。采用DRM求解静力增量步时：在虚拟时间域采用显式时间积分算式进行a的向量式更新，从而保留EQSA的优点；可采用结点虚拟质量缓解EQSA的第一个局限性，采用动态阻尼规避第二个局限性。此外，DRM的另一优点为：若SHDM结构兼具应变软化性质(如混凝土结构)，它能解决“应变软化结构因静力定解问题失去适定性”而导致的迭代不收敛，相关论证可见作者的先前工作[16]。Rezaieepajand等[17]进行了桁架、框架结构静力问题在施加不同虚拟阻尼时的DRM分析，结果表明kDRM(动态阻尼DRM)的计算成本最少。DRM的结构分析应用始于20世纪五六十年代[18－19]。而后应用范围逐渐扩大，详见文献[20－27]。但是，SHDM结构的DRM分析报道较少见诸于文献。
本文采用kDRM求解SHDM结构的静力平衡路径。迭代求解增量步时，为了避免虚假应变历史效应，构造了高斯点应力更新的非线性弹性增量算法(NEIA)——尽管材料一般不是非线性弹性材料，但在一个增量步内，将材料视为非线性弹性材料。全文首先给出问题的描述，然后详述问题的求解，接着描述三个算例及其计算结果，其中两个SHDM结构模型兼具应变软化特征，最后给出讨论、结论。
1 问题的描述
本文讨论的依赖应变历史材料(SHDM)是指材料点应变加卸载时应力响应不同的材料。在一维情况下，SHDM满足 [1]：
式中：σ为应力标量；ε为应变；ε0为材料点在应变加卸载过程中任意时刻的应变。在多轴应力情况下，只要应力、应变的某一分量(或等效应力、等效应变)满足式(1)时，即为本文所指的SHDM。土木工程材料中常用的塑性本构模型、损伤本构模型为SHDM的本构模型。
结构有限元分析时，自由度分为两类：其一，施加外部结点荷载的加载自由度；其二，未加载自由度。本文在进行结构非线性有限元静力分析时：结构为边界受到充分约束而不具有刚体位移的几何不变体系；主要采用比例加荷载方式。由此，基于全部自由度建立的待解非线性平衡方程组为：
式中：a为结点位移向量；p为内部结点荷载向量；λ为外部结点荷载系数；fref为外部结点荷载参照向量；ψ为结点残值荷载向量；a、p、fref均为r维向量(有限元模型的自由度数)。p按式
计算[1]，其中B为nσ×r应变矩阵，σ 为nσ×1应力向量，nσ为σ 内的元素数，Ω 为固体所占空间。
本文采用加荷载方式或加位移方式求解式(2)所描述的问题(以下简称问题2)。求解方式为加荷载时：已知量为λ；待求量为a。求解方式为加位移时：已知量为被选取的加载位移aI(下标表示第I个自由度)；待求量为aII——a中除aI外的其他元素，以及λ；记a=[aI, aII]T。
结构有限元模型内的所有高斯点的本构状态为集合
其中s为本构关系状态向量，σ 为应力，ε 为应变。向量v ={v1, v2,…, vn}的无穷范数
2 问题的求解
本文中：① 问题2总体上用荷载或位移增量法求解；② 增量步用kDRM求解——结点位移静态解借助对虚拟动力系统施加足够次的动态阻尼而获得；③ 通过高斯点应力更新的NEIA使得静态解符合静力解特征——高斯点应力由真实应变历史更新。本文中，施加动态阻尼的含义为：置零“有限元离散系统结点总动能极大值时刻的”结点速度。
2.1 静力问题的增量解法
总体上，问题2用荷载或位移增量法求解。求解时，整个平衡路径连续地分解成序列{Rn}，其中Rn(n=1,2,…)为第n个增量步。Rn的解为an、状态量Cn、λn。an、Cn中的εn、σn，以及λn可写为：
式中：Δan为结点位移增量；Δεn为应变增量；Δσn为应力增量；Δλn为外部结点荷载系数增量。
加荷载增量或者加位移增量求解Rn。默认先加荷载增量求解R1；R1施加前的初值λ0、a0、C0为已知量，λ0和C0内的σ0满足式(2)。Rn求解要点如下：
1) 已知量：① Rn-1的解λn-1、an-1、Cn-1；② 加荷载时施加的固定荷载增量系数Δλ，或加位移时对第I个自由度施加的固定位移增量ΔaI，即：
2) 未知量：① 加荷载时为an、Cn，其中an为基本未知量；② 加位移时为an, II、Δλn、Cn，其中an, II为基本未知量。
3) 未知量用kDRM迭代求解，详见2.2节。
4) 记i为kDRM求解的迭代步序号，对于求得的
① 若满足收敛判据
式中，c1为限值，则得问题的解为
加位移时，若此解符合增量变换规则：
则表明此解处在平衡路径硬化段，则后面使用加荷载方式求解Rn+1。文献[22]在采用DRM求解张拉结构平衡路径时，即采用式(3)所示的收敛判据，并取c1=0.001，此判据和限值符合文献[1]的建议。
② 加荷载时，若
满足发散判据
式中，c2为限值，则说明λn fref超越或临近平衡路径的峰值点，此时更换为加位移方式后重新求解Rn。
2.2 增量步的kDRM解法
2.2.1 kDRM的控制方程
kDRM求解静力问题的策略为[28]：在虚拟时间域监测无阻尼质点系统的动能，在各结点静力平衡前，循环执行“在动能极大值时刻先置零结点速度，再重启无阻尼运动”。kDRM求解Rn时，记第k个施加动态阻尼循环为Sk，则其内结点运动问题写为：
式中：
分别为结点的速度、加速度；Mn为对角质量矩阵；tk-1、tk分别为Sk的起始时间、结束时间；
为初始条件值。式(9)所述问题(以下简称问题9)中各参数、变量的确定方法如下(an(t)的确定方法见2.2.2)：
1) tk为在t ＞tk-1上的首个系统动能的极大值时间；tk也为循环Sk+1的起始时间。
2) tk-1为循环Sk-1的结束时间；当k=1时，tk-1=0，此时即为增量步Rn的加荷载(位移)时间。
3) Mn为元素是已知常量的矩阵，2.2.2小节给出了Mn的具体构造方法。
4) Cn(t)为元素是未知变量εn(t)、σn(t)、sn(t)的集合。εn(t)=Ban(t)；sn(t)的自变量为εn(t)；σn(t)的自变量为sn(t)及εn(t)。所以，Cn(t )是an(t)的函数。
5) pn(t)为元素是未知变量的向量，算式为[8]：
式中：e表示单元编号；ne表示系统内的单元数量；p{e}、B{e}、σ{e}、Ω{e}分别为单元e的内部结点荷载、应变矩阵、应力向量、所占空间；p{e}、B{e}、σ{e}分别为
的矩阵(
为单元e的结点自由度数)；T{e}为单元e的转换矩阵，尺寸为
(r为系统全部自由度数)。因为C(t )是n an(t )的函数，所以pn(t)也为an(t)的函数。
6) 对于λn(t)：① 加荷载时为已知的常量
② 加位移时为未知的变量
式中：fref,I、pn,I为fref、pn(t)中与第I个自由度(加位移自由度)对应的元素值；λn(t )是an(t )的函数。
7)
取Sk-1结束时刻的值，但当k=1时：① 若加位移，
由Cn-1、an-1和
算得，
由式(12)算得；② 若加荷载，
由本小节分析可知：kDRM求解由问题9所描述的Rn的关键是，给出基本未知量an(t)的数学算法。因为问题9是一个强非线性动力问题，所以求解an(t)时，宜选用显式时间积分算法[29]。
2.2.2 控制方程的中心差分算法求解
本文选用中心差分算法(CFDM)求解问题9；求解时，ih时刻
的迭代式为[30]：
式中：ik-1为Sk的起始时刻序号；h为离散时间域的固定时间增量；
的计算式为：
其中，
按式(11)或式(12)确定。用式(13)迭代求解问题9中的基本未知量an(t)时，需解决三个问题：
1) 确定Mn及h。一般采用CFDM的稳定条件构造Mn、h。Barnes给出的CFDM稳定条件[22]为：
式中：下标表示结点序号；mj为结点质量；kjmax为结点沿各坐标正负方向位移刚度最大值。本文估算kjmax的方式为：首先给结点施加微小位移，然后计算相应结点荷载，最后算出kjmax。本文取
式中：c3为系数，值大于1；mmin是结点质量下限值，其作用是避免结点负刚度而构造出负的结点质量。本文按式(16)逐结点算得mj后构造Mn。
2) 规避高斯点的虚假应变历史。本构状态集合
、a、C计算，即采用本文提出的NEIA n-1n-1更新C中
NEIA的进一步说明及与现有增量算法——非线性增量算法(NIA)的比较见2.3节。
3) 确定Sk的结束时刻。时刻(i-1/2)h、(i-3/2)h的动能分别为
若满足：
则认为(i-1)h时离散质点系统出现了动能极大值，Sk的结束时间以及Sk+1的起始时间取为(i-1)h，以
为初始条件启动Sk+1。
至此，基于式(13)能解出问题9中的基本未知量an(t)，其中：① 加荷载时，由式(11)知，
由迭代式(13)直接解出；② 加位移时，由式(12)、式(13)知，
保持不变，
由迭代式(13)直接解出，
由式(12)解出。
2.2.3 完整的kDRM数值计算流程
依据2.1节、2.2节的讨论，Rn的完整kDRM求解流程，见表1。
2.3 高斯点应力更新：非线性弹性增量算法
依赖应变历史材料(SHDM)在应变加载、卸载时具有不同材料响应，故对结构有限元模型进行受力分析时，s需记录高斯点加载、卸载时应力-应变关系的状态。
结构静力分析时，结果须符合静力解特征[1]：高斯点应力的更新必须基于增量步起始状态，即基于
而不能基于增量步中间状态，即不能基于
为此，采用DRM求解增量步的静力解时，本文构造了高斯点应力更新的非线性弹性增量算法——在当前增量步内，材料视作非线性弹性，从而ih时刻的
由ih时刻的
及增量步起点处的εn-1、sn-1更新，即：
表1 Rn的完整kDRM求解流程
Table 1 kDRM algorithm for Rn
步骤 操作1 Rn的初始化：已知Rn-1的解λn-1、an-1、Cn-1；已知Rn的加载(位移、荷载)方式；已知Δλ、fref、aI、c1、c2、c3、mmin的值；时刻序号i=0，时间增量h；由式(16)构造Mm 2 S1的计算…… ……k+1 Sk的计算Sk的初始化：已知Sk-1结束时刻序号ik-1的k1 n a、k1 i-iλ-、n 1)C；Sk的起始时刻序号为ik-1；2) 对于迭代步i ＞ ik-1，依序执行：1 knia满足式(8)，则重启Rn且转换为加位移求解；② 采用NEIA由ain、an-1、Cn-1算得Cin；③ 由式(11)或式(12)确定inλ；由式(13)算得ia；若i①n n若满足式(6)：取i④nλ、ain、Cin为Rn的解，进入Rn+1；若Rn为加位移，inλ、ain同时满足(7)，进入Rn+1且加荷载求解；若i ＞ ik-1+2且ain、1 n a、2 iλ为初始条件启动Sk+1 k+2 Sk+1的计算i-i-n⑤a满足式(17)，则：a) 判定时间(i-1)h时系统出现动能极大值；b) Sk结束时间、Sk+1起始时间的序号取为ik=i-1；c) 以k a、ik i n C、k n n…… ……
式中，fσ为材料点应力计算函数。按此算法，仅当满足收敛判据时才更新C。NEIA不同于结构动力分析时(包括EQSA)的NIA——在当前增量步内，ih时刻的
由ih时刻的
及(i-1)h时刻的
更新(注：
更新)，即：
按此算法，按时间序列{ih}更新C。
按上述，图1和图2示意了依赖应变历史的、一维非弹性本构模型的四个高斯点(1、2、3、4)，在增量步的DRM求解过程中，采用NEIA、NIA时的本构状态递变历程。显然，有限元离散系统静态时，前者为具备真实应变历史的真实静力平衡状态，后者可能为具备虚假应变历史的虚假静力平衡状态。
图1 非线性弹性增量算法下的本构递变
Fig.1 A illustration of constitutive-state-changing under nonlinear elastic increment-algorithm
图2 非线性增量算法下的本构递变
Fig.2 A illustration of constitutive-state-changing under nonlinear increment-algorithm
3 算例
下面采用结点位移更新的kDRM方法，以及分别采用高斯点应力更新的NEIA、NIA算法对两个有限元模型进行了静力平衡路径计算，计算结果证实了NEIA的优越性。
3.1 算例I
模型I为一个拉压杆组成的平面屋架(图3)，编号1、2、3杆的截面面积A均为1500 mm2，编号4、5、6、7杆的A均为16000 mm2，其余杆的A均为900 mm2。分析时，每根杆子视为一个二结点线单元，其形函数采用一次多项式；单元轴向应变ε=(l- lini)/lini，式中l、lini分别为已变形的、未变形的轴向长度；由ε 算得单元轴向应力σ 后，可通过(结点的)直接平衡法算得p{e}。
模型I使用《混凝土结构设计规范》(GB 50010―2010)C.1节定义的依赖应变历史的、应变硬/软化的、一维弹塑性本构模型σ =σ(ε)，其为奇函数；当ε ≥ 0时，σ表达式为：
式中：E为弹性模量；εy为屈服应变；fy为屈服应力；k为硬/软化段斜率(当k＞0时，为应变硬化；当k=0时，为理想弹塑性，当k＜0时，为应变软化)；εuy为硬/软化段起点应变；εu为峰值应变。各拉压杆的本构模型参数取值见表2。应力-应变关系的示意图见图3。
图3 模型I的几何信息、离散和本构关系示意
Fig.3 The model I: geometry，meshing and constitutive relation illustration
表2 模型I 杆的本构关系参数及截面面积
Table 2 The model I: parameters of constitutive model and cross-section area of bar
杆号 E/(N•mm-2) k/(N•mm-2) εy εuy εu A/mm2 1～3 200000 2048 0.0019 0.0078 0.0751500 4, 7 18000 1000 0.0016 0.0032 0.01516000 5, 6 18000 -1500 0.0019 0.0035 0.01516000 8～11 200000 2048 0.0019 0.0078 0.075900
模型I的加载结点编号、加载自由度见图3。分析时：mmin=0.3 kg，c1=2×10-4，c3=4.0，h=0.1 s，加荷载时Δλ=1；外部结点荷载参照向量fref、c2、加位移控制时加位移自由度上施加的固定位移增量ΔaI的配置见表3。计算中，考虑了几何非线性，但是未考虑欧拉失稳；当任一单元内出现应变大于εu情形时，终止增量步分析。
3.2 算例II
模型II为一个悬臂梁(图4)，厚度为10 mm，自由端承受集中力。此悬臂梁采用形函数为二次多项式的四结点面单元进行离散，采用4个高斯点近似计算p{e}，详见图4。
表3 模型I、II、III的分析参数
Table 3 Analysis parameters for model I, II and III
注：1) 下标A、B、C、D表示采用NEIA分析，下标a、b、c表示采用NIA分析；2) IIID采用3.3节中论述的解法2求解。
模型号 外部结点荷载参照向量fref /kN 位移增量ΔaI /mmc2/mm IA、Ia f1y=6，其余为0 Δa1y=3 1000 IB、Ib f1y=30，其余为0 Δa1y=10 1000 IC、Ic f1y=50，其余为0 Δa1y=20 1000 ID f1y=60，其余为0 Δa1y=30 1000 IIA、IIaf1y=4，其余为0 Δa1y=2 360 IIB、IIbf1y=12，其余为0 Δa1y=6 360 IIC、IIcf1y=20，其余为0 Δa1y=10 360 IIIA、IIIaf2y= f3y= f4y=4，f1y= f5y=2，其余为0 Δa3y=0.013 IIIB、IIIbf2y= f3y= f4y=12，f1y= f5y=6，其余为0 Δa3y=0.033 IIIC、IIIcf2y= f3y= f4y=32，f1y= f5y=16，其余为0 Δa3y=0.063 IIID f2y= f3y= f4y=4，f1y= f5y=2，其余为0 - 3
模型II采用描述理想塑性材料的Prandtl-Reuss本构模型[31](图4中三个灰色单元除外，原因见本节下文)，相应的增量本构关系为：
式中：dσ为应力增量；dε为应变增量；D为材料刚度矩阵。D为：
式中：De、Dp分别为材料弹性、塑性刚度矩阵，具体由弹性模量E、泊松比ν、增量前应力σ0构造，详见文献[32]；
为由试探应力σtry =σ0+Dedε算得的偏应力张量第二不变量，平面应力问题中，
β为各向同性参数。本算例中，材料参数取值为：ν =0.3，E=200000 N·mm-2，β =300 N·mm-2；按此材料参数计算，轴向拉压情形下，材料屈服强度
集中力(位置见图4)作用附近可能存在复杂的高应力，此高应力可能超越给定材料参数的Prandtl-Reuss本构模型描述范围，从而出现分析无法进行的情况；为了避免此情况，本算例中，集中力作用附近的三个单元(图4中的灰色单元)采用各向同性弹性本构模型，构造此模型的ν =0.3，E=200000 N·mm-2。
图4 模型II的几何信息和离散
Fig.4 The model II: geometry and meshing
模型II的加载结点编号、加载自由度见图4。分析时：c1=2×10-4，c3=4.0，h=0.1s，mmin=0.3 kg，加荷载时Δλ=1；fref、c2、ΔaI的配置见表3。当1号结点的y向位移a1y ＞80 mm时，终止增量步分析。
3.3 算例III
模型III是一个素混凝土构件(图5)，厚度为150 mm，承受轴心压力。有限元分析时，采用与模型Ⅱ相同的单元类型、单元形函数、p{e}的计算方式。
图5 模型III的几何信息、离散和本构关系示意
Fig.5 The model III: geometry，meshing and constitutive relation illustration
模型III使用Mazars提出的单标量弹性损伤本构模型[32－33]：
式中：Dini是由弹性模量E、泊松比ν 构造的材料初始刚度矩阵(各向同性)；标量ω 是损伤变量。ω 的计算公式为：
式中：ωmax是材料点所达到过的最大损伤值；α+ω++α-ω-是此材料点的即时损伤值，其中α+、α-是依赖于ε 的系数，ω-、ω+是受压、受拉损伤变量。α+、α-的算式见文献[33]。ω+、ω-的算式为：
式中：Kini、A+、B+、A-、B-均为材料常数，其值通过实验测定；
是值依赖于ε 的等效应变，计算式参见文献[33]。本算例中，材料参数为：ν =0.2，E=25500 N·mm-2； A+ =0.95，A-=1.38，B+=11500，B-=2000，Kini= 9.8×10-5；利用上述参数可给出材料的单轴受压应力应变关系(图5为示意图)，其中峰值强度fc=24.4 N·mm-2，与fc对应的应变εc=1.75×10-3。
模型III的加载结点编号、加载自由度见图5。模型III是一个受轴心压力的混凝土构件，使用两种方法分析：第一种方法(解法1)，使用比例加荷载增量分析完整平衡路径；第二种方法(解法2)，使用比例加荷载增量分析峰值点前平衡路径，而使用比例加等位移增量分析峰值点后平衡路径。比例加等位移增量时，即将相等位移增量施加于加载自由度，此时平衡方程组为ψ(fn, an)=fn-p(an)=0，式中fn是外部结点荷载。求解此方程组时：非加载自由度上的位移增量Δaα为未知量，而荷载增量Δfα=0；加载自由度上的荷载增量位移增量Δaβ为已知量，而Δfβ为未知量；由kDRM求解Δaα、Δfβ的计算流程可见文献[34]。无论是解法1，还是解法2，本文都采用kDRM求解增量步，且kDRM分析参数中的c1、c3、h、mmin、Δλ取值相同，即c1=2×10-4，c3=4.0，h=0.1 s，mmin=0.3 kg，Δλ=1。采用解法1时，kDRM分析时的其余参数配置见表3。采用解法2时，仅对一组参数配置下的模型IIID进行分析，分析时：c2及求解平衡路径硬化段时的fref见表3；求解平衡路径软化段时，对加载自由度施加的结点位移增量满足Δa1y=Δa2y=Δa3y=Δa4y=Δa5y=0.015 mm。当迭代解时，终止计算。由计算后的数据分析发现：解法2获得的解答同样满足表3给出的比例加荷载要求。
3.4 计算结果及分析
1) 模型I：图6为1号结点y向位移-荷载关系曲线(f1y-a1y曲线)，其标识了模型I加载的硬化阶段、屈服阶段(限模型号IA、IB和IC)、软化阶段；图7为IA情况下4号、5号杆的ε -a1y的关系曲线。模型II：图8为f1y-a1y关系曲线，其标识了模型II加载的硬化阶段、屈服阶段；图9为IIA情况下的1号单元、2号单元的1号高斯点的
关系曲线和
关系曲线。模型III：图10为f3y-a3y关系曲线，其标识了模型III加载的硬化阶段、软化阶段；图11为IIIB、IIID情况下1号、2号单元的εy-a3y的关系曲线。对于图10：图形★标识IIIA的最终一个收敛解(图中未标识临近★的IIIB、IIIC的最终一个收敛解)；图形◆标识IIIa、IIIb、IIIc的最终一个收敛解；图形▼标识“增量步迭代终止条件(式(26))满足、但收敛条件(式(6))不满足”的迭代解；图例IIID标示采用解法2和NEIA算得的平衡路径。对于图11：图形●标识最终一个收敛解；图形▲或▼标识“增量步迭代终止条件(式(26))满足、但收敛条件(式(6))不满足”的迭代解。
2) 从图6、图8、图10可知：① 采用NEIA时，不同的增量幅度(IA、IB、IC，IIA、IIB、IIC，IIIA、IIIB、IIIC)给出相同的平衡路径；② 采用NIA时，不同的增量幅度(Ia、Ib、Ic，IIa、IIb、IIc，IIIa、IIIb、IIIc)给出不同的平衡路径，增量幅度越小，算得的平衡路径越接近NEIA算得的平衡路径。
图6 模型I的f1y-a1y曲线
Fig.6 The model I: f1y-a1ycurve
图7 模型I的ε-a1y曲线 (IA)
Fig.7 The model I: ε-a1ycurve (IA)
图8 模型II的f1y-a1y曲线
Fig.8 The model II f1y-a1ycurve
3) 对于模型I，由图7可知(结合表2)：4号杆、5号杆的材料性能阶段与相应的平衡路径阶段(图6)相关，如结构的软化阶段始于5号单元的应变软化，且同时伴随着4号杆的应变卸载，此现象即为形变局部化。
4) 对于模型II，对于x=121 mm处的梁横截面S，假定S应变保持平面，则由材料力学公式得：① 当S仅有最外边缘纤维应力等于fy时(状态e)，f1y=59.7 kN；② 当S的所有纤维应力等于fy时(状态p)，f1y=89.6 kN。依据图8、图9，并结合图4给出的单元1的1号高斯点、单元2的2号高斯点的y坐标，得：① 当a1y≈24 mm时，S刚越过状态e，此时结构刚度开始降低，f1y≈68.0 kN，此与理论值接近；② 当a1y=49 mm时，S临近状态p，此时结构刚度临近零，f1y ≈92.0 kN，此与理论值接近。
图10 模型III的f3y-a3y曲线
Fig.10 The model III: f3y-a3ycurve
图11 模型III的εy-a3y曲线(IIIB、IIID)
Fig.11 The model III: εy-a3ycurve(IIIB、IIID)
4 讨论
1) 当加载结点的增量幅度超过一定范围时，NEIA可能会给出不真实的平衡路径，如算例I的ID的平衡路径不同于IA、IB、IC的平衡路径。从计算所得的应变数据可知：ID的5号杆跨越了应变屈服阶段。此导致f1y-a1y曲线不存在屈服阶段、算得的平衡路径偏离真实平衡路径。由此可知，即使采用NEIA求解平衡路径，也应避免采用过大的增量幅度。
2) 尽管算例Ⅰ采用了虚构的本构关系，算例Ⅱ采用了较为符合实际的钢材本构关系，算例III采用的单标量弹性损伤本构模型不能很好反映混凝土多轴应力状态下的本构关系[35]，但因为DRM所构造的虚拟动力过程均符合热力学第一定律，动能的耗散必将使得动力解答逼近具有真实应变历史的静力解答。所以，可以认为在应用上，本文方法不受有限单元类型和材料本构模型特征的限制，可进一步研究此方法在由实验标定的、更能真实反映材料复杂受力情况的本构模型的常见结构静力分析中的应用。
3) 在时间域进行结构受力问题分析时，问题具有适定性[36]。可尝试应用DRM、NEIA对依赖应变历史的、应变软化的、变形局部化的钢筋混凝土结构进行静力分析。
4) 若仅从求解静力增量解的迭代次数衡量，DRM算法远大于ISA算法。但DRM的优越性在于：其一，它规避了ISA求解静力解时的总体刚度矩阵K的组装和求逆，此运算涉及较多计算成本，且此计算成本与计算模型尺寸成指数比例关系；其二，它仅涉及向量运算，计算成本与计算模型尺寸成线性比例关系。此外，也可构造一些降低DRM迭代次数的技巧，如线性外推策略。
5 结论
本文提出采用非线性弹性增量算法(NEIA)，更新SHDM结构(依赖应变历史材料制成的结构)非线性有限元静力增量步动力松弛法(DRM)迭代求解时的高斯点应力，以此来规避迭代解中的虚假应变历史。NEIA不同于非线性增量算法(NIA)：使用NEIA时，一个增量步内高斯点的应变加卸载路径被施加该增量步增量前的高斯点本构状态所固定；而使用NIA时，高斯点的应变加卸载路径随增量步内虚拟时间而变化。本文从总体上的增量步解、增量步中的迭代解、迭代步中的应力更新步骤三个层次阐述了结构非线性有限元静力问题的DRM解法。最后，分别采用NEIA、NIA进行了三个算例的DRM求解。据此，本文的结论如下：
(1) 采用基于虚拟动力过程构造的DRM求解增量步时，原静力问题适定可解，且此问题的解等于上述虚拟动力过程的稳态解。采用NIA时，此稳态解为包含虚假应变历史效应的近似解；而采用NEIA时，此稳态解为消除虚假应变历史效应意义上的正确解，从而NEIA能给出具有真实应变历史的结构静力平衡路径。
(2) 在保持本构关系正确的前提下，NEIA给出的平衡路径不依赖于增量大小；而NIA给出的平衡路径依赖于增量大小，当增量趋小时，平衡路径趋于NEIA给出的平衡路径。
(3) 增量步的DRM解法通过动态阻尼、收敛发散判据、荷载位移增量转换规则进行静力平衡路径硬化段、极值点、软化段的全过程分析。
(4) 增量步的DRM解法无需隐式静力算法的整体刚度矩阵K的组装、求逆。
本文方法解决了SHDM结构非线性有限元全过程静力计算中的虚假应变历史问题。
[1]Zienkiewicz O C, Taylor R L.The finite element method:for solid and structural mechanics[M].6th edition.Beijing: Beijing World Publishing Corporation, 2009:22―84.
[2]Dassault Systemes Simulia Corporation.Abaqus 6.14 analysis user’s guide: Volume II[M].Providence:Dassault Systemes Simulia Corporation, 2014: 6.1―6.3.
[3]Day A S.An introduction to dynamic relaxation[J].The Engineer, 1965, 219: 218―221.
[4]苏小卒.预应力混凝土框架抗震性能研究[M].上海:上海科学技术出版社, 1998: 58―79.Su Xiaozu.Seismic behavior of prestressed concrete frame[M].Shanghai: Shanghai Scientific ＆ Technical Publishers, 1998: 58―79.(in Chinese)
[5]Chung W J, Cho J W, Belytschko T.On the dynamic effects of explicit FEM in sheet metal forming analysis[J].Engineering Computations, 1998, 15(6): 750―776.
[6]Dassault Systemes Simulia Corporation.Abaqus 6.14 theory guide[M].Providence: Dassault Systemes Simulia Corporation, 2014: 2.2, 2.4.5―1.
[7]卡德斯图赛 H, 诺里 D H, 布雷齐 F, 等.有限元法手册[M].诸德超, 傅子智, 龚尧南, 等译.北京: 科学出版社, 1996: 1290―1294.Kardestuncer H, Norrie D H, Brezzi F, et al.Finite element handbook[M].Translated by Zhu Dechao, Fu Zizhi, Gong Yaonan, et al.Beijing: Science Press, 1996:1290―1294.(in Chinese)
[8]Borst R D, Crisfield M A, Remmers J J C, et al.Non-linear finite element analysis of solids and structures[M].Chester: John Wiley ＆ Sons, 2012.
[9]何政, 欧进萍.钢筋混凝土结构非线性分析[M].哈尔滨: 哈尔滨工业大学出版社, 2007: 184―188.He Zheng, Ou Jinping.No-linear finite element analysis of reinforced concrete structures[M].Harbin: Harbin Institute of Technology Press, 2007: 184―188.(in Chinese)
[10]段云岭.非线性方程组的解法: 局部弧长法[J].力学学报, 1997, 29(1): 116―121.Duan Yunling.Local arc-length method-A solution procedure for non-linear finite element method[J].Chinese Journal of Theoretical and Applied Mechanics,1997, 29(1): 116―121.(in Chinese)
[11]May I M, Duan Y.A local arc-length procedure for strain softening[J].Computers ＆ Structures, 1997, 64(1-4):297―303.
[12]Yang Z, Chen J.Fully automatic modelling of cohesive discrete crack propagation in concrete beams using local arc-length methods[J].International Journal of Solids ＆Structures, 2004, 41(3/4): 801―826.
[13]Alfano G, Crisfield M A.Solution strategies for the delamination analysis based on a combination of localcontrol arc-length and line searches[J].International Journal for Numerical Methods in Engineering, 2003,58(7): 999―1048.
[14]Natario P, Silvestre N, Camotim D.Web crippling failure using quasi-static FE models[J].Thin-Walled Structures,2014, 84: 34―49.
[15]Cundall P A.Explicit finite-difference methods in geomechanics[J].Numerical Methods in Geomechanics,1976.
[16]王伟, 苏小卒.动力松弛法在应变软化类结构有限元静力分析中的应用[J].计算力学学报, 2018, 35(2):230―237.Wang Wei, Su Xiaozu.Application of dynamic relaxation method in finite element static-solution of strainsoftening-type structure[J].Chinese Journal of Computational Mechanics, 2018, 35(2): 230―237.(in Chinese)
[17]Rezaiee-pajand M, Sarafrazi S R, Rezaiee H.Efficiency of dynamic relaxation methods in nonlinear analysis of truss and frame structures[J].Computers ＆ Structures,2012, s112/113(s112/113): 295―310.
[18]Otter J R H, Cassell A C, Hobbs R E.Dynamic relaxation[J].ICE Proceedings, 1966, 35(4): 633―656.
[19]Rushton K R.Large deflexion of variable-thickness plates[J].International Journal of Mechanical Sciences,1968, 10(9): 723―735.
[20]Frieze P A, Hobbs R E, Dowling P J.Application of dynamic relaxation to the large deflection elasto-plastic analysis of plates[J].Computers ＆ Structures, 1978,8(2): 301―310.
[21]Ramesh G, Krishnamoorthy C S.Geometrically noninear analysis of plates and shallow shells by dynamic relaxation[J].Computer Methods in Applied Mechanics and Engineering, 1995, 123(1): 15―32.
[22]叶继红, 李爱群, 刘先明.动力松弛法在索网结构形状确定中的应用[J].土木工程学报, 2002, 35(6): 14―19.Ye Jihong, Li Aiqun, Liu Xianming.Form finding of cable nets by modified dynamic relaxation[J].China Civil Engineering Journal, 2002, 35(6): 14―19.(in Chinese)
[23]阚远, 叶继红.动力松弛法在索穹顶结构形状确定中的应用[J].工程力学, 2007, 24(9): 50―55.Kan Yuan, Ye Jihong.Form finding of cable domes by modified dynamic relaxation[J].Engineering Mechanics,2007, 24(9): 50―55.(in Chinese)
[24]张志宏, 李志强, 董石麟.杂交空间结构形状确定问题的探讨[J].工程力学, 2010, 27(11): 56―63.Zhang Zhihong, Li Zhihong, Dong Shilin.Discussions on shape determination of hybrid spatial structures[J].Engineering Mechanics, 2010, 27(11): 56―63.
[25]Ali N B H, Rhode-barbarigos L, Smith I F C.Analysis of clustered tensegrity structures using a modified dynamic relaxation algorithm[J].International Journal of Solids ＆Structures, 2011, 48(5): 637―647.
[26]Yu R C, Saucedo L, Ruiz G.Finite-element study of the diagonal-tension failure in reinforced concrete beams[J].International Journal of Fracture, 2011, 169(2): 169―182.
[27]Garcia J R, Rio G, Cadou J M.Numerical study of dynamic relaxation methods: Contribution to the modeling of inflatable lifejackets[M].Saarbrücken:LAP LAMBERT Academic Publishing Gmbh ＆ Co.KG, 2012.
[28]Baverel O, Nabaei S S, Weinand Y.Mechanical form-finding of the timber fabric structures with dynamic relaxation method[J].International Journal of Space Structures, 2013, 28(3/4): 197―214.
[29]Oliver J, Huespe A E, Cante J C.An implicit/explicit integration scheme to increase computability of non-linear material and contact/friction problems[J].Computer Methods in Applied Mechanics ＆ Engineering,2008, 197(21): 1865―1889.
[30]Zienkiewicz O C, Taylor R L, Zhu J Z.The finite element method: Its basis and fundamentals[M].6th ed.Beijing:Beijing World Publishing Corporation, 2008: 615―618.
[31]陈惠发, 萨里 A F.弹性与塑性力学[M].余天庆, 王勋文, 刘再华, 译.北京: 中国建筑工业出版社, 2004:179―188, 225―234.Chen W F, Saleeb A F.Elasticity and plasticity[M].Translated by Yu Tianqing, Wang Xunwen, Liu Zaihua.Beijing: China Architecture and Building Press, 2004:179―188, 225―234.(in Chinese)
[32]Mazars J.A description of micro- and macroscale damage of concrete structures[J].Engineering Fracture Mechanics, 1986, 25(5/6): 729―737.
[33]Mazars J, Pijaudier J.Continuum damage theoryapplication to concrete[J].Journal of Engineering Mechanics, 1989, 115(2): 345―365.
[34]彭芳乐, 李福林, 白晓宇.针对砂土应变软化强非线性问题的动态松弛有限元法研究[J].岩土力学, 2012,33(2): 590―596.Peng Fangle, Li Fulin, Bai Xiaoyun.A dynamic relaxation-finite element method for strong nonlinearity caused by post-peak strain softening of sands[J].Rock ＆ Soil Mechanics, 2012, 33(2): 590―596.(in Chinese)
[35]李杰, 任晓丹.混凝土静力与动力损伤本构模型研究进展述评[J].力学进展, 2010, 40(3): 284―297.Li Jie, Ren Xiaodan.A review on the constitutive model for static and dynamic damage of concrete[J].Advances in Mechanics, 2010, 40(3): 284―297.(in Chinese)
[36]维诺格拉多夫 I M.数学百科全书: 第3卷[M].《数学百科全书》编译委员会, 译.北京: 科学出版社,1997: 6―9.Vinogradov I M.Encyclopaedia of mathematics: Vol.3[M].Translated by Mathematics encyclopedia compilation committee.Beijing: Science Press, 1997: 6―9.(in Chinese)