模拟退火算法在悬臂浇筑拱桥索力计算中的应用研究

孙波， 王祺顺， 田仲初， 成魁

(1.湖南交通国际经济工程合作有限公司， 湖南 长沙 410005； 2.长沙理工大学 土木工程学院； 3.湖南省交通科学研究院有限公司)

悬臂浇筑拱桥施工一般采用斜拉扣挂系统，以扣锚索作为拱圈节段临时约束，通过调整扣锚索索力值以确保主拱圈截面上下缘拉应力不超限，因此该类桥梁施工过程中索力值调整成了整桥施工及监测的重点与难点。目前，中国多名学者采用不同的方法对悬浇拱索力优化进行了研究，如应力平衡法、零弯矩法、索力可行域法等，均取得了较好效果，但是以上算法均需要对拱圈节段进行正装迭代验算，计算繁琐，耗费机时，且影响矩阵庞大，随着跨径增大、拱圈节段数量增多，其弊端更将凸显。另有学者提出运用多目标算法进行索力调整，但是多目标优化方法参数多，且某些参数之间的相互关系很难量化确定，在实际运用中存在一定困难，因此只能简化处理，影响了计算结果的精度及准确性。该文在借鉴前人研究成果的基础上，基于模拟退火算法原理，提出一种索力优化调整的新方法，研究成果可为悬臂浇筑拱桥索力计算提供一种新的思路。

1 模拟退火算法原理

模拟退火算法源于固体退火原理，是一种基于蒙特卡洛法的迭代随机寻优算法，其出发点是固体物质的退火过程与一般组合优化的一致性。模拟退火算法从某一较高初温出发，伴随温度参数的不断下降,结合概率突跳特性在解空间中随机寻找目标函数的全局最优解，即在局部最优解能概率性地跳出并最终趋于全局最优，该法可有效避免陷入局部最小死循环，是一种有效全局搜素最优化解法。

模拟退火算法主要收敛判断准则为Metropolis准则，其一般数学表达式为：

(1)

式中：exp为趋于平衡状态的概率；其中E为达到某温度时的内能；ΔE为其改变量。

将p与[0,1]之间的一个随机数α作比较，若p<α，则接受新状态xnew为当前状态,取代原状态xold，否则舍弃。不断重复上述过程，产生新状态，直到系统达到能量最低的平衡状态为止，固体状态的分布概率趋于吉布斯正则分布。

2 工程概况

该文以某在建悬臂浇筑拱桥为研究对象，该桥为单箱双室钢筋混凝土拱桥，主拱圈净跨径为240 m，矢跨比为1/6，拱轴系数为1.85，采用斜拉扣挂悬臂浇筑施工工艺，沿拱圈方向共分为37个节段，其中包括34个悬臂浇筑节段、2个拱脚支架施工节段及1个合龙段，临时扣塔采用Q345a钢材，扣锚索采用预应力钢绞线。桥型布置图及斜拉扣挂系统示意图如图1、2所示。

图1 桥型布置图(单位：cm)

图2 斜拉扣挂系统示意图

3 优化模型建立

3.1 参数化有限元模型建立

为便于后续优化程序的实现，首先需建立参数化有限元模型，该文使用Ansys APDL建立全桥实体、梁、杆空间混合有限元模型命令流文件，其中混凝土拱圈、交界墩使用Solid65 8节点混凝土实体单元模拟，扣塔主要承重钢管、横撑、斜撑、腹杆及竖杆使用Beam189梁单元模拟，两拱脚及交界墩底采用固结约束方式，锚索在锚锭处模拟为固结，扣塔立柱与横向联系间设置虚拟刚性梁。扣索扣点与主拱圈实体单元共节点。利用Ansys生死单元技术模拟主拱圈实际施工阶段。进行各施工阶段计算分析前，提前在前处理器中生成所有单元。进行施工阶段分析时，先用EKILL命令“杀死”所有单元，某施工阶段达到强度或安装完成的单元用“ELIVE”命令激活。采用集中力和集中弯矩模拟拱圈节段混凝土湿重及挂篮模板重的施加，因主拱圈截面为实体单元，不便施加集中力和集中弯矩，在每个节段浇筑分界面处设置质量节点MASS21作为加载节点，然后将节段浇筑分界面连同质量节点建立刚性约束，将集中力和集中弯矩施加于加载节点。

