基于锈损单元畸变控制的结构3D荷载路径拓扑搜寻方法

Topology search method for structural 3D load paths based on distortion control of corrosion-damaged elements

摘要

锈蚀损伤易引起局部单元畸变，导致结构非线性拓扑优化数值失稳。针对现有抑制方法主要用于大变形弹性结构，存在效率低、适用性差等问题，难以满足锈损结构3D荷载路径搜寻的不足，提出了基于锈损单元畸变控制的结构3D荷载路径拓扑搜寻方法。首先，该方法采用局部刚度自适应缩放的思想来控制锈损单元畸变，通过定义的锈损单元畸变度对其局部刚度按不同比例进行缩放，同时抑制单元畸变与优化误差。其次，给出了锈损结构荷载路径搜寻的非线性拓扑优化数学表达式，考虑材料劣化和粘结退化的影响，基于伴随法推导了锈损单元灵敏度计算公式。最后，开展数值算例验证，表明了本文方法能合理地生成锈蚀损伤结构3D荷载路径，揭示了锈蚀RC梁荷载路径发展规律。

Abstract

Corrosion-damage can easily cause local elements distortion， leading to numerical instability in the nonlinear topology optimization of structures. Existing suppression methods are mainly applied to elastic structures undergoing large deformation， but these methods suffer from the issues such as low efficiency and limited applicability， making them difficult to meet the requirements of 3D load paths search for corrosion-damaged structures. Therefore， a method for structural 3D load paths topology search based on distortion control of corrosion-damaged elements is proposed in this paper. Firstly， this method employs the concept of adaptive scaling of local stiffness to control the distortion of corrosion-damaged elements. Their local stiffness is scaled at different proportions based on a defined distortion degree of corrosion-damaged elements， while restraining element distortion and optimization errors. Secondly， a mathematical expression of nonlinear topology optimization of load paths search in corrosion-damaged structures is given. By taking into account the effects of material deterioration and bond degradation， a sensitivity calculation formula for corrosion-damaged elements is derived based on the adjoint method. Finally， numerical examples are conducted for verification， demonstrating that the proposed method in this paper can reasonably generate 3D load paths for corrosion-damaged structures， revealing the development law of load paths in corroded RC beams.

Graphical abstract

关键词

土木工程 / 锈损单元畸变控制 / 数值失稳 / 局部刚度自适应缩放 / 荷载路径 / 拓扑搜寻

Key words

civil engineering / distortion control of corrosion-damaged elements / numerical instability / local stiffness adaptive scaling / load paths / topological search

0 引　言

钢筋锈蚀是造成钢筋混凝土（Reinforce concret， RC）结构破坏的主要因素之一^［1，2］。锈蚀会导致钢筋与混凝土材料力学性能损失、粘结退化和混凝土开裂等损伤，进而影响结构安全性与服役性能^［3-5］。当RC结构遭受上述损伤时，RC结构D区的内力轨迹形式（即荷载路径）会发生改变。锈蚀损伤RC结构荷载路径对于理解结构传力机理和指导结构维修加固具有重要意义。

近年来，已有学者采用线弹性拓扑优化方法得到了RC结构的荷载路径^［6-8］，但该类方法主要侧重于搜寻素混凝土结构弹性荷载路径，忽略了混凝土材料的拉-压不对称响应与埋置钢筋对结构荷载路径的影响^{［7，9，10］}。随后，一些学者通过考虑钢筋与混凝土不同的拉、压力学性能，提出了基于双相材料拓扑优化的荷载路径搜寻方法^{［6，8，11］}。该类方法得到了RC结构抗压混凝土与抗拉钢筋的最优布局，可指导RC结构的初始配筋设计。然而，上述方法忽略了钢筋锈蚀导致的系列损伤的影响，难以适用于锈蚀RC结构实际荷载路径轨迹的搜寻。

钢筋锈蚀会影响RC结构的非线性响应，进而改变RC结构的荷载路径轨迹。首先，锈蚀引起的系列损伤会降低局部区域结构刚度，使得相同外部荷载下结构的变形增大，导致结构应力或应变能分布发生变化，在荷载路径搜寻过程中低效混凝土单元（即应变能较小的单元）会被逐渐删除，进而影响RC结构荷载路径轨迹的演化。此外，拓扑优化方法中最小密度单元刚度通常被人为地降低了10^-9倍^［12-14］，加之锈蚀损伤的影响，使得刚度无穷小的锈损单元易发生畸变，进而导致锈蚀RC结构荷载路径拓扑搜寻过程出现数值失稳问题^{［15，16］}。目前，已有学者给出了一些非线性拓扑优化过程中的数值失稳抑制方法，如单元删除法^［17］、收敛准则松弛法^［18-20］、能量插值法^［21］、单元连通性参数化法^［15］和添加超弹性材料技术^［22］。然而，这些方法主要针对大变形弹性结构的几何非线性拓扑优化数值失稳问题，其数值实现过程较为复杂且求解速率缓慢，难以适用于锈损RC结构的荷载路径搜寻。此外，Yuan等^［23］提出了一种基于刚度自动调节的锈蚀RC结构荷载路径生成方法。然而，该方法主要针对二维锈蚀RC梁的荷载路径生成问题，难以适用于荷载路径更为复杂的三维（Three-dimensions， 3D）锈损RC结构，相关问题还有待进一步研究。

为此，本文提出了一种基于锈损单元畸变控制的结构3D荷载路径拓扑搜寻方法。首先，给出了一种局部区域锈损单元畸变控制方法来克服荷载路径搜寻过程中锈损单元畸变引发的数值失稳问题。其次，通过考虑材料性能损失和粘结退化的影响，基于伴随法推导了锈损单元灵敏度计算公式。最后，通过数值算例说明了本文所提方法的有效性与适用性，并讨论了不同锈蚀水平下RC结构的3D荷载路径发展规律。

1 局部区域锈损单元畸变控制方法

本文提出了一种局部区域锈损单元畸变控制方法来克服荷载路径搜寻过程中的数值失稳问题。该方法采用局部刚度自适应缩放的思想来控制锈损单元畸变，结合锈损单元畸变度对单元刚度进行多比例缩放来抑制数值失稳问题。该方法可控制局部区域单元刚度增长与衰减的速度，从而有效地抑制单元畸变与优化误差。

在双向渐进结构优化（Bi-directional evolutionary structural optimization， BESO）方法中，单元的删除是通过将材料相对密度从1转换到{\tilde{x}}_{\mathrm{m}}来实现的^［24］。各单元刚度矩阵{\mathit{k}}_i与材料相对密度{\tilde{x}}_i的关系为：

