基于小波有限元法的連續(xù)梁移動荷載識別

修建權(quán)，尤 瓊，史治宇

(1.中航工業(yè)東安發(fā)動機(集團)有限公司，哈爾濱 150066;2.南京航空航天大學(xué) 機械結(jié)構(gòu)力學(xué)與控制國家重點實驗室，南京 210016)

傳統(tǒng)有限元法(Traditional Finite Element Method，TFEM)理論成熟，原理簡單，并且有強大的商業(yè)軟件支持，在工程問題的數(shù)值模擬中占據(jù)著重要地位。近些年來人們又提出了一些新型的有限元方法，包括小波有限元法(Wavelet Finite Element Method，WFEM)。Ko［1］最早正式提出了小波有限元，構(gòu)造的WFEM單元主要用于求解簡單的數(shù)學(xué)方程。臺灣學(xué)者Chen［2-3］則利用樣條小波單元，分析了工程結(jié)構(gòu)的振動問題。西安交大向家偉構(gòu)造了區(qū)間B樣條小波(B-spline Wavelet on the Interval，BSWI)有限元［4］，并對各種結(jié)構(gòu)進(jìn)行振動分析和簡單的損傷診斷［5］。但總體而言，基于這種方法的模型能否用于復(fù)雜模型，如連續(xù)梁結(jié)構(gòu)，以便真正應(yīng)用到實際工程中去，仍需要進(jìn)一步研究。

針對上述原因，本文借助BSWI小波，采用了基于WFEM的連續(xù)梁車橋系統(tǒng)模型。在實際工程問題中，加速度傳感器使用范圍最廣，價格相對低廉且穩(wěn)定性高。因此，為提高文中方法的實用性，反問題分析時只采用仿真的測得加速度信號。

移動荷載識別方法繁多［6-10］，小波除了可以與有限元結(jié)合進(jìn)行建模，還可以用于分析識別問題［9-10］。文中僅采用測點的加速度動響應(yīng)數(shù)據(jù)，并借助動態(tài)規(guī)劃技術(shù)［11］與一階 Tikhonov正則化法［12］實現(xiàn)移動荷載識別。仿真結(jié)果表明，基于WFEM的移動荷載識別精度與相同條件下基于TFEM的類似，且所需單元數(shù)較少;一階Tikhonov正則化法不僅可以解決荷載通過約束位置時的振蕩問題，且具有一定的平滑噪聲能力。

1 系統(tǒng)建模

移動荷載作用對象為簡支Bernoulli-Euler連續(xù)梁模型，該模型只考慮縱向動態(tài)響應(yīng)與轉(zhuǎn)角作為梁單元的節(jié)點自由度。以兩個固定間距，且以相同速度從左向右移動的時變力向量表示移動荷載。這種系統(tǒng)模型可以用于分析車橋系統(tǒng)的特性［7，13］。

2 基于小波的連續(xù)梁車橋系統(tǒng)模型

2.1 轉(zhuǎn)換矩陣的推導(dǎo)

在TFEM中，多采用插值多項式和樣條函數(shù)為插值函數(shù)。本文中的WFEM模型采用BSWI作為插值函數(shù)來構(gòu)造小波空間的小波單元［4］。為滿足邊界的連續(xù)性及兼容性，以及方便引入邊界條件，單元剛度矩陣和質(zhì)量矩陣必須由小波空間轉(zhuǎn)換到物理空間，相應(yīng)的單元自由度從小波系數(shù)轉(zhuǎn)換為未知場函數(shù)。因此，首先引入在區(qū)間B樣條小波單元構(gòu)造中起關(guān)鍵作用的轉(zhuǎn)換矩陣。

圖1 離散求解域ΩiFig.1 Discretization of a solution domain Ωi

因此，單元每個節(jié)點實際坐標(biāo)值可表示為xh∈［x1，xn+1］，1≤h≤n+1。若定義:

上式將x映射成標(biāo)準(zhǔn)求解區(qū)間［0，1］，每個節(jié)點xh的映射值ξh為:

最終，在小波空間張成的未知場函數(shù)可表示為:

由經(jīng)典彎曲梁理論，梁的轉(zhuǎn)角等于梁橫向位移的一階導(dǎo)數(shù)。因此一個小波單元有n+3個自由度，定義如下:

將式(3)中不同節(jié)點的位移w(ξi)分別代入式(4)可得:

式中，矩陣［Re］定義為:

將式(5)代入式(3)，［Re］-1可得:

式中Ne為梁模型小波單元形函數(shù):

矩陣［Re］-1改寫為［Te］=［Re］-1，顯見小波空間與物理空間的關(guān)系為:

上式通過轉(zhuǎn)換矩陣［Te］建立了小波空間的小波尺度函數(shù)Φ與物理空間的小波單元形函數(shù)之間的關(guān)系。

2.2 基于小波的有限元連續(xù)梁及移動荷載構(gòu)造

對于一根單元長度le=b-a的Bernoulli-Euler梁，根據(jù)梁彎曲單元勢能泛函，并由變分原理，可得到單元求解方程:

式中，單元剛度矩陣、分布荷載作用下、集中荷載作用下及集中彎矩作用下的荷載列陣分別如下

其中，EI為抗彎剛度;q(x)為分布荷載;Pj為集中荷載;Mk為集中彎矩;

根據(jù)Bernoulli-Euler梁自由振動問題的勢能泛函及單元剛度矩陣，同樣經(jīng)由變分后，可得單元一致質(zhì)量矩陣為:

式中，ρ為材料密度，A為梁橫截面面積。

上式與式(9)分別給出了梁單元基于小波空間BSWI的各物理參數(shù)，則整個梁單元中的物理參數(shù)以及作用荷載的構(gòu)造方式與TFEM一致

其中s表示總體單元數(shù);［C］為Rayleigh阻尼矩陣，α和β由前幾階固有頻率決定;{f}為真實的移動荷載列向量;［Y］為移動荷載的時變位置矩陣;［M］，［C］，［K］分別為總體質(zhì)量矩陣、阻尼矩陣與剛度矩陣。

3 移動荷載識別

本文的研究目標(biāo)在于識別出與實測力情況最吻合的力向量{f}。由于實際中，只有某一部分的加速度響應(yīng)能被測得，因此只有對應(yīng)部分的{X}能被求出。此處，對應(yīng)部分{X}的仿真值可定義為{Z}2ns×1=［Q］2ns×2nn{X}2nn×1。式中，ns個測得點數(shù)通常遠(yuǎn)小于系統(tǒng)的自由度數(shù)。這樣，上述問題可通過非線性最小二乘解法，即Tikhonov正則化解決。

其中，(x，y)表示兩向量x與y的內(nèi)積;［A］，［B］同為權(quán)矩陣。［A］通常為一ns×ns的單位矩陣，［B］為一nf×nf的對角矩陣，該矩陣含有最優(yōu)正則化參數(shù)λ，可平滑待識別量。}表示數(shù)值計算中含有縱向位移與速度的狀態(tài)向量，N表示時間步總長。

一階Tikhonov正則化將求解待識別量{f}的問題轉(zhuǎn)換為求解待識別量的一階導(dǎo)數(shù){}。則式(12)可變換為:

4 仿真算例分析

算例中，將連續(xù)橋梁模型分別分為12個TFEM單元，3個WFEM(m=4，j=3)單元進(jìn)行移動荷載識別，模型如圖2所示。已知連續(xù)梁模型總長L=60 m，等分為三跨，并采用鉸支支撐;連續(xù)梁模型抗彎剛度為EI=1.274 9×1011N·m2;沿梁軸向的線密度為ρA=12 000 kg/m;前三階阻尼比均為ξi=0.02，且前三階頻率分別為f1=12.8 Hz，f2=16.4 Hz，f3=24.0 Hz;采用 Reyleigh阻尼比，根據(jù)阻尼比與頻率，可獲得Reyleigh阻尼系數(shù)。

圖2 移動荷載作用下的連續(xù)梁模型Fig.2 Vehicle-Bridge model under moving forces

仿真移動荷載速度為v=32 m/s;兩個移動荷載的固定間距為Vd=4 m;并假定待識別的移動荷載形式如下，由一恒載部分與兩高頻部分組成:

分析中，測得的加速度動態(tài)響應(yīng)數(shù)據(jù)包括有噪聲及無噪聲情況，由Newmark法仿真算出。沿連續(xù)梁模型，在邊緣跨的1/4與3/4跨長處各布置兩個測點，中間跨則等間距布置3個測點，共計7個測點用于采集加速度信號;時間步數(shù)N取決于作用在連續(xù)橋梁模型上的時間總長T及數(shù)據(jù)的采樣頻率。

