郭 力，李兆霞，韓曉林

(1.東南大學(xué) 土木工程學(xué)院，南京 210096;2.江蘇省工程力學(xué)分析重點實驗室(東南大學(xué))，南京 210096)

自20世紀(jì)60年代以來，國內(nèi)外廣大的工程師一直在追求建立精確的有限元模型來研究結(jié)構(gòu)的動力響應(yīng)問題。但是，由于在有限元建模過程中對結(jié)構(gòu)的材料特性、構(gòu)件聯(lián)結(jié)形式以及外部作用進行了理想化，使得依據(jù)設(shè)計圖紙而建立的有限元模型與實際結(jié)構(gòu)在構(gòu)件尺寸甚至結(jié)構(gòu)形態(tài)上存在一定的偏差。為了提高模型的計算精度，必須對初始模型進行修正［1］。

模型修正研究工作始于1965年Guyan［2］發(fā)表的模型縮聚方面的文章，經(jīng)過幾十年的發(fā)展，現(xiàn)有的模型修正理論包括參考基準(zhǔn)法、混合矩陣法、特征結(jié)構(gòu)匹配法、特征值取逆法、誤差矩陣法，貝葉斯估計等。由于修正目標(biāo)的不同，各種模型修正方法間差異顯著，但就具體的模型修正過程來看，模型修正方法總體上可以分為矩陣型和參數(shù)型兩類。關(guān)于這些經(jīng)典的模型修正方法，文獻［3］中有詳細的評述。

根據(jù)修正過程中依據(jù)的基本測量信息，模型修正又可分為動力模型修正和靜力模型修正［4］。動力模型主要依據(jù)結(jié)構(gòu)的速度、加速度、頻率和振型等動態(tài)測試信息進行修正，由于這些信息能較好地反映結(jié)構(gòu)的整體力學(xué)行為，因此修正大型工程結(jié)構(gòu)模型時較多采用動態(tài)測量信息。但是，動態(tài)測量信息在數(shù)據(jù)采集、傳輸和儲存過程中容易受到噪聲污染，修正結(jié)果的精度需要進一步考證。靜力模型修正采用的基礎(chǔ)數(shù)據(jù)為結(jié)構(gòu)在靜力載荷下的位移、應(yīng)力等信息。由于結(jié)構(gòu)靜力響應(yīng)測量簡單，且精度很高，因此，一般認為靜力模型修正過程比動力模型修正過程精度高。

為了實現(xiàn)基于動態(tài)測量信息的結(jié)構(gòu)模型修正，同時保持靜力模型修正過程精度高的優(yōu)點，本文研究如何把動力模型等效為靜力模型進行修正。同時考察有限測量信息的擴展問題，分析測量誤差對模型修正精度的影響程度，以建立依據(jù)動態(tài)測量信息和靜力模型修正技術(shù)的擬靜力模型修正方法。下面先介紹靜力模型修正過程及基于復(fù)變量求導(dǎo)法的參數(shù)靈敏度分析方法，再研究動力模型的擬靜力修正過程，進一步對修正結(jié)果的可靠性進行統(tǒng)計分析，最后對本文提出的方法進行數(shù)值驗證，并對實際橋塔的模型進行修正。

1 靜力模型修正及靈敏度分析

結(jié)構(gòu)的靜力響應(yīng)是結(jié)構(gòu)或材料參數(shù)的函數(shù)，即:

這里g為結(jié)構(gòu)的靜力響應(yīng)，如位移、應(yīng)力或應(yīng)變等;x為計算模型中的參數(shù)，如材料的彈性模量、泊松比、密度和外載荷等。關(guān)于計算模型中單元間聯(lián)結(jié)條件的影響，這里暫不考慮。

給定設(shè)計參數(shù)x0后，利用有限元模型由式(1)可以計算出結(jié)構(gòu)的響應(yīng)g0=f(x0)。另一方面，結(jié)構(gòu)響應(yīng)可以通過測試得到，記響應(yīng)的測試值為g*。一般情況下，由于初始有限元模型存在一些待修正的誤差，結(jié)構(gòu)響應(yīng)的計算值和測量值間有一定的偏差，記為:

