模擬退火算法與Python-ABAQUS結合的有限元模型修正方法研究

成 魁, 王祺順,, 張祖軍, 田仲初, 吳 欣

(1.廣西交通科學研究院有限公司, 廣西 南寧 530005; 2.湖南省交通科學研究院有限公司, 湖南 長沙 410015; 3.長沙理工大學, 湖南 長沙 410114)

0 引言

在對中低速磁浮軌道梁等復雜結構模型的仿真分析中，通常有限元模型單純按設計圖紙所給尺寸和材料參數等建立，往往得出的是理想的狀態結果，考慮到結構實際的運行工作情況，此模型結果難以反映出軌道梁結構的實際工作狀態，故需利用靜動力試驗分析結果，如撓度、應變、頻率、振型、阻尼比等，在保證模型符合實際尺寸和精度的狀態下，調整或修正部分由于施工質量、環境等因素影響的設計參數值，使計算結果與試驗結果盡可能統一，從而更加逼近于實際狀態。

1 工程概況

本文以長沙中低速磁浮交通工程(25+35+25)m連續梁橋為工程背景(見圖1)，對其進行有限元模型修正研究。該橋單線軌道梁截面形式為單箱單室等高箱梁截面，單線軌道梁之間通過端橫梁及中橫梁連接，箱梁對應橫梁位置設置隔板。梁高2.4m，箱梁頂板寬1.3m，底板寬1.4m，線間距4.4m，兩線箱梁總寬5.7m。支座中心至梁端0.35m，箱梁橫橋向每線下設一個支座，邊、中支座橫橋向中心距均為4.4m。

圖1 (25+35+25)m雙線連續軌道梁實拍現場圖

2 初始有限元模型建立

為使所建立的模型更接近于磁浮軌道梁的現場實際運營的狀態，本文利用ABAQUS軟件對(25+35+25)m雙線連續軌道梁主體結構和軌道結構均采用實體單元精細化模擬，橋墩與箱梁之間采用運動耦合形式，約束部分自由度來模擬實際的邊界條件。

本文采用經驗法與ABAQUS敏感度分析方法相結合來確定有限元模型修正的設計參數和設計響應。由于中低速磁浮軌道梁上部結構比較復雜，影響其動態特性的參數很多，如箱梁主體的質量與剛度、軌道結構的質量與剛度、橋墩的質量與剛度等，因此，在進行有限元模型之前確定適當的設計參數尤為重要。

各個參數影響不一，對于那些對軌道梁動力特性的影響較大的結構參數，才是有限元模型修正應該著重加以控制的設計參數，才能盡可能接近實際狀態，確定實際橋梁的承載能力和工作性能，因此對眾多的影響結構動力特性的參數進行敏感性分析就顯得很有必要。

具體截面幾何尺寸和材料參數主要依據工程施工圖紙確定。

a.建模初始參數(見表1)。

表1 建模初始參數表Table 1 Modeling initial parameter list材料性質彈性模量/(×104MPa)泊松比密度/(kg·m-3)箱梁3.550.22 500橫梁3.550.22 500承軌臺3.400.22 400橋墩3.400.22 400H型鋼軌枕21.000.37 850F軌21.000.37 850

b.箱梁實際截面尺寸(見圖2、圖3)。

圖2 軌道梁箱梁跨中斷面的實際尺寸圖(單位： cm)

圖3 軌道梁箱梁支點斷面的實際尺寸圖(單位： cm)

c.軌道結構實際尺寸圖(見圖4)。

(a) 軌道梁軌道結構部分橫斷面圖

d.建立的初始有限元模型(見圖5)。

圖5 ABAQUS初始有限元模型示意圖

3 模擬退火優化算法

退火算法合理依據了物體退火結晶原理，當固體溫度很高時，內部粒子處于快速且無序的運動狀態，內部能量較大，當慢慢開始降低固體溫度時，粒子的運動狀態慢慢變得規則有序，其內部能量也開始慢慢減少，直到當固體的溫度降至常溫時，其內部能量可以達到最小，此時粒子的運動最為規則和穩定。該算法首先從指定的一個較高溫度出發開始做迭代計算，隨機生成一個初始解，并在該解的領域范圍內做隨機擾動產生一個新解，根據Metropolis，確認新解是否可接受，并且如不能滿足程序終止條件，則以某種方式降溫，并開始新的迭代。

exp(-ΔE/(kT))用于表示粒子在溫度T下趨于平衡的概率，其中E為溫度T時的內能，ΔE為其改變數；k為Boltzmann常數。Metropolis準則常表示為：

p=

將p與隨機數α∈[0,1]進行比較，p<α，則接受xnew，用來替代xold，否則舍棄。該過程在不斷循環重復中可產生新狀態，直至系統能量最低的平衡狀態，終止此過程。

