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

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

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

0 引言

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

1 工程概況

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

圖1 (25+35+25)m雙線連續(xù)軌道梁實(shí)拍現(xiàn)場圖

2 初始有限元模型建立

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

本文采用經(jīng)驗(yàn)法與ABAQUS敏感度分析方法相結(jié)合來確定有限元模型修正的設(shè)計(jì)參數(shù)和設(shè)計(jì)響應(yīng)。由于中低速磁浮軌道梁上部結(jié)構(gòu)比較復(fù)雜，影響其動(dòng)態(tài)特性的參數(shù)很多，如箱梁主體的質(zhì)量與剛度、軌道結(jié)構(gòu)的質(zhì)量與剛度、橋墩的質(zhì)量與剛度等，因此，在進(jìn)行有限元模型之前確定適當(dāng)?shù)脑O(shè)計(jì)參數(shù)尤為重要。

各個(gè)參數(shù)影響不一，對于那些對軌道梁動(dòng)力特性的影響較大的結(jié)構(gòu)參數(shù)，才是有限元模型修正應(yīng)該著重加以控制的設(shè)計(jì)參數(shù)，才能盡可能接近實(shí)際狀態(tài)，確定實(shí)際橋梁的承載能力和工作性能，因此對眾多的影響結(jié)構(gòu)動(dòng)力特性的參數(shù)進(jìn)行敏感性分析就顯得很有必要。

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

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

表1 建模初始參數(shù)表Table 1 Modeling initial parameter list材料性質(zhì)彈性模量/(×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.箱梁實(shí)際截面尺寸(見圖2、圖3)。

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

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

c.軌道結(jié)構(gòu)實(shí)際尺寸圖(見圖4)。

(a) 軌道梁軌道結(jié)構(gòu)部分橫斷面圖

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

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

3 模擬退火優(yōu)化算法

退火算法合理依據(jù)了物體退火結(jié)晶原理，當(dāng)固體溫度很高時(shí)，內(nèi)部粒子處于快速且無序的運(yùn)動(dòng)狀態(tài)，內(nèi)部能量較大，當(dāng)慢慢開始降低固體溫度時(shí)，粒子的運(yùn)動(dòng)狀態(tài)慢慢變得規(guī)則有序，其內(nèi)部能量也開始慢慢減少，直到當(dāng)固體的溫度降至常溫時(shí)，其內(nèi)部能量可以達(dá)到最小，此時(shí)粒子的運(yùn)動(dòng)最為規(guī)則和穩(wěn)定。該算法首先從指定的一個(gè)較高溫度出發(fā)開始做迭代計(jì)算，隨機(jī)生成一個(gè)初始解，并在該解的領(lǐng)域范圍內(nèi)做隨機(jī)擾動(dòng)產(chǎn)生一個(gè)新解，根據(jù)Metropolis，確認(rèn)新解是否可接受，并且如不能滿足程序終止條件，則以某種方式降溫，并開始新的迭代。

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

p=

將p與隨機(jī)數(shù)α∈[0,1]進(jìn)行比較，p<α，則接受xnew，用來替代xold，否則舍棄。該過程在不斷循環(huán)重復(fù)中可產(chǎn)生新狀態(tài)，直至系統(tǒng)能量最低的平衡狀態(tài)，終止此過程。

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

模擬退火算法總共分為兩層循環(huán)，其中一層循環(huán)是在初始解鄰域范圍內(nèi)隨機(jī)擾動(dòng)產(chǎn)生一個(gè)新解并計(jì)算和比較兩者目標(biāo)函數(shù)值，決定是否被接受；而另一層則是由于算法給定的溫度初始值較高，需經(jīng)過一定數(shù)量的迭代才能緩慢降溫，最終收斂得到全局最優(yōu)解。具體步驟如下：

a.隨機(jī)給出初始模型參數(shù)的初始解，并計(jì)算其目標(biāo)函數(shù)值E(m0)。

b.設(shè)置溫度初始值T0和迭代計(jì)數(shù)器t。

c.在初始解鄰域范圍內(nèi)隨機(jī)進(jìn)行一個(gè)擾動(dòng)，產(chǎn)生新的模型參數(shù)值m1=m0+Δm，并計(jì)算新的目標(biāo)函數(shù)值E(m1)。

d.計(jì)算兩種情況下的目標(biāo)函數(shù)值之差ΔE，若ΔE≤0，接受新模型參數(shù)值，否則根據(jù)Metropolis標(biāo)準(zhǔn)以一定的概率可否接受新值。

e.若計(jì)算迭代次數(shù)t小于設(shè)定的最大迭代次數(shù)，則t=t+1，返回至c.，繼續(xù)計(jì)算循環(huán)迭代過程。

f.如果已到達(dá)設(shè)定的最大迭代次數(shù)，則判斷是否達(dá)到溫度終止條件，若未至冷卻狀態(tài)，則通過指定的降溫方式冷卻，并轉(zhuǎn)至c.；若達(dá)到終止條件，則停止算法，輸出的當(dāng)前模型參數(shù)值為全局最優(yōu)解。

4 Python-ABAQUS的模擬退火算法設(shè)計(jì)

Python作為計(jì)算機(jī)編程語言，擁有非常強(qiáng)大的面向?qū)ο蟮哪芰?，在不同的平臺上，用戶可二次開發(fā)和編寫腳本。

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