依據(jù)設(shè)計圖紙建立的初始有限元模型，其中的計算參數(shù)x0可以看作是結(jié)構(gòu)真實參數(shù)x*的近似，由式(1)將f(x)在x0處泰勒展開，可得:

由于g*=f(x*)，在上式中取一階近似展開，則式(2)可變?yōu)?

式中 Δx=x*－x0。

令 x={x1，x2，…，xm}T，其中 xi為待修正的參數(shù)，設(shè)測量點共有n個，由式(4)可得:

將式(5)簡記為:

由于測點處的計算值和測量值已知，即Δg已知，利用式(6)求得靈敏度矩陣S后，由最小二乘法可得:

進一步可得:

如果采用疊代求解，則相應(yīng)的疊代模式為:

利用式(10)即可實現(xiàn)結(jié)構(gòu)模型的參數(shù)修正。

在疊代過程中，通過殘差比來判斷計算過程是否收斂，殘差比定義為:

式(5)～式(11)的求解過程即為經(jīng)典的基于靈敏度分析的靜力模型修正過程，其中的關(guān)鍵是如何計算靈敏度矩陣S。由式(6)可知，靈敏度矩陣中的元素為結(jié)構(gòu)響應(yīng)對待修正參數(shù)的一階偏導(dǎo)數(shù)。常規(guī)的基于差分法的靈敏度計算方法精度較低，且計算量較大，當(dāng)結(jié)構(gòu)響應(yīng)對修正參數(shù)不敏感時，采用差分法將難以計算靈敏度。為了解決高精度快速計算靈敏度的問題，文獻［5，6］利用復(fù)變量求導(dǎo)法，建立了復(fù)域靈敏度分析方法，并成功地用于結(jié)構(gòu)非線性參數(shù)反演。因此，這里也采用該方法進行靈敏度分析，下面簡要地給出分析過程。

對復(fù)變量函數(shù)f(x+ih)，當(dāng)h很小時將其泰勒展開可得［7］:

比較上式兩邊的虛部，可得:

當(dāng)h很小時(如10－20)，由上式計算得到的函數(shù)一階偏導(dǎo)數(shù)的精度為10－40，遠高于計算機的機器精度，計算量是常規(guī)差分法的一半以下。

綜合式(13)和式(6)，可得靈敏度矩陣為:

利用式(14)進行靈敏度分析時，只要把結(jié)構(gòu)計算的有限元程序通過變量聲明(IMPLICIT)的辦法很方便地改為復(fù)變量的程序，將待計算靈敏度的參數(shù)加上很小的虛部h，提取計算的結(jié)構(gòu)響應(yīng)的虛部，除以h即得到結(jié)構(gòu)響應(yīng)對該參數(shù)的靈敏度。因此，復(fù)域靈敏度分析過程可以在結(jié)構(gòu)響應(yīng)計算時同時進行，對結(jié)構(gòu)響應(yīng)計算的程序改動量小，分析過程容易實現(xiàn)，在下面的模型修正過程中采用復(fù)變量求導(dǎo)法進行參數(shù)靈敏度分析。

2 擬靜力模型修正過程

從上面的分析可以看出，靜力模型修正過程簡單，且精度較高。因此，如果可以得到結(jié)構(gòu)的靜力響應(yīng)信息，盡可能采用靜力模型修正方法來修正模型。但是，當(dāng)結(jié)構(gòu)的尺寸較大時，難以施加靜載荷，且進行靜力響應(yīng)測量很困難，一般對大型工程結(jié)構(gòu)多采用動力響應(yīng)(如結(jié)構(gòu)的頻率和模態(tài))測量來了解結(jié)構(gòu)的力學(xué)特性。常規(guī)基于動力測量信息的模型修正過程需要構(gòu)造關(guān)于結(jié)構(gòu)頻率或模態(tài)的目標(biāo)函數(shù)，修正過程中必須計算結(jié)構(gòu)的頻率和模態(tài)。計算頻率或模態(tài)將涉及到特征值計算問題，對大型矩陣的特征值求解過程較復(fù)雜，且存在一些固有的誤差，因此，編寫動力模型修正過程的程序較煩瑣。如果能把測量的結(jié)構(gòu)動力響應(yīng)信息經(jīng)過一定的轉(zhuǎn)換，使得可以利用靜力模型修正過程來修正模型，將顯著降低模型修正過程的復(fù)雜程度，這正是本文提出擬靜力模型修正方法的宗旨。下面給出轉(zhuǎn)換過程。

