基于混合萤火虫-模拟退火算法的斜拉桥纵向振动控制与黏滞阻尼器优化设计

2020-12-28龙关旭辛公锋徐兴伟任士朴苏文明

公路交通科技 2020年12期

龙关旭，辛公锋，徐兴伟，任士朴，苏文明

(1.山东高速工程检测有限公司，山东济南 250002；2.山东华鉴工程检测有限公司，山东济南 250100；3.桥梁结构大数据与性能诊治提升重点实验室，山东济南 250002)

0 引言

大跨径桥梁结构轻柔、阻尼较低，对地震、风、车等外部荷载引起的振动更为敏感，影响着桥梁的结构运营安全和车辆驾驶舒适度。基于被动控制理论而设置的黏滞阻尼器技术发展趋于成熟，在不改变桥梁刚度的前提下可有效提高结构阻尼，减小结构振动响应，并对伸缩缝起到一定的保护作用，已被广泛应用到桥梁结构振动控制中[1-2]。

20世纪90年代国内外就已通过模型试验与工程应用对黏滞阻尼器的减振消能作用展开研究[3-4]。叶爱君[5]研究得出了黏滞阻尼器对斜拉桥在地震作用下的结构响应有着良好的控制作用，可明显减小桥梁纵向位移和塔底弯矩。焦驰宇等[6]研究了黏滞阻尼器在市政桥梁加固上的应用，得出了该加固措施在工程造价、交通影响和可实施性等方面都有着明显优势。丁幼亮等[7]通过非线性时域分析，研究了黏滞阻尼器对多塔斜拉桥风致振动的控制效果，黏滞阻尼器的设置可使得各桥塔和主梁的风振响应幅值趋于一致。Guo[8]以润扬大桥为依托，得出了跑车工况下黏滞阻尼器可有效控制梁端纵向位移。同时，很多学者对阻尼器参数的优化做了相关研究[9-10]，王波等[11]基于斜拉桥的随机地震响应，采用非线性规划法对黏滞阻尼器关键参数进行了优化取值。目前，不论是对黏滞阻尼器的振动控制研究，还是对黏滞阻尼器关键参数的优化研究，都是考虑以地震为主的单一荷载作用，且参数优化多采用传统标准算法来实现。

为此，以某大跨斜拉桥为例，对地震、风、随机车流作用下黏滞阻尼器的振动控制进行研究，并基于标准萤火虫算法和模拟退火算法提出混合萤火虫-模拟退火算法(HFSA)，对黏滞阻尼器在3种荷载的共同作用下的关键参数进行优化设计。

1 黏滞阻尼器及斜拉桥模型的建立

黏滞阻尼器是一种无刚度、速度型阻尼器，通常由缸体、导杆、黏滞流体材料、活塞和阻尼孔等组成，构造如图1所示。

图1 黏滞阻尼器构造

黏滞阻尼器的力学模型理论计算公式为：

Fd=Cvαsgn(v)，

(1)

式中，Fd为阻尼力；C为阻尼系数；v为相对速度；α为速度指数，其范围为0.1～1。

以某大跨双塔斜拉桥为依托，跨径布置为 (100+100+300+1 088+300+100+100) m，主梁采用扁平流线形钢箱梁，索塔采用倒Y形混凝土塔，全线桥涵设计荷载采用汽车-超20级，挂车-120级。

图2 斜拉桥ANSYS模型

采用有限元软件ANSYS建立桥梁空间三维有限元模型，如图2所示。模型采用Beam4单元来模拟结构的主塔、主梁、盖板和墩柱等构件，斜拉索使用Link10单元进行模拟，利用Ernst公式进行弹性模量修正，过渡墩和辅助墩墩顶压重质量使用Mass21模拟，只考虑单元的质量不考虑刚度。全桥共870个节点，1 134个单元，将斜拉索单元设置初始应变参数，考虑大变形非线性求解，激活应力刚度效应。塔梁之间采用节点耦合的方式模拟支座约束，质量单元和主梁之间采用共节点的方式连接，索梁之间建立刚臂连接。梁端约束根据设计资料确定，墩底采用固结约束。

