張松涵，高芳清，張偉偉

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

0 引言

近年來，基于響應(yīng)面的有限元模型修正方法在工程中已得到廣泛應(yīng)用，其方法在于使用顯式函數(shù)關(guān)系逼近設(shè)計(jì)參數(shù)與響應(yīng)值之間的隱式復(fù)雜關(guān)系，與傳統(tǒng)有限元模型修正相比，有效避免了迭代過程中重復(fù)調(diào)用有限元計(jì)算的缺點(diǎn)，使迭代速度大大加快［1］。

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

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

為驗(yàn)證2階式-高斯組合徑向基函數(shù)響應(yīng)面在實(shí)際工程中應(yīng)用的有效性，筆者通過模型修改獲得合理的設(shè)計(jì)參數(shù)，在重復(fù)試驗(yàn)的基礎(chǔ)上建立了模態(tài)頻率響應(yīng)面函數(shù)，構(gòu)造殘差型目標(biāo)函數(shù)，利用粒子群優(yōu)化算法求解，得到設(shè)計(jì)參數(shù)最優(yōu)解。在保證求解精度不受影響的前提下，避免了求解過程中反復(fù)調(diào)用有限元分析與靈敏度矩陣計(jì)算，適用于結(jié)構(gòu)形式復(fù)雜的橋梁有限元模型修正。

1 基于響應(yīng)面的有限元模型修正方法

基于響應(yīng)面的有限元模型修正通過重復(fù)試驗(yàn)，將復(fù)雜結(jié)構(gòu)的有限元模型轉(zhuǎn)化為簡單的響應(yīng)面數(shù)學(xué)模型［7］，結(jié)合實(shí)測數(shù)據(jù)，采用優(yōu)化算法尋求到響應(yīng)面模型中的最優(yōu)解，再代入有限元模型進(jìn)行驗(yàn)證，從而達(dá)到模型修正的目的，其具體流程見圖1。

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

1.1 試驗(yàn)設(shè)計(jì)

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

1.2 組合函數(shù)響應(yīng)面

徑向基函數(shù)是一種徑向?qū)ΨQ的函數(shù)，形式簡單、不依賴于空間維數(shù)，且具有良好的非線性擬合能力。構(gòu)造RBF的線形組合模型形式如式(1):

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

將工程中普遍常用的2階多項(xiàng)式與高斯徑向基函數(shù)結(jié)合，得到改進(jìn)后的響應(yīng)面函數(shù)如式(2):

式中:核寬度參數(shù)σ一般取設(shè)計(jì)空間半徑［4］。

1.3 經(jīng)典最小二乘法回歸響應(yīng)面

將樣本參數(shù)取值視為輸入，計(jì)算樣本頻率視為輸出，建立關(guān)于待定系數(shù)向量的方程:

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

要使殘差平方和

最小，對θ求偏導(dǎo)，令其等于0，即:

得到:

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

1.4 基于自然選擇的混合粒子群優(yōu)化算法

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

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

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

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

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

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

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

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

2 懸臂梁模型實(shí)驗(yàn)

采用Q235工字鋼，將其一端固接形成懸臂梁，在距離固定螺栓中心點(diǎn)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 動力特性測試

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

圖3 懸臂梁振動實(shí)驗(yàn)Fig.3 Vibration test of cantilever beam

實(shí)驗(yàn)測得加速度時間歷程曲線如圖4。

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

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

2.2 有限元模型改進(jìn)

結(jié)合懸臂梁實(shí)際狀況，發(fā)現(xiàn)有以下3個非理想因素:

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

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

3)模型計(jì)算長度為螺栓中心點(diǎn)到自由端距離，而實(shí)際懸臂梁有效長度偏短，具體數(shù)值難以準(zhǔn)確測定。

考慮以上因素改進(jìn)計(jì)算模型如圖5。

圖5 計(jì)算模型Fig.5 Calculation model

根據(jù)計(jì)算模型，使用ANSYS建立改進(jìn)后的有限元模型如圖6。

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

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

圖7 設(shè)計(jì)參數(shù)對頻率的影響Fig.7 Influence of the design parameters on frequency

2.3 模型修正

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

表1 設(shè)計(jì)參數(shù)修改范圍Table 1 Value range of the design parameters

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

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

各階響應(yīng)面有效性檢驗(yàn)結(jié)果如表3。

表3 響應(yīng)面有效性檢驗(yàn)Table 3 Validity tests of the response surfaces

結(jié)合各階響應(yīng)面與實(shí)測結(jié)果構(gòu)造殘差型目標(biāo)函數(shù)如式(7):

采用基于自然選擇的混合粒子群算法對目標(biāo)函數(shù)進(jìn)行求解，迭代過程如圖8。

圖8 優(yōu)化求解迭代曲線Fig.8 Iteration curves of the optimization

將所選設(shè)計(jì)參數(shù)解除歸一化，得到優(yōu)化求解后的參數(shù)取值如表4。

表4 優(yōu)化求解結(jié)果Table 4 Optimization results

優(yōu)化后的參數(shù)值均在合理范圍內(nèi)，將其代回有限元模型，再次計(jì)算結(jié)構(gòu)模態(tài)頻率，與實(shí)測結(jié)果比較，相對誤差見表5。

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

