混凝土结构疲劳损伤全过程模拟的加速算法

摘要

基于混凝土随机疲劳损伤本构模型改进了循环跳跃加速算法。将外推变量从损伤变量修改为累积耗能变量，研究了跳跃阈值和外推方式对加速算法计算结果的影响。当选取合适的跳跃阈值时，可以利用较少的逐步计算结果获得更为精确的外推变量信息，使得混凝土疲劳损伤全过程模拟具有更高的精度和更小的计算成本。最后，采用改进加速算法对受压素混凝土梁进行疲劳损伤全过程模拟。结果表明，改进加速算法可以准确高效地模拟混凝土结构的疲劳破坏全过程。

关键词

混凝土; 损伤本构模型; 疲劳损伤; 加速算法

混凝土疲劳破坏的本质是材料内部微裂纹在疲劳荷载作用下不断扩展、聚合，最终形成宏观不稳定裂纹的过程。连续介质损伤力学理论^［1-4］引入损伤内变量来描述材料性能的劣化，为混凝土结构疲劳全过程分析提供了基础。将混凝土疲劳损伤本构模型与有限元分析方法相结合，能够再现混凝土结构在疲劳荷载作用下的非线性发展过程。混凝土疲劳损伤的演化过程往往达数年时间，因此采取逐步循环计算方式实现结构全寿命疲劳损伤演化分析是不可取的。目前疲劳加速算法主要分为2类：时域多尺度算法^［5-7］和循环跳跃算法^［8-17］。前者借鉴均匀化理论中空间材料均匀化的思想，对疲劳损伤的演化进行时间尺度的分解，通过耦合求解时间域上的双尺度控制方程实现疲劳计算的加速。后者以有限次循环加载计算结果合理外推损伤变量的方式达到加速效果。时间双尺度算法需要在不同时间尺度上传递响应信息，所以本构关系不同时需要重新推导控制方程。循环加速算法虽然简单实用，但是跳跃阈值和外推方式决定了加速解是否能收敛到逐步计算解，如果选取不当，就有可能与逐步计算解发生较大偏离。本文将循环跳跃加速算法作为研究对象。目前，已有众多学者针对不同的疲劳损伤本构关系开发了相应的循环跳跃算法。Kiewel等^［9］基于一类黏塑性模型发展了循环跳跃加速算法，并考察了多种外推函数包括多项式和样条函数的影响。Sun等^［10-11］发展了自适应区块跳跃算法，将许多循环加载打包为一个区块，当加载区块所包含的周期数合适时可以达到较好的加速效果。针对计算中各积分点演化速率不同的特点，Van Paepegem等^［12-13］研究了结构计算中跳跃阈值选取问题。Sally等^［14］将循环跳跃算法推广至准随机荷载。

从循环跳跃加速算法的基本思想可以发现，外推变量和外推方式的选取会直接影响计算结果的准确性。目前疲劳损伤计算中采用的外推变量大多是损伤变量^［15-19］。当损伤演化曲线较为光滑时，该方法配合一定的精度控制策略可以获得较好的精度和较快的计算效率。当损伤变量的演化呈现阶梯状等不光滑增长形式时，少量循环计算往往无法获得准确的损伤增长信息，导致计算精度差甚至出现错误，而通过大量循环计算获取当前状态的损伤增长信息，又与加速算法节约计算成本的初衷相悖。

由微‒细观随机断裂模型发展而来的随机疲劳损伤本构模型在数值实现中不可避免地会遇到疲劳损伤增长不光滑的问题。为此，本文尝试将外推变量从损伤变量改为累积耗能变量。结果表明，改进加速算法利用更少的计算量就可得到更精确的外推结果，从而提高疲劳损伤全过程模拟的计算精度和效率。

1 混凝土疲劳损伤本构关系

基于随机介质和随机损伤的观点，李杰等^［1，20-21］发展了混凝土损伤分析的微‒细观随机断裂模型。该模型采用微弹簧的随机断裂表征混凝土内部由微裂纹随机萌生、扩展导致的随机损伤演化。以单轴受拉为例，将混凝土代表性体积单元（RVE）抽象为一系列由两端刚性板约束的并联弹簧，如图1所示。图中，$\sigma$、$\epsilon$分别为单轴应力和应变，$N_{\mathrm{s}}$为微弹簧数量。假定微弹簧的断裂应变为一随机场，根据随机介质假定^［22-23］，随机损伤的演化可表示为

