基于混合螢火蟲-模擬退火算法的斜拉橋縱向振動(dòng)控制與黏滯阻尼器優(yōu)化設(shè)計(jì)

龍關(guān)旭，辛公鋒，徐興偉，任士樸，蘇文明

(1.山東高速工程檢測(cè)有限公司，山東 濟(jì)南 250002；2.山東華鑒工程檢測(cè)有限公司，山東 濟(jì)南 250100；3.橋梁結(jié)構(gòu)大數(shù)據(jù)與性能診治提升重點(diǎn)實(shí)驗(yàn)室，山東 濟(jì)南 250002)

0 引言

大跨徑橋梁結(jié)構(gòu)輕柔、阻尼較低，對(duì)地震、風(fēng)、車等外部荷載引起的振動(dòng)更為敏感，影響著橋梁的結(jié)構(gòu)運(yùn)營安全和車輛駕駛舒適度。基于被動(dòng)控制理論而設(shè)置的黏滯阻尼器技術(shù)發(fā)展趨于成熟，在不改變橋梁剛度的前提下可有效提高結(jié)構(gòu)阻尼，減小結(jié)構(gòu)振動(dòng)響應(yīng)，并對(duì)伸縮縫起到一定的保護(hù)作用，已被廣泛應(yīng)用到橋梁結(jié)構(gòu)振動(dòng)控制中[1-2]。

20世紀(jì)90年代國內(nèi)外就已通過模型試驗(yàn)與工程應(yīng)用對(duì)黏滯阻尼器的減振消能作用展開研究[3-4]。葉愛君[5]研究得出了黏滯阻尼器對(duì)斜拉橋在地震作用下的結(jié)構(gòu)響應(yīng)有著良好的控制作用，可明顯減小橋梁縱向位移和塔底彎矩。焦馳宇等[6]研究了黏滯阻尼器在市政橋梁加固上的應(yīng)用，得出了該加固措施在工程造價(jià)、交通影響和可實(shí)施性等方面都有著明顯優(yōu)勢(shì)。丁幼亮等[7]通過非線性時(shí)域分析，研究了黏滯阻尼器對(duì)多塔斜拉橋風(fēng)致振動(dòng)的控制效果，黏滯阻尼器的設(shè)置可使得各橋塔和主梁的風(fēng)振響應(yīng)幅值趨于一致。Guo[8]以潤揚(yáng)大橋?yàn)橐劳校贸隽伺苘嚬r下黏滯阻尼器可有效控制梁端縱向位移。同時(shí)，很多學(xué)者對(duì)阻尼器參數(shù)的優(yōu)化做了相關(guān)研究[9-10]，王波等[11]基于斜拉橋的隨機(jī)地震響應(yīng)，采用非線性規(guī)劃法對(duì)黏滯阻尼器關(guān)鍵參數(shù)進(jìn)行了優(yōu)化取值。目前，不論是對(duì)黏滯阻尼器的振動(dòng)控制研究，還是對(duì)黏滯阻尼器關(guān)鍵參數(shù)的優(yōu)化研究，都是考慮以地震為主的單一荷載作用，且參數(shù)優(yōu)化多采用傳統(tǒng)標(biāo)準(zhǔn)算法來實(shí)現(xiàn)。

為此，以某大跨斜拉橋?yàn)槔瑢?duì)地震、風(fēng)、隨機(jī)車流作用下黏滯阻尼器的振動(dòng)控制進(jìn)行研究，并基于標(biāo)準(zhǔn)螢火蟲算法和模擬退火算法提出混合螢火蟲-模擬退火算法(HFSA)，對(duì)黏滯阻尼器在3種荷載的共同作用下的關(guān)鍵參數(shù)進(jìn)行優(yōu)化設(shè)計(jì)。

1 黏滯阻尼器及斜拉橋模型的建立

黏滯阻尼器是一種無剛度、速度型阻尼器，通常由缸體、導(dǎo)桿、黏滯流體材料、活塞和阻尼孔等組成，構(gòu)造如圖1所示。

圖1 黏滯阻尼器構(gòu)造

黏滯阻尼器的力學(xué)模型理論計(jì)算公式為：

Fd=Cvαsgn(v)，

(1)