將目標函數值f作為內能E模，控制參數t作為溫度T，便可衍生出解組合優化問題的模擬退火演算法：從控制參數的初始值t與初始解i開始，重復當前解“① 產生新解；② 計算目標函數差；③ 接受或丟棄”，該循環在逐漸衰減t值的同時，最優解可進行近似處理。利用Cooling Schedule(冷卻進度表)控制退火過程，惡化解可以在一定的概率下依據準則被接受。

模擬退火算法總共分為兩層循環，其中一層循環是在初始解鄰域范圍內隨機擾動產生一個新解并計算和比較兩者目標函數值，決定是否被接受；而另一層則是由于算法給定的溫度初始值較高，需經過一定數量的迭代才能緩慢降溫，最終收斂得到全局最優解。具體步驟如下：

a.隨機給出初始模型參數的初始解，并計算其目標函數值E(m0)。

b.設置溫度初始值T0和迭代計數器t。

c.在初始解鄰域范圍內隨機進行一個擾動，產生新的模型參數值m1=m0+Δm，并計算新的目標函數值E(m1)。

d.計算兩種情況下的目標函數值之差ΔE，若ΔE≤0，接受新模型參數值，否則根據Metropolis標準以一定的概率可否接受新值。

e.若計算迭代次數t小于設定的最大迭代次數，則t=t+1，返回至c.，繼續計算循環迭代過程。

f.如果已到達設定的最大迭代次數，則判斷是否達到溫度終止條件，若未至冷卻狀態，則通過指定的降溫方式冷卻，并轉至c.；若達到終止條件，則停止算法，輸出的當前模型參數值為全局最優解。

4 Python-ABAQUS的模擬退火算法設計

Python作為計算機編程語言，擁有非常強大的面向對象的能力，在不同的平臺上，用戶可二次開發和編寫腳本。

ABAQUS軟件有效地使用Python腳本語言為使用者提供了許多庫函數，實現有限元建模、單元網格劃分、指定材料屬性、提交作業和可視化后處理等操作。為了充分利用ABAQUS的二次開發功能，采用Python語言編寫了(25+35+25)m雙線連續軌道梁的全橋代碼。模擬退火算法優化過程見圖6。

圖6 模擬退火算法優化過程示意圖解

a.確定設計參數。

利用ABAQUS設計敏感性分析程序計算得到箱梁主體、橋墩、承軌臺、H型鋼軌枕和F軌的相關材料參數對雙線連續軌道梁前10階固有頻率的歸一化敏感度，軌道梁動力特性的主要影響因素是箱梁和墩的質量，以及F軌的剛度，其他為次要因素。本次模型修正選用了5個設計參數，分別是箱梁、橋墩的密度和彈性模量，以及F軌的彈性模量，選取結果匯總見表2。

表2 設計參數表Table 2 Design Parameter Table材料性質彈性模量/MPa密度/(kg·m-3)箱梁3.55×1042 500橋墩3.40×1042 400F軌21.00×104—

b.選取狀態變量。

對(25+35+25)m雙線連續軌道梁，根據動力特性測試試驗所得結果，取橫向前3階頻率和豎向前3階頻率作為狀態變量，匯總見表3。

表3 狀態變量表Table 3 State Variable Table頻率階次實測值/Hz15.365豎向頻率28.95039.85912.450橫向頻率23.33734.370

c.目標函數構造。

本文取計算頻率和實測頻率的相對誤差平方和作為目標函數，使其取得最小值時所對應的設計參數值即為最優狀態值，頻率相對誤差平方和可以表示為：

式中：λai，λti分別代表第i階計算頻率和實測頻率。

d.算法設計。

① 初始溫度。

在模擬退火過程中，要想較大可能性獲得全局最優值，則可設定足夠高的初始溫度。然而，過高的初始溫度往往導致數字較大的CPU時間。需要注意的是：模擬退火算法受啟發于金屬材料退火過程，其初始溫度、退火溫度、終止溫度作為模擬退火的重要參數，并不像金屬退火過程具有℃或K等單位。實際上，應根據以下條件選擇初始溫度：根據Metropolis判斷規則初始溫度時新的狀態能量接受率(狀態能量發生轉移的粒子數與新的狀態能量總數粒子總數)近乎100%；接受率隨退火溫度的降低逐漸降低，取為100 ℃。

② 溫度衰減。

溫度衰減對于模擬退火算法收斂性和CPU時間對于具有決定性作用。如果降溫速率太慢，CPU時間數字較大，太快則易于得到局部最優。本文采用Kirkpatrick的幾何溫度衰減函數，溫度衰減系數可在0.80～0.99范圍內進行選擇，取為0.98。

③ 終止規則。

模擬退火算法經常使用結束溫度作為其終止準則，即當溫度達到結束溫度時退火過程終止。隨著模擬退火，退火溫度呈指數下降。實際上，溫度會降到一個無限接近于0的值。此時，更小的狀態能量難以出現，根據Metropolis判斷規則目標函數難以獲得更優值。以溫度衰減系數為0.9時為例，溫度從200降到0.0001將發生138次轉移，足以獲得全局最優方案，故本文中結束溫度皆設為一個定值0.0001。

