一种基于荷载路径的锈蚀损伤混凝土结构加固设计方法

1.本发明涉及混凝土结构非线性拓扑优化技术领域，尤其涉及一种基于荷载路径的锈蚀损伤混凝土结构加固设计方法。


背景技术：

2.混凝土桥梁在我国现役桥梁中占有较大比重。混凝土结构在遭受锈蚀损伤时易出现钢筋截面积减小、保护层锈胀开裂以及混凝土与力筋间粘结退化等病害，使得混凝土结构内部应力分布变得紊乱，不再符合平截面假定，该区域即为混凝土结构应力扰动区(d区)。混凝土结构d区是关键的受力构件与部位，其抗裂性能及承载力涉及到整个结构的耐久性和安全性。目前，复杂边界条件下锈蚀损伤混凝土结构d区的传力机理仍难以准确地表征。拉压杆模型被认为是混凝土结构d区设计最为有效的方法。在采用拉压杆模型对结构d区进行设计时，最为关键的问题是如何确定拉压杆模型的构型。目前，应力迹线法、荷载路径法和拓扑优化法已被成功应用于混凝土结构d区的拉压杆模型设计。然而，基于应力迹线法与荷载路径法得到的拉压杆模型并不是唯一的，其往往依赖于设计者的直觉和经验，具有较大的主观性。
3.为克服基于应力迹线法与荷载路径法的主观性，一些学者借助拓扑优化方法对混凝土结构的拉压杆模型进行了初步研究。现有方法侧重于健全的弹性混凝土结构荷载路径搜寻，进而确定相应的拉压杆模型。然而，锈蚀损伤混凝土结构荷载路径搜寻涉及材料劣化与粘结退化等损伤，且为一个非线性拓扑优化问题，其相较于健全的弹性混凝土结构更加复杂。首先，锈蚀损伤混凝土结构荷载路径搜寻过程中材料密度的转换会引起最小密度单元发生畸变，导致荷载路径搜寻过程的数值失稳问题。然后，锈蚀导致的材料劣化、粘结退化及混凝土开裂等损伤及其耦合作用大大增加了荷载路径搜寻过程的复杂程度。目前，锈蚀损伤混凝土结构荷载路径搜寻研究尚未见报道。锈蚀损伤混凝土结构荷载路径对于描述其d区传力机理、评估结构性能与指导维修加固决策具有重要意义。
4.现有的混凝土结构加固方法有粘贴钢板、增大截面、体外预应力、纤维复合材料加固、新型材料加固及组合加固等，但目前这些加固方法多是全截面布设，存在自重大、成本高与周期长等问题。理想的加固设计应遵循结构荷载传递路径，充分发挥加固材料性能，提高其利用效率。结构损伤程度与外部荷载类型等都会改变结构荷载路径与加固材料布设轨迹。如何基于荷载路径对锈蚀损伤混凝土结构加固材料进行布设也有待进一步研究。因此，亟需一种基于荷载路径的锈蚀损伤混凝土结构加固设计方法，既能合理地生成锈蚀损伤混凝土结构荷载路径，又能有效地克服由最小密度单元畸变引起的数值失稳问题，还可以基于损伤混凝土结构荷载路径指导加固材料的布设。


技术实现要素：