式中，F(xiàn)d為阻尼力；C為阻尼系數(shù)；v為相對(duì)速度；α為速度指數(shù)，其范圍為0.1～1。

以某大跨雙塔斜拉橋?yàn)橐劳校鐝讲贾脼?(100+100+300+1 088+300+100+100) m，主梁采用扁平流線形鋼箱梁，索塔采用倒Y形混凝土塔，全線橋涵設(shè)計(jì)荷載采用汽車-超20級(jí)，掛車-120級(jí)。

圖2 斜拉橋ANSYS模型

采用有限元軟件ANSYS建立橋梁空間三維有限元模型，如圖2所示。模型采用Beam4單元來模擬結(jié)構(gòu)的主塔、主梁、蓋板和墩柱等構(gòu)件，斜拉索使用Link10單元進(jìn)行模擬，利用Ernst公式進(jìn)行彈性模量修正，過渡墩和輔助墩墩頂壓重質(zhì)量使用Mass21模擬，只考慮單元的質(zhì)量不考慮剛度。全橋共870個(gè)節(jié)點(diǎn)，1 134個(gè)單元，將斜拉索單元設(shè)置初始應(yīng)變參數(shù)，考慮大變形非線性求解，激活應(yīng)力剛度效應(yīng)。塔梁之間采用節(jié)點(diǎn)耦合的方式模擬支座約束，質(zhì)量單元和主梁之間采用共節(jié)點(diǎn)的方式連接，索梁之間建立剛臂連接。梁端約束根據(jù)設(shè)計(jì)資料確定，墩底采用固結(jié)約束。

在兩個(gè)主塔和主梁間分別設(shè)置兩個(gè)縱向阻尼器，全橋共4個(gè)，阻尼器單元采用非線性彈簧單元(combin37)模擬，由兩個(gè)活動(dòng)節(jié)點(diǎn)(I，J)和兩個(gè)可選控制節(jié)點(diǎn)(K，L)組成，如圖3所示。在計(jì)算過程中，可以通過選取控制節(jié)點(diǎn)得出單元內(nèi)部的相對(duì)位移對(duì)單元實(shí)常數(shù)進(jìn)行設(shè)置，阻尼器單元與主梁之間的作用采用剛臂進(jìn)行連接。

圖3 combin37單元示意圖

2 外荷載作用下縱向黏滯阻尼器振動(dòng)控制

分別選用地震作用、風(fēng)荷載、隨機(jī)車流荷載進(jìn)行荷載輸入：根據(jù)橋址和《公路工程抗震規(guī)范》(JTG-B02—2013)[12]確定地震設(shè)計(jì)反應(yīng)譜，進(jìn)而作得3組加速度時(shí)程曲線，選用最大的一組進(jìn)行縱向和豎向地震加速度輸入，如圖4所示。根據(jù)依托斜拉橋相關(guān)技術(shù)資料，采用高雷諾數(shù)主梁靜三分力系數(shù)作為靜風(fēng)力輸入，如圖5所示。

采用諧波合成法生成三維脈動(dòng)風(fēng)場[13-14]，同時(shí)根據(jù)《公路橋梁抗風(fēng)設(shè)計(jì)規(guī)范》(JTG/T D60-01—2004)[15]，分別采用Lumley譜和Simiu譜對(duì)豎向和橫向脈動(dòng)風(fēng)速時(shí)程進(jìn)行模擬輸入，采用MATLAB軟件編制脈動(dòng)風(fēng)模擬程序。圖6給出了全橋關(guān)鍵風(fēng)速模擬示意點(diǎn)，以1#節(jié)點(diǎn)為例，15 m/s風(fēng)速時(shí)程曲線如圖7所示。

圖4 地面加速度時(shí)程曲線

圖5 主梁及橋塔三分力系數(shù)

圖6 關(guān)鍵風(fēng)速模擬點(diǎn)

圖7 1#節(jié)點(diǎn)脈動(dòng)風(fēng)速時(shí)程曲線

圖8 單主梁模型車輪荷載分配示意圖

