修建權,尤 瓊,史治宇
(1.中航工業(yè)東安發(fā)動機(集團)有限公司,哈爾濱 150066;2.南京航空航天大學 機械結(jié)構力學與控制國家重點實驗室,南京 210016)
傳統(tǒng)有限元法(Traditional Finite Element Method,TFEM)理論成熟,原理簡單,并且有強大的商業(yè)軟件支持,在工程問題的數(shù)值模擬中占據(jù)著重要地位。近些年來人們又提出了一些新型的有限元方法,包括小波有限元法(Wavelet Finite Element Method,WFEM)。Ko[1]最早正式提出了小波有限元,構造的WFEM單元主要用于求解簡單的數(shù)學方程。臺灣學者Chen[2-3]則利用樣條小波單元,分析了工程結(jié)構的振動問題。西安交大向家偉構造了區(qū)間B樣條小波(B-spline Wavelet on the Interval,BSWI)有限元[4],并對各種結(jié)構進行振動分析和簡單的損傷診斷[5]。但總體而言,基于這種方法的模型能否用于復雜模型,如連續(xù)梁結(jié)構,以便真正應用到實際工程中去,仍需要進一步研究。
針對上述原因,本文借助BSWI小波,采用了基于WFEM的連續(xù)梁車橋系統(tǒng)模型。在實際工程問題中,加速度傳感器使用范圍最廣,價格相對低廉且穩(wěn)定性高。因此,為提高文中方法的實用性,反問題分析時只采用仿真的測得加速度信號。
移動荷載識別方法繁多[6-10],小波除了可以與有限元結(jié)合進行建模,還可以用于分析識別問題[9-10]。文中僅采用測點的加速度動響應數(shù)據(jù),并借助動態(tài)規(guī)劃技術[11]與一階 Tikhonov正則化法[12]實現(xiàn)移動荷載識別。仿真結(jié)果表明,基于WFEM的移動荷載識別精度與相同條件下基于TFEM的類似,且所需單元數(shù)較少;一階Tikhonov正則化法不僅可以解決荷載通過約束位置時的振蕩問題,且具有一定的平滑噪聲能力。
移動荷載作用對象為簡支Bernoulli-Euler連續(xù)梁模型,該模型只考慮縱向動態(tài)響應與轉(zhuǎn)角作為梁單元的節(jié)點自由度。以兩個固定間距,且以相同速度從左向右移動的時變力向量表示移動荷載。這種系統(tǒng)模型可以用于分析車橋系統(tǒng)的特性[7,13]。
在TFEM中,多采用插值多項式和樣條函數(shù)為插值函數(shù)。本文中的WFEM模型采用BSWI作為插值函數(shù)來構造小波空間的小波單元[4]。為滿足邊界的連續(xù)性及兼容性,以及方便引入邊界條件,單元剛度矩陣和質(zhì)量矩陣必須由小波空間轉(zhuǎn)換到物理空間,相應的單元自由度從小波系數(shù)轉(zhuǎn)換為未知場函數(shù)。因此,首先引入在區(qū)間B樣條小波單元構造中起關鍵作用的轉(zhuǎn)換矩陣。
圖1 離散求解域ΩiFig.1 Discretization of a solution domain Ωi
因此,單元每個節(jié)點實際坐標值可表示為xh∈[x1,xn+1],1≤h≤n+1。若定義:
上式將x映射成標準求解區(qū)間[0,1],每個節(jié)點xh的映射值ξh為:
最終,在小波空間張成的未知場函數(shù)可表示為:
由經(jīng)典彎曲梁理論,梁的轉(zhuǎn)角等于梁橫向位移的一階導數(shù)。因此一個小波單元有n+3個自由度,定義如下:
將式(3)中不同節(jié)點的位移w(ξi)分別代入式(4)可得:
式中,矩陣[Re]定義為:
將式(5)代入式(3),[Re]-1可得:
式中Ne為梁模型小波單元形函數(shù):
矩陣[Re]-1改寫為[Te]=[Re]-1,顯見小波空間與物理空間的關系為:
上式通過轉(zhuǎn)換矩陣[Te]建立了小波空間的小波尺度函數(shù)Φ與物理空間的小波單元形函數(shù)之間的關系。
對于一根單元長度le=b-a的Bernoulli-Euler梁,根據(jù)梁彎曲單元勢能泛函,并由變分原理,可得到單元求解方程:
式中,單元剛度矩陣、分布荷載作用下、集中荷載作用下及集中彎矩作用下的荷載列陣分別如下
其中,EI為抗彎剛度;q(x)為分布荷載;Pj為集中荷載;Mk為集中彎矩;
根據(jù)Bernoulli-Euler梁自由振動問題的勢能泛函及單元剛度矩陣,同樣經(jīng)由變分后,可得單元一致質(zhì)量矩陣為:
式中,ρ為材料密度,A為梁橫截面面積。
上式與式(9)分別給出了梁單元基于小波空間BSWI的各物理參數(shù),則整個梁單元中的物理參數(shù)以及作用荷載的構造方式與TFEM一致
其中s表示總體單元數(shù);[C]為Rayleigh阻尼矩陣,α和β由前幾階固有頻率決定;{f}為真實的移動荷載列向量;[Y]為移動荷載的時變位置矩陣;[M],[C],[K]分別為總體質(zhì)量矩陣、阻尼矩陣與剛度矩陣。
本文的研究目標在于識別出與實測力情況最吻合的力向量{f}。由于實際中,只有某一部分的加速度響應能被測得,因此只有對應部分的{X}能被求出。此處,對應部分{X}的仿真值可定義為{Z}2ns×1=[Q]2ns×2nn{X}2nn×1。式中,ns個測得點數(shù)通常遠小于系統(tǒng)的自由度數(shù)。這樣,上述問題可通過非線性最小二乘解法,即Tikhonov正則化解決。
其中,(x,y)表示兩向量x與y的內(nèi)積;[A],[B]同為權矩陣。[A]通常為一ns×ns的單位矩陣,[B]為一nf×nf的對角矩陣,該矩陣含有最優(yōu)正則化參數(shù)λ,可平滑待識別量。}表示數(shù)值計算中含有縱向位移與速度的狀態(tài)向量,N表示時間步總長。
一階Tikhonov正則化將求解待識別量{f}的問題轉(zhuǎn)換為求解待識別量的一階導數(shù){}。則式(12)可變換為:
算例中,將連續(xù)橋梁模型分別分為12個TFEM單元,3個WFEM(m=4,j=3)單元進行移動荷載識別,模型如圖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)響應數(shù)據(jù)包括有噪聲及無噪聲情況,由Newmark法仿真算出。沿連續(xù)梁模型,在邊緣跨的1/4與3/4跨長處各布置兩個測點,中間跨則等間距布置3個測點,共計7個測點用于采集加速度信號;時間步數(shù)N取決于作用在連續(xù)橋梁模型上的時間總長T及數(shù)據(jù)的采樣頻率。
反問題求解中,求逆導致的數(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é)果對噪聲都相當敏感,當噪聲由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
本文采用測點的加速度動響應數(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]何正嘉,陳雪峰,李 兵,等.小波有限元理論及其工程應用[M].北京:科學出版社,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é)構動態(tài)載荷識別研究[J].振動與沖擊,2010,29(10):1-3.
[9]黃 林,袁向榮.小波分析在橋上移動荷載識別中的應用[J].鐵道學報,2003,25(4):97-101.
[10]吳紹慶,史治宇.由有限元-Wavelet-Galerkin法識別橋面移動載荷[J].振動工程學報,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.