为简化计算，该文基于线形徐变理论对悬臂浇筑过程中混凝土徐变效应进行分析，在Ansys中选取6号徐变方程作为徐变准则，同时采用应变强化准则以考虑应力随时间的变化情况，不考虑温度对徐变的影响。同时分节段激活拱圈各节段，施加索力及其他外荷载，进入PRESS模块进行结构静力分析，在每个时间步长内将材料徐变参数重新赋值更新，进行正装迭代求解，依次类推，直至全桥合龙。

为保证计算结果精度满足要求，采用映射分网技术，同时拱圈实体使用六面扫略分网，单元尺寸控制在50 cm，为简化模型，挂篮荷载用等效节点荷载代替。

3.2 优化模型建立

模拟退火算法优化求解程序分为两层：第1层为在初始解附件随机生成某解，与目标值进行比较，不断全局搜索最终被接受；第2层为设置某较大初始解，通过多次迭代后缓慢降温最终收敛，具体实施步骤如下：

(1) 设立目标函数，随机设置初始解。

(2) 设置迭代计数器。

(3) 在初始解邻域范围内随机产生某扰动，生成新的模型参数值m1=m0+Δm，并计算新的目标函数值E(m1)。

(4) 对比两种情况下的目标函数值之差ΔE，若ΔE≤0，则接受新的模型参数值，否则按照Metropolis准则以一定的概率判断是否接受新值。

(5) 若计算迭代次数t小于设定的最大迭代次数，则t=t+1，转到步骤(3)，继续计算迭代过程。

(6) 若已到达设定的最大迭代次数，则判断是否达到温度终止条件，若未达到冷却状态，则按指定的降温方式冷却，并转到(3)；若达到终止条件，则停止算法，输出的当前模型参数值为全局最优解。

根据以上求解步骤，设置以下算法参数：

(1) 设计参数变量：该文选取扣索初拉力值作为设计修正变量。

(2) 状态变量：对于采用预应力钢绞线的扣索，采用规范要求2.5倍安全系数下的索力值为状态变量之一，即扣索峰值应力应小于1 860/2.5=744 MPa，同时，该桥主拱圈采用C60高强混凝土，在施工过程中各节段各截面顶底板最大拉应力峰值均应小于C60混凝土抗拉强度设计值2.04 MPa。

(3) 目标函数：选取调整后索力作用下拱圈各节段最大拉应力与设计一次成拱状态下拱圈各节段最大拉应力之差作为目标函数，其差的平方和可表述为：

(2)

模拟退火算法优化程序需设定一个较高的初始温度(即初始解)方能获得较为精确的全局最优解，但初始温度设置过高，迭代次数会大幅增加，根据Metropolis判断规则，取初始温度为100 ℃，同时为避免程序陷入局部最优死循环，根据Kirkpatrick几何温度衰减函数，温度衰减系数取0.9。

(4) 终止规则：考虑到退火算法迭代过程中温度呈指数函数的形式降低，最终无限趋于0，若不设置终止温度，程序将无限迭代陷入死循环。参考以往研究经验，将终止温度设为0.000 1 ℃。

4 修正后结果及对比

考虑到该有限元模型规模庞大，节点单元数较多，为提高计算效率，在Matlab中求解得到模拟退火算法数学解，将其作为索力初始值，以实常数的形式施加于Ansys有限元参数模型中，同时，调用Batch批处理，对施工阶段进行多点启动求解，在后处理模块中查看有限元结果，如此反复，直至收敛。该法将Ansys作为求解器使用，同时借助于Matlab强大的矩阵运算能力，可极大缩短机时。

4.1 索力修正前后对比

通过Ansys及Matlab联合求解，调取程序收敛于全局最优解后的设计参数修正值如表1所示。

表1 扣索力修正结果

续表1

表1计算结果表明：

(1) 原设计索力偏小，由模拟退火算法优化调整后的各节段扣索索力值均有不同幅度增长，其中增幅最大为东西岸1#节段扣索索力值，增量最大为东西岸16#扣索。