根據(jù)《橋梁用黏滯流體阻尼器》(JTT 926—2014)[18]和已有黏滯阻尼器在大跨橋梁上的應(yīng)用，選定速度指數(shù)α取值0.1～0.6和1，阻尼系數(shù)C取值1 000，2 500，5 000，7 000，10 000，15 000，20 000 kN/(m/s)α，共設(shè)定49組(α，C)參數(shù)組合，分別選用主梁跨中彎矩、梁端縱向位移、梁端縱向加速度和塔根彎矩進(jìn)行結(jié)構(gòu)響應(yīng)分析。以地震作用為例，選用阻尼系數(shù)C=5 000 kN/(m/s)α，結(jié)構(gòu)在地震作用下的主要響應(yīng)時(shí)程曲線如圖9所示，對(duì)比3種荷載引起的結(jié)構(gòu)主要響應(yīng)極大值變化，如圖10所示。

圖9 地震作用下主要響應(yīng)時(shí)程曲線

圖10 斜拉橋主要響應(yīng)極大值變化

由圖9和圖10可看出，黏滯阻尼器的設(shè)置可有效減小結(jié)構(gòu)振動(dòng)響應(yīng)，且基本隨著速度指數(shù)的減小、阻尼系數(shù)的增大，振動(dòng)控制效果更為明顯。同時(shí)，地震引起的結(jié)構(gòu)主要響應(yīng)較為顯著，隨機(jī)車流次之，風(fēng)荷載較小。

將縱向黏滯阻尼器的振動(dòng)控制作用進(jìn)行直觀化，定義控制效率的概念，具體表達(dá)式為：

(2)

式中，θ為阻尼器的控制效率；A1為設(shè)置阻尼器后結(jié)構(gòu)的響應(yīng)；A2為無阻尼器的結(jié)構(gòu)響應(yīng)。由于黏滯阻尼器關(guān)鍵參數(shù)組合較多，表1只對(duì)部分參數(shù)組合進(jìn)行了表述。

3 混合螢火蟲-模擬退火優(yōu)化算法

3.1 標(biāo)準(zhǔn)算法3.1.1 螢火蟲算法

螢火蟲算法(Firefly Algorithm，F(xiàn)A)是一種仿生智能優(yōu)化算法，由劍橋?qū)W者Yang[19]受啟發(fā)于螢火蟲群體尾部亮光行為而首次提出。螢火蟲算法模型包括兩種關(guān)鍵因子：亮度和吸引度。亮度即目標(biāo)函數(shù)，決定著螢火蟲所處的位置優(yōu)劣進(jìn)而決定單只螢火蟲的移動(dòng)方向；吸引度決定螢火蟲移動(dòng)的距離，吸引度越大移動(dòng)距離就越長。在迭代循環(huán)過程中，螢火蟲亮度和吸引度不斷更新信息，從而對(duì)目標(biāo)值進(jìn)行優(yōu)化。

(1)螢火蟲相對(duì)熒光亮度

(3)

式中，I0為螢火蟲自身(r=0)的最大熒光亮度；γ為光吸收系數(shù)；rij為螢火蟲i和j之間的距離。

(2)螢火蟲吸引度

(4)

式中β0為光源處(r=0)最大吸引度，其余與式(3)含義相同。

(3)螢火蟲移動(dòng)位置更新

(5)

3.1.2模擬退火算法

模擬退火算法(Simulated Annealing， SA)是一種貪心算法，最早是由Metropolis等根據(jù)統(tǒng)計(jì)力學(xué)中的物質(zhì)退火方法而提出，即固體加熱至高溫再慢慢冷卻的過程[20]。

定義固體熱能函數(shù)(目標(biāo)函數(shù))為f(x)，當(dāng)前解為x，新解為x′，則熱能增量為Δf=f(x′)-f(x)，根據(jù)Metropolis準(zhǔn)則接受新解的概率為：

(6)

3.2 混合螢火蟲-模擬退火算法

螢火蟲算法對(duì)于全局優(yōu)化能力更強(qiáng)，局部搜索能力較差，而模擬退火算法具有較強(qiáng)的局部最優(yōu)，對(duì)全局最優(yōu)搜索效率較低。可將螢火蟲算法和模擬退火算法相互結(jié)合，取長補(bǔ)短，設(shè)計(jì)提出混合螢火蟲-模擬退火算法(HFSA)。