{\mathit{k}}_i={\tilde{x}}_i^p{\int }_{v_{\mathrm{c}\mathrm{o}}}{\mathit{B}}^{\mathrm{T}}{\mathit{D}}_{\mathrm{C}}\mathit{B}\mathrm{d}{v}_i={\tilde{x}}_i^p{\mathit{k}}_{\mathrm{c}},{\tilde{x}}_i={\tilde{x}}_{\mathrm{m}}\mathrm{o}\mathrm{r}\mathrm{1}（1）

式中：\mathit{B}为应变-位移矩阵；{\mathit{k}}_{\mathrm{c}}为混凝土单元刚度矩阵；v_{\mathrm{c}\mathrm{o}}为混凝土单元的体积；{\mathit{D}}_{\mathrm{C}}为弹塑性本构矩阵；p为惩罚因子，通常取p=\mathrm{3}^［25］。单元密度{\tilde{x}}_i=x_{\mathrm{m}}表明第i个单元为“空”单元，其单元刚度矩阵为{\tilde{x}}_{\mathrm{m}}^p{\mathit{k}}_{\mathrm{c}}。材料最小密度通常设置为{\tilde{x}}_{\mathrm{m}}=\mathrm{0.001}。因此，单元的刚度被人为地降低了10^-9倍。此外，加之钢筋锈蚀导致的系列损伤的影响，该无穷小刚度的锈损单元会发生畸变，从而导致荷载路径搜寻过程出现数值失稳问题。

一种简单的方法是通过增大局部区域锈损单元的刚度以缓解其畸变问题，如图1所示。在拓扑优化过程中，无穷小刚度的锈损单元对原始结构非线性响应的影响可以忽略不计。然而，过度地增加锈损单元刚度会影响原始结构的非线性响应，从而引起拓扑优化中的误差，影响结构荷载路径的搜寻^［16］。

锈损单元的刚度极小，其本构关系可以表示为：