5.为解决上述技术问题，本发明提出一种基于荷载路径的锈蚀损伤混凝土结构加固设计方法，能够合理地生成锈蚀损伤混凝土结构荷载路径，又能够有效地克服由最小密度
单元畸变引起的数值失稳问题，还可以基于损伤钢筋混凝土(rc)结构荷载路径指导加固材料的布设。
6.为实现上述目的，本发明提供一种基于荷载路径的锈蚀损伤混凝土(rc)结构加固设计方法，其主要步骤包括：
7.s1：确定锈蚀损伤rc结构荷载路径搜寻问题的非线性拓扑优化公式；
8.s2：计算锈蚀材料性能、粘结强度与弹簧刚度；
9.s3：建立考虑材料性能损失与粘结退化的锈蚀损伤rc结构数值模型；
10.s4：对锈蚀损伤rc结构进行非线性拓扑优化分析：首先，对步骤s3建立的锈蚀损伤rc结构数值模型进行有限元分析；然后，获取各单元的弹、塑性应变能，基于所给出的锈蚀损伤单元应变能形式的灵敏度公式计算各单元灵敏度；最后，采用模糊滤波方案更新单元灵敏度信息，确定下一个优化步的目标体积，更新设计变量和单元类型；
11.s5：基于所提出的最小密度单元刚度自动调节方法调整最小密度单元刚度；
12.s6：计算锈蚀损伤rc结构荷载路径搜寻过程中的收敛误差；
13.s7：生成锈蚀损伤rc结构荷载路径，并验证结构荷载路径的合理性；
14.s8：基于锈蚀损伤rc结构荷载路径指导加固材料布设。
15.进一步的，所述步骤s1中，损伤结构荷载路径搜寻问题的非线性拓扑优化公式如下：
16.锈蚀损伤rc结构荷载路径搜寻问题的优化目标是在体积约束v
*
下生成结构柔度最小(即刚度最大)的荷载路径，将混凝土单元材料密度作为设计变量，锈蚀损伤rc结构非线性优化问题的数学表达式为：
[0017][0018][0019]
其中，为目标函数，也就是结构柔度；为单元设计变量的向量；为第i个单元的材料密度，在双向渐进结构优化方法(beso方法)中可以分别取为1或者f
ext
为外部荷载向量；u
η
为结构位移向量；vi和v
*
分别为第i个单元的体积和规定材料的体积分数。n为设计域中的混凝土单元数量；为单元最小密度；为非线性锈蚀损伤rc结构的残余力向量。
[0020]
进一步的，所述步骤s2中，锈蚀材料性能、粘结强度与弹簧刚度的计算方法如下：
[0021]
1)锈蚀材料性能包括锈蚀混凝土抗压强度计算和锈蚀钢筋极限抗拉强度，锈蚀混凝土的抗压强度f
′
cm
可以表示为：
[0022][0023]
其中，f
cm
为未锈蚀损伤混凝土的抗压强度；kr为经验系数；ε
cu
为最大压缩应变；ε1为开裂混凝土的平均受拉应变，可表示为：
[0024]
ε1＝n
bs
·wcr
/bwꢀꢀꢀꢀ
(4)
[0025]
其中，bw为未锈蚀损伤rc梁的初始宽度；n
bs
为rc梁受压区钢筋数量；w
cr
为锈蚀损伤
rc梁总裂缝宽度，可计算为：
[0026]wcr
＝2πx(u
rs-1)
ꢀꢀꢀꢀ
(5)
[0027]
其中，u
rs
为体积膨胀率；x为钢筋锈蚀深度，可表示为：
[0028][0029]
其中，η和rs分别为钢筋锈蚀率和钢筋半径。
[0030]
2)锈蚀钢筋极限抗拉强度可表示为：
[0031]fy，c
＝(1-αyη)fyꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(7)
[0032]fu，c
＝(1-αuη)fuꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(8)
[0033]
ε
u，c
＝(1-α1η)εuꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(9)
[0034]
其中，f
y，c
，f
u，c
与ε
u，c
分别为锈蚀钢筋的屈服强度，极限抗拉强度和极限应变；fy，fu与εu分别为未锈蚀钢筋的屈服强度、极限强度和极限应变；αy，αu与α1为相关经验系数；
[0035]
3)锈蚀钢筋与混凝土间粘结强度τ
a，η
可表示为：
[0036][0037]
其中，cc为混凝土保护层厚度；d
l
为纵向钢筋的直径；a
s，η
为锈蚀箍筋截面积；λ，ζ，k为经验系数；s是箍筋间距；f
ct
为混凝土的抗拉强度，f
cm
为混凝土的抗压强度；f
sy，η
为锈蚀箍筋的屈服强度，可表示为：
[0038][0039]
其中，ηs为箍筋平均锈蚀率；f
sy
为未锈蚀箍筋的屈服强度。弹簧单元的刚度可表示为：
[0040]ksp
＝π
·d·
lr·
τ
a，η
/s1ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(12)
[0041]
其中，lr为混凝土单元的长度；s1为钢筋最大粘结应力处对应的最大滑移值。
[0042]
进一步的，所述步骤s4中，给出的锈蚀损伤单元应变能形式的灵敏度公式推导过程如下：
[0043]
单元灵敏度为目标函数对设计变量的偏导数，其表达式如下所示：
[0044][0045]
其中，k
t，η
为结构平衡状态下的切线刚度矩阵；为结构平衡状态下的内力。公式(13)中的项在锈蚀损伤rc结构非线性有限元分析中可近似地表示为：
[0046][0047]
因此，锈蚀损伤rc结构末端柔度可以近似地表示为：
[0048][0049]
锈蚀损伤rc结构处于平衡状态时内力f
int
等于外部荷载f
ext
。此项可以被看作是外荷载所做的功(即应变能)。灵敏度公式可以近似地重新写为：
[0050][0051]
其中，w
i，η
为每个迭代步中单元i的应变能。锈蚀损伤rc结构由混凝土材料、弱弹性材料和钢筋组成。结构内部荷载可以写成：
[0052][0053]
其中，f
co，η
，f
we
和f
cr，η
分别是混凝土单元、弱弹性材料单元和钢筋单元的内部荷载。钢筋单元的内部荷载与材料虚拟密度无关。因此，在式(17)中此项等于零。rc结构的灵敏度公式可以重新写成：
[0054][0055]
其中，w
co，η
和w
we
分别为混凝土和弱弹性材料单元的应变能。
[0056]
在非线性拓扑优化过程中，锈蚀损伤rc结构在平衡状态下的应变能更易于操作，可通过非线性有限元分析直接提取得到各单元的应变能。设计域中单元i的应变能w
i，η
为：
[0057][0058]
锈蚀损伤rc结构灵敏度公式可表示为：
[0059][0060]
其中，λo是一个非常小的常数，可以设为λo＝10-10
。λ加到右边第二项的分母上，以避免在时出现歧义。
[0061]
进一步的，所述步骤s5中，提出的最小密度单元刚度自动调节方法，其步骤包括：
[0062]
首先，提取各优化步中最小密度单元的最大应变然后，判断最大应变与应变阈值ε
*
的大小；如果使用刚度多比例增长策略增大下一个迭代步中最小密度单元刚度，以抑制荷载路径搜寻过程中最小密度单元畸变引起的数值失稳问题；如果则减小下一个迭代步中最小密度单元刚度，以抑制因单元刚度增加而引起的优化误差；
[0063]
所述由最小密度单元畸变引起的数值失稳问题如下：
[0064]
在beso方法中，单元的删除是通过将相对密度从1.0转换到来实现的，而不是完全删除。各单元刚度矩阵ki与相对密度的关系为：
[0065][0066]
其中，b为单元应变矩阵；kc为混凝土单元的刚度矩阵；p为惩罚因子。设计变量表明第i个单元为“空心”单元，也就是最小密度单元，其刚度矩阵为最小密度通常设置为惩罚因子取为p＝3。因此，最小密度单元的刚度被人为地降低了10-9
倍。这个无限小的刚度会引起最小密度单元的畸变，从而导致荷载路径搜寻过程的数值失稳问题。
[0067]
所述最小密度单元刚度多比例增长策略如下：
[0068]
最小密度单元的刚度极小，可被视为一种弱弹性材料，其本构关系可以表示为：
[0069][0070]
其中，de为弱弹性材料的本构矩阵；e
min
与ve分别为弹性模量和泊松比。因此，可通过改变材料弹性模量来调节单元刚度。
[0071]
在非线性拓扑优化过程中，弹性模量e
min
值过小会导致数值失稳问题，但过度的增大弹性模量ee会带来优化误差。因此，给出了最小密度单元刚度多比例增长策略，在每个优化迭代步中可根据最小密度单元的最大应变增量以不同比例调节其单元刚度，既可以克服由最小密度单元畸变导致的数值失稳问题，还可以抑制由增大单元刚度导致的优化误差。在每次优化迭代中，通过追踪最小密度单元的最大应变自动调节其弹性模量在第k 1个优化步骤中，可使用以下准则进行自动调节：
[0072][0073]
其中，与分别为最小密度单元在第k和k 1优化步中的弹性模量；是第k个优化步中最小密度单元的最大应变；α1，α2，α3和α4是用来控制单元刚度增长速率的系数；在每个优化步骤中将阈值ε
*
更新为f(x)为sigmoid函数，可以表示为：
[0074][0075]
其中，采用sigmoid函数可以逐步地调整弹性模量如图3所示。
[0076]
在非线性拓扑优化过程中，如果则意味着结构中部分最小密度单元可能会发生畸变。在下一个优化迭代步中，应增加最小密度单元的弹性模量以增强其单元刚度，进而防止最小密度单元畸变。如果则应在下一个优化迭代步中降低以抑制由增加最小密度单元刚度引起的优化误差。在荷载路径搜寻过程中，弱弹性材料的弹性模量最终接近于阈值e
max
，如此的设置可以克服锈蚀损伤rc结构荷载路径搜寻过程中数值失稳问题。
[0077]
进一步的，所述步骤s4中，各单元的总应变能wi可表示为：
[0078][0079]
其中，与分别为各单元的弹、塑性应变能。
[0080]
进一步的，所述步骤s4中，模糊滤波方案与下一个优化步的目标体积的确定方法
如下：
[0081]
采用模糊滤波方案来得到网格无关性的拓扑优化解，在该滤波方案中需对原始单元灵敏度进行进一步处理，该滤波方案可表示为：
[0082][0083]
w(r
ij
)＝max(0，r
min-r
ij
)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(27)
[0084]
其中，r
ij
为单元i和j中心之间的距离；w为平均原始灵敏度的权重函数；ηj为权重因子；r
min
为过滤半径。当前迭代步的灵敏度与前一个迭代步的灵敏度进一步平均，以获得收敛的解决方案。单元灵敏度可以表示为：
[0085][0086]
根据当前迭代步的体积vk和进化比er，可以得到下一迭代步的目标体积v
k 1
。
[0087]vk 1
＝vk(1
±
er)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(29)
[0088]
进一步的，所述步骤s4中，设计变量和单元类型更新方法如下：
[0089]
如果混凝土单元的灵敏度小于灵敏度阈值则单元密度(也就是设计变量)由1切换为对于弱弹性体单元，如果单元密度从切换为1；然后根据单元密度对混凝土和弱弹性材料截面属性进行重新分配。
[0090]
进一步的，所述步骤s4中，锈蚀损伤rc结构荷载路径搜寻过程中的收敛准则为：
[0091][0092]
其中，erk为锈蚀损伤rc结构的收敛误差；∈为允许的收敛误差；n为迭代步数量。
[0093]
进一步的，对步骤s4的锈蚀损伤rc结构进行拓扑优化分析前，还需要定义初始优化参数：如过滤半径r
min
，惩罚因子p，目标体积分数v
*
，进化率er，初始应变阈值ε
*
，材料最小密度弹性模量阈值e
max
和初始弹性模量
[0094]
初始化设计变量与混凝土截面属性：设置“实体”与“空心”材料截面属性，“实体”材料即为混凝土材料；“空心”材料为一种弹性材料，其泊松比为0.3，弹性模量为其中ec为混凝土弹性模量；将所有的初始设计变量设置为设计域内所有单元截面属性设置为“实体”材料截面属性，后续将根据每个优化步中材料密度取值对截面属性进行变换。
[0095]
进一步的，所述步骤s8中，基于锈蚀损伤rc结构荷载路径指导加固材料布设方法如下：
[0096]
理想的加固材料布设应遵循结构荷载传递路径，充分发挥加固材料性能，提高其利用效率。基于损伤结构的荷载路径分析其中每个构件受力特性，对于受拉构件则在相应位置布设抗拉强度较高的加固材料，对于受压构件则需布设抗压强度较高的加固材料，加固材料的布设应遵循锈蚀损伤rc结构荷载路径。
[0097]
由于本发明采用了以上技术方案，使本发明具备的有益效果在于：
[0098]
本发明提供一种基于荷载路径的锈蚀损伤混凝土结构加固设计方法，考虑了锈蚀损伤下材料性能损失与钢筋粘结退化的影响，提出了一种最小密度单元刚度自动调节方
法，给出了锈蚀损伤单元应变能形式的设计灵敏度公式，可克服由最小密度单元畸变引起的数值失稳问题，能合理地生成损伤结构最优荷载路径，可指导加固材料布设。本发明生成的锈蚀损伤rc结构荷载路径可为描述锈蚀损伤rc结构d区传力机理、评估结构性能与指导维修加固决策提供了理论指导。
附图说明
[0099]
为了更清楚地说明本发明实施方式或现有技术中的技术方案，下面将对实施方式或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施方式，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图示出的结构获得其他的附图。
[0100]
图1为本发明中混凝土cdp模型的应力-应变曲线；
[0101]
图2为本发明中锈蚀钢筋双线性本构模型示意图；
[0102]
图3为本发明中基于sigmoid函数的刚度调节示意图；
[0103]
图4为本发明中锈蚀损伤rc结构荷载路径搜寻流程图；
[0104]
图5为本发明中钢筋混凝土梁优化模型示意图；
[0105]
图6为本发明中锈蚀率η＝5％的rc梁荷载路径；
[0106]
图7为本发明中不同荷载工况下rc梁的荷载路径：(a)beso方法；(b)f＝5n；(c)f＝50n；(d)f＝100n；(e)f＝150n；(f)f＝200n；(g)f＝250n和(h)f＝300n；
[0107]
图8为本发明中外部载荷为f＝150n时不同锈蚀率下锈蚀损伤rc梁的荷载路径：(a)η＝0％；(b)η＝1.5％；(c)η＝3％；(d)η＝5％；(e)η＝10％；(f)η＝15％和(g)η＝20％。
具体实施方式
[0108]
下面将结合本发明附图，对本发明实施技术方案进行清楚、完整地描述，显然，所描述的实施方式仅仅是本发明的一部分实施方式，而不是全部的实施方式。基于本发明中的实施方式，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施方式，都属于本发明保护的范围。
[0109]
本发明各个实施方式之间的技术方案可以相互结合，但是必须是以本领域普通技术人员能够实现为基础，当技术方案的结合出现相互矛盾或无法实现时应当认为这种技术方案的结合不存在，也不在本发明要求的保护范围之内。
[0110]
参照图1-图8，本发明提供一种基于荷载路径的锈蚀损伤混凝土结构加固设计方法，其步骤包括：
[0111]
s1：确定锈蚀损伤rc结构荷载路径搜寻问题非线性拓扑优化公式；
[0112]
锈蚀损伤rc结构荷载路径搜寻问题的优化目标是在体积约束v
*
下生成结构柔度最小(即刚度最大)的荷载路径。将混凝土单元材料密度作为设计变量，锈蚀损伤rc结构非线性优化问题的数学表达式为：
[0113]
[0114][0115]
其中，为目标函数，也就是结构柔度；为单元设计变量的向量；为第i个单元的材料密度，在双向渐进结构优化方法(beso方法)中可以分别取为1或者f
ext
为外部荷载向量；u
η
为结构位移向量；vi和v
*
分别为第i个单元的体积和规定材料的体积分数。n为设计域中的混凝土单元数量；为单元最小密度；为非线性锈蚀损伤rc结构的残余力向量。
[0116]
s2：计算锈蚀材料性能、粘结强度与弹簧刚度；
[0117]
1)钢筋锈蚀会降低混凝土的抗压强度，假设纵向钢筋受均匀锈蚀，锈蚀下混凝土的抗压强度f
′
cm
可以表示为：
[0118][0119]
其中，f
cm
为未锈蚀混凝土的抗压强度；kr为经验系数；ε
cu
为最大压缩应变；ε1为开裂混凝土的平均受拉应变，可表示：
[0120]
ε1＝n
bs
·wcr
/bwꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(4)
[0121]
其中，bw为未锈蚀损伤rc梁的初始宽度；n
bs
为rc梁受压区钢筋数量；w
cr
为锈蚀损伤rc梁总裂缝宽度，可计算为：
[0122]wcr
＝2πx(u
rs-1)
ꢀꢀꢀꢀ
(5)
[0123]
其中，u
rs
为体积膨胀率；x为钢筋锈蚀深度，可表示为：
[0124][0125]
其中，η和rs分别为钢筋锈蚀率和钢筋半径。
[0126]
2)锈蚀钢筋极限抗拉强度可表示为：
[0127]fy，c
＝(1-αyη)fyꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(7)
[0128]fu，c
＝(1-αuη)fuꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(8)
[0129]
ε
u，c
＝(1-α1η)εuꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(9)
[0130]
其中，f
y，c
，f
u，c
与ε
u，c
分别为锈蚀钢筋的屈服强度，极限抗拉强度和极限应变；fy，fu与εu分别为未锈蚀钢筋的屈服强度、极限强度和极限应变；αy，αu与α1为相关经验系数。
[0131]
3)锈蚀钢筋与混凝土间的粘结滑移行为可以采用弹簧单元进行模拟，弹簧单元的刚度可由钢筋的粘结强度和相应的滑移量来确定。每对重合钢筋与混凝土节点之间采用长度为零的弹簧单元连接来模拟两者之间的粘结滑移行为。锈蚀钢筋与混凝土间粘结强度τ
a，η
可表示为：
[0132][0133]
其中，cc为混凝土保护层厚度；d
l
为纵向钢筋的直径；a
s，η
为锈蚀箍筋截面积；λ，ζ，k为经验系数；s是箍筋间距；混凝土的抗拉强度f
ct
与其抗压强度f
cm
有关，可以表示为有关，可以表示为f
sy，η
为锈蚀箍筋的屈服强度，可表示为：
[0134][0135]
其中，ηs为箍筋平均锈蚀率；f
sy
为未锈蚀箍筋的屈服强度。弹簧单元的刚度可表示为：
[0136]ksp
＝π
·d·
lr·
τ
a，η
/s1ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(12)
[0137]
其中，lr为混凝土单元的长度；s1为钢筋最大粘结应力处对应的最大滑移值。
[0138]
s3：建立锈蚀损伤混凝土结构有限元模型
[0139]
1)建立几何模型：根据实际工程结构，确定混凝土结构中各构件的几何尺寸，在本实施例中采用abaqus软件part模块构建混凝土、钢筋与垫块等构件的几何模型，采用assembly模块对各构件进行组装。
[0140]
2)定义锈蚀材料本构模型：混凝土受拉和受压时的非对称响应采用abaqus中混凝土损伤塑性模型(concrete damaged plasticity)模拟。采用property模块对锈蚀混凝土与钢筋材料进行定义。采用中国规范gb50010中的混凝土本构关系对混凝土硬化与软化行为进行模拟。混凝土受拉和受压力学响应，如图1所示。采用双线性本构模型来表征锈蚀钢筋的力学性能，如图2所示。
[0141]
3)定义网格划分与加载方式：采用mesh模块对混凝土与钢筋等材料进行网格划分，尽量采用规则四边形单元，钢筋与混凝土单元尺寸尽量保持一致。采用四节点cps4r单元和两节点线性二维桁架t2d2单元分别对混凝土与锈蚀钢筋进行数值模拟。
[0142]
4)采用load模块对混凝土结构定义边界条件。集中力和约束施加于刚性垫板上来防止局部混凝土单元过度变形，边界条件可根据工程实际确定。
[0143]
5)设置混凝土与锈蚀钢筋间的弹簧单元：采用interaction模块中special工具及其中的springs/dashpots工具，在锈蚀钢筋与混凝土单元的重合节点间采用两根长度为零的弹簧单元连接。在两根弹簧中，竖向弹簧变形相对于纵向的变形可忽略不计，因此，可以将竖向弹簧的刚度系数k取为无穷大。钢筋与混凝土间的粘结滑移主要通过纵向弹簧模拟，弹簧单元刚度由式(12)中k
sp
确定。假定箍筋和上部纵向钢筋分别与周围的混凝土完全粘结。采用interaction模块中嵌入约束(embedded equation)将混凝土单元节点分别与上部纵向钢筋和箍筋单元节点进行连接。采用绑定约束(tie)将刚性钢板与设计域的接触面进行连接。
[0144]
6)设置求解器并定义有限元输出：在step模块中定义有限元模块的输出，即单元应力、应变、弹性与塑性应变能等。为提高计算效率，在job模块中采用多个处理器与gpu加速，其他设置保持默认。
[0145]
s4：对锈蚀损伤rc结构进行拓扑优化分析，具体步骤如下：
[0146]
1)定义初始优化参数：过滤半径r
min
，惩罚因子p，目标体积分数v
*
，进化率er，初始应变阈值ε
*
，材料最小密度弹性模量阈值e
max
和初始弹性模量
[0147]
2)初始化设计变量与混凝土截面属性：设置“实体”与“空心”材料截面属性，“实体”材料即为混凝土材料；“空心”材料为一种弹性材料，其泊松比为0.3，弹性模量为其中ec为混凝土弹性模量。将所有的初始设计变量设置为设计域内所有单元截面属性设置为“实体”材料截面属性。
[0148]
3)基于abaqus进行非线性有限元分析。提取该优化迭代步中各单元的弹、塑性应
变能，单元i的总应变能wi可表示为：
[0149][0150]
其中，与分别为各单元的弹、塑性应变能。
[0151]
4)计算锈蚀损伤rc结构单元灵敏度，对结构设计域中单元灵敏度进行计算；
[0152]
锈蚀损伤rc结构灵敏度公式可表示为：
[0153][0154]
其中，λo是一个非常小的常数，可以设为λo＝10-10
。λ加到右边第二项的分母上，以避免在时出现歧义。
[0155]
5)更新单元灵敏度信息：采用模糊滤波方案得到网格无关性的拓扑优化解，对原始单元灵敏度进行进一步处理，单元灵敏度可表示为：
[0156][0157]
w(r
ij
)＝max(0，r
min-r
ij
)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(27)
[0158]
其中，r
ij
为单元i和j中心之间的距离；w为平均原始灵敏度的权重函数；ηj为权重因子；r
min
为过滤半径；当前迭代步的灵敏度与前一个迭代步的灵敏度进行平均；单元灵敏度可以表示为：
[0159][0160]
6)确定下一个优化步的目标体积：根据当前迭代步的体积vk和进化比er，可以得到下一迭代步的目标体积v
k 1
。
[0161]vk 1
＝vk(1
±
er)
ꢀꢀꢀꢀ
(29)
[0162]
7)更新设计变量和单元类型：如果混凝土单元的灵敏度小于灵敏度阈值则单元密度由1切换为对于弱弹性体单元，如果单元密度从切换为1。再根据单元密度对混凝土和弱弹性材料截面属性进行重新分配。
[0163]
s5：对最小密度单元刚度自动调节：首先提取最小密度单元的最大应变然后，判断该迭代步中最大应变与应变阈值ε
*
的大小。如果使用刚度多比例增长策略增大下一个迭代步中最小密度单元的弹性模量以抑制荷载路径搜寻过程中最小密度单元畸变引起的数值失稳问题；如果则减小下一个迭代步中的以抑制因单元刚度增加而引起的优化误差，见公式(23)。
[0164]
s6：计算锈蚀损伤rc梁荷载路径搜寻过程中的收敛误差erk。重复优化迭代过程，直到满足目标体积和收敛准则为止，引入收敛准则来终止荷载路径搜寻过程，该收敛准则为：
[0165][0166]
其中，∈为允许的收敛误差；n为迭代步数量。当结构柔度保持稳定时，该优化过程终止，生成锈蚀损伤rc结构最优荷载路径。
[0167]
s7：验证锈蚀损伤rc结构荷载路径的合理性。通过锈蚀损伤rc梁的数值算例阐明所生成的荷载路径合理性。
[0168]
s8：基于锈蚀损伤rc结构荷载路径指导加固材料布设。理想的加固材料布设应遵循结构荷载传递路径，充分发挥加固材料性能，提高其利用效率。基于损伤结构的荷载路径分析其中每个构件受力特性，对于受拉构件则在相应位置布设抗拉强度较高的加固材料，对于受压构件则需布设抗压强度较高的加固材料。加固材料的布设应遵循锈蚀损伤rc结构荷载路径。
[0169]
为了进一步的说明本发明提供的一种锈蚀损伤混凝土结构荷载路径拓扑搜寻方法的操作步骤，下面对如附图5所示的锈蚀损伤rc梁进行荷载路径搜寻与加固材料布设研究，试件的几何尺寸为600
×
100mm，垫块尺寸为20
×
10mm，相关rc梁几何参数如表1所示，混凝土和钢筋材料的力学性能如表2所示。
[0170]
表1 rc梁几何参数
[0171][0172]
表2混凝土与钢筋力学性能
[0173][0174]
本算例的基于荷载路径的锈蚀损伤混凝土结构加固设计方法具体过程如下：
[0175]
s1：确定锈蚀损伤rc梁荷载路径搜寻问题非线性拓扑优化公式；
[0176]
锈蚀损伤rc梁荷载路径搜寻问题的优化目标是在体积约束v
*
下生成结构柔度最小的荷载路径，其数学表达式为：
[0177][0178][0179]
s2：计算锈蚀材料性能、粘结强度与弹簧刚度。锈蚀损伤rc梁中钢筋锈蚀率为5％，根据文献“failure analysis of corroded rc beams subjected to shear-flexural actions，antonino recupero，nino spinella，francesco tondolo，engineering failure analysis，2018，93：26-37”可知：经验系数kr＝0.1，体积膨胀率u
rs
＝2，受压区钢筋数量n
bs
＝2，经验系数αy，αu与α1可以分别取为：0.012，0.011与0.03；基于公式(3)-(9)可求得：钢筋锈蚀深度x为0.304mm，总裂缝宽度w
cr
为1.909mm，平均受拉应变ε1应小于等于混凝土极限应变ε
s，u
，ε1可取为0.0033；因此，锈蚀下混凝土的抗压强度f
′
cm
为36.36mpa，锈蚀钢筋屈服强度为367mpa，锈蚀钢筋极限强度为548.1mpa，锈蚀钢筋极限应变为0.425。
[0180]
根据文献“shear and flexural strength prediction of corroded r.c.beams，g.campione，f.cannella，l.cavaleri，construction and building materials，2017，149：395-405”可知：经验系数ζ，k，λ分别可取为0.1，0.16，0.4，保护层厚度cc为20mm，箍筋间距s为50mm。混凝土单元的长度lr为5mm，最大滑移值s1为1mm。基于公式(10)-(12)可求得：锈蚀混凝土抗拉强度f
ct
为3.54mpa，箍筋屈服强度f
sy，η
为300mpa，箍筋截面积a
s，η
为50.24mm2，锈蚀钢筋与混凝土间粘结强度τ
a，η
为3.83mpa，弹簧单元刚度k
sp
为1443.14kn/mm。
[0181]
s3：建立锈蚀损伤混凝土结构有限元模型。
[0182]
1)建立几何模型。根据rc梁尺寸与钢筋尺寸分别创建混凝土、钢筋与垫块部件。采用assembly模块对各部件进行组装。
[0183]
2)定义锈蚀材料本构模型。采用property模块＞mechanical＞plasticity＞concrete damage plasticity对混凝土与钢筋材料进行定义。cdp模型中参数膨胀角ψ＝36
°
、流动势偏心∈＝0.1、双轴抗压强度与单轴抗压强度之比静水压力轴线与压缩和拉伸子午线之间距离的比值k
ζ
＝0.6667和黏度参数μ＝0.0005。具体有关混凝土cdp模型使用说明参照abaqus说明书。
[0184]
3)定义网格划分与加载方式。采用四边形单元网格将混凝土设计域离散为120
×
20个网格，混凝土单元类型为cps4r，单元尺寸为5mm。钢筋材料单元类型为t2d2，单元长度为5mm。采用load模块中concentrated force选项将集中荷载f＝150n竖向施加于rc梁中部刚性垫板节点上，将左端刚性垫板节点位移与旋转displacement/rotation选项全部约束，将右端刚性垫板节点displacement/rotation选项竖向位移约束。
[0185]
4)设置混凝土与锈蚀钢筋间的弹簧单元。采用interaction模块＞special＞springs/dashpots工具，在锈蚀钢筋与混凝土单元的重合节点间采用两根长度为零的弹簧单元连接，设置弹簧单元刚度k
sp
为1443.14kn/mm。采用interaction模块中嵌入约束(embedded equation)将钢筋部件嵌入混凝土部件。采用绑定约束(tie)将刚性钢板与设计域的接触面的单元节点进行连接。
[0186]
5)设置求解器并定义有限元输出。在step模块＞static，general定义分析步，在此案例中不考虑几何非线性，关闭nlgeom选项，在incrementation中设置最大迭代步maximum number of increments为10000，初始荷载步大小initial为0.001，最小荷载步大小minimum为10-8
，最大荷载步大小maximum为1，打开非对称求解unsymmetric选项。定义有限元输出create field output与create history output，打开stresses中s，stress components and invariants、strains，e，total strain components、energy中ener，all energy magnitudes与energy中allen，all energy totals选项。在job模块中创建工作步，打开parallelization选项，设置use multiple processors为8，设置use gpgpu acceleration为6，其他设置保持默认，锈蚀损伤rc梁优化模型如图5所示。
[0187]
s4：对锈蚀损伤rc梁进行拓扑优化分析。采用beso方法对锈蚀损伤rc梁进行拓扑优化，有关beso方法的python程序，可参考文献“a simple and compact python code for complex 3d topology optimization，zhi hao zuo，yi min xie，advances in engineering software，2015，85：1-11”。此外，最小密度单元刚度自动调节方法采用
python语言编写。
[0188]
1)定义初始优化参数：过滤半径r
min
为20mm，惩罚因子p为3，目标体积分数v
*
为50％，进化率er为0.02，初始应变阈值ε
*
为0.0015，材料最小密度为0.001，弹性模量阈值e
max
为10mpa，初始弹性模量为10-4
ec。
[0189]
2)初始化设计变量与混凝土截面属性。设置“实体”与“空心”材料截面属性，“实体”材料即为混凝土材料；“空心”材料为一种弹性材料，其泊松比为0.3，弹性模量为将所有的设计变量设置为设计域内所有单元截面属性设置为“实体”材料截面属性。
[0190]
3)对锈蚀损伤rc梁进行非线性有限元分析。分别提取每个优化迭代步中各单元的弹、塑性应变能，单元i的总应变能wi可由公式(25)计算所得。
[0191]
4)基于公式(20)计算锈蚀损伤rc结构各单元灵敏度。
[0192]
5)基于公式(26)-(28)更新单元灵敏度信息。
[0193]
6)基于公式(29)确定下一个优化步的目标体积。
[0194]
7)更新设计变量和单元类型。采用二分法确定该迭代步灵敏度阈值如果混凝土单元的灵敏度小于灵敏度阈值则单元密度由1切换为对于弱弹性体单元，如果单元密度从切换为1。再根据单元密度对混凝土和弱弹性材料的截面属性进行重新分配。
[0195]
s5：最小密度单元刚度自动调节。单元刚度自动调节方法中系数α1、α2、α3和α4可分别取为1.2、0.5、1和0.9。首先，提取该迭代步中最小密度单元的最大应变然后，判断最大应变与应变阈值ε
*
的大小；如果基于公式(23)增大最小密度单元的弹性模量如果则基于公式(23)减小最小密度单元的弹性模量
[0196]
s6：基于公式(30)计算锈蚀损伤rc梁荷载路径搜寻过程中的收敛准则erk。根据文献“bi-directional evolutionary topology optimization of continuum structures with one or multiple materials，huang x，xie ym，computational mechanics，2009，43(3)：393-401.”可知：收敛误差∈为0.01％，迭代步数量为n＝5。重复优化迭代过程，直到满足目标体积和收敛准则为止。当结构柔度在至少连续10次迭代步中保持稳定时，该优化过程终止，生成锈蚀损伤rc结构最优荷载路径。
[0197]
s7：验证锈蚀损伤rc结构最优荷载路径的合理性。基于发明方法生成了外部载荷为f＝150n且锈蚀率η＝5％的rc梁的荷载路径，如图6所示。荷载路径从加载点与支撑点向梁1/4处延伸，形成了两条几乎平行的荷载路径轨迹，体现了rc梁弯剪段的剪切行为，可合理地描述锈蚀损伤rc梁的剪力传递机理。
[0198]
进一步验证本发明分析方法的合理性与正确性，采用本发明方法分别研究了不同荷载与锈蚀率下锈蚀损伤rc梁的荷载路径发展规律。图7给出了不同荷载下(beso方法；f＝5n；f＝50n；f＝100n；f＝150n；f＝200n：f＝250n和f＝300n)rc梁的荷载路径发展规律。图7(a)为基于beso方法的弹性rc梁荷载路径，在最小荷载f＝5n的情况下，rc梁的荷载路径与已有文献中弹性的荷载路径几乎是一致的。当载荷f增加到100n时，rc梁的载荷路径开始发生变化，荷载路径从加载点直接延伸到支撑点，与已有文献中经典的拉压模型内力轨迹一
致，如图7(c)所示。当外部荷载进一步增大时，支撑点附近的荷载路径末端逐渐向rc梁跨中移动，如图7(c)-(f)所示。此后，随着外部荷载进一步增大至接近极限状态，荷载路径从加载点与支撑点向梁1/4处延伸，形成了两条几乎平行的荷载路径轨迹，体现了rc梁的剪切行为，如图7(g)-(h)所示。如上所述，本发明方法有效地解决了非线性拓扑优化过程中数值失稳问题，成功地生成了rc梁的荷载路径。可基于损伤混凝土结构荷载路径指导加固材料的布设。
[0199]
图8给出了不同锈蚀率下(0～20％)锈蚀损伤rc梁在外部荷载f＝150n时的荷载路径变化。锈蚀损伤rc梁的荷载路径随着锈蚀率的增加逐渐发生变化。不同锈蚀水平与不同荷载工况下锈蚀损伤rc梁的荷载路径的变化规律一致。锈蚀会加速rc梁的荷载路径由非极限状态向极限状态转变。
[0200]
s8：基于锈蚀损伤rc结构荷载路径指导加固材料布设。基于损伤结构的荷载路径分析其中每个构件受力特性，对于受拉构件则在相应位置布设抗拉强度较高的加固材料，对于受压构件则需布设抗压强度较高的加固材料。加固材料的布设应遵循锈蚀损伤rc结构荷载路径。
[0201]
以上仅为本发明的优选实施例，并非因此限制本发明的专利范围，凡是在本发明的发明构思下，利用本发明说明书及附图内容所作的等效结构变换，或直接/间接运用在其他相关的技术领域均包括在本发明的专利保护范围。