基于組合函數響應面的橋梁有限元模型修正

張松涵，高芳清，張偉偉

(1.西南交通大學土木工程學院，四川成都610031;2.西南交通大學 力學與工程學院，四川成都610031)

0 引言

近年來，基于響應面的有限元模型修正方法在工程中已得到廣泛應用，其方法在于使用顯式函數關系逼近設計參數與響應值之間的隱式復雜關系，與傳統有限元模型修正相比，有效避免了迭代過程中重復調用有限元計算的缺點，使迭代速度大大加快［1］。

響應面模型能否準確地反映出有限元模型中的設計參數-響應值關系，是模型修正成敗的關鍵。目前，橋梁有限元模型修正中大多采用2階多項式或3階多項式作為響應面函數［2］，而隨著函數階次的提高和設計參數的增多，自變量數目急劇增加，在顯著性檢驗、響應面擬合及優化求解中都會造成計算量過大，故高階多項式函數響應面在有限元模型修正中難以普及。尋求綜合能力更好、更易于通用的響應面函數是當今有待解決的重要問題。

周林仁等［3］采用高斯徑向基函數響應面完成了大跨度斜拉橋有限元模型修正，證實了其擬合精度和抗噪能力滿足工程應用需要;孔憲仁等［4］和秦玉靈等［5］采用1階式-高斯組合徑向基函數響應面完成蜂窩板和機翼有限元模型修正，證實其擬合能力較1階多項式函數和高斯徑向基函數有所提高;張松涵等［6］將1階式-高斯組合徑向基函數、2階多項式、2階式-高斯組合徑向基函數、3階多項式用于簡支梁數值模型修正，經對比發現，在含有5%噪聲的條件下，2階式-高斯組合徑向基函數能保持較強的擬合能力和較高的求解精度。

為驗證2階式-高斯組合徑向基函數響應面在實際工程中應用的有效性，筆者通過模型修改獲得合理的設計參數，在重復試驗的基礎上建立了模態頻率響應面函數，構造殘差型目標函數，利用粒子群優化算法求解，得到設計參數最優解。在保證求解精度不受影響的前提下，避免了求解過程中反復調用有限元分析與靈敏度矩陣計算，適用于結構形式復雜的橋梁有限元模型修正。

1 基于響應面的有限元模型修正方法

基于響應面的有限元模型修正通過重復試驗，將復雜結構的有限元模型轉化為簡單的響應面數學模型［7］，結合實測數據，采用優化算法尋求到響應面模型中的最優解，再代入有限元模型進行驗證，從而達到模型修正的目的，其具體流程見圖1。

圖1 基于響應面的有限元模型修正流程Fig.1 The updating process of finite element model based on response surfaces

1.1 試驗設計

常見的試驗設計方法有均勻試驗設計、正交設計、中心復合設計、D-最優試驗設計等。試驗樣本數量應保證樣本矩陣具有一定超定次數，避免求解過程產生奇異。考慮試驗因素、水平數量，筆者采用U30(303)均勻試驗設計表［8］，共生成30個試驗樣本。

1.2 組合函數響應面

徑向基函數是一種徑向對稱的函數，形式簡單、不依賴于空間維數，且具有良好的非線性擬合能力。構造RBF的線形組合模型形式如式(1):

式中:[a1，a2，a3，…，an]為回歸系數向量;n為修正參數的個數;x為修正參數值;xi為第i個參數的樣本均值;x－x( )i表示修正參數取值到對應樣本均值的Euclidean范數;σ為高斯核寬度參數。

將工程中普遍常用的2階多項式與高斯徑向基函數結合，得到改進后的響應面函數如式(2):

式中:核寬度參數σ一般取設計空間半徑［4］。

1.3 經典最小二乘法回歸響應面

將樣本參數取值視為輸入，計算樣本頻率視為輸出，建立關于待定系數向量的方程:

式中:H為輸入矩陣;Z為輸出向量;θ為待定系數向量;e為殘差向量。

要使殘差平方和

最小，對θ求偏導，令其等于0，即:

得到:

即為待定系數向量的經典最小二乘解，將其帶回原組合函數，即可得到擬合后的回歸響應面方程。

1.4 基于自然選擇的混合粒子群優化算法