該算法以標(biāo)準(zhǔn)螢火蟲算法為總體框架，從空間內(nèi)隨機(jī)產(chǎn)生的初始解中進(jìn)行全局最優(yōu)搜索，然后對(duì)產(chǎn)生的較優(yōu)全局解分別進(jìn)行模擬退火操作，并更新到下一個(gè)螢火蟲群體中。整個(gè)過程反復(fù)迭代更新判斷，直到滿足收斂條件為止。這樣兩種算法的優(yōu)點(diǎn)充分結(jié)合了起來，大大提高了算法的效率，可迅速找出全局最優(yōu)解，增強(qiáng)了算法跳出局部最優(yōu)解的能力。具體流程如下。

表1 縱向黏滯阻尼器控制效率(單位：%)

(1)初始化算法參數(shù)

①螢火蟲算法信息：種群規(guī)模(螢火蟲個(gè)數(shù)m)，最大迭代次數(shù)(螢火蟲進(jìn)化次數(shù)Mmax)，相對(duì)吸引度β0，熒光吸收系數(shù)γ，步長因子δ，目標(biāo)函數(shù)f(x)。

②模擬退火算法信息：初始溫度T0，退火因子D，馬爾科夫鏈長L，進(jìn)行模擬退火算法的螢火蟲個(gè)數(shù)k。

(2)空間隨機(jī)初始化螢火蟲群體，更新確定螢火蟲初始位置狀態(tài)xi(i=1, 2, 3,…,m)，計(jì)算螢火蟲亮度Ii，即初始目標(biāo)函數(shù)值。

(3)螢火蟲種群更新迭代

①i=1，根據(jù)式(4)吸引度大小判斷移動(dòng)半徑和式(5)更新螢火蟲空間位置。

②根據(jù)式(3)更新計(jì)算螢火蟲相對(duì)亮度。

③重新賦值i，j。

④判斷i，j是否大于m，否則循環(huán)②～④。

(4)以收斂種群中較優(yōu)的k個(gè)解模擬退火操作，進(jìn)行局部最優(yōu)搜索。

①迭代次數(shù)it=0，產(chǎn)生螢火蟲個(gè)體新解x′k。

②計(jì)算目標(biāo)函數(shù)差值Δfk=f(x′k)-f(xk)。

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

④循環(huán)L次，然后結(jié)束局部更新搜索。

(5)進(jìn)行退火冷卻操作：it=it+1，溫度T=DT。

(6)判斷it>Mmax，算法結(jié)束，否則循環(huán)步驟(3)～(5)。

具體流程如圖11所示。

圖11 混合螢火蟲-模擬退火算法流程

3.3 HFSA算法算例驗(yàn)證

選用4個(gè)標(biāo)準(zhǔn)測(cè)試函數(shù)來進(jìn)行HFSA優(yōu)化算法收斂性能和尋優(yōu)效果的驗(yàn)證，并與標(biāo)準(zhǔn)螢火蟲算法和模擬退火算法進(jìn)行對(duì)比分析，具體參數(shù)取值如表2所示。

表2 測(cè)試函數(shù)參數(shù)

3.3.1Ackley函數(shù)

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

2.2.2 概率性敏感性分析 根據(jù)前文的假設(shè)，成本服從正態(tài)分布，轉(zhuǎn)移概率服從β分布。成本數(shù)據(jù)的狀態(tài)分布參考相關(guān)文獻(xiàn)[17-18，24，27]報(bào)道的參數(shù)的 95%CI值，計(jì)算其標(biāo)準(zhǔn)差；轉(zhuǎn)移概率的狀態(tài)分布則根據(jù)其均值和標(biāo)準(zhǔn)差計(jì)算其α、β值。采用TreeAge Pro 2011軟件對(duì)模型結(jié)果進(jìn)行概率性敏感性分析，對(duì)模型進(jìn)行次數(shù)為1 000的蒙特卡洛模擬，相關(guān)參數(shù)見表8。對(duì)模擬結(jié)果繪制成本-效用可接受曲線，詳見圖3。

(7)

3.3.2ShafferF6函數(shù)

(8)

3.3.3Sphere函數(shù)

(9)

3.3.4Rastrigrin函數(shù)

(10)

圖12 測(cè)試函數(shù)示意圖