在两个主塔和主梁间分别设置两个纵向阻尼器，全桥共4个，阻尼器单元采用非线性弹簧单元(combin37)模拟，由两个活动节点(I，J)和两个可选控制节点(K，L)组成，如图3所示。在计算过程中，可以通过选取控制节点得出单元内部的相对位移对单元实常数进行设置，阻尼器单元与主梁之间的作用采用刚臂进行连接。

图3 combin37单元示意图

2 外荷载作用下纵向黏滞阻尼器振动控制

分别选用地震作用、风荷载、随机车流荷载进行荷载输入：根据桥址和《公路工程抗震规范》(JTG-B02—2013)[12]确定地震设计反应谱，进而作得3组加速度时程曲线，选用最大的一组进行纵向和竖向地震加速度输入，如图4所示。根据依托斜拉桥相关技术资料，采用高雷诺数主梁静三分力系数作为静风力输入，如图5所示。

采用谐波合成法生成三维脉动风场[13-14]，同时根据《公路桥梁抗风设计规范》(JTG/T D60-01—2004)[15]，分别采用Lumley谱和Simiu谱对竖向和横向脉动风速时程进行模拟输入，采用MATLAB软件编制脉动风模拟程序。图6给出了全桥关键风速模拟示意点，以1#节点为例，15 m/s风速时程曲线如图7所示。

图4 地面加速度时程曲线

图5 主梁及桥塔三分力系数

图6 关键风速模拟点

图7 1#节点脉动风速时程曲线

图8 单主梁模型车轮荷载分配示意图

根据《桥梁用黏滞流体阻尼器》(JTT 926—2014)[18]和已有黏滞阻尼器在大跨桥梁上的应用，选定速度指数α取值0.1～0.6和1，阻尼系数C取值1 000，2 500，5 000，7 000，10 000，15 000，20 000 kN/(m/s)α，共设定49组(α，C)参数组合，分别选用主梁跨中弯矩、梁端纵向位移、梁端纵向加速度和塔根弯矩进行结构响应分析。以地震作用为例，选用阻尼系数C=5 000 kN/(m/s)α，结构在地震作用下的主要响应时程曲线如图9所示，对比3种荷载引起的结构主要响应极大值变化，如图10所示。

图9 地震作用下主要响应时程曲线

图10 斜拉桥主要响应极大值变化

由图9和图10可看出，黏滞阻尼器的设置可有效减小结构振动响应，且基本随着速度指数的减小、阻尼系数的增大，振动控制效果更为明显。同时，地震引起的结构主要响应较为显著，随机车流次之，风荷载较小。

将纵向黏滞阻尼器的振动控制作用进行直观化，定义控制效率的概念，具体表达式为：

(2)

式中，θ为阻尼器的控制效率；A1为设置阻尼器后结构的响应；A2为无阻尼器的结构响应。由于黏滞阻尼器关键参数组合较多，表1只对部分参数组合进行了表述。

3 混合萤火虫-模拟退火优化算法

3.1 标准算法3.1.1 萤火虫算法

萤火虫算法(Firefly Algorithm，FA)是一种仿生智能优化算法，由剑桥学者Yang[19]受启发于萤火虫群体尾部亮光行为而首次提出。萤火虫算法模型包括两种关键因子：亮度和吸引度。亮度即目标函数，决定着萤火虫所处的位置优劣进而决定单只萤火虫的移动方向；吸引度决定萤火虫移动的距离，吸引度越大移动距离就越长。在迭代循环过程中，萤火虫亮度和吸引度不断更新信息，从而对目标值进行优化。

(1)萤火虫相对荧光亮度

(3)

式中，I0为萤火虫自身(r=0)的最大荧光亮度；γ为光吸收系数；rij为萤火虫i和j之间的距离。

(2)萤火虫吸引度

(4)

式中β0为光源处(r=0)最大吸引度，其余与式(3)含义相同。

(3)萤火虫移动位置更新

(5)

3.1.2模拟退火算法

模拟退火算法(Simulated Annealing， SA)是一种贪心算法，最早是由Metropolis等根据统计力学中的物质退火方法而提出，即固体加热至高温再慢慢冷却的过程[20]。

定义固体热能函数(目标函数)为f(x)，当前解为x，新解为x′，则热能增量为Δf=f(x′)-f(x)，根据Metropolis准则接受新解的概率为：