基于自然選擇的混合粒子群算法將基本粒子群算法與自然選擇機理相結合，其主要手段在于將整個群體在迭代過程中按照適應值排序，用群體中最優的一半粒子替換最差的一半，同時將每個粒子所經歷過的歷史最優點全部保留，能進行全局搜索，有效避免陷入局部極小值，且能夠避開反復的靈敏度矩陣計算［9-10］。文中的目標函數優化求解均采用此算法。算法步驟如下:

1)將各個微粒按照隨機原則賦予初始位置和初始速度;

2)計算各粒子的適應值并對其進行評價，把每個微粒所對應位置和適應度儲存到個體極值pbest中，在所有pbest中挑出適應值最優微粒的位置xi和適應值pbest，儲存到總體極值gbest中;

3)更新每個粒子的速度和位置;

4)將各個微粒的適應值與其途徑過的最佳位置作比較，如果當前位置更好，則將其替換到當前最佳位置;

5)比較所有微粒的pbest和全局gbest值，更新gbest;

6)按照適應值高低將整個粒子群排序，將群體中最差的一半微粒以最好的一半替換，并保持pbest和gbest不改變;

7)達到迭代停止條件，則終止搜索，輸出結果，否則返回3)更新繼續搜索。

2 懸臂梁模型實驗

采用Q235工字鋼，將其一端固接形成懸臂梁，在距離固定螺栓中心點200 mm處有一損傷。已知彈性模量 E=2.1 ×105MPa，I=0.137 ×10－5m4，泊松比μ =0.3，密度 ρ=7 800 kg/m3，幾何尺寸如圖2。

圖2 懸臂梁模型Fig.2 Model of cantilever beam

2.1 動力特性測試

在懸臂端設置加速度傳感器，采用錘擊法測試該懸臂梁的橫向、豎向前兩階固有頻率，見圖3。

圖3 懸臂梁振動實驗Fig.3 Vibration test of cantilever beam

實驗測得加速度時間歷程曲線如圖4。

圖4 加速度-時間曲線Fig.4 Acceleration-time curves

經FFT變換可識別出結構的豎向前兩階固有頻率。建立理想懸臂梁有限元模型，計算豎向前兩階固有頻率，并與實測結果相對比:1階固有頻率的實測值 =166.26 Hz，計算值 =191.36 Hz;2 階固有頻率的實測值 =966.50 Hz，計算值 =993.80 Hz。

2.2 有限元模型改進

結合懸臂梁實際狀況，發現有以下3個非理想因素:

1)豎向激振時，約束端隨梁段共同運動，可認為固定端的豎向轉動約束非理想，可用扭轉彈簧約束代替;

2)距離螺栓中心點20 cm處存在損傷，采用單元抗彎剛度降低予以考慮;

3)模型計算長度為螺栓中心點到自由端距離，而實際懸臂梁有效長度偏短，具體數值難以準確測定。

考慮以上因素改進計算模型如圖5。

圖5 計算模型Fig.5 Calculation model

根據計算模型，使用ANSYS建立改進后的有限元模型如圖6。

圖6 有限元模型Fig.6 Finite element model

選取扭轉彈簧剛度、損傷單元彈性模量、計算長度為設計參數，在可能的取值范圍內進行多次試算，得到其變化對各階頻率的影響規律如圖7。

圖7 設計參數對頻率的影響Fig.7 Influence of the design parameters on frequency

2.3 模型修正

選用2階式-高斯組合徑向基函數為響應面函數，根據2.2章節試算所得各因素對頻率的影響程度，選定設計參數的合理取值范圍如表1。

表1 設計參數修改范圍Table 1 Value range of the design parameters

選用U30(303)均勻設計表形成3因素30水平均勻試驗，共計30個試驗樣本［8］。并利用逐步消去法篩除顯著性低的自變量，擬合得到橫向、豎向前兩階固有頻率回歸響應面，各階響應面待定系數回歸值如表2。

表2 待定系數最小二乘解Table 2 Least square solution of parameters

各階響應面有效性檢驗結果如表3。

表3 響應面有效性檢驗Table 3 Validity tests of the response surfaces

結合各階響應面與實測結果構造殘差型目標函數如式(7):

采用基于自然選擇的混合粒子群算法對目標函數進行求解，迭代過程如圖8。