結(jié)構(gòu)自由振動時，有如下方程［8］:

其中ωi為結(jié)構(gòu)的第i階頻率，K為結(jié)構(gòu)整體剛度矩陣，M為結(jié)構(gòu)整體質(zhì)量矩陣，Φi為結(jié)構(gòu)的第i階模態(tài)。

如果在結(jié)構(gòu)上施加一組載荷為:

這里Pi稱為第i階振型的特征載荷。

設(shè)Pi引起的結(jié)構(gòu)位移為δi，則有:

比較式(15)和式(17)可得δi=Φi，即在特征載荷作用下，結(jié)構(gòu)產(chǎn)生的位移等于該振型的特征向量(模態(tài))。這樣，模型的特征值計算過程可以等效為在特征載荷作用下的結(jié)構(gòu)靜力響應(yīng)計算。

在建立結(jié)構(gòu)的初始有限元模型時，結(jié)構(gòu)和材料參數(shù)的誤差往往引起剛度矩陣K有較大的誤差，對質(zhì)量矩陣M的影響較小。因此，在模型修正過程中，認為系統(tǒng)的質(zhì)量矩陣沒有誤差，只對剛度矩陣K進行修正。

當(dāng)結(jié)構(gòu)的頻率和模態(tài)通過測量得到后，可以求得相應(yīng)的特征載荷，即式(17)的右邊已知，進一步可計算出結(jié)構(gòu)的靜位移δi。由于初始模型的剛度矩陣K有一定的誤差，計算得到的位移一般不等于測量的模態(tài)Φi，令:

利用式(5)～式(11)的修正過程，借助上述的復(fù)域靈敏度分析過程，實現(xiàn)結(jié)構(gòu)參數(shù)的修正。

上述模型修正的目標(biāo)是實現(xiàn)計算模態(tài)與測試模態(tài)的高度相關(guān)，即δi=Φi。當(dāng)該目標(biāo)實現(xiàn)后，由式(17)可知，依據(jù)修正后模型計算的頻率值一定高精度逼近其測量值。因此，上述修正過程最終可以實現(xiàn)模態(tài)和頻率的雙重修正目標(biāo)。這樣，基于動力測量信息的結(jié)構(gòu)模型修正可以通過擬靜力模型修正過程來實現(xiàn)。

采用式(16)計算特征載荷時一般需要完備的測量信息，但是，實際測量時結(jié)構(gòu)的某些自由度(如轉(zhuǎn)動自由度)信息很難測量。為了避免轉(zhuǎn)動自由度缺失對結(jié)果的影響，在式(16)中采用集中質(zhì)量矩陣進行計算。這樣，在計算特征載荷時轉(zhuǎn)動自由度的影響自動消除。

對于測量信息與計算信息間不匹配的問題，可以通過模態(tài)擴展的方法進行解決，文獻［1］中關(guān)于此問題有詳細地討論，本文采用樣條插值的辦法進行模態(tài)擴展，模型修正過程中不需要進行特征值和特征向量的求解，有效降低了修正過程中的計算量。

3 誤差分析

模型修正過程中依據(jù)的測量信息在測量和存儲過程中往往會受到噪聲污染，同時由于測量系統(tǒng)自身的缺陷以及外界偶然因素的影響，難免存在一定的測量誤差。因此，分析測量誤差對修正結(jié)果的影響是必要的。

由式(18)和式(7)可得:

假設(shè)測量中的誤差是均值為0的隨機誤差，記為X，把式(19)兩邊關(guān)于X求偏導(dǎo)，可得:

由于隨機誤差來源于測量過程，而δi及S由計算得到，因此有:

將上式代入式(20)可得:

由最小二乘法可得:

依據(jù)此式可以采用統(tǒng)計分析的方法來考察測量誤差對修正結(jié)果的影響［9－10］。

由于這里假設(shè)測量誤差的均值為0，對式(8)兩邊取數(shù)學(xué)期望，有:

進一步可得疊代過程中修正參數(shù)的期望為:

疊代過程中修正參數(shù)的方差為:

由文獻［11］，可知:

在結(jié)構(gòu)測量過程中，頻率的測量精度很高，這里認為頻率的測量值是精確的，僅考慮模態(tài)具有測量誤差的情形，因此有:

將上式結(jié)合式(23)代入式(27)和(28)，進一步再代入式(26)，即可求得疊代修正過程中修正參數(shù)的方差。

4 算例

為了考察上面提出的擬靜力模型修正過程是否有效，考察如圖1所示懸臂梁，其長度為10 m，橫截面為0.1 m×0.1 m。假設(shè)梁只發(fā)生平面內(nèi)彎曲，采用兩節(jié)點平面梁單元進行劃分，各節(jié)點的自由度為3，待修正的模型包括10個單元(圖中用①，②，…，⑩表示)，節(jié)點數(shù)為11個(圖中用1，2，…，11表示)。懸臂粱的材料為非均勻的，其中5～7節(jié)點間材料的彈性模量為E2=180 GPa，其余部分的彈性模量為 E1=240 GPa，泊松比均為0.3，質(zhì)量密度為 2.8 ×103kg/m3。

首先檢驗基于樣條插值的模態(tài)擴展方法是否可靠。取測量自由度遠低于模型的完備自由度，缺失自由度信息通過樣條插值的方法產(chǎn)生。插值的結(jié)果如圖2所示，圖中的振幅為歸一化振幅，將懸臂粱均勻劃分為100個單元后計算的模態(tài)近似地看作真實模態(tài)(“exact curve”)。圖2(a)中的測量值為一階振型中4、7、11三個節(jié)點處的橫向振動信息，圖2(b)中的測量值為二階振型中3、6、9、11四個節(jié)點處的橫向振動信息?？梢钥闯觯逯档玫降哪B(tài)與真實模態(tài)吻合很好，因此可以認為基于樣條插值的模態(tài)擴展方法是可行的。

圖1 懸臂梁及其有限元模型Fig.1 Cantilever and its FE model

圖2 樣條插值模態(tài)與測量模態(tài)Fig.2 Spline interpolated modals and measured ones

利用樣條插值擴展得到模態(tài)后，可進一步利用式(16)計算特征載荷。但是，這里的模態(tài)為歸一化模態(tài)，必須轉(zhuǎn)化為正則化模態(tài)后才可用于特征載荷計算。轉(zhuǎn)化過程中質(zhì)量矩陣取為有限元模型的集中質(zhì)量陣，由于過程簡單，這里不再具體列出。

初始有限元模型中，取E1=E2=200 GPa，其余參數(shù)取為真實值，利用上述的擬靜力模型修正過程對該模型進行修正，修正時僅依據(jù)一階模態(tài)中測點處的信息進行修正，修正過程中彈性模量的變化情況如圖3中“E1”和“E2”所示，可以看出，相應(yīng)參數(shù)均收斂于真實值。

為了考察測量噪聲對修正結(jié)果的影響，把測量信息中加入強度為2%的白噪聲，相應(yīng)的彈性模量的修正結(jié)果如圖3所示，其中的“ER1”和“ER2”即為包含測量噪聲的彈性模量E1和E2修正值。由于測量噪聲的影響，彈性模量無法精確修正到真實值，但總體來說，修正結(jié)果仍然很好地逼近真實值。E1和E2修正后與真實值的相對誤差分別為0.37%和3.8%，其中E2的相對誤差相對高一些的原因是模型中僅有兩個單元的彈性模量為E2，結(jié)構(gòu)響應(yīng)對參數(shù)E2的變化不敏感，因此測量噪聲對E2修正結(jié)果的影響更大些。

修正過程中參數(shù)E1和E2的方差由式(26)計算可得，其變化情況如圖4所示，其中“cov1”和“cov2”分別代表彈性模量E1和E2的方差，可以看出它們的方差均很低，因此，可以認為修正結(jié)果的可靠度很高。

圖3 修正過程中參數(shù)變化情況Fig.3 Model parameters history in updating process