3種算法采用相同的參數(shù)設(shè)置，如表3所示，分別運(yùn)行20次，表4給出了4種算法對(duì)4個(gè)測(cè)試函數(shù)的優(yōu)化結(jié)果匯總。

表3 算法主要參數(shù)

由表4可以看出，3種方法對(duì)二維函數(shù)優(yōu)化都較好，迭代次數(shù)較少，但是FA算法最優(yōu)值方差較大，SA和HFSA算法的魯棒性較好；而對(duì)多維函數(shù)，SA算法體現(xiàn)了較差的全局優(yōu)化性能，尤其對(duì)于多峰函數(shù)容易對(duì)局部產(chǎn)生過長搜索(5 000次迭代)。優(yōu)化改進(jìn)的HFSA算法，不論在單峰函數(shù)還是多峰函數(shù)都表現(xiàn)出了較好的魯棒性和適應(yīng)度，收斂速度更快，方差也較小。

表4 優(yōu)化結(jié)果

4 縱向黏滯阻尼器關(guān)鍵參數(shù)優(yōu)化

4.1 HFSA優(yōu)化參數(shù)的確定

設(shè)計(jì)變量選為黏滯阻尼器的關(guān)鍵參數(shù)：速度指數(shù)α∈[0.1,0.6]和阻尼系數(shù)C∈[1 000, 20 000]。

以經(jīng)濟(jì)效益為核心，綜合考慮地震-風(fēng)-隨機(jī)車流的外荷載作用下的阻尼力之和，構(gòu)建目標(biāo)函數(shù)：

Z=min(aFe+bFw+cFv)，

(11)

約束條件選定為斜拉橋設(shè)計(jì)規(guī)范中對(duì)結(jié)構(gòu)變形的要求，即控制主梁跨中豎向位移，并同時(shí)考慮梁端縱向位移變化：

w≤l/500，

(12)

Δ≤50 cm，

(13)

式中，l為橋梁主跨跨徑；w為斜拉橋跨中豎向撓度；Δ為梁端縱向位移。

4.2 基于響應(yīng)面法的目標(biāo)函數(shù)擬合

采用Design-Expert軟件的BBD法進(jìn)行試驗(yàn)設(shè)計(jì)，根據(jù)第2節(jié)的有限元分析計(jì)算結(jié)果進(jìn)行輸入，基于最小二乘法構(gòu)建多項(xiàng)式系數(shù)，從而得到響應(yīng)面模型函數(shù)。以地震作用下的最大阻尼力F為算例，然后采用響應(yīng)面法擬合函數(shù)，擬合回歸方程為：

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

(14)

對(duì)擬合的目標(biāo)函數(shù)進(jìn)行精度檢驗(yàn)，如表5所示。F統(tǒng)計(jì)值為3 493.67，說明模型具有較好的統(tǒng)計(jì)意義(大于95%的默認(rèn)置信區(qū)間)。變異系數(shù)1.58<10，表明試驗(yàn)的可信度和精確度較高，相關(guān)系數(shù)R2為0.999 8，表明擬合的函數(shù)具有較高的精度；0.996 4的預(yù)測(cè)R2與0.999 5的調(diào)整值具有合理的一致性。205.197的精度度量表明了模型具有較高的精密度。

表5 精度檢驗(yàn)指標(biāo)值

4.3 權(quán)重分配及工況設(shè)置

考慮到縱向黏滯阻尼器的設(shè)計(jì)初衷，以地震作用的權(quán)重系數(shù)為主因素，具體加權(quán)系數(shù)取值如表6所示。

考慮4種地震動(dòng)加速度峰值的抗震設(shè)防烈度區(qū)域、3種平均風(fēng)速荷載和1種隨機(jī)車流工況，分別進(jìn)行縱向黏滯阻尼器的優(yōu)化設(shè)計(jì)研究。具體荷載工況設(shè)置如表7所示。

表6 加權(quán)系數(shù)組合

表7 荷載工況組合

4.4 抗震設(shè)防烈度影響