(6)

3.2 混合萤火虫-模拟退火算法

萤火虫算法对于全局优化能力更强，局部搜索能力较差，而模拟退火算法具有较强的局部最优，对全局最优搜索效率较低。可将萤火虫算法和模拟退火算法相互结合，取长补短，设计提出混合萤火虫-模拟退火算法(HFSA)。

该算法以标准萤火虫算法为总体框架，从空间内随机产生的初始解中进行全局最优搜索，然后对产生的较优全局解分别进行模拟退火操作，并更新到下一个萤火虫群体中。整个过程反复迭代更新判断，直到满足收敛条件为止。这样两种算法的优点充分结合了起来，大大提高了算法的效率，可迅速找出全局最优解，增强了算法跳出局部最优解的能力。具体流程如下。

表1 纵向黏滞阻尼器控制效率(单位：%)

(1)初始化算法参数

①萤火虫算法信息：种群规模(萤火虫个数m)，最大迭代次数(萤火虫进化次数Mmax)，相对吸引度β0，荧光吸收系数γ，步长因子δ，目标函数f(x)。

②模拟退火算法信息：初始温度T0，退火因子D，马尔科夫链长L，进行模拟退火算法的萤火虫个数k。

(2)空间随机初始化萤火虫群体，更新确定萤火虫初始位置状态xi(i=1, 2, 3,…,m)，计算萤火虫亮度Ii，即初始目标函数值。

(3)萤火虫种群更新迭代

①i=1，根据式(4)吸引度大小判断移动半径和式(5)更新萤火虫空间位置。

②根据式(3)更新计算萤火虫相对亮度。

③重新赋值i，j。

④判断i，j是否大于m，否则循环②～④。

(4)以收敛种群中较优的k个解模拟退火操作，进行局部最优搜索。

①迭代次数it=0，产生萤火虫个体新解x′k。

②计算目标函数差值Δfk=f(x′k)-f(xk)。

③根据式(6)按一定概率接受新解x′k。

④循环L次，然后结束局部更新搜索。

(5)进行退火冷却操作：it=it+1，温度T=DT。

(6)判断it>Mmax，算法结束，否则循环步骤(3)～(5)。

具体流程如图11所示。

图11 混合萤火虫-模拟退火算法流程

3.3 HFSA算法算例验证

选用4个标准测试函数来进行HFSA优化算法收敛性能和寻优效果的验证，并与标准萤火虫算法和模拟退火算法进行对比分析，具体参数取值如表2所示。

表2 测试函数参数

3.3.1Ackley函数

e0.5[cos(2πx)+cos(2πy)]+e。

2.2.2 概率性敏感性分析根据前文的假设，成本服从正态分布，转移概率服从β分布。成本数据的状态分布参考相关文献[17-18，24，27]报道的参数的 95%CI值，计算其标准差；转移概率的状态分布则根据其均值和标准差计算其α、β值。采用TreeAge Pro 2011软件对模型结果进行概率性敏感性分析，对模型进行次数为1 000的蒙特卡洛模拟，相关参数见表8。对模拟结果绘制成本-效用可接受曲线，详见图3。

(7)

3.3.2ShafferF6函数

(8)

3.3.3Sphere函数

(9)

3.3.4Rastrigrin函数

(10)

图12 测试函数示意图

3种算法采用相同的参数设置，如表3所示，分别运行20次，表4给出了4种算法对4个测试函数的优化结果汇总。

表3 算法主要参数

由表4可以看出，3种方法对二维函数优化都较好，迭代次数较少，但是FA算法最优值方差较大，SA和HFSA算法的鲁棒性较好；而对多维函数，SA算法体现了较差的全局优化性能，尤其对于多峰函数容易对局部产生过长搜索(5 000次迭代)。优化改进的HFSA算法，不论在单峰函数还是多峰函数都表现出了较好的鲁棒性和适应度，收敛速度更快，方差也较小。

表4 优化结果

4 纵向黏滞阻尼器关键参数优化

4.1 HFSA优化参数的确定

设计变量选为黏滞阻尼器的关键参数：速度指数α∈[0.1,0.6]和阻尼系数C∈[1 000, 20 000]。