圖4 修正過程中參數(shù)方差Fig.4 Parameters variance in updating process

從上面的算例可以看出，提出的擬靜力模型修正方法能有效地實現(xiàn)基于動力測量數(shù)據(jù)的結(jié)構(gòu)模型參數(shù)修正。用該方法來修正潤揚長江大橋北汊斜拉橋的橋塔有限元模型。潤揚大橋斜拉橋橋塔凈高143 m［圖5(a)］，對于高聳的大型結(jié)構(gòu)無法用常規(guī)的人工激勵振動的方法進行測試，而采用環(huán)境激勵振動測試的方法。測試系統(tǒng)選用13只941B型伺服式超低頻加速度傳感器，其中一只傳感器作為參考點，其余12只傳感器按圖5b所示方式布置，測量方向為水平方向，測試過程按測量方向分為順橋向和橫橋向兩個部分進行。信號采集與分析采用AZ－Cras振動與動態(tài)信號采集系統(tǒng)，采用環(huán)境激勵下的模態(tài)識別方法進行參數(shù)識別。

圖5 潤揚大橋斜拉橋及橋塔傳感器位置Fig.5 Runyang cable-stayed bridge tower and its sensors location

在各測點信號與參考點信號之間互譜分析的基礎(chǔ)上識別橋塔結(jié)構(gòu)的固有頻率和振型，模態(tài)振型(前4階)如圖6所示。

圖6 識別的橋塔振動模態(tài)Fig.6 Identified modals of bridge tower

在建立潤揚斜拉橋橋塔有限元模型時，考慮到橋塔屬于高聳結(jié)構(gòu)，截面尺寸相對于自身縱、橫向尺寸較小，因此采用了空間二節(jié)點等參梁單元，共86個單元。

由于橋塔是鋼筋混凝土結(jié)構(gòu)，各部位的鋼筋配比率以及施工過程中的養(yǎng)護條件有所差異，因此各部位相應(yīng)的材料參數(shù)可能差別較大。在建立有限元模型時，對上、中、下塔柱以及三個橫梁采用不同的材料參數(shù)，因此待修正的參數(shù)共有6組。由于施工過程中鋼筋的配比量及混凝土的標(biāo)號已知，相應(yīng)區(qū)段的平均質(zhì)量密度可以較準(zhǔn)確地確定，這樣，在建模過程中各區(qū)段的材料密度按計算的平均密度給定，模型修正過程中不考慮密度參數(shù)的誤差，僅修正各區(qū)段的抗彎剛度，共12個參數(shù)。

在樣條插值擴展模態(tài)時，為了降低插值運算量，對塔柱的上、中、下三段分別進行插值，這樣，由于插值點與測量點在一條直線上，插值計算量得以降低。橫梁上的點插值時，由于各個橫梁的兩個端點均為測點，容易進行插值。

由于橋墩墩臺處的測量信號較弱，修正時沒有采用墩臺上2個測點處的信息，其余11個測點處的信息均作為響應(yīng)的測量值來修正模型。由于待修正的參數(shù)共12個，任一階模態(tài)的11個測點共有22個測量值，完全可以滿足參數(shù)的修正需要。為提高修正模型的精度，這里采用了前兩階的模態(tài)測量值進行修正，修正后各區(qū)段的參數(shù)如表1所示。表中的“EI11”和“EI22”分別為構(gòu)件橫截面繞形心主慣性軸的抗彎剛度。

表1 模型參數(shù)的修正結(jié)果Tab.1 Results of updated model parameters

為了檢驗修正結(jié)果的精度，根據(jù)修正得到的模型參數(shù)，利用ABAQUS軟件選取工程粱單元建立了橋塔結(jié)構(gòu)的有限元模型，網(wǎng)格劃分與修正模型的一致，計算的頻率值及模態(tài)如表2所示，表中“MAC”表示計算模態(tài)與測量模態(tài)的相關(guān)程度［4］，可以看出修正后模型的計算值與測試值吻合良好，其中第1、2階模態(tài)的計算值和測試值相關(guān)性很高，這主要是因為修正過程中直接采用了該兩階模態(tài)的測量值進行修正。雖然第3、4階模態(tài)的測量值在模型修正過程中沒有直接采用，但可以看出計算值與測量值仍然很接近，這從另一方面也證明了修正結(jié)果的正確性。