$$D={\int }_{\mathrm{0}}^{\mathrm{1}}H(\epsilon -\Delta (x))\mathrm{d}x$$

（1）

式中：$\Delta (x)$为微弹簧的断裂应变随机场，一维概率密度函数服从对数正态分布^［1，21］；$x$为微弹簧的空间坐标；$H(·)$为Heaviside函数。令$Z(x)=\mathrm{l}\mathrm{n}\Delta (x)$，$Z(x)$的均值和标准差分别为$\lambda$和$\zeta$。

图1  受拉微‒细观随机断裂模型

Fig.1  Tensile micro-meso stochastic fracture model

考虑疲劳荷载作用下疲劳损伤的累积效应，Ding等^［24］假定每个微弹簧存在内部结构并将微弹簧视为一个微观能量耗散单元。当微弹簧的累积耗能$E_{\mathrm{f}}$超过弹性应变能$E_{\mathrm{s}}$时，微观单元破坏，则可将式（1）改写为

$$D={\int }_{\mathrm{0}}^{\mathrm{1}}H\left(E_{\mathrm{f}}-E_{\mathrm{s}}\right)\mathrm{d}x,E_{\mathrm{s}}=\frac{\mathrm{1}}{\mathrm{2}}{E}_{\mathrm{0}}{\Delta }^{\mathrm{2}}(x)$$

（2）

式中，$E_{\mathrm{0}}$为混凝土弹性模量。

Ding等^［24］基于速率过程理论^［25］和裂纹分形假设^［26］分析了微弹簧的纳微观耗能过程，并建立了$E_{\mathrm{f}}$的表达式：

$$E_{\mathrm{f}}=$$

$${\int }_{\mathrm{0}}^t{C}_{\mathrm{0}}\mathrm{e}\mathrm{x}\mathrm{p}(-\kappa DY)(Y-\gamma (\vartheta ))\left(\frac{Y^p}{\Gamma \mathrm{e}\mathrm{x}\mathrm{p}(-\beta \vartheta )}\right)\mathrm{d}t$$

（3）

式中：$Y$为损伤能释放率；$\Gamma$、$\gamma$分别为宏观和微观均匀化表面能的代表值；$C_{\mathrm{0}}$、$\kappa$、$\vartheta$、$\beta$、$p$为模型参数。

为了便于工程应用，在$E_{\mathrm{f}}$的建模过程中忽略表面能的影响。考虑裂纹间的相互作用，进一步引入文献［27］中提出的损伤愈合函数和损伤扩展函数修正耗能表达式。$E_{\mathrm{f}}$表达式修正为

$$E_{\mathrm{f}}=$$

$${\int }_{\mathrm{0}}^t{C}_{\mathrm{1}}(\mathrm{e}\mathrm{x}\mathrm{p}(-{\theta }_{\mathrm{1}}D)+\mathrm{e}\mathrm{x}\mathrm{p}(-{\theta }_{\mathrm{2}}(\mathrm{1}-D)))Y^{p+\mathrm{2}}\mathrm{d}t$$

（4）

式中：$C_{\mathrm{1}}$为指前因子；$p+\mathrm{2}$为尺度参数，表示从纳观到宏观的裂纹尺度数；${\theta }_{\mathrm{1}}$、${\theta }_{\mathrm{2}}$分别为控制疲劳损伤第一阶段和第三阶段发展速率的模型参数。

经过系统严密的推导，Wu等^［4］给出的受拉和受压损伤能释放率显式表达式为

$$Y^{+}=\sqrt[]{E_{\mathrm{0}}\left({\overline{\mathit{\sigma }}}^{+}:{\mathit{C}}_{\mathrm{0}}:\overline{\mathit{\sigma }}\right)},Y^{-}=\alpha {\overline{I}}_{\mathrm{1}}+\sqrt[]{\mathrm{3}{\overline{J}}_{\mathrm{2}}}$$

（5）

式中：${\mathit{C}}_{\mathrm{0}}$为四阶柔度张量；$\overline{\mathit{\sigma }}$和${\overline{\mathit{\sigma }}}^{+}$分别为有效应力张量及其正分量；${\overline{I}}_{\mathrm{1}}$、${\overline{J}}_{\mathrm{2}}$分别为有效应力张量的第一和第二不变量；$\alpha$为与双轴强度提高有关参数，通常取0.121 2。