選定風(fēng)速為10 m/s，利用HFSA算法進(jìn)行優(yōu)化設(shè)計(jì)，得到5種加權(quán)系數(shù)組合情況下斜拉橋各工況的最優(yōu)黏滯阻尼器參數(shù)組合，如表8所示。可以看出，組合1和組合2的結(jié)果較為接近，說明地震加權(quán)系數(shù)為1和0.8時(shí)對(duì)結(jié)果無明顯影響；組合4和組合5優(yōu)化結(jié)果較為接近，說明不考慮地震作用和較小考慮地震作用對(duì)振動(dòng)控制影響并不明顯。當(dāng)抗震烈度增加時(shí)，阻尼系數(shù)C有著明顯的變大趨勢(shì)，而速度指數(shù)α一直維持在0.2～0.4之間。因此，阻尼系數(shù)C對(duì)振動(dòng)控制起到?jīng)Q定作用。

表8 不同抗震設(shè)防烈度下參數(shù)優(yōu)化組合(α, C)

圖13 不同抗震設(shè)防烈度下結(jié)構(gòu)響應(yīng)控制效果

圖13給出了不同加權(quán)系數(shù)組合下主梁跨中豎向位移w和梁端縱向位移Δ的控制效果對(duì)比，以下只基于組合1、組合3和組合5進(jìn)行分析。從組合1至組合5，隨著地震作用所占權(quán)重降低，控制效率也逐漸降低，說明在地震烈度高的地區(qū)，應(yīng)加大地震作用權(quán)重，又不能完全忽略風(fēng)和車輛荷載，可獲得對(duì)結(jié)構(gòu)豎向變形更優(yōu)的控制效果。

隨著抗震烈度的增加，主梁跨中豎向位移逐漸增大，控制效率逐漸降低。抗震烈度增加，梁端縱向位移逐漸增大，而控制效率維持在65%上下，對(duì)梁端縱向位移有著明顯的控制效果。

4.5 風(fēng)速影響

選定抗震設(shè)防烈度為6(加速度峰值0.05g)的E1地震動(dòng)輸入，選取平均風(fēng)速為10，20，30 m/s的風(fēng)荷載，得到5種加權(quán)系數(shù)組合下黏滯阻尼器最優(yōu)參數(shù)組合，如表9所示。可以看出，從組合1～組合5，阻尼系數(shù)C和速度指數(shù)α都有著逐漸增大的趨勢(shì)，阻尼系數(shù)C對(duì)振動(dòng)控制起到?jīng)Q定作用。

圖14給出了斜拉橋在最優(yōu)阻尼器參數(shù)組合下, 主梁跨中最大豎向位移和梁端最大縱向位移數(shù)據(jù)及控制效果對(duì)比。加權(quán)系數(shù)組合1工況下，風(fēng)速越大，控制效率越低。組合3和組合5隨著平均風(fēng)速的增加，梁端縱向位移逐漸增大，控制效率逐漸變大；隨著地震作用所占權(quán)重的降低，跨中豎向位移控制效率整體降低，而梁端縱向位移控制效率整體變大，說明地震作用對(duì)豎向振動(dòng)起決定作用。

表9 不同風(fēng)速下參數(shù)優(yōu)化組合(α, C)

5 結(jié)論

(1)隨著黏滯阻尼器速度指數(shù)α的減小、阻尼系數(shù)的增大，控制效率也逐漸增大。地震引起的結(jié)構(gòu)響應(yīng)較為顯著，隨機(jī)車流次之，風(fēng)荷載較小。

(2)提出的混合螢火蟲-模擬退火算法(HFSA)具有收斂速度快、優(yōu)化精度高的特點(diǎn)。

(3)綜合3種荷載組合下黏滯阻尼器優(yōu)化結(jié)果，速度指數(shù)α在[0.2, 0.4]之間，阻尼系數(shù)C在[4 000, 5 000]kN/(m/s)α?xí)r控制效果達(dá)到最優(yōu)。

圖14 不同風(fēng)速下結(jié)構(gòu)響應(yīng)控制效果

(4)在地震多發(fā)、抗震設(shè)防烈度高的地區(qū)，建議選取加權(quán)系數(shù)組合2來進(jìn)行黏滯阻尼器設(shè)計(jì)考慮；在沿海、山谷等風(fēng)速較大區(qū)域，建議選取加權(quán)系數(shù)組合4來進(jìn)行黏滯阻尼器設(shè)計(jì)考慮。