反問題求解中，求逆導(dǎo)致的數(shù)值計算病態(tài)問題一般采用Tikhonov正則化方法求解。表1所示為TFEM 12個單元與WFEM 3個單元模型下，加速度信號在有無噪聲影響時，移動荷載的識別結(jié)果。

表1 Tikhonov正則化下的移動荷載識別結(jié)果Tab.1 Moving force identification with Tikhonov regularization

由表1可知，WFEM 3個單元下的識別結(jié)果與TFEM 12個單元下的識別情況非常接近。但是兩種識別結(jié)果對噪聲都相當(dāng)敏感，當(dāng)噪聲由1%增加至5%時，識別的移動荷載RPE值就超出了可接受范圍，如圖3(a)中所示，且識別結(jié)果在支撐處有較大振蕩。針對此問題，文中采用一階Tikhonov正則化法來削弱噪聲信號對移動荷載識別效果的影響，如圖3(b)所示。顯而易見:一階Tikhonov正則化法下，在移動荷載通過支撐處的瞬時時刻，識別結(jié)果并無振蕩現(xiàn)象發(fā)生;并且噪聲沿整個時程的影響明顯被平滑了。

由表2可見，與TFEM相比，在WFEM下需要相對稍長的計算時間獲得類似的識別精度;一階Tikhonov正則化法比Tikhonov正則化法需消耗稍多的計算時間。但總體而言，消耗的計算時間都在可接受范圍內(nèi)。

圖3 不同正則化法下的移動荷載識別結(jié)果Fig.3 Moving force identification with different regularization methods

表2 各模型及正則化下的計算時間Tab.2 Computation time required for different regularization methods and models

5 結(jié)論

本文采用測點的加速度動響應(yīng)數(shù)據(jù)，驗證了WFEM模型用于移動荷載識別的可行性。同時，采用一階Tikhonov正則化法識別移動荷載。仿真結(jié)果表明，與TFEM相比，WFEM只需要較少的單元即可獲得類似的識別精度;一階Tikhonov正則化法對噪聲具有較好的平滑性能，且可避免時程分析中移動荷載在支撐處的振蕩現(xiàn)象，但消耗的計算時間相對較多。

［1］ Ko J，Kurdila A J，Pilant M S.A class of finite element methods based on orthonormal，compactly supported wavelets［J］.Computational Mechanics，1995，16(4):235-244.

［2］ Chen W H，Wu C W.Spline wavelets element method for frame structures vibration［J］.Computational Mechanics，1995，16(1):11-21.

［3］Chen W H，Wu C W.Extension of spline wavelets element method to membrane vibration analysis［J］.Computational Mechanics，1996，18(1):46-54.

［4］何正嘉，陳雪峰，李 兵，等.小波有限元理論及其工程應(yīng)用［M］.北京:科學(xué)出版社，2006.

［5］Xiang J W，Chen X F，Mo Q Y，et al.Identification of crack in a rotor system based on wavelet finite element method［J］.Finite elements in analysis and design，2007，43(14):1068-1081.

［6］ Yu L.Accounting for bridge dynamic loads using Moving Force Identification System(MFIS)［D］.Hong Kong:The Hong Kong Polytechnic University，2002.

［7］ Yu L，Chan T H T.Recent research on identification of moving loads on bridges［J］.Journal of Sound and Vibration，2007，305(1-2):3-21.

［8］毛玉明，郭杏林，趙 巖，等.基于靈敏度分析的結(jié)構(gòu)動態(tài)載荷識別研究［J］.振動與沖擊，2010，29(10):1-3.

［9］黃 林，袁向榮.小波分析在橋上移動荷載識別中的應(yīng)用［J］.鐵道學(xué)報，2003，25(4):97-101.

［10］吳紹慶，史治宇.由有限元-Wavelet-Galerkin法識別橋面移動載荷［J］.振動工程學(xué)報，2006，19(4):494-498.

［11］Law S S，F(xiàn)ang Y L.Moving force identification:optimal state estimation approach［J］.Journal of Sound and Vibration，2001，239(2):233-254.

［12］ Busby H R，Trujillo D M.Optimal regularization of an inverse dynamic problem［J］.Computers and Structures，1997，63(2):243-248.

［13］ Chan T H T，Ashebo D B.Theoretical study of moving force identification on continuous bridges［J］.Journal of Sound and Vibration，2006，295(3-5):870-883.