修正后的有限元模型豎向模態(tài)頻率計(jì)算值與實(shí)測結(jié)果更加接近，能夠更準(zhǔn)確地反映實(shí)際結(jié)構(gòu)的真實(shí)狀態(tài)，表明該模型修正方法在實(shí)際應(yīng)用中具有較好可行性。

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

3.1 工程概況

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

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

3.2 環(huán)境振動測試

采用環(huán)境隨機(jī)激勵法測試橋梁的自振特性，拾振器采用4只DSP0.5-2-H(V)型超低頻位移傳感器，拾振范圍≥0.1 Hz;動態(tài)數(shù)據(jù)采集使用美國產(chǎn)WaveBook216a型信號采集分析系統(tǒng)，橋跨結(jié)構(gòu)典型橫向和豎向脈動響應(yīng)頻譜圖見圖10。

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

圖10 實(shí)測頻譜Fig.10 Measured spectrogram

3.3 有限元模型改進(jìn)

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

結(jié)合橋梁實(shí)際狀態(tài)，發(fā)現(xiàn)兩端支點(diǎn)與橋臺間采用橡膠支座，與理想有限元模型中的固定約束不符，造成計(jì)算頻率偏大?？紤]改用COMBIN14單元模擬支座橫向、豎向剛度，將剛度系數(shù)視為設(shè)計(jì)參數(shù)，計(jì)算模型如圖11。

圖11 計(jì)算模型Fig.11 Calculation model

根據(jù)以上計(jì)算模型建立如圖12的有限元模型。

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

3.4 模型修正

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

根據(jù)試驗(yàn)樣本數(shù)據(jù)，采用最小二乘法回歸所得各階響應(yīng)面系數(shù)回歸值如表6。

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

各階響應(yīng)面有效性檢驗(yàn)值如表7。

表7 響應(yīng)面有效性檢驗(yàn)Table 7 Validity tests of the response surfaces

根據(jù)式(8)建立目標(biāo)函數(shù):

通過基于自然選擇的混合粒子群算法對目標(biāo)函數(shù)優(yōu)化求解，迭代過程如圖13。

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

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

4.5 動力響應(yīng)驗(yàn)證

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

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

在修正后有限元模型的基礎(chǔ)上，采用ANSYS瞬態(tài)動力學(xué)分析對橋梁動力響應(yīng)進(jìn)行仿真計(jì)算，使用與現(xiàn)場試驗(yàn)相同的荷載，通過橋梁后，保持端部荷載不撤除，模擬現(xiàn)場實(shí)測過程，得到跨中動撓度曲線如圖15。

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

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

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

表8表明，經(jīng)修正后的橋梁有限元模型的動力響應(yīng)與實(shí)際結(jié)構(gòu)相近，表明其能較準(zhǔn)確地反映出結(jié)構(gòu)的真實(shí)動力特性，可用于日后的優(yōu)化設(shè)計(jì)、健康監(jiān)測、結(jié)構(gòu)動力響應(yīng)再分析等。

4 結(jié)論

1)2階式-高斯組合徑向基函數(shù)響應(yīng)面具有較強(qiáng)的非線性適應(yīng)性，求解過程中，能準(zhǔn)確再現(xiàn)設(shè)計(jì)參數(shù)-響應(yīng)值之間的復(fù)雜關(guān)系，利用其可避免有限元模型的反復(fù)調(diào)用。

2)設(shè)計(jì)參數(shù)的選取應(yīng)根據(jù)工程實(shí)際情況而定，必要時可對有限元模型進(jìn)行局部改進(jìn)，使其更接近于實(shí)際結(jié)構(gòu)，從而產(chǎn)生新的設(shè)計(jì)參數(shù)作為修正對象。

3)文中殘差型目標(biāo)函數(shù)對各自變量平等對待，優(yōu)化過程中自動優(yōu)先考慮殘差絕對值較大的自變量，后續(xù)對于重要自變量未被優(yōu)先考慮的情形，可設(shè)置加權(quán)值予以改善。

4)文中方法修正所得有限元模型能有效再現(xiàn)原結(jié)構(gòu)的動力特性，在初始模型與設(shè)計(jì)參數(shù)合理的前提下，既使實(shí)測數(shù)據(jù)有限，也可得到準(zhǔn)確的有限元模型。

［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ìn)萍.基于徑向基函數(shù)響應(yīng)面方法的大跨度斜拉橋有限元模型修正［J］.中國鐵道科學(xué)，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ìn)高斯徑向基函數(shù)響應(yīng)面方法的蜂窩板模型修正［J］.復(fù)合材料學(xué)報(bào)，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］ 秦玉靈，孔憲仁，羅文波.基于徑向基函數(shù)響應(yīng)面的機(jī)翼有限元模型修正［J］.北京航空航天大學(xué)學(xué)報(bào)，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］ 張松涵.基于組合函數(shù)動力響應(yīng)面的橋梁有限元模型修正［D］.成都:西南交通大學(xué)，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］ 方開泰.均勻設(shè)計(jì)與均勻設(shè)計(jì)表［M］.北京:科學(xué)出版社，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最優(yōu)化計(jì)算［M］.2版.北京:電子工業(yè)出版社，2012.Gong Chun，Wang Zhenglin.Proficient in Optimization Calculation of MATLAB［M］.2nd ed.Beijing:Publishing House of Electronics Industry，2012.