表2 修正模型的計算結(jié)果和測量值Tab.2 Computed results with updated model and the measured ones

對于修正結(jié)果的統(tǒng)計分析，由于測量過程中的誤差信息無法得到，這里假設(shè)測量過程中有方差強度為2%均值為0的白噪聲。利用式(26)計算參數(shù)修正過程中的方差情況，方差較大的為下塔柱的抗彎剛度EI11，疊代修正結(jié)束時其方差為5.7%，表明在修正后的下塔柱參數(shù)中尚有部分未確定因素沒有消除。下塔柱實際上是變截面構(gòu)件，而有限元模型中近似處理為等截面構(gòu)件，這一近似過程使得修正結(jié)果具有一定的未確定性。要進一步降低相關(guān)的未確定因素，可以把下塔柱細分為結(jié)構(gòu)參數(shù)不同的幾段，通過進一步地修正來降低其未確定性，這樣處理會增加待修正參數(shù)的數(shù)目，造成修正過程的計算量較大，在工程精度允許的前提下，可以適當(dāng)損失一定的精度而采用簡單模型。關(guān)于參數(shù)細分問題的研究這里不做進一步展開。

5 結(jié)論

本文提出了擬靜力模型修正方法，借助特征載荷概念，將基于動力特性測量信息的較為復(fù)雜的動力模型修正轉(zhuǎn)化為靜力模型修正，避免了模型修正過程中特征值和特征向量的計算問題，降低了修正過程的復(fù)雜程度。對修正過程中測量信息和計算信息間的匹配問題，通過樣條插值的辦法加以實現(xiàn)。采用復(fù)域靈敏度分析法，解決了修正過程中參數(shù)靈敏度高效計算問題。對參數(shù)修正過程中的誤差進行了統(tǒng)計分析，為評估修正結(jié)果的可靠性提供了基礎(chǔ)。數(shù)值算例表明，擬靜力模型修正方法進行參數(shù)修正時需要的測量信息有限，利于實際工程結(jié)構(gòu)模型的修正。對潤揚大橋北汊斜拉橋橋塔模型進行了修正，與測量值對比表明修正結(jié)果較為精確可靠。

［1］Mottershead J E，F(xiàn)riswell M I.Model updating in structural dynamics:a survey［J］.Journal of Sound and Vibration，1993，167(2):347－375.

［2］Guyan R J.Reduction of stiffness and mass matrices［J］.American Institute of Aeronautics and Astronautics Journal，1965，3(2):380.

［3］朱安文，曲廣吉，高耀南，等.結(jié)構(gòu)動力模型修正技術(shù)的發(fā)展［J］.力學(xué)進展，2002，32(3):337－348.

［4］郭 力，李兆霞，陳鴻天.基于子結(jié)構(gòu)分析的多重子步模型修正方法［J］.中國工程科學(xué)，2006，8(9):42－48.

［5］Gao X W，Liu D D，Chen P C.Internal stresses in inelastic BEM using complex-variable differentiation ［J］.Computational Mechanics，2002，28(1):40 －46.

［6］郭 力，高效偉.復(fù)變量求導(dǎo)法靈敏度分析及彈塑性參數(shù)反演［J］.東南大學(xué)學(xué)報，2008，38(1):141－145.

［7］Lyness J N，Moler C B.Numerical differentiation of analytic functions［J］.SIAM Journal of Numerical Analysis，1967，4(2):202－210.

［8］朱伯芳.有限單元法原理與應(yīng)用，第二版［M］.北京:中國水利水電出版社，1998.

［9］Li X Y，Law S S.Damage identification of structures including system uncertainties and measurement noise［J］.AIAA Journal，2008，46(1):263 －276.

［10］Law S S，Li X Y，Lu Z R.Structural damage detection from wavelet coefficient sensitivity with model errors［J］.Journal of Engineering Mechanics，2006，132(10):1077－1087.

［11］Papadopoulas L，Garcia E.Structural damage identification:a probabilistic approach［J］.AIAA Journal，1998，36(11):2137－2145.