(2) 扣索索力调整前后结果表明：东西岸14#～18#拱圈节段索力增量普遍较大，分析原因为：随着拱圈悬臂浇筑过程的进行，拱圈节段与水平线夹角呈减小趋势，即扣索的水平倾角逐渐减小，故作用于节段上的竖直分力减小，而在斜拉扣挂体系中，起主要作用的为扣索的竖直分力。为满足状态变量的约束条件，大编号的扣索索力增量均较大。

(3) 修正后的索力与设计索力变化趋势相同，均为拱脚和拱顶较小，拱腰位置最大，索力修正未改变索力矩阵的排列方式。

4.2 索力修正后拱圈应力结果

拱圈截面应力值为悬臂浇筑拱桥施工过程中结构安全重要控制性指标，关系到拱圈合龙后应力分布。将修正后的索力导入Ansys APDL参数化模型中进行正装迭代计算，提取修正索力下拱圈截面应力结果。

计算结果表明：拱圈节段顶板应力峰值出现位置及工况为下一节段浇筑时相邻两节段拱圈连接处位置，底板应力峰值出现工况为该阶段扣索张拉时靠近拱圈端头位置。提取各节段该位置截面在各施工阶段下的峰值应力如表2所示。

表2表明：

(1) 在原设计索力下，由于索力偏小，在拱圈节段自重及其他外荷载作用下，拱圈节段在悬臂浇筑过程中拉应力水平普遍较高，其峰值达到4.73 MPa，超出C60混凝土抗拉强度设计值，混凝土存在开裂现象，而在扣锚索张拉时，底板拉应力始终位于较低水平。较小的扣索力未充分挖掘截面的承载能力，顶底板应力差较大，应力分布不均，施工过程中存在一定风险性。

表2 拱圈截面峰值应力 MPa

(2) 索力修正后，拱圈顶底板应力幅及应力分布均发生一定变化，拱圈在浇筑过程中最大峰值拉应力为1.72 MPa，底板拉应力峰值为1.62 MPa，较原设计索力下拉应力具有大幅下降，且均满足C60混凝土抗拉强度设计限值。

(3) 索力修正后，改善了之前应力分布不均的现象，顶底板应力在整个施工过程中交替变化且均匀，有效利用了截面强度，发挥了材料性能。

4.3 施工过程中扣塔局部分析

索力修正后，张拉力较原设计普遍较大，为了保持扣塔偏位满足规范要求，相应锚索力不可避免地随之增大。扣塔虽为临时性结构，但是在施工过程中决定了结构的安全及稳定，增大的扣锚索索力值对扣塔受力有不利影响，为确保扣塔安全性，提取各施工阶段最不利工况下扣塔应力计算结果。

扣塔立柱及横撑斜撑在各施工阶段最不利工况下应力计算结果如图3所示。

图3 立柱及横撑Von Mises应力云图(单位：MPa)

图3表明：在索力修正后，最不利工况下，扣塔立柱Von Mises应力为180 MPa，横撑及斜撑等局部构件Von Mises应力为177 MPa，扣塔除锚固点有轻微应力集中现象外，其他大部分区域应力分布较为均匀，Von Mises应力值小于屈服应力限值。

5 结论

以某悬臂浇筑拱桥为研究背景，基于模拟退火全局最优搜索法，求解了全局最优索力解，并建立该桥三维实体、梁、杆空间混合有限元模型，对求解的索力进行验证，得到以下结论：

(1) 基于模拟退火算法，设置一定的约束条件，可有效对扣索力进行全局最优修正，修正后的索力较设计索力均有不同幅度的增长，最大索力增幅为37%，最大增量为563.8 kN。

(2) 调索后拱圈应力有明显改善，顶板最大拉应力由4.73 MPa降至1.72 MPa，底板最大拉应力由0.54 MPa增至1.62 MPa，索力调整后，拱圈顶板应力大幅下降，底板拉应力有一定上升，保证了截面上应力流随拱圈混凝土浇筑-扣索张拉循环工况下缓和变化，应力分布均匀，改善了顶底板应力差值大的情况。

(3) 索力调整后，扣塔最大Von Mises应力值为180 MPa，且大部分区域应力分布流畅，应力无明显跳点及集中现象。