圖8 優化求解迭代曲線Fig.8 Iteration curves of the optimization

將所選設計參數解除歸一化，得到優化求解后的參數取值如表4。

表4 優化求解結果Table 4 Optimization results

優化后的參數值均在合理范圍內，將其代回有限元模型，再次計算結構模態頻率，與實測結果比較，相對誤差見表5。

表5 修正后的豎向模態頻率Table 5 Vertical modal frequencies of updated model /Hz

修正后的有限元模型豎向模態頻率計算值與實測結果更加接近，能夠更準確地反映實際結構的真實狀態，表明該模型修正方法在實際應用中具有較好可行性。

3 貝雷梁鋼橋有限元模型修正

3.1 工程概況

四川省某水電站建設工程臨時用橋梁主要用于現場公路臨時改線后造成的臨時公路運輸。該鋼便橋設計使用年限＜8年，設計荷載為汽-20。全橋為平橋設計，單車道，全橋長 45.3 m，寬 3.6 m，右岸有26 m長的引堤。鋼橋上部結構采用型鋼結構，由縱橫向分配梁I 14、I 25a和10 mm厚鋼板組成。I 14間距35 cm，I 25 a間距1.5 m。主縱梁選用“321”型貝雷梁、45 m主跨、采用三排雙層加強貝雷梁組合型式構成，兩端簡支。橋梁總體布置如圖9。

圖9 總體布置Fig.9 General arrangement drawing

3.2 環境振動測試

采用環境隨機激勵法測試橋梁的自振特性，拾振器采用4只DSP0.5-2-H(V)型超低頻位移傳感器，拾振范圍≥0.1 Hz;動態數據采集使用美國產WaveBook216a型信號采集分析系統，橋跨結構典型橫向和豎向脈動響應頻譜圖見圖10。

該橋豎向1階頻率為3.36 Hz，橫向前兩階頻率分別為 2.17，3.69 Hz。

圖10 實測頻譜Fig.10 Measured spectrogram

3.3 有限元模型改進

依據圖紙建立理想有限元模型，經動力分析得其豎向1階頻率為5.02 Hz，橫向前兩階頻率分別為1.83，5.71 Hz。

結合橋梁實際狀態，發現兩端支點與橋臺間采用橡膠支座，與理想有限元模型中的固定約束不符，造成計算頻率偏大。考慮改用COMBIN14單元模擬支座橫向、豎向剛度，將剛度系數視為設計參數，計算模型如圖11。

圖11 計算模型Fig.11 Calculation model

根據以上計算模型建立如圖12的有限元模型。

圖12 有限元模型Fig.12 Finite element model

3.4 模型修正

將支座的橫向、豎向剛度作為設計參數，根據多次試算所得規律與工程經驗，確定參數取值范圍:橫向剛度的最小值=1 720 N/mm，最大值=5 200 N/mm;豎向剛度的最小值=1 720 N/mm，最大值=5 200 N/mm。選用U30(303)均勻設計表形成3因素30水平均勻試驗，共計30個試驗樣本［8］。

根據試驗樣本數據，采用最小二乘法回歸所得各階響應面系數回歸值如表6。

表6 待定系數最小二乘解Table 6 Least square solution of parameters

各階響應面有效性檢驗值如表7。

表7 響應面有效性檢驗Table 7 Validity tests of the response surfaces

根據式(8)建立目標函數:

通過基于自然選擇的混合粒子群算法對目標函數優化求解，迭代過程如圖13。

圖13 迭代收斂曲線Fig.13 Iteration curves of the optimization

優化求解經歷迭代150次左右達到收斂，得到設計參數支座橫向剛度2 547.11 N/mm，支座豎向剛度3 746.12 N/mm。優化值均在樣本范圍內，將其代回有限元模型，通過模態分析得到豎向一階模態頻率為3.37 Hz，與實測值相差僅0.30%，表明修正后有限元模型的整體動力特性與實際橋梁更加符合。

4.5 動力響應驗證

現場使用25T載重汽車分別以5，10，20 km/h速度通過橋梁。受現場條件限制，車輛到達甘洛端后就地停留，未駛出橋面。實測跨中測點動撓度曲線如圖14。

圖14 跨中實測動撓度曲線Fig.14 The measured dynamic deflection curves of the mid-span