5 修正前后軌道梁結構動力特性結果對比分析

利用ABAQUS基于Python腳本接口的二次開發功能和設計敏感性分析的基本原理與方法，運用經驗法選取了5個與結構動力特性直接相關且能直觀反映橋梁實際狀態的結構參數，并對其進行了敏感性分析。

通過比較分析箱梁主體、橋墩、承軌臺、H型鋼軌枕和F軌的相關材料參數對雙線連續軌道梁前10階固有頻率的歸一化敏感度，得出了軌道梁動力特性的主要影響因素為箱梁和橋墩的質量和剛度、F軌的剛度，對于軌道梁動力特性的影響，可以看出其他部分均為次要因素。

利用模擬退火算法與Python-ABAQUS相結合的有限元模型修正方法對影響(25+35+25)m雙線連續軌道梁動力特性的主要設計參數在一定的可行域內進行了修正，經多次迭代過程，得到了最佳設計參數來，此時選定的目標函數取得最小值，修正前后的計算結果對比見表4、表5，部分振型見圖7～圖14。

圖7 實測第1階振型(f=2.450 Hz，橫向)

圖8 修正后第1階振型(f=2.433 Hz，橫向)

圖9 實測第2階振型(f=3.337 Hz，橫向)

圖10 修正后第2階振型(f=3.294 Hz，橫向)

圖11 實測第4階振型(f=5.365 Hz，豎向)

圖12 修正后第4階振型(f=5.270 Hz，豎向)

圖13 實測第7階振型(f=8.950 Hz，豎向)

圖14 修正后第7階振型(f=8.811 Hz，豎向)

由上述分析結果可知，通過選取對軌道梁結構動力特性影響較大的5個主要設計參數，以計算頻率和實測頻率的相對誤差平方和作為目標函數，利用模擬退火算法與Python-ABAQUS二次開發結合的有限元模型修正方法，成功地對(25+35+25)m雙線連續軌道梁初始有限元模型進行了修正，其中參數變化值最大的為箱梁彈性模量，比初始有限元模型增加了15.2%，并與動力試驗的實測頻率和振型的對比結果表明，修正后的(25+35+25)m雙線連續軌道梁模型能正確地反映結構的實際動力特性。

表4 修正前后各主要設計參數的變化對比Table 4 Comparison of changes in main design parameters before and after correction主要設計參數箱梁質量密度/(kg·m-3)箱梁彈性模量/MPa橋墩質量密度/(kg·m-3)橋墩彈性模量/MPaF軌彈性模量MPa修正前2 500.0003.550 0×1042 400.0003.400 0×10421.000×104修正后2 635.4524.089 6×1042 557.8253.693 1×10421.018×104參數變化5.42%15.2%6.58%8.62%0.09%

表5 修正前后前10階實測頻率與計算頻率的對比Table 5 Comparison of the measured frequency and the calculated frequency of the first 10 steps before and after correction階次實測頻率/Hz修正前修正后計算頻率/Hz誤差/%計算頻率/Hz誤差/%振型描述12.4502.307-5.842.433-0.68橫向彎曲一階23.3372.996-10.223.294-1.28橫向彎曲二階34.3703.685-15.684.3870.39橫向彎曲三階45.3654.764-6.015.270-1.80豎向彎曲一階55.4325.181-4.625.5251.72橫向彎曲四階66.7026.257-6.646.604-1.47橫向彎曲五階78.9508.249-7.838.811-1.56豎向彎曲二階89.5368.888-6.809.422-1.20橫向彎曲六階99.8599.539-3.2510.0842.28豎向彎曲三階1010.9509.918-9.4210.593-3.30扭轉一階

6 結論

本文提出了一種利用Python-ABAQUS二次開發結合模擬退火算法的有限元模型修正方法，并將這種方法應用于長沙磁浮交通工程(25+35+25)m雙線軌道梁的空間有限元模型修正中，得到了較為理想的修正模型。主要研究結論如下：

a.通過敏感度分析，選取了對軌道梁結構動力特性影響較大的5個主要設計參數，選取計算頻率和實測頻率的相對誤差平方和作為目標函數，利用模擬退火算法與Python語言結合的有限元模型修正方法，得到了使目標函數取得最小值時的最優設計參數，其中參數變化值最大的為箱梁彈性模量，比初始有限元模型增加了15.2%，變化最小的為F軌彈性模量，僅增加了0.09%。

b.通過修正前后的軌道梁動力特性與動力試驗的實測頻率和振型的對比結果表明，利用模擬退火算法與Python語言結合的有限元模型修正方法修正后的(25+35+25)m雙線連續軌道梁模型能正確地反映結構的實際動力特性。