结合弹塑性损伤力学框架^［1］，混凝土随机疲劳损伤本构关系为

$$\mathit{\sigma }=\left(\mathrm{1}-D\right):{\mathit{E}}_{\mathrm{0}}:(\mathit{\epsilon }-{\mathit{\epsilon }}^{\mathrm{p}})$$

（6）

式中：$\mathit{\sigma }$为Cauchy应力张量；${\mathit{E}}_{\mathrm{0}}$为四阶刚度张量；$\mathit{\epsilon }$、${\mathit{\epsilon }}^{\mathrm{p}}$分别为应变张量和塑性应变张量。

值得指出的是，在单调加载中常用的损伤准则为Kuhn-Tucker条件，该条件规定，当损伤能释放率超过历史最大值$r$时引发损伤：

$$\dot{d}=\left\{\begin{array}{l}G(Y),Y>r\\ \mathrm{0},Y\le r\end{array}\right.$$

（7）

这一条件无法反映损伤面内的损伤。因此，在疲劳损伤分析中，应将Kuhn-Tucker条件替换为疲劳加卸载不可逆条件^［28］：

$$\dot{d}=\left\{\begin{array}{l}G(Y),\dot{Y}>\mathrm{0}\\ \mathrm{0},\mathrm{其他}\end{array}\right.$$

（8）

2 循环跳跃加速算法

将损伤变量作为外推变量进行疲劳分析。根据循环跳跃算法的基本思想^［15，18-19］，先对少量循环加载进行精细有限元计算，得出$t$时刻损伤变量$d_t$的变化趋势。根据设定的跳跃阈值$\delta d_{\mathrm{t}\mathrm{h}}$，将损伤外推至$(t+\mathrm{\Delta }t)$时刻得到的损伤$d_{t+\mathrm{\Delta }t}$，然后以此时的损伤状态作为初始条件代回原模型进行精细化计算。如此重复外推、回代直至单元破坏。基于疲劳损伤本构关系的循环跳跃加速算法基本原理如图2所示，${\sigma }_{\mathrm{m}\mathrm{a}\mathrm{x}}$、${\sigma }_{\mathrm{m}\mathrm{i}\mathrm{n}}$分别表示疲劳加载的最大和最小应力。在$(t+\mathrm{\Delta }t)$时刻的损伤变量$d_{t+\mathrm{\Delta }t}$可用$t$时刻损伤变量$d_t$的泰勒展开表示：

$$d_{t+\mathrm{\Delta }t}=d_t+{\dot{d}}_t\mathrm{\Delta }t+\frac{\mathrm{1}}{\mathrm{2}}{\ddot{d}}_t\mathrm{\Delta }{t}^{\mathrm{2}}+O(\mathrm{\Delta }{t}^{\mathrm{2}})$$

（9）

图2  循环跳跃加速算法原理

Fig.2  Principle of cycle jump acceleration algorithm

考虑到疲劳损伤的演化趋势并不复杂，所以只考虑变化的二阶信息。由式（9）可知，在外推过程中需要给出外推时间$\mathrm{\Delta }t$以及损伤变量的变化速度${\dot{d}}_t$和加速度${\ddot{d}}_t$。由逐步计算结果可知，需要至少计算1个完整循环以获得${\dot{d}}_t$，至少计算2个完整循环以获得${\ddot{d}}_t$。

考虑2个加载循环，令$\delta t$为一个加载周期，可得

$${\dot{d}}_t=\frac{d_{t_{\mathrm{3}}}-d_{t_{\mathrm{1}}}}{t_{\mathrm{3}}-t_{\mathrm{1}}}=\frac{\delta d_{\mathrm{13}}}{\mathrm{2}\delta t}$$

（10）

$${\ddot{d}}_t=\frac{\delta d_{\mathrm{23}}-\delta d_{\mathrm{12}}}{\delta t}$$

（11）

外推时间$\mathrm{\Delta }t$可以根据设定的跳跃阈值$\delta d_{\mathrm{t}\mathrm{h}}$确定，一般考虑线性外推，则

$$\mathrm{\Delta }t=\frac{\delta d_{\mathrm{t}\mathrm{h}}}{{\dot{d}}_t}$$

（12）

若考虑二阶外推，则有

$$\frac{\mathrm{1}}{\mathrm{2}}{\ddot{d}}_t\mathrm{\Delta }{t}^{\mathrm{2}}+{\dot{d}}_t\mathrm{\Delta }t-\delta d_{\mathrm{t}\mathrm{h}}=\mathrm{0}$$

（13）

此时，$\mathrm{\Delta }t$的解与判别式$({\dot{d}}_t^{\mathrm{2}}+\mathrm{2}{\ddot{d}}_t\delta d_{\mathrm{t}\mathrm{h}})$的正负有关，需要分类讨论。

${\ddot{d}}_t<\mathrm{0}$对应损伤发展的第一阶段，方程存在2个正根，对应图3中的2个解。由二次函数的性质可知，较小的正根是合理的外推结果，故跳跃步长取为

$$\mathrm{\Delta }t=\frac{-{\dot{d}}_t+\sqrt[]{{\dot{d}}_t^{\mathrm{2}}+\mathrm{2}{\ddot{d}}_t\delta d_{\mathrm{t}\mathrm{h}}}}{{\ddot{d}}_t},{\ddot{d}}_t<\mathrm{0}$$

（14）

图3  外推时间的求解

Fig.3  Solution of extrapolation time

${\ddot{d}}_t\ge \mathrm{0}$对应损伤发展第三阶段，方程存在2个正根，对应图3中的2个解。由二次函数的性质可知，较小的正根是合理的外推结果，故跳跃步长取为

$$\mathrm{\Delta }t=\frac{-{\dot{d}}_t+\sqrt[]{{\dot{d}}_t^{\mathrm{2}}+\mathrm{2}{\ddot{d}}_t\delta d_{\mathrm{t}\mathrm{h}}}}{{\ddot{d}}_t},{\ddot{d}}_t\ge \mathrm{0}$$

（15）

可以发现，式（14）与式（15）在最终表达上一致。

当判别式小于零时，方程无解，对应图3中的目标线2，说明此时跳跃阈值过大，应利用二分法将阈值减小，直至方程有解。

值得说明的是，文献［18-19］中均没有求解式（13），而是通过假定外推前后损伤演化加速度不变给出二阶外推时间。因此，一般需要引入额外的精度控制机制，才能使外推结果与逐步计算结果有较好的吻合度。

研究发现，直接将损伤变量作为外推变量存在如下问题：

（1） 需要较多循环数的精细化计算才能得到较为准确的损伤变量速度和加速度，而少量循环无法得到有效的速度和加速度。

（2） 由于数值分析中的微弹簧样本是离散的，因此无法确定外推后的累积耗能值。累积耗能值有可能是2个邻近微弹簧应变能之间的任意值，更新不当会使整个程序陷入死循环。

（3） 除了损伤变量外，模型中还涉及其他内变量如塑性应变的更新。如果同样采取外推的方式，就有可能出现外推后损伤变量和塑性应变不匹配的问题。

以一个样本的疲劳损伤为例，疲劳损伤演化曲线如图4所示。该样本对应的疲劳损伤速度和加速度演化曲线如图5所示。可以看出，少量循环计算获得的损伤发展速度和加速度严重依赖循环次数的选取，循环次数选取不当，可能导致外推结果与真实值相差巨大。如果精确计算多个循环以获得当前损伤状态的整体发展趋势，计算成本就会上升，达不到预想的加速计算效果。

图4  样本疲劳损伤演化曲线

Fig.4  Fatigue damage evolution curve of the sample

图5  样本疲劳损伤速度和加速度演化曲线

Fig.5  Fatigue damage velocity and acceleration evolution curves of the sample

注意到，累积耗能变量的变化是较为光滑的。对同一样本，计算累积耗能随循环次数变化的速度和加速度，结果如图6所示。可以发现，累积耗能随循环次数的变化较为平滑，更适合作为外推变量。

图6  样本累积耗能速度和加速度演化曲线

Fig.6  Energy dissipation velocity and acceleration evolution curves of the sample

将累积耗能作为控制变量，并进行泰勒展开，有

$$E_{\mathrm{f},t+\mathrm{\Delta }t}=E_{\mathrm{f},t}+{\dot{E}}_{\mathrm{f},t}\mathrm{\Delta }t+\frac{\mathrm{1}}{\mathrm{2}}{\ddot{E}}_{\mathrm{f},t}\mathrm{\Delta }{t}^{\mathrm{2}}+O(\mathrm{\Delta }{t}^{\mathrm{2}})$$

（16）

线性外推时间为

$$\mathrm{\Delta }t=\frac{E_{\mathrm{f},\mathrm{t}\mathrm{h}}}{{\dot{E}}_{\mathrm{f},t}}$$

（17）

二阶外推时间为

$$\mathrm{\Delta }t=\frac{-{\dot{E}}_{\mathrm{f},t}\pm \sqrt[]{{\dot{E}}_{\mathrm{f},t}^{\mathrm{2}}+\mathrm{2}{\ddot{E}}_{\mathrm{f},t}{E}_{\mathrm{f},\mathrm{t}\mathrm{h}}}}{{\ddot{E}}_{\mathrm{f},t}}$$

（18）

式中：${\dot{E}}_{\mathrm{f},t}$和${\ddot{E}}_{\mathrm{f},t}$分别为累积耗能的增长速度和加速度；$E_{\mathrm{f},\mathrm{t}\mathrm{h}}$为累积耗能跳跃阈值。

外推回代过程中塑性应变与损伤变量不匹配的问题可以通过经验塑性模型避免。经验塑性模型的塑性应力增量记为^［29］

$${\dot{\mathit{\sigma }}}^{\mathrm{p}\pm }={\eta }^{\pm }{\dot{\overline{\mathit{\sigma }}}}^{\pm }$$

（19）

式中，${\eta }^{\pm }$为与损伤相关的塑性因子。${\eta }^{\pm }$的表达式为

$${\eta }^{\pm }={\left(\frac{{\xi }^{\pm }{d}^{\pm }}{\mathrm{1}-d^{\pm }}\right)}^{n_{\mathrm{p}}^{\pm }}$$

（20）

式中，${\xi }^{\pm }$、$n_{\mathrm{p}}^{\pm }$为塑性模型参数。

在外推前后将式（19）进行数值离散，可得

$${\mathit{\sigma }}_{t+\mathrm{\Delta }t}^{\mathrm{p}\pm }-{\mathit{\sigma }}_t^{\mathrm{p}\pm }={\eta }_{t+\mathrm{\Delta }t}^{\pm }({\overline{\mathit{\sigma }}}_{t+\mathrm{\Delta }t}^{\pm }-{\overline{\mathit{\sigma }}}_t^{\pm })$$

（21）

根据每个循环疲劳荷载最大Cauchy应力相等的条件，可以求得外推后的塑性应力为

$${\mathit{\sigma }}_{t+\mathrm{\Delta }t}^{\mathrm{p}\pm }={\mathit{\sigma }}_t^{\mathrm{p}\pm }+{\eta }_{t+\mathrm{\Delta }t}^{\pm }\left(\frac{\mathrm{1}-d_t}{\mathrm{1}-d_{t+\mathrm{\Delta }t}}-\mathrm{1}\right){\mathit{\sigma }}_t^{\pm }$$

（22）

塑性应变与塑性应力的转换关系为

$${\mathit{\epsilon }}^{\mathrm{p}}={\mathit{C}}_{\mathrm{0}}:{\mathit{\sigma }}^{\mathrm{p}}$$

（23）

在实际计算中，需要对结构或构件进行网格划分，并为网格确定单元类型，不同的单元内部会有一个或多个积分点。一般情况下，不同积分点对应的状态变量是不同的。以Abaqus中四边形平面应力单元CPS4为例，该单元内部有4个高斯积分点，计算时需要对单元每个积分点提取计算结果。通常，不同积分点对应的跳跃步长不一致，跳跃步长可以取为所有积分点跳跃步长的最小值。Van Paepegem等^［12-13］建议取为所有跳跃步长的某一个分位值。

3 实例分析

3.1　混凝土单轴受压

以一个混凝土试块的单轴受压疲劳为例，对单元积分点处的代表性体积单元取微弹簧数为1 000，采用隐式分析，计算单元类型为CPS4。计算中，首先只考虑线性外推，计算2个完整循环，并以2个循环的平均累积耗能速度外推。不同跳跃阈值下第1个积分点的计算结果如图7所示。

图7  不同累积耗能跳跃阈值的线性外推计算结果

Fig.7  Linear extrapolation calculation results of different cumulative energy dissipation jump thresholds

同时采用逐步精细化分析的方式获取分析结果。图7中参考曲线即为逐步计算结果，可以认为该计算结果为精确结果。逐步计算求解耗时约为45 min。线性加速情况下，累积耗能跳跃阈值$E_{\mathrm{t}\mathrm{h}}$为0.976 5、1.953 1、3.906 2、7.812 4 kJ·m^-3时对应的计算时间分别为12、6、4、2 min。第1个累积耗能跳跃阈值对应的损伤跳跃阈值$d_{\mathrm{t}\mathrm{h}}=\mathrm{0.01}$。与预期相同，随着累积耗能跳跃阈值的增大，加速效果明显提升。累积耗能跳跃阈值反映出损伤发展精度逐步下降的趋势，如果累积耗能跳跃阈值过大，相比参考曲线就会出现较大偏差。

考虑二阶效应的外推，同样只计算2个完整循环，累积耗能跳跃阈值$E_{\mathrm{t}\mathrm{h}}$为3.906 2、7.812 4、66.312 0、191.210 0 kJ·m^-3时的计算结果如图8所示，对应的计算耗时分别为25、15、6、4 min。

图8  不同累积耗能跳跃阈值的二阶外推计算结果

Fig.8  Second-order extrapolation calculation results of different cumulative energy dissipation jump thresholds

与线性外推类似，随着累积耗能跳跃阈值的增大，计算效率提升，但随之而来的是精度下降，累积耗能跳跃阈值过大会出现较大误差。与线性外推不同的是，对于相同累积耗能跳跃阈值的情况（见图7c、d和图8a、b），二阶外推的步长远小于线性外推，精度也随之上升，这是由于引入了二阶项修正的结果。对比图7d和图8d可以发现，相比于线性外推，即使二阶外推采用很大的跳跃阈值，二阶修正效应也会给精度带来较大的提升。

由于疲劳损伤第二阶段发展缓慢且数值上微弹簧样本存在最小的损伤分辨率，因此2个完整循环的损伤增量有可能小于微弹簧计算损伤的分辨率，而从宏观来看损伤变量并未增加。这就导致传统基于损伤变量的外推方法在损伤分辨率较大的情况下不能适用，所以本文没有对比2种外推方法的计算效率。通过提高样本中微弹簧数可以提高损伤分辨率，但在提高损伤分辨率的同时带来更大的计算量。

为了将循环跳跃算法的误差定量化，引入循环跳跃算法和逐步计算精细化结果的1范数来定义相对误差。2种外推方式下不同累积耗能跳跃阈值的相对误差（e）见表1，表中的a~d对应图7和图8中4种工况的分析结果。相对误差计算式为

$$e=\frac{{‖d_{\mathrm{s}\mathrm{t}\mathrm{e}\mathrm{p}\mathrm{‑}\mathrm{b}\mathrm{y}\mathrm{‑}\mathrm{s}\mathrm{t}\mathrm{e}\mathrm{p}}-d_{\mathrm{c}\mathrm{y}\mathrm{c}\mathrm{l}\mathrm{e}\mathrm{j}\mathrm{u}\mathrm{m}\mathrm{p}}‖}_{\mathrm{1}}}{{‖d_{\mathrm{s}\mathrm{t}\mathrm{e}\mathrm{p}\mathrm{‑}\mathrm{b}\mathrm{y}\mathrm{‑}\mathrm{s}\mathrm{t}\mathrm{e}\mathrm{p}}‖}_{\mathrm{1}}}$$

（24）

表1  2种外推方式的相对误差

Tab.1  Relative errors of two extrapolation methods

外推方式	不同加载工况相对误差
	a	b	c	d
线性外推	0.034 1	0.052 0	0.077 5	0.138 8
二阶外推	0.020 2	0.022 9	0.050 6	0.107 3

以第1个累积耗能跳跃阈值作为基准阈值，其他累积耗能跳跃阈值除以基准阈值得到正则化跳跃阈值。2种外推方式的相对误差和计算时间随正则化跳跃阈值的变化如图9所示。可以发现，两者构成了一个典型效率和精度的平衡问题。显然，2种外推格式对应的最优累积耗能跳跃阈值不同。

图9  外推计算时间和相对误差曲线

Fig.9  Extrapolation calculation time and relative error curve

3.2　素混凝土梁疲劳损伤全过程模拟

选取文献［23］中的素混凝土三点弯梁进行数值模拟。梁截面尺寸信息如图10所示。对C50强度等级混凝土的均值参数进行模拟，具体参数取值为^［24］：${\lambda }^{+}=\mathrm{5.202}$，${\zeta }^{+}=\mathrm{0.463}$，${\lambda }^{-}=\mathrm{7.49}$，${\zeta }^{-}=\mathrm{0.19}$。计算得到抗压强度$f_{\mathrm{c}}=\mathrm{43.4}$ MPa，抗拉强度$f_{\mathrm{t}}=\mathrm{3.0}$ MPa，弹性模量$E_{\mathrm{0}}={\mathrm{30}}^{}$ kN·mm^-2。

图10  素混凝土梁尺寸 （单位：mm）

Fig.10  Dimension of plain concrete beam (unit: mm)

素混凝土梁在单调加载下峰值荷载$P_{\mathrm{m}\mathrm{a}\mathrm{x}}=\mathrm{19.2}$ kN。疲劳数值模拟中选取了4个加载工况，4个工况对应疲劳荷载最大值分别为0.9、0.8、0.7、0.6倍的峰值荷载，最小疲劳荷载为0.1倍的峰值荷载，加载频率为5 Hz。

基于Abaqus平台建立有限元模型，采用CPS4平面四边形单元，单元总数为216，采用的模型参数如表2所示。拉压累积耗能跳跃阈值分别为 $E_{\mathrm{f},\mathrm{t}\mathrm{h}}^{+}=$25.856 0 J·m^-3，$E_{\mathrm{f},\mathrm{t}\mathrm{h}}^{-}=$6.574 7 J·m^-3，对应的损伤跳跃阈值均为$\mathrm{0.05}$。采用线性外推策略，结构层面跳跃步长取所有积分点跳跃步长的最小值。

表2  模型参数取值

Tab.2  Parameter values of the model

受力状态	$C_1$	$p$	${\theta }_1$	${\theta }_2$	$\xi$	$n_{\mathrm{p}}$
受拉	$2.25\times {10}^{-5}$	5	30	10	0.5	3
受压	$3.75\times {10}^{-18}$	12	30	10	0.5	5

不同疲劳加载工况下计算的疲劳寿命和计算时长如表3所示。可以看出，采用加速算法后，最大疲劳荷载水平为0.8至0.6的3个加载工况的计算时长不会随着疲劳寿命的增加而大幅增加。这是因为跳跃步长会随最大疲劳荷载水平的降低而增大，计算时长仅与加速算法跳跃次数相关。

表3  素混凝土梁疲劳分析结果

Tab.3  Fatigue analysis results of plain concrete beam

计算工况	疲劳寿命	计算时长/s	跳跃次数	最大跳跃步长
$P/P_{\mathrm{m}\mathrm{a}\mathrm{x}}=0.9$	4 430	5 389	95	848
$P/P_{\mathrm{m}\mathrm{a}\mathrm{x}}=0.8$	16 513	11 391	179	3276
$P/P_{\mathrm{m}\mathrm{a}\mathrm{x}}=0.7$	33 532	11 542	189	6 270
$P/P_{\mathrm{m}\mathrm{a}\mathrm{x}}=0.6$	60 809	11 990	220	10 841

图11给出了$P/P_{\mathrm{m}\mathrm{a}\mathrm{x}}=\mathrm{0.8}$工况下受拉损伤（SDV15）随循环次数增加的演化情况。其中，$N/N_{\mathrm{f}}$为疲劳加载次数与疲劳寿命之比。图12对比了不同加载下梁底部跨中位置x方向应变随循环次数的变化。可以发现，本文提出的模型能较好地反映疲劳荷载作用下素混凝土梁裂纹的扩展过程及变形特点。

图11  素混凝土梁受拉损伤云图

Fig.11  Tensile damage contour of plain concrete beam

图12  梁底部跨中x方向应变随循环次数的变化

Fig.12  Variation of strain in x direction at the bottom of beam span with cycles

4 结论

（1）使用循环跳跃加速算法时，外推变量应尽量光滑，才能做到使用少量精细化的计算信息获得准确的外推量信息。

（2）在进行结构层次的分析计算之前，应针对一个单元进行试算，以选择容许误差范围内的最大跳跃阈值。

（3）采用加速算法的计算时间与构件的疲劳寿命没有直接关系，跳跃步长会随疲劳寿命的增大而增大，故有望使用该算法进行超高周疲劳的全过程计算。

作者贡献声明

虢成功：模型提出，数值模拟，论文撰写。

李 杰：提出研究方向，审阅、修改论文。