圖6 模擬退火算法優(yōu)化過程示意圖解

a.確定設(shè)計(jì)參數(shù)。

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

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

b.選取狀態(tài)變量。

對(25+35+25)m雙線連續(xù)軌道梁，根據(jù)動(dòng)力特性測試試驗(yàn)所得結(jié)果，取橫向前3階頻率和豎向前3階頻率作為狀態(tài)變量，匯總見表3。

表3 狀態(tài)變量表Table 3 State Variable Table頻率階次實(shí)測值/Hz15.365豎向頻率28.95039.85912.450橫向頻率23.33734.370

c.目標(biāo)函數(shù)構(gòu)造。

本文取計(jì)算頻率和實(shí)測頻率的相對誤差平方和作為目標(biāo)函數(shù)，使其取得最小值時(shí)所對應(yīng)的設(shè)計(jì)參數(shù)值即為最優(yōu)狀態(tài)值，頻率相對誤差平方和可以表示為：

式中：λai，λti分別代表第i階計(jì)算頻率和實(shí)測頻率。

d.算法設(shè)計(jì)。

① 初始溫度。

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

② 溫度衰減。

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

③ 終止規(guī)則。

模擬退火算法經(jīng)常使用結(jié)束溫度作為其終止準(zhǔn)則，即當(dāng)溫度達(dá)到結(jié)束溫度時(shí)退火過程終止。隨著模擬退火，退火溫度呈指數(shù)下降。實(shí)際上，溫度會降到一個(gè)無限接近于0的值。此時(shí)，更小的狀態(tài)能量難以出現(xiàn)，根據(jù)Metropolis判斷規(guī)則目標(biāo)函數(shù)難以獲得更優(yōu)值。以溫度衰減系數(shù)為0.9時(shí)為例，溫度從200降到0.0001將發(fā)生138次轉(zhuǎn)移，足以獲得全局最優(yōu)方案，故本文中結(jié)束溫度皆設(shè)為一個(gè)定值0.0001。

5 修正前后軌道梁結(jié)構(gòu)動(dòng)力特性結(jié)果對比分析

利用ABAQUS基于Python腳本接口的二次開發(fā)功能和設(shè)計(jì)敏感性分析的基本原理與方法，運(yùn)用經(jīng)驗(yàn)法選取了5個(gè)與結(jié)構(gòu)動(dòng)力特性直接相關(guān)且能直觀反映橋梁實(shí)際狀態(tài)的結(jié)構(gòu)參數(shù)，并對其進(jìn)行了敏感性分析。

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

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

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

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

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

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

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

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

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

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

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

表4 修正前后各主要設(shè)計(jì)參數(shù)的變化對比Table 4 Comparison of changes in main design parameters before and after correction主要設(shè)計(jì)參數(shù)箱梁質(zhì)量密度/(kg·m-3)箱梁彈性模量/MPa橋墩質(zhì)量密度/(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參數(shù)變化5.42%15.2%6.58%8.62%0.09%

表5 修正前后前10階實(shí)測頻率與計(jì)算頻率的對比Table 5 Comparison of the measured frequency and the calculated frequency of the first 10 steps before and after correction階次實(shí)測頻率/Hz修正前修正后計(jì)算頻率/Hz誤差/%計(jì)算頻率/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扭轉(zhuǎn)一階

6 結(jié)論

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

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

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