在修正后有限元模型的基礎上，采用ANSYS瞬態動力學分析對橋梁動力響應進行仿真計算，使用與現場試驗相同的荷載，通過橋梁后，保持端部荷載不撤除，模擬現場實測過程，得到跨中動撓度曲線如圖15。

圖15 跨中仿真動撓度曲線Fig.15 The simulation dynamic deflection curves of the mid-span

分別從修正前、修正后的有限元仿真分析結果中提取跨中動撓度最大值，與實測值相比較，相對誤差見表8。

表8 最大動撓度對比Table 8 Comparison of the maximum dynamic deflection

表8表明，經修正后的橋梁有限元模型的動力響應與實際結構相近，表明其能較準確地反映出結構的真實動力特性，可用于日后的優化設計、健康監測、結構動力響應再分析等。

4 結論

1)2階式-高斯組合徑向基函數響應面具有較強的非線性適應性，求解過程中，能準確再現設計參數-響應值之間的復雜關系，利用其可避免有限元模型的反復調用。

2)設計參數的選取應根據工程實際情況而定，必要時可對有限元模型進行局部改進，使其更接近于實際結構，從而產生新的設計參數作為修正對象。

3)文中殘差型目標函數對各自變量平等對待，優化過程中自動優先考慮殘差絕對值較大的自變量，后續對于重要自變量未被優先考慮的情形，可設置加權值予以改善。

4)文中方法修正所得有限元模型能有效再現原結構的動力特性，在初始模型與設計參數合理的前提下，既使實測數據有限，也可得到準確的有限元模型。

［1］ Wei Xinren，Hua Bingchen.Finite element model updating in structural dynamics by using the response surface method［J］.Engineering Structures，2010，32(4):2455-2465.

［2］ Lu Deng，Cai C S.Bridge model updating using response surface method and genetic algorithm ［J］.Journal of Bridge Engineering，2010，15(5):553-564.

［3］ 周林仁，歐進萍.基于徑向基函數響應面方法的大跨度斜拉橋有限元模型修正［J］.中國鐵道科學，2012，33(3):8-15.Zhou Linren，Ou Jinping.Finite element model updating of longspan cable-stayed bridge based on the response surface method of radial basis function ［J］.China Academy of Railway Sciences，2012，33(3):8-15.

［4］ 孔憲仁，秦玉靈，羅文波.基于改進高斯徑向基函數響應面方法的蜂窩板模型修正［J］.復合材料學報，2011，28(5):220-227.Kong Xianren，Qin Yuling，Luo Wenbo.Improved Gaussian RBF response surface method-based model updating for the honeycomb sandwich panel［J］.Acta Materiae Compositae Sinica，2011，28(5):220-227.

［5］ 秦玉靈，孔憲仁，羅文波.基于徑向基函數響應面的機翼有限元模型修正［J］.北京航空航天大學學報，2011，37(11):1465-1470.Qin Yuling，Kong Xianren，Luo Wenbo.Finite element model updating of airplane wing based on Gaussian radial basis function response surface［J］.Journal of Beijing University of Aeronautics and Astronautics，2011，37(11):1465-1470.

［6］ 張松涵.基于組合函數動力響應面的橋梁有限元模型修正［D］.成都:西南交通大學，2013.Zhang Songhan.Bridge Finite Element Model Updating Based on Dynamic Response Surface of Combination Function［D］.Chengdu:Southwest Jiaotong University，2013.

［7］ Tshilidzi Marwala.Finite-Element-Model Updating Using Computational Intelligence Techniques［M］.New York:Springer，2010.

［8］ 方開泰.均勻設計與均勻設計表［M］.北京:科學出版社，1994.Fang Kaitai.Uniform Design and Uniform Design Table［M］.Beijing:Science Press，1994.

［9］ Anula Khare，Saroj Rangnekar.Particle swarm optimization:A review ［J］.Applied Soft Computing Journal，2012(6):467-484.

［10］ 龔純，王正林.精通MATLAB最優化計算［M］.2版.北京:電子工業出版社，2012.Gong Chun，Wang Zhenglin.Proficient in Optimization Calculation of MATLAB［M］.2nd ed.Beijing:Publishing House of Electronics Industry，2012.