基于混合螢火蟲-模擬退火算法的斜拉橋縱向振動控制與黏滯阻尼器優化設計

龍關旭，辛公鋒，徐興偉，任士樸，蘇文明

(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隨著平均風速的增加，梁端縱向位移逐漸增大，控制效率逐漸變大；隨著地震作用所占權重的降低，跨中豎向位移控制效率整體降低，而梁端縱向位移控制效率整體變大，說明地震作用對豎向振動起決定作用。

表9 不同風速下參數優化組合(α, C)

5 結論

(1)隨著黏滯阻尼器速度指數α的減小、阻尼系數的增大，控制效率也逐漸增大。地震引起的結構響應較為顯著，隨機車流次之，風荷載較小。

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

(3)綜合3種荷載組合下黏滯阻尼器優化結果，速度指數α在[0.2, 0.4]之間，阻尼系數C在[4 000, 5 000]kN/(m/s)α時控制效果達到最優。

圖14 不同風速下結構響應控制效果

(4)在地震多發、抗震設防烈度高的地區，建議選取加權系數組合2來進行黏滯阻尼器設計考慮；在沿海、山谷等風速較大區域，建議選取加權系數組合4來進行黏滯阻尼器設計考慮。