以经济效益为核心，综合考虑地震-风-随机车流的外荷载作用下的阻尼力之和，构建目标函数：

Z=min(aFe+bFw+cFv)，

(11)

约束条件选定为斜拉桥设计规范中对结构变形的要求，即控制主梁跨中竖向位移，并同时考虑梁端纵向位移变化：

w≤l/500，

(12)

Δ≤50 cm，

(13)

式中，l为桥梁主跨跨径；w为斜拉桥跨中竖向挠度；Δ为梁端纵向位移。

4.2 基于响应面法的目标函数拟合

采用Design-Expert软件的BBD法进行试验设计，根据第2节的有限元分析计算结果进行输入，基于最小二乘法构建多项式系数，从而得到响应面模型函数。以地震作用下的最大阻尼力F为算例，然后采用响应面法拟合函数，拟合回归方程为：

F=6 551.19α2-3.43×10-6C2-4 360.93α-0.97C+553.92。

(14)

对拟合的目标函数进行精度检验，如表5所示。F统计值为3 493.67，说明模型具有较好的统计意义(大于95%的默认置信区间)。变异系数1.58<10，表明试验的可信度和精确度较高，相关系数R2为0.999 8，表明拟合的函数具有较高的精度；0.996 4的预测R2与0.999 5的调整值具有合理的一致性。205.197的精度度量表明了模型具有较高的精密度。

表5 精度检验指标值

4.3 权重分配及工况设置

考虑到纵向黏滞阻尼器的设计初衷，以地震作用的权重系数为主因素，具体加权系数取值如表6所示。

考虑4种地震动加速度峰值的抗震设防烈度区域、3种平均风速荷载和1种随机车流工况，分别进行纵向黏滞阻尼器的优化设计研究。具体荷载工况设置如表7所示。

表6 加权系数组合

表7 荷载工况组合

4.4 抗震设防烈度影响

选定风速为10 m/s，利用HFSA算法进行优化设计，得到5种加权系数组合情况下斜拉桥各工况的最优黏滞阻尼器参数组合，如表8所示。可以看出，组合1和组合2的结果较为接近，说明地震加权系数为1和0.8时对结果无明显影响；组合4和组合5优化结果较为接近，说明不考虑地震作用和较小考虑地震作用对振动控制影响并不明显。当抗震烈度增加时，阻尼系数C有着明显的变大趋势，而速度指数α一直维持在0.2～0.4之间。因此，阻尼系数C对振动控制起到决定作用。

表8 不同抗震设防烈度下参数优化组合(α, C)

图13 不同抗震设防烈度下结构响应控制效果

图13给出了不同加权系数组合下主梁跨中竖向位移w和梁端纵向位移Δ的控制效果对比，以下只基于组合1、组合3和组合5进行分析。从组合1至组合5，随着地震作用所占权重降低，控制效率也逐渐降低，说明在地震烈度高的地区，应加大地震作用权重，又不能完全忽略风和车辆荷载，可获得对结构竖向变形更优的控制效果。

随着抗震烈度的增加，主梁跨中竖向位移逐渐增大，控制效率逐渐降低。抗震烈度增加，梁端纵向位移逐渐增大，而控制效率维持在65%上下，对梁端纵向位移有着明显的控制效果。

4.5 风速影响

选定抗震设防烈度为6(加速度峰值0.05g)的E1地震动输入，选取平均风速为10，20，30 m/s的风荷载，得到5种加权系数组合下黏滞阻尼器最优参数组合，如表9所示。可以看出，从组合1～组合5，阻尼系数C和速度指数α都有着逐渐增大的趋势，阻尼系数C对振动控制起到决定作用。

图14给出了斜拉桥在最优阻尼器参数组合下, 主梁跨中最大竖向位移和梁端最大纵向位移数据及控制效果对比。加权系数组合1工况下，风速越大，控制效率越低。组合3和组合5随着平均风速的增加，梁端纵向位移逐渐增大，控制效率逐渐变大；随着地震作用所占权重的降低，跨中竖向位移控制效率整体降低，而梁端纵向位移控制效率整体变大，说明地震作用对竖向振动起决定作用。