\left\{\begin{array}{l}{\mathit{k}}_{\mathrm{e},\eta }={\int }_{v_i}{\mathit{B}}^{\mathrm{T}}{\mathit{D}}_{\mathrm{e},\eta }\mathit{B}\mathrm{d}{v}_i,\\ {\mathit{D}}_{\mathrm{e},\eta }=\frac{E_{\mathrm{m},\eta }}{\mathrm{1}-v_{\mathrm{e}}^{\mathrm{2}}}\left[\begin{array}{l}\mathrm{1}{v}_{\mathrm{e}}\mathrm{0}\\ v_{\mathrm{e}}\mathrm{1}\mathrm{0}\\ \mathrm{0}\mathrm{0}\frac{\mathrm{1}-v_{\mathrm{e}}}{\mathrm{2}}\end{array}\right]\end{array}\right.（2）

式中：{\mathit{k}}_{\mathrm{e},\eta }为锈损单元刚度矩阵；{\mathit{D}}_{\mathrm{e},\eta }为本构矩阵；v_i为第i个单元的体积；v_{\mathrm{e}}和E_{\mathrm{m},\eta }分别为泊松比和弹性模量。

由式（2）可知，锈损单元刚度与v_{\mathrm{e}}和E_{\mathrm{m},\eta }有关。因此，在荷载路径搜寻过程中，可通过改变锈损单元弹性模量来调节单元刚度。

针对上述问题，本文给出了一种局部刚度自适应缩放策略。在每次优化迭代中，通过追踪局部区域锈损单元最大畸变度，自动调节该区域锈损单元刚度k_{\mathrm{m},\eta }^k。为避免在荷载路径搜寻过程中E_{\mathrm{m},\eta }值过小导致的数值失稳问题，将初始弹性模量设为一个较大的值E_{\mathrm{m},\eta }^{\mathrm{0}}={\mathrm{10}}^{-\mathrm{4}}{E}_{\mathrm{c}}。在第k+\mathrm{1}优化迭代步中，当锈损单元的最大应变{\epsilon }_{\mathrm{m}\mathrm{a}\mathrm{x}}^k≥应变阈值{\epsilon }^{\mathrm{*}}时，本文所提的局部刚度自适应缩放策略可表示为：

（3）

当{\epsilon }_{\mathrm{m}\mathrm{a}\mathrm{x}}^k\le{\epsilon }^{\mathrm{*}}时，局部刚度自适应缩放策略可表示为：

E_{\mathrm{m},\eta }^{k+\mathrm{1}}=\mathrm{m}\mathrm{a}\mathrm{x}(f(x)+\mathrm{1}-f(\mathrm{1}),{\alpha }_{\mathrm{4}})\cdot E_{\mathrm{m},\eta }^k（4）

式中：E_{m,\eta }^k和E_{m,\eta }^{k+\mathrm{1}}分别为第k步和第k+1步锈损单元的弹性模量；{\alpha }_{\mathrm{1}}、{\alpha }_{\mathrm{2}}、{\alpha }_{\mathrm{3}}和{\alpha }_{\mathrm{4}}为控制单元刚度增长速率的系数^［25］。初始应变阈值{\epsilon }^{\mathrm{*}}等于原数值模型中混凝土单元的最大应变，在每个优化步骤中将应变阈值{\epsilon }^{\mathrm{*}}更新为{\epsilon }_{\mathrm{m}\mathrm{a}\mathrm{x}}^k=({\epsilon }_{\mathrm{m}\mathrm{a}\mathrm{x}}^k+{\epsilon }^{\mathrm{*}})/\mathrm{2}。f(x)为sigmoid函数，可以表示为：

f(x)=\frac{\mathrm{1}}{\mathrm{1}+{\mathrm{e}}^{-x}},x>\mathrm{0}（5）

式中：x={\epsilon }_{max}^k/{\epsilon }^{*}为锈损单元的畸变程度。

采用sigmoid函数可以逐步地调整弹性模量E_{\mathrm{m},\eta }^k。

在结构荷载路径搜寻过程中，通过对比锈损单元最大应变与应变阈值的大小，基于式（3）和式（4），对局部区域锈损单元刚度进行调节。若{\epsilon }_{\mathrm{m}\mathrm{a}\mathrm{x}}^k>{\epsilon }^{\mathrm{*}}，则增大下一个迭代步中局部区域单元刚度{\mathit{k}}_{\mathrm{m},\eta }^k，克服单元畸变导致的数值失稳问题。若{\epsilon }_{\mathrm{m}\mathrm{a}\mathrm{x}}^k\le {\epsilon }^{\mathrm{*}}，则减小下一个迭代步中局部区域单元刚度{\mathit{k}}_{\mathrm{m},\eta }^k，抑制因单元刚度增加而引起的优化误差。

2 锈蚀RC结构荷载路径拓扑搜寻优化公式与灵敏度分析

锈蚀导致的材料性能损失和粘结退化会影响结构位移 U，进而改变锈蚀RC结构的荷载路径轨迹。本文给出了锈蚀RC结构荷载路径搜寻问题的非线性拓扑优化公式，基于伴随变量法，推导了锈损单元灵敏度计算公式，展示了锈蚀RC结构荷载路径搜寻数值实现方案。

2.1　材料性能损失与粘结退化模拟

锈蚀RC结构荷载路径与钢筋和混凝土材料非线性行为以及两者间的粘结行为有关。混凝土材料力学行为采用混凝土塑性损伤（Concrete damage plasticity， CDP）模型进行模拟^［19］。锈蚀会导致钢筋与混凝土的力学性能以及两者间的粘结性能退化，也会降低钢筋与混凝土的抗压强度。

钢筋锈蚀产物膨胀也会影响结构非线性响应，但荷载路径搜寻往往需要进行几十至数百次迭代和非线性有限元分析，直接模拟钢筋锈蚀产物膨胀会大大增加荷载路径搜寻的复杂程度与计算成本。因此，本文通过考虑钢筋锈蚀产物膨胀对混凝土材料性能的影响来简化RC梁数值模型。

此外，结构荷载路径也受混凝土裂缝的影响，混凝土开裂会将RC梁分割为多个混凝土柱，结构通过这些混凝土柱与钢筋来传递内部应力，荷载路径搜寻过程中会逐渐删除低效的混凝土单元，使得一些低效的混凝土柱逐渐消失，优化迭代过程中混凝土会进行应力重分布，最终保留最有效的混凝土柱（即荷载路径），数值模型中混凝土开裂行为由CDP模型表征。

假设纵向钢筋受均匀锈蚀，锈蚀混凝土抗压强度f_{\mathrm{c}\mathrm{m}}^{\mathrm{\text{'}}}可表示为^［26-28］：

f_{\mathrm{c}\mathrm{m}}^{\mathrm{\text{'}}}=\frac{f_{\mathrm{c}\mathrm{m}}}{\mathrm{1}+k_{\mathrm{r}}\cdot {\epsilon }_{\mathrm{1}}/{\epsilon }_{\mathrm{c}\mathrm{u}}}（6）

式中：f_{\mathrm{c}\mathrm{m}}为未锈蚀混凝土抗压强度；k_{\mathrm{r}}=\mathrm{0.1}为经验系数；{\epsilon }_{\mathrm{c}\mathrm{u}}为最大压缩应变；{\epsilon }_{\mathrm{1}}为开裂混凝土平均受拉应变。

{\epsilon }_{\mathrm{1}}的表达式为^{［28，29］}：

{\epsilon }_{\mathrm{1}}=n_{\mathrm{b}\mathrm{s}}\cdot w_{\mathrm{c}\mathrm{r}}/b_{\mathrm{w}}（7）

式中：b_{\mathrm{w}}为未锈蚀RC梁初始宽度；n_{\mathrm{b}\mathrm{s}}为RC梁受压区钢筋数量；w_{\mathrm{c}\mathrm{r}}为锈蚀RC梁总裂缝宽度。

w_{\mathrm{c}\mathrm{r}}的表达式为^{［27，30］}：

w_{\mathrm{c}\mathrm{r}}=\mathrm{2}\mathrm{\pi }X(u_{\mathrm{r}\mathrm{s}}-\mathrm{1})（8）

式中：u_{\mathrm{r}\mathrm{s}}=\mathrm{2}，为体积膨胀率；X为钢筋锈蚀深度。

X的表达式为^［30］：

X=(\mathrm{1}-\sqrt[]{\mathrm{1}-{\eta }_{\mathrm{l}}})r_{\mathrm{1}}（9）

式中：r_{\mathrm{1}}为钢筋半径；{\eta }_{\mathrm{l}}为钢筋锈蚀率。

{\eta }_{\mathrm{l}}的表达式为^［31］：

{\eta }_{\mathrm{l}}=\frac{\mathrm{\Delta }{m}_{\mathrm{l},\eta }}{m_{\mathrm{l}}}=\frac{\mathrm{\Delta }{A}_{\mathrm{l},\eta }}{A_{\mathrm{l}}}\cdot \mathrm{100}\mathrm{\%}（10）

式中：\mathrm{\Delta }{m}_{\mathrm{l},\eta }为锈蚀钢筋的质量损失；m_{\mathrm{l}}为未锈蚀钢筋质量；\mathrm{\Delta }{A}_{\mathrm{l},\eta }为锈蚀钢筋截面面积损失；A_{\mathrm{l}}为未锈蚀钢筋截面面积。

锈蚀钢筋的屈服强度和极限抗拉强度可分别表示为^［32］：

f_{\mathrm{l}\mathrm{y},\eta }=(\mathrm{1}-{\alpha }_{\mathrm{y}}{\eta }_{\mathrm{l}})f_{\mathrm{l}\mathrm{y}}（11）

f_{\mathrm{l}\mathrm{u},\eta }=(\mathrm{1}-{\alpha }_{\mathrm{u}}{\eta }_{\mathrm{l}})f_{\mathrm{l}\mathrm{u}}（12）

式中：f_{\mathrm{l}\mathrm{y},\eta }和f_{\mathrm{l}\mathrm{u},\eta }分别为锈蚀钢筋的屈服强度与极限抗拉强度；f_{\mathrm{l}\mathrm{y}}和f_{\mathrm{l}\mathrm{u}}分别为未锈蚀钢筋的屈服强度与抗压极限强度；{\alpha }_{\mathrm{y}}=\mathrm{0.015}，{\alpha }_{\mathrm{u}}=\mathrm{0.013}，均为相关经验系数^［32］。

锈蚀钢筋与混凝土间的粘结滑移行为采用弹簧单元进行模拟。假定箍筋和上部纵向钢筋与混凝土完全粘结。锈蚀钢筋与混凝土间粘结强度{\tau }_{\mathrm{a},\eta }可表示为^{［23，26，31］}：

{\tau }_{\mathrm{a},\eta }=\mathrm{0.60}\left(\mathrm{0.5}+\frac{C_{\mathrm{c}}}{D_{\mathrm{l}}}\right)f_{\mathrm{c}\mathrm{t}}(\mathrm{1}-\lambda X^{\zeta })+\frac{k{A}_{\mathrm{s},\eta }{f}_{\mathrm{s}\mathrm{y},\eta }}{s{D}_{\mathrm{l}}}（13）

式中：C_{\mathrm{c}}为混凝土保护层厚度；D_{\mathrm{l}}为纵向钢筋直径；A_{\mathrm{s},\eta }为锈蚀箍筋截面面积；\lambda =\mathrm{0.4}，\zeta =\mathrm{0.1}，k=\mathrm{0.16}，均为经验系数^{［33，34］}；s为箍筋间距；f_ct为混凝土的抗拉强度，与其抗压强度f_{\mathrm{c}\mathrm{m}}有关，可表示为f_{\mathrm{c}\mathrm{t}}=\mathrm{0.56}\sqrt[]{f_{\mathrm{c}\mathrm{m}}}^［35］；f_{\mathrm{s}\mathrm{y},\eta }为锈蚀箍筋屈服强度^［36］。

f_{\mathrm{s}\mathrm{y},\eta }=\frac{\mathrm{1}-\mathrm{1.1}\cdot {\eta }_{\mathrm{s}}}{\mathrm{1}-{\eta }_{\mathrm{s}}}{f}_{\mathrm{s}\mathrm{y}}（14）

式中：{\eta }_{\mathrm{s}}为箍筋平均锈蚀率；f_{\mathrm{s}\mathrm{y}}为未锈蚀箍筋屈服强度。

弹簧单元刚度由钢筋粘结强度与相应的滑移量确定，表达式为^［37］：

K_{\mathrm{s}\mathrm{p}}=\mathrm{\pi }{D}_{\mathrm{l}}\cdot l_{\mathrm{r}}\cdot {\tau }_{\mathrm{a},\eta }/S_{\mathrm{1}}（15）

式中：l_{\mathrm{r}}为混凝土单元长度；S_{\mathrm{1}}为钢筋最大粘结应力处对应的最大滑移值。

2.2　锈蚀RC结构非线性拓扑优化公式与灵敏度分析

锈蚀RC结构荷载路径搜寻问题的优化目标是在体积约束V^{\mathrm{*}}下生成结构柔度最小（即刚度最大）的拓扑构型。锈蚀RC结构荷载路径搜寻为一个非线性拓扑优化问题，其数学表达式为：

Minimize：C(\tilde{\mathit{x}},{\mathit{U}}_{\eta \mathrm{l}})={\mathit{f}}_{\mathrm{e}\mathrm{x}\mathrm{t}}^{\mathrm{T}}{\mathit{U}}_{\eta \mathrm{l}}（16）

Subject to：\left\{\begin{array}{l}R(\tilde{\mathit{x}},{\mathit{U}}_{\eta \mathrm{l}})=\mathrm{0}\\ {\displaystyle \sum _{i=\mathrm{1}}^n}{\tilde{x}}_i{v}_i\le V^{\mathrm{*}}\\ \tilde{\mathit{x}}={\left[{\tilde{x}}_{\mathrm{1}}{\tilde{x}}_{\mathrm{2}}\dots {\tilde{x}}_n\right]}^{\mathrm{T}}\\ {\tilde{x}}_i={\tilde{x}}_{\mathrm{m}}\mathrm{o}\mathrm{r}\mathrm{1}(i=\mathrm{1,2},\dots ,n)\end{array}\right.（17）

式中：C(\tilde{\mathit{x}},{\mathit{U}}_{\eta \mathrm{l}})为目标函数，即结构柔度；\tilde{\mathit{x}}为单元设计变量的向量；{\tilde{x}}_i为第i个单元设计变量（即材料相对密度），可以取1或者{\tilde{x}}_{\mathrm{m}}；{\mathit{f}}_{\mathrm{e}\mathrm{x}\mathrm{t}}为外部荷载向量；{\mathit{U}}_{\eta \mathrm{l}}为结构位移向量；v_i和V^{\mathrm{*}}分别为第i个单元体积和目标体积分数；n为设计域中的混凝土单元数量；{\tilde{x}}_{\mathrm{m}}为单元相对密度的下限；R(\tilde{\mathit{x}},{\mathit{U}}_{\eta \mathrm{l}})为锈蚀RC结构残余力。

锈蚀RC结构荷载路径搜寻的灵敏度公式为目标函数C(\tilde{\mathit{x}},{\mathit{U}}_{\eta \mathrm{l}})对设计变量{\tilde{x}}_i的偏导数，表达式为^{［15，16，21］}：

\frac{\partial C(\tilde{\mathit{x}},{\mathit{U}}_{\eta \mathrm{l}})}{\partial {\tilde{x}}_i}=-({\mathit{f}}_{\mathrm{e}\mathrm{x}\mathrm{t}}{\mathit{K}}_{\mathrm{t},\eta \mathrm{l}}^{-\mathrm{1}}{)}^{\mathrm{T}}\frac{\partial {\mathit{f}}_{\mathrm{i}\mathrm{n}\mathrm{t}}(\tilde{\mathit{x}},{\mathit{U}}_{\eta \mathrm{l}})}{\partial {\tilde{x}}_i}（18）

式中：{\mathit{K}}_{\mathrm{t},\eta \mathrm{l}}为结构平衡状态下的切线刚度矩阵。

{\mathit{f}}_{\mathrm{e}\mathrm{x}\mathrm{t}}{\mathit{K}}_{\mathrm{t},\eta \mathrm{l}}^{-\mathrm{1}}项可近似地表示为^［15］：

{\mathit{f}}_{\mathrm{e}\mathrm{x}\mathrm{t}}{\mathit{K}}_{\mathrm{t},\eta \mathrm{l}}^{-\mathrm{1}}\approx {\mathit{U}}_{\eta \mathrm{l}}（19）

因此，锈蚀RC结构末端柔度可以近似地表示为：

\partial C(\tilde{\mathit{x}},{\mathit{U}}_{\mathrm{\eta }\mathrm{l}})\approx {\mathit{U}}_{\mathrm{\eta }\mathrm{l}}^{\mathrm{T}}{\mathit{f}}_{\mathrm{i}\mathrm{n}\mathrm{t}}(\tilde{\mathit{x}},{\mathit{U}}_{\mathrm{\eta }\mathrm{l}})（20）

当锈蚀RC结构处于平衡状态时，内力{\mathit{f}}_{\mathrm{i}\mathrm{n}\mathrm{t}}等于外部荷载{\mathit{f}}_{\mathrm{e}\mathrm{x}\mathrm{t}}，则{\mathit{U}}_{\eta \mathrm{l}}^{\mathrm{T}}{\mathit{f}}_{\mathrm{i}\mathrm{n}\mathrm{t}}(\tilde{\mathit{x}},{\mathit{U}}_{\eta \mathrm{l}})可以被看作是外荷载所做的功（即应变能）。上述灵敏度公式可以近似地重新写为：

\frac{\partial C(\tilde{\mathit{x}},{\mathit{U}}_{\eta \mathrm{l}})}{\partial {\tilde{x}}_i}=-\frac{w_{i,\eta \mathrm{l}}}{\partial {\tilde{x}}_i}（21）

式中：w_{i,\eta \mathrm{l}}为每个迭代步中单元i的应变能。

锈蚀RC结构由混凝土材料、锈损材料和钢筋组成，结构内部荷载{\mathit{f}}_{\mathrm{i}\mathrm{n}\mathrm{t}}(\tilde{\mathit{x}},{\mathit{U}}_{\eta \mathrm{l}})可表示为：

{\mathit{f}}_{\mathrm{i}\mathrm{n}\mathrm{t}}(\tilde{\mathit{x}},{\mathit{U}}_{\eta \mathrm{l}})={\mathit{f}}_{\mathrm{c}\mathrm{o},\eta \mathrm{l}}+{\mathit{f}}_{\mathrm{w}\mathrm{e}}+{\mathit{f}}_{\mathrm{c}\mathrm{r},\eta \mathrm{l}}（22）

式中：{\mathit{f}}_{\mathrm{c}\mathrm{o},\eta \mathrm{l}}、{\mathit{f}}_{\mathrm{w}\mathrm{e}}和{\mathit{f}}_{\mathrm{c}\mathrm{r},\eta \mathrm{l}}分别为混凝土单元、锈损单元和钢筋单元的内部荷载。

钢筋单元的内部荷载与材料密度无关，因此\frac{\partial {\mathit{f}}_{\mathrm{c}\mathrm{r},\eta \mathrm{l}}}{\partial {\tilde{x}}_i}=0。RC结构的灵敏度公式\frac{\partial C(\tilde{\mathit{x}},{\mathit{U}}_{\eta \mathrm{l}})}{\partial {\tilde{x}}_i}可以写成：

\frac{\partial C(\tilde{\mathit{x}},{\mathit{U}}_{\eta \mathrm{l}})}{\partial {\tilde{x}}_i}\approx -\frac{\partial w_{\mathrm{c}\mathrm{o},\eta \mathrm{l}}}{\partial {\tilde{x}}_i}-\frac{\partial w_{\mathrm{w}\mathrm{e}}}{\partial {\tilde{x}}_i}（23）

式中：w_{\mathrm{c}\mathrm{o},\eta \mathrm{l}}和w_{\mathrm{w}\mathrm{e}}分别为混凝土单元和锈损单元的应变能。

在锈蚀RC结构荷载路径搜寻过程中，单元应变能可通过非线性有限元分析直接提取得到。设计域中单元i的应变能w_{i,\eta }可表示为：

w_{i,\eta }={\tilde{x}}_i^p{w}_{\mathrm{c}\mathrm{o},\eta \mathrm{l}}+(\mathrm{1}-{\tilde{x}}_i^p)w_{\mathrm{w}\mathrm{e}}（24）

当锈损单元畸变控制方法被激活时，锈蚀RC结构灵敏度公式可表示为：

\begin{array}{l}\frac{\partial C(\tilde{\mathit{x}},{\mathit{U}}_{\eta \mathrm{l}})}{\partial {\tilde{x}}_i}=\frac{p}{{\tilde{x}}_i}{\tilde{x}}_i^p{w}_{\mathrm{c}\mathrm{o},\eta \mathrm{l}}-\\ \frac{p{\tilde{x}}_i^{p-\mathrm{1}}}{\mathrm{1}-{\tilde{x}}_i^p+{\lambda }_{\mathrm{0}}}(\mathrm{1}-{\tilde{x}}_i^p)w_{\mathrm{w}\mathrm{e}}\end{array}（25）

式中：{\lambda }_{\mathrm{0}}为非常小的常数^［16］，可以设为{\lambda }_{\mathrm{0}}={\mathrm{10}}^{-\mathrm{10}}。

2.3　荷载路径搜寻数值实现方案

锈蚀RC结构荷载路径搜寻过程包括3个阶段，即预处理阶段、拓扑搜寻阶段和局部区域锈损单元畸变控制阶段。图2为基于局部区域锈损单元畸变控制方法的RC结构荷载路径拓扑搜寻流程图，概括如下：

在预处理阶段，计算材料参数（f_{\mathrm{c}\mathrm{m}}^{\mathrm{\text{'}}}，f_{\mathrm{l}\mathrm{y},\eta }，f_{\mathrm{l}\mathrm{u},\eta }，{\tau }_{\mathrm{a},\mathrm{\eta }}），定义混凝土和钢筋的材料本构模型。建立锈蚀RC结构有限元模型，定义有限元分析输出（即单元应变能、应变等）。

在拓扑搜寻阶段，拓扑优化迭代过程与BESO方法相似。

步骤1：定义初始的优化参数（r_{\mathrm{m}\mathrm{i}\mathrm{n}}，V^{\mathrm{*}}，\mathrm{E}\mathrm{R}，{\epsilon }^{\mathrm{*}}，E_{\mathrm{m}\mathrm{a}\mathrm{x}}和E_{\mathrm{m}\mathrm{a}\mathrm{x}}^{\mathrm{0}}）。

步骤2：基于ABAQUS进行非线性有限元分析，提取各单元的弹、塑性应变能。单元i的应变能w_i可表示为^［38］：

w_i=w_i^{\mathrm{e}}+w_i^{\mathrm{p}}（26）

式中：w_i^{\mathrm{e}}和w_i^{\mathrm{p}}分别为第i单元的弹性应变能和塑性应变能。

步骤3：用式（25）计算单元灵敏度。

步骤4：通过过滤和平均历史信息来更新单元灵敏度。采用Yang等^［24］提出的模糊滤波方案得到网格无关性的荷载路径。该滤波方案可表示为：

{\widehat{\alpha }}_i=\frac{\partial C(\tilde{\mathit{x}},{\mathit{U}}_{\eta \mathrm{l}})}{\partial {\tilde{x}}_i}={\sum }_j\left(\frac{w(r_{ij})}{{\sum }_{j}w(r_{ij})}{\alpha }_j\right)={\sum }_j{\eta }_j{\alpha }_j（27）

w(r_{ij})=\mathrm{m}\mathrm{a}\mathrm{x}(\mathrm{0},r_{\mathrm{m}\mathrm{i}\mathrm{n}}-r_{ij})（28）

式中：r_{ij}为单元i和单元j中心之间的距离；w为平均原始灵敏度权重函数；{\eta }_j为权重因子；r_{\mathrm{m}\mathrm{i}\mathrm{n}}为过滤半径。当前迭代步的敏感度{\widehat{\alpha }}_i与前一个迭代步的敏感度需要进一步平均^［25］，以获得收敛的解决方案。单元敏感度{\widehat{\alpha }}_i可以表示为：

{\widehat{\alpha }}_i=\frac{{\widehat{\alpha }}_i^k+{\widehat{\alpha }}_i^{k+\mathrm{1}}}{\mathrm{2}}（29）

步骤5：确定下一个优化迭代步的目标体积。下一迭代步的目标体积V^{k+\mathrm{1}}可以表示为^［39］：

V^{k+\mathrm{1}}=V^k(\mathrm{1}\pm \mathrm{E}\mathrm{R})（30）

式中：V^k为当前迭代步的体积；\mathrm{E}\mathrm{R}为进化率。

步骤6：更新设计变量和单元类型。若混凝土单元灵敏度{\widehat{\alpha }}_i≤灵敏度阈值{\widehat{\alpha }}_{\mathrm{t}\mathrm{h}}，则单元密度{\tilde{x}}_i由1切换为{\tilde{x}}_{\mathrm{m}}。对于最小密度单元，若{\widehat{\alpha }}_i>{\widehat{\alpha }}_{\mathrm{t}\mathrm{h}}，则单元密度{\tilde{x}}_i从{\tilde{x}}_{\mathrm{m}}切换为1。再根据单元密度{\tilde{x}}_i更新混凝土和最小密度单元类型。

在局部区域锈损单元畸变控制阶段，首先，提取锈损区单元集合{\Omega }_{\mathrm{v}}内的最大应变{\epsilon }_{\mathrm{m}\mathrm{a}\mathrm{x}}^k，并确定其单元编号。其次，明确畸变度最大单元过滤半径内的单元编号，并去除过滤半径内实体单元，形成需要进行刚度缩放的单元集合{\Omega }_{\mathrm{n}\mathrm{m}}。最后，比较锈损单元的最大应变{\epsilon }_{\mathrm{m}\mathrm{a}\mathrm{x}}^k与应变阈值{\epsilon }^{\mathrm{*}}的大小，如果{\epsilon }_{\mathrm{m}\mathrm{a}\mathrm{x}}^k>{\epsilon }^{\mathrm{*}}，增大下一个迭代步中局部区域锈损单元刚度{\mathit{k}}_{\mathrm{m},\eta }^k。如果{\epsilon }_{\mathrm{m}\mathrm{a}\mathrm{x}}^k\le {\epsilon }^{\mathrm{*}}，减小下一个迭代步中的{\mathit{k}}_{\mathrm{m},\eta }^k，以抑制因单元刚度增加而引起的优化误差。

重复上述优化迭代过程，直到满足目标体积和收敛准则为止。锈蚀RC结构荷载路径搜寻终止的收敛准则为^［30］：

E_{\mathrm{r}}^k=\frac{\left|{\displaystyle \sum _{i=\mathrm{1}}^N}(C_{k-i+\mathrm{1}}-C_{k-N+i+\mathrm{1}})\right|}{{\displaystyle \sum _{i=\mathrm{1}}^N}{C}_{k-i+\mathrm{1}}}\le \theta（31）

式中：\theta为允许的收敛误差，通常取\theta =\mathrm{0.01}\mathrm{\%}；N为迭代步数量，可取N=\mathrm{5}。

当结构柔度在至少连续10次迭代步中保持稳定时，荷载路径搜寻过程终止。

3 验证与讨论

本文首先通过未锈蚀RC梁的数值算例验证了所提方法的有效性与适用性。其次，研究了不同锈蚀水平下RC梁的荷载路径。最后，讨论了锈蚀钢筋的抗拉强度损失和粘结退化对荷载路径的影响。

3.1　验证

（1） 未锈蚀三维RC梁荷载路径搜寻

采用RC梁荷载路径搜寻来说明本文所提方法的有效性与适用性，如图3所示。由于RC梁结构的对称性，取该结构的1/4进行荷载路径搜寻以节约计算成本。1/4 RC梁尺寸为10 mm×15 mm×60 mm，采用六面体单元网格将设计域离散为\mathrm{10}\times \mathrm{15}\times \mathrm{60}个网格。集中力和约束施加于刚性垫板节点上，以防止局部混凝土单元过度变形。采用嵌入约束将混凝土单元节点分别与上部纵向钢筋和箍筋单元节点进行连接。采用弹簧单元连接下部纵向钢筋节点与混凝土单元节点。

在本算例中，惩罚因子P、过滤器半径r_{\mathrm{m}\mathrm{i}\mathrm{n}}、目标体积分数V^{\mathrm{*}}和进化率\mathrm{E}\mathrm{R}分别设置为3、4 mm、15%和0.04。局部刚度自适应缩放策略中系数{\alpha }_{\mathrm{1}}、{\alpha }_{\mathrm{2}}、{\alpha }_{\mathrm{3}}和{\alpha }_{\mathrm{4}}分别取1.2、0.5、1和0.9。锈损材料的初始弹性模量和泊松比分别设置为E_{\mathrm{m},\eta }^{\mathrm{0}}={\mathrm{10}}^{-\mathrm{4}}{E}_{\mathrm{c}}和v_{\mathrm{0}}=\mathrm{0.4}。不同工况下锈损材料的弹性模量阈值E_{\mathrm{m}\mathrm{a}\mathrm{x}}分别设置为1、10、25、50、100、150 MPa。混凝土的力学性能如下：弹性模量为32 000 MPa，泊松比为0.2，抗压强度为40 MPa，极限应变为0.003 0，滑移植为1 mm。钢筋材料的力学性能如下：弹性模量为190 000 MPa，泊松比为0.3，钢筋直径为3 mm，屈服强度为400 MPa，极限强度为580 MPa。

图4给出了基于BESO方法的未锈蚀RC梁3D最优荷载路径。图4（a）为无埋置钢筋的弹性RC梁3D荷载路径，其中混凝土被模拟为各向同性弹性材料，该荷载路径中存在受拉混凝土轨迹，而实际结构中混凝土材料难以承受较大的拉应力，因此该荷载路径与结构传力机理存在一定差异。当考虑埋置钢筋与材料非线性等因素的影响时，从图4（b）中可以看出RC梁底部受拉的混凝土轨迹逐渐被纵筋替代，但由于出现了由单元畸变引起的数值失稳问题，导致了荷载路径搜寻过程终止。在最后的迭代步中，局部区域钢筋出现了应力集中，使得周边具有无限小刚度的最小密度单元发生畸变。因此，局部区域单元畸变会导致结构荷载路径搜寻过程中出现数值失稳问题。

图5给出了基于本文方法的不同荷载下未锈蚀RC梁3D荷载路径。在最小荷载F=\mathrm{1}\mathrm{N}的情况下，由于混凝土材料拉-压性能不对称特征，图4（a）中底部受拉混凝土轨迹在优化过程中被去除，其拉应力由图5（a）中底部钢筋承受。此外，两者在剪跨段的传力轨迹也基本一致，外部荷载主要由加载点传递至底部钢筋，一部分外力由底部钢筋传递至支撑点，另一部分由倾斜混凝土杆传递至支撑点，如图5（a）所示。在较小外部荷载作用下结构仍处于弹性阶段，因此图5（a）与图4（a）中的荷载路径除顶部水平受压混凝土轨迹外基本一致。

当载荷F从100 N增加到800 N时，RC梁的荷载路径从加载点直接延伸至支撑点，如图5（b）~（d）所示。在此过程中，RC梁底部逐渐开裂，从跨中竖向裂缝向剪跨段倾斜裂缝发展。RC梁被分割为上部较大的拱形混凝土块与下部较小的混凝土柱，低效的混凝土单元会被逐渐去除，外部荷载主要由上部的拱形混凝土块直接传递至支撑点，RC梁底部拉力由纵筋承受。

当荷载F增加到1 000 N时，RC梁的荷载路径从加载点向梁1/3处延伸，将一部分外部荷载传递至底部钢筋，上部仍为拱形的荷载路径，在剪跨段两条荷载路径轨迹几乎平行，体现了RC梁的剪切行为，如图5（e）所示。RC梁的荷载路径与其传力机理基本一致。随着荷载增大到1 000 N，混凝土裂缝宽度逐渐变大，外部荷载一部分由上部拱形混凝土块传递至支撑点。另一部分外部荷载从加载点由下部混凝土块向剪跨段底部纵筋传递，再由纵筋传递至支撑点。其中，部分荷载由箍筋传递至上部拱形混凝土块，再传递至支撑点。随着荷载进一步增大到1 500 N，RC梁荷载路径轨迹与上一工况下的荷载路径基本保持一致，如图5（f）所示。

综上所述，本文方法可合理地生成RC梁荷载路径，且无数值失稳问题。实际工程中RC构件多受弯、剪、扭等复合作用，复合受力作用会改变RC构件的应变能分布，导致结构荷载路径发生变化，本文所提方法是基于BESO方法改进所得，其适用性并不会受边界条件变化的影响^［40］。

（2）锈蚀三维RC梁荷载路径搜寻

图6为不同锈蚀水平下（0%~20%）锈蚀RC梁在外部荷载F=\mathrm{800}\mathrm{N}时的荷载路径变化规律。锈蚀RC梁的荷载路径随着锈蚀率的增加逐渐发生变化。钢筋锈蚀率增加与荷载增大下锈蚀RC梁的荷载路径变化规律基本一致。结果表明，锈蚀会加速RC梁的荷载路径由非极限状态向极限状态转变。这一现象可能是由于锈蚀RC梁的极限承载力随着锈蚀率的增加而降低所导致的。随着钢筋锈蚀率的增加，外部荷载F=\mathrm{800}\mathrm{N}可能成为RC梁的极限荷载。

3.2　讨论

图7为不同粘结强度下锈蚀RC梁的荷载路径。根据式（13）可确定锈蚀RC梁的粘结强度，锈蚀率分别取0%、3%、5%、10%、15%和20%。随着锈蚀钢筋与混凝土间粘结强度的降低，锈蚀RC梁荷载路径逐渐发生变化。由于钢筋与混凝土间粘结性能的退化，在荷载F=\mathrm{800}\mathrm{N}作用下的荷载路径逐渐转变为极限状态下荷载路径。

图8为不同抗拉强度下锈蚀RC梁的荷载路径。采用与上文相同的锈蚀率（即0%、3%、5%、10%、15%和20%），可根据式（11）和式（12）确定锈蚀钢筋的抗拉强度。不同抗拉强度下锈蚀RC梁的荷载路径基本相同，除了图8（f）中的荷载路径。这是因为在锈蚀钢筋力学性能损失下RC梁未达到极限状态。结果表明，在相同锈蚀水平下，荷载路径对粘结退化比对钢筋材料性能损失更为敏感。

4 结　论

（1）本文方法能够合理地生成锈蚀RC梁的3D荷载路径。锈蚀RC梁的3D荷载路径随着锈蚀程度与外部荷载的增加逐渐发生变化。锈蚀损伤会加速锈蚀RC梁荷载路径由非极限状态向极限状态转换。

（2）本文提出的局部区域锈损单元畸变控制方法可有效地克服由锈损单元畸变引起的数值失稳问题，给出的局部刚度自适应缩放策略可合理地根据单元畸变度对单元刚度进行多比例缩放，有效地抑制优化误差。

（3）锈损RC梁的3D荷载路径会随粘结强度与钢筋抗拉强度的退化而发生变化。在相同锈蚀水平下，锈损RC结构的荷载路径发展规律对粘结退化比对锈蚀钢筋的力学性能损失更为敏感。

参考文献

[1] 陈华鹏, 刘昌雨, 江钰, 等. 钢筋混凝土结构锈蚀开裂模型研究[J]. 铁道学报, 2023, 45(10): 180-188.

[2] Chen Hua-peng, Liu Chang-yu, Jiang Yu,et al. Cracking model of reinforced concrete structure due to rebar corrosion[J]. Journal of the China Railway Society, 2023, 45(10): 180-188.

[3] 戴理朝, 周亮, 杨晓文, 等. 基于Connector单元的锈蚀RC梁界面粘结性能细观数值模拟[J]. 吉林大学学报: 工学版, 2023, 53(10): 2886-2896.

[4] Dai Li-zhao, Zhou Liang, Yang Xiao-wen, et al. Meso-scale numerical simulation of interfacial bond behavior of corroded RC beams based on Connector element[J]. Journal of Jilin University(Engineering and Technology Edition), 2023, 53(10): 2886-2896.

[5] 张艺博, 郑山锁, 董立国, 等. 往复荷载作用下锈蚀钢筋与混凝土黏结性能试验研究[J]. 建筑结构学报, 2024, 45(2): 237-248.

[6] Zhang Yi-bo, Zheng Shan-suo, Dong Li-guo, et al. Experimental study on bond behavior of corroded reinforcements in concrete under reversed cyclic loading[J]. Journal of Building Structures, 2024, 45(2): 237-248.

[7] Ye Z, Zhang W, Gu X. Deterioration of shear behavior of corroded reinforced concrete beams[J]. Engineering Structures, 2018, 168: 708-720.

[8] 张仁波, 杨鸿申, 金浏, 等. 锈蚀钢筋混凝土梁动态抗弯性能细观数值研究[J]. 振动与冲击, 2023, 42(9): 77-85.

[9] Zhang Ren-bo, Yang Hong-shen, Jin Liu, et al. Meso numerical study on dynamic anti-bending performance of corroded RC beam[J]. Journal of Vibration and Shock, 2023, 42(9): 77-85.

[10] 刘霞, 易伟建. 开孔深梁压杆-拉杆模型构造[J]. 工程力学, 2012, 29(12): 141-146.

[11] Liu Xia, Yi Wei-jian. Strut-and-tie model construction for deep beams with openings[J]. Engineering Mechanics, 2012, 29(12): 141-146.

[12] Bruggi M. Generating strut-and-tie patterns for reinforced concrete structures using topology optimization[J]. Computers & Structures, 2009, 87(23/24): 1483-1495.

[13] Luo Y, Kang Z. Layout design of reinforced concrete structures using two-material topology optimization with Drucker-Prager yield constraints[J]. Structural and Multidisciplinary Optimization, 2013, 47(1): 95-110.

[14] Leu L J, Huang C W, Chen C S, et al. Strut-and-tie design methodology for three-dimensional reinforced concrete structures[J]. Journal of Structural Engineering, 2006, 132(6): 929-938.

[15] Liang Q Q, Xie Y M, Steven G P. Topology optimization of strut-and-tie models in reinforced concrete structures using an evolutionary procedure[J]. ACI Structural Journal, 2000, 97(2): 322-332.

[16] Amir O. A topology optimization procedure for reinforced concrete structures[J]. Computers and Structures, 2013, 114: 46-58.

[17] Bendsøe M P. Optimal shape design as a material distribution problem[J]. Structural Optimization, 1989, 1(4): 193-202.

[18] Zhou M, Rozvany G I N. The COC algorithm, Part II: Topological, geometrical and generalized shape optimization[J]. Computer Methods in Applied Mechanics and Engineering, 1991, 89(1-3): 309-336.

[19] Xie Y M, Steven G P. A simple evolutionary procedure for structural optimization[J]. Computers & Structures, 1993, 49(5): 885-896.

[20] Chen Q, Zhang X, Zhu B. A 213-line topology optimization code for geometrically nonlinear structures[J]. Structural and Multidisciplinary Optimization, 2019, 59(5): 1863-1879.

[21] Luo Y, Wang M Y, Kang Z. Topology optimization of geometrically nonlinear structures based on an additive hyperelasticity technique[J]. Computer Methods in Applied Mechanics and Engineering, 2015, 286: 422-441.

[22] Rad M M, Habashneh M, Lógó J. Reliability based bi-directional evolutionary topology optimization of geometric and material nonlinear analysis with imperfections[J]. Computers & Structures, 2023, 287: No.107120.

[23] Buhl T, Pedersen C B W, Sigmund O, et al. Stiffness design of geometrically nonlinear structures using topology optimization[J]. Structural and Multidisciplinary Optimization, 2000, 19(2): 93-104.

[24] Pedersen C B W, Buhl T, Sigmund O. Topology synthesis of large‐displacement compliant mechanisms[J]. International Journal for Numerical Methods in Engineering, 2001, 50(12): 2683-2705.

[25] Sigmund O. Design of multiphysics actuators using topology optimization—part I: one-material structures[J]. Computer Methods in Applied Mechanics and Engineering, 2001, 190(49/50): 6577-6604.

[26] Wang F, Lazarov B S, Sigmund O, et al. Interpolation scheme for fictitious domain techniques and topology optimization of finite strain elastic problems[J]. Computer Methods in Applied Mechanics and Engineering, 2014, 276: 453-472.

[27] Yoon G H, Kim Y Y. Element connectivity parameterization for topology optimization of geometrically nonlinear structures[J]. International Journal of Solids and Structures, 2005, 42(7): 1983-2009.

[28] Yuan P, Wang L, Floyd R W. Generation of optimal load paths for corroded reinforced concrete beams—part I: automatic stiffness adjustment technique[J]. ACI Structural Journal, 2023, 120(4): 103-114.

[29] Yang X Y, Xie Y M, Steven G P, et al. Bidirectional evolutionary method for stiffness optimization[J]. AIAA Journal, 1999, 37(11): 1483-1488.

[30] Huang X, Xie Y M. Convergent and mesh-independent solutions for the bi-directional evolutionary structural optimization method[J]. Finite Elements in Analysis and Design, 2007, 43(14): 1039-1049.

[31] Kmiecik P, Kamiński M. Modelling of reinforced concrete structures and composite structures with concrete strength degradation taken into consideration[J]. Archives of Civil and Mechanical Engineering, 2011, 11(3): 623-636.

[32] Recupero A, Spinella N, Tondolo F. Failure analysis of corroded RC beams subjected to shear-flexural actions[J]. Engineering Failure Analysis, 2018, 93: 26-37.

[33] Antoniou M, Mantakas A, Nikitas N, et al. A numerical case study on the long-term seismic assessment of reinforced concrete tunnels in corrosive environments[J]. Journal of Rock Mechanics and Geotechnical Engineering, 2023, 15(3): 551-572.

[34] Molina F J, Alonso C, Andrade C, et al. Cover cracking as a function of rebar corrosion: part 2—Numerical model[J]. Materials and Structures, 1993, 26(9): 532-548.

[35] Dekoster M, Buyle-Bodin F, Maurel O, et al. Modelling of the flexural behaviour of RC beams subjected to localised and uniform corrosion[J]. Engineering Structures, 2003, 25(10): 1333-1341.

[36] Chen H P, Nepal J. Analytical model for residual bond strength of corroded reinforcement in concrete structures[J]. Journal of Engineering Mechanics, 2016, 142(2): No.04015079.

[37] Cairns J, Plizzari G A, Du Y, et al. Mechanical properties of corrosion-damaged reinforcement[J]. ACI Materials Journal, 2005, 102(4): 256-264.

[38] Campione G, Cannella F, Cavaleri L. Shear and flexural strength prediction of corroded RC beams[J]. Construction and Building Materials, 2017, 149: 395-405.

[39] Rodriguez J, Ortega L, Garcia A. Corrosion of reinforcing bars and service life of RC structures: corrosion and bond deterioration[J]. Concrete Across Borders, 1994, 2: 315-326.

[40] ACI Committee. Building Code Requirements for Structural Concrete(ACI 318-11) and Commentary[M]. Farmington Hills: American Concrete Institute, 2011.

[41] Lin H, Zhao Y, Yang J Q, et al. Effects of the corrosion of main bar and stirrups on the bond behavior of reinforcing steel bar[J]. Construction and Building Materials, 2019, 225: 13-28.

[42] Khalaf J, Huang Z. Analysis of the bond behaviour between prestressed strands and concrete in fire[J]. Construction and Building Materials, 2016, 128: 12-23.

[43] Huang X, Xie Y M. Bidirectional evolutionary topology optimization for structures with geometrical and material nonlinearities[J]. AIAA Journal, 2007, 45(1): 308-313.

[44] Bi M, Tran P, Xia L, et al. Topology optimization for 3D concrete printing with various manufacturing constraints[J]. Additive Manufacturing, 2022, 57: No.102982.

[45] Huang X, Xie Y M. A further review of ESO type methods for topology optimization[J]. Structural and Multidisciplinary Optimization, 2010, 41: 671-683.

基金资助

国家重点研发计划项目(2021YFB2600900)

国家自然科学基金项目(52278140)

长沙理工大学研究生科研创新项目(CXCLY2022035)