楊秋偉, 周衛(wèi)東, 梁超鋒

(1.紹興文理學院 土木工程系,浙江 紹興 312000;2.華匯工程設(shè)計集團股份有限公司,浙江 紹興 312000)

重大工程結(jié)構(gòu)在服役期間，由于外力碰撞、環(huán)境腐蝕、材料老化等因素的影響，將不可避免的出現(xiàn)損傷。結(jié)構(gòu)的局部損傷將可能導致結(jié)構(gòu)整體的迅速破壞而釀成重大的工程事故。為了確保結(jié)構(gòu)的安全性能，必須對結(jié)構(gòu)損傷的發(fā)生、損傷的位置和程度作出及時的判斷。近年來，關(guān)于結(jié)構(gòu)損傷識別方面的研究已取得了長足的進展[1-4]。目前大部分損傷識別方法均是基于結(jié)構(gòu)有限元模型的模型修正方法，其原理是通過不斷修正結(jié)構(gòu)的有限元模型使其和測量所得的反應(yīng)數(shù)據(jù)相匹配，模型的修正量即為結(jié)構(gòu)的損傷量。模型修正方法通?？梢苑譃樗念悾壕仃噧?yōu)化方法[5-6]，靈敏度方法[7-9]，特征結(jié)構(gòu)分派方法[10-11]，最小秩擾動方法[12-16]。一般情況下，結(jié)構(gòu)損傷均發(fā)生于結(jié)構(gòu)中的少數(shù)區(qū)域，而最小秩方法中對更新矩陣的秩的限制最符合這種實際情況，因此，最小秩方法有著其獨特的優(yōu)越性[16]。

本文研究了結(jié)構(gòu)損傷與系統(tǒng)矩陣秩變化之間的對應(yīng)關(guān)系，并聯(lián)合利用最小秩和柔度擾動新方法提出了一種很精確的損傷識別方法。和已有的最小秩方法相比，本文方法的創(chuàng)新性在于：

(1) 已有的最小秩方法均屬于動力方法，都利用了結(jié)構(gòu)的振動響應(yīng)參數(shù)；而本文方法則屬于靜力方法，利用了結(jié)構(gòu)在靜力荷載作用下的反應(yīng)參數(shù)。動力方法和靜力方法各有優(yōu)缺點，Wang等[17]指出了動力方法中所存在的幾個必須解決的問題：①結(jié)構(gòu)的動力反應(yīng)數(shù)據(jù)除與結(jié)構(gòu)的剛度有關(guān)以外，還與結(jié)構(gòu)的質(zhì)量和阻尼有關(guān)。而結(jié)構(gòu)的損傷多數(shù)情況下只造成結(jié)構(gòu)剛度的損失，所以，許多動力方法都忽略了結(jié)構(gòu)阻尼和質(zhì)量變化對動力測試數(shù)據(jù)的影響，這顯然與實際工況有偏差；②對于巨大的土木工程結(jié)構(gòu)，難以測量出精確的振動數(shù)據(jù)；③高階的模態(tài)往往對損傷更為敏感，然而實踐中卻只能測量出低價模態(tài)數(shù)據(jù)。相比而言，結(jié)構(gòu)的靜力測試數(shù)據(jù)只和結(jié)構(gòu)的剛度有關(guān)，且一般可以很精確的測量得到，所以靜力方法在土木工程損傷識別領(lǐng)域有著可觀的應(yīng)用前景。當然，靜力方法也有其相應(yīng)的缺點：比如需要進行額外的靜力加載工作，某些情況下可能會對結(jié)構(gòu)造成二次損傷，不利于在線測量數(shù)據(jù)等等。因此，靜力方法和動力方法都各有優(yōu)缺點，應(yīng)根據(jù)具體的工程實際情況來選用。

(2) 已有的最小秩方法均是求解出剛度矩陣擾動的最小秩解，而本文方法則是直接求出柔度矩陣擾動的限定秩解，它是一種理論上的精確解。本文的研究還表明：損傷前后結(jié)構(gòu)的靜力位移差向量之間必然存在著某種線性相關(guān)性，所以所施加的靜力荷載組數(shù)并非越多越好，而是只要大于損傷前后柔度變化矩陣的秩即可。文中以一個桁架結(jié)構(gòu)為例對所提方法進行了驗證，結(jié)果表明了所提方法的可行性。

1 基于秩分析的損傷識別

1.1 柔度改變量的秩分析

本節(jié)首先分析由于結(jié)構(gòu)損傷所導致的柔度改變量的秩變化。結(jié)構(gòu)損傷前后，其剛度和柔度矩陣必然滿足以下關(guān)系式：

FK=FdKd=I

(1)

式中F和K是完好結(jié)構(gòu)的柔度矩陣和剛度矩陣(n×n維);Fd和Kd是結(jié)構(gòu)損傷后的柔度矩陣和剛度矩陣;I是n×n維單位矩陣。

一般而言，結(jié)構(gòu)損傷將導致結(jié)構(gòu)剛度減小而柔度增大。令ΔF和ΔK為結(jié)構(gòu)損傷前后柔度和剛度的變化量，即有：

Fd=F+ΔF

(2)

Kd=K-ΔK

(3)

將方程(2)和(3)代入(1)，整理可得

ΔFKd=FΔK

(4)

顯然，方程(4)中矩陣Kd和F均為對稱且滿秩的矩陣。因此，根據(jù)矩陣理論，由方程(4)必然有ΔF的秩和ΔK的秩相等，即

rank(ΔF)=rank(ΔK)

(5)

由有限元理論，結(jié)構(gòu)損傷前后剛度變化矩陣ΔK為

(0≤αi≤1)

(6)

式中;Ki是第i個單元剛度矩陣，αi為其相應(yīng)的損傷參數(shù)，N為單元總數(shù)。

一般情況下，結(jié)構(gòu)發(fā)生損傷僅限于局部少數(shù)單元，由方程(6)可知ΔK一般為虧損矩陣，且當發(fā)生損傷的單元數(shù)目改變時，ΔK的秩也將相應(yīng)改變(例如：對于桁架結(jié)構(gòu)(可參考算例)，若1個單元發(fā)生損傷，則對應(yīng)rank(ΔK)=1；若2個單元發(fā)生損傷，則對應(yīng)rank(ΔK)=2，其它具體的工程結(jié)構(gòu)亦可類似分析)。再根據(jù)方程(5)可知，ΔF一般也是虧損矩陣，且ΔF的秩也將隨著損傷單元的數(shù)目的變化而變化(仍以桁架結(jié)構(gòu)為例，若1個單元發(fā)生損傷，則對應(yīng)rank(ΔF)=1；若2個單元發(fā)生損傷，則對應(yīng)rank(ΔF)=2，其它具體結(jié)構(gòu)亦可類似分析)。

1.2 柔度擾動量的限定秩解

本節(jié)從結(jié)構(gòu)的靜力響應(yīng)方程出發(fā)，推導損傷前后結(jié)構(gòu)柔度擾動的限定秩解。對于未損傷結(jié)構(gòu)，在已知的外荷載li(n維列向量)的作用下，可以測量或者通過計算得到結(jié)構(gòu)相應(yīng)的靜力位移ui，即

Kui=li

(7)

方程(7)可以改寫為

ui=K-1li=Fli

(8)

(9)

方程(9)減(7)可得損傷前后靜力位移差Δui為

Δui=ΔFli

(10)

若在結(jié)構(gòu)上分別加載了n個線性無關(guān)的荷載向量li(i=1～n)，則可得n個方程，可以將其組合為如下的矩陣方程

ΔU=ΔFL

(11)

式中:矩陣ΔU=[Δu1,Δu2,…,Δun]，矩陣L=[l1,l2,…,ln]。顯然，由于所加的荷載向量是線性無關(guān)的，故矩陣L是滿秩矩陣，根據(jù)矩陣理論，有

rank(ΔU)=rank(ΔF)

(12)

方程(12)表明，即使測量了很多組靜力位移向量，所得到的靜力位移差向量Δui中只有少數(shù)幾個是獨立的(仍以桁架結(jié)構(gòu)為例，若僅一個單元發(fā)生損傷(rank(ΔF)=1)，則即使施加了n組荷載向量并測量得到n組位移數(shù)據(jù)，最終所得的位移差矩陣ΔU的秩仍然為1，即n個Δui中只有一個是獨立的，其它的均可用其線性表示。因此，應(yīng)用靜力方法進行結(jié)構(gòu)損傷識別時，施加靜力荷載的組數(shù)可以不用很多，只要大于損傷前后剛度或者柔度變化矩陣的秩即可。工程實踐中，由于結(jié)構(gòu)損傷前后的剛度或柔度改變量的秩事先是未知的，我們可以采取嘗試的辦法，首先施加一定組數(shù)的靜力荷載并計算相應(yīng)的ΔU的秩(采用計算ΔU的奇異值來確定其秩的方法)，如果計算結(jié)果中出現(xiàn)了近似為零的奇異值，則表明所施加的靜力荷載組數(shù)是足夠了的。反之則繼續(xù)增加靜力荷載組數(shù)，直到出現(xiàn)近似為零的奇異值為止即可。由于大部分情況下結(jié)構(gòu)損傷所引起的秩變化都是比較小的，因此通常所施加的靜力荷載組數(shù)都不會很多，通過上述嘗試過程很快就可以知道所施加的靜力荷載組數(shù)是否夠用。

計算柔度改變矩陣限定秩解的過程如下：假設(shè)ΔU的秩為r，則可以從ΔU中取出某個最大線性無關(guān)組ΔUr=[Δu1,Δu2,…,Δur](其中r=rank(ΔU))，該組對應(yīng)的荷載向量矩陣為Lr=[l1,l2,…,lr]，則方程(11)簡化為

ΔUr=ΔFLr

(13)

類似于最小秩理論[12-16]，由方程(13)可得ΔF的限定秩解為

(14)

1.3 損傷識別

當ΔF的限定秩解計算出來以后，接下來可用文獻[9]所提的柔度擾動新方法來計算各單元損傷參數(shù)αi(i=1～N)，據(jù)此便可判斷結(jié)構(gòu)中哪個單元損傷及其損傷程度。該方法的突出優(yōu)點在于計算過程不需要迭代運算或高階靈敏度分析，計算量小且精度高，且對于靜定結(jié)構(gòu)而言，該方法是一種精確方法。該方法的主要公式簡述如下(詳細推導過程見文獻[9])：

首先利用單元剛度矩陣的特征值分解與重新組合，可得結(jié)構(gòu)損傷前后的總剛度矩陣的分解形式：

K=CPCT

(15)

Kd=CPdCT

(16)

根據(jù)柔度矩陣和剛度矩陣互逆可得

F=(C+)TP-1C+

(17)

(18)

方程(18)減去方程(17)可得

ΔF=EΔBET

(19)

式中:矩陣E稱為柔度聯(lián)系矩陣，其計算公式為：

E=(C+)T=(CCT)-1C

(20)

對角矩陣ΔB為

ΔB=diag(β1,β2,…,βN)

(21)

式中:βi為第i個單元的柔度擾動參數(shù)，它和剛度擾動參數(shù)之間的關(guān)系為

(22)

綜上所述，計算損傷參數(shù)(即剛度擾動參數(shù))的主要步驟為：首先，根據(jù)方程(14)計算出ΔF的限定秩解，然后根據(jù)方程(19)計算出各柔度擾動參數(shù)ΔB，最后利用方程(22)計算出各剛度擾動參數(shù)αi(i=1～N)，根據(jù)所得結(jié)果便可對損傷情況作出判斷。

2 算例

以圖1所示桁架結(jié)構(gòu)為例，驗證本文所提的損傷識別方法。該結(jié)構(gòu)基本參數(shù)為：彈性模量E=200 GPa, 密度ρ=7.8×103kg/m3,單元長度L=1 m和橫截面面積A=7.85×10-5m2?？赡艿撵o力加載點也見圖1中。考慮3種加載方式，方式1：F1=10 kN，F(xiàn)2=0，F(xiàn)3=0；方式2：F1=10 kN，F(xiàn)2=10 kN，F(xiàn)3=0；方式3：F1=10 kN，F(xiàn)2=10 kN，F(xiàn)3=10 kN。顯然，這3種加載方式是線性無關(guān)的。假設(shè)3種損傷情況：① 單個損傷：單元17剛度損傷20%；② 多個小損傷：單元10和15剛度損傷10%和15%；③ 多個大損傷：單元10和15剛度損傷30%和40%。

對于第1種損傷情況，不考慮測量誤差時三種加載方式下所得的靜力位移差向量列于表1中(限于篇幅，表中只給出了節(jié)點7-12所對應(yīng)的數(shù)據(jù)，即對應(yīng)于自由度12-23處的數(shù)據(jù))。為說明所得三組位移差向量的線性相關(guān)性，各自由度所對應(yīng)的位移差比值也同時列于表1中。由表1可見，Δu1，Δu2和Δu3成比例，說明這3個位移差向量是線性相關(guān)的，只有一個是獨立的，由這3個向量所組成的位移差矩陣ΔU的秩為1。因此，可以任意取一個位移差向量Δui(i=1,2或3)，用方程(14)來計算ΔF的限定秩解，然后應(yīng)用方程(19)來計算各柔度擾動參數(shù)βi(i=1～N)，最后用方程(22)計算出各單元損傷參數(shù)αi(i=1～N)，結(jié)果見于圖2中。圖2中單元17的損傷程度計算值為α17=0.2，和假設(shè)值完全一致。這說明了在不考慮測量誤差的情況下，ΔF的限定秩解是一種理論上的精確解，而柔度擾動新方法對于靜定結(jié)構(gòu)而言，亦是一種精確方法。

圖1 平面桁架結(jié)構(gòu)及其靜力加載

表1 靜力位移差向量 (損傷情況1)

表2 靜力位移差向量 (損傷情況2)

圖2 不考慮測量誤差時的損傷識別結(jié)果(損傷情況1)

對于第2種損傷情況，不考慮測量誤差時三種加載方式下所得的靜力位移差向量列于表2中。由表2可見，Δu1，Δu2和Δu3并不完全成比例，為了更好的判斷這組位移差向量的線性相關(guān)性，我們通過計算這3個向量所組成的位移差矩陣ΔU的奇異值，根據(jù)所得非零奇異值的個數(shù)來判斷。表2中Δu1，Δu2和Δu3組成的矩陣ΔU的奇異值為：{1.080 4,0.027 3,0.000 0}。顯然，ΔU的秩為2，這說明Δu1，Δu2和Δu3三個向量中只有兩個是獨立的，這也說明僅兩個單元發(fā)生損傷。可從ΔU中取出兩個向量用方程(14)計算出ΔF的限定秩解，然后應(yīng)用柔度擾動新方法計算出各單元損傷參數(shù)，結(jié)果示于圖3中，其中單元10和15損傷程度計算值為α10=0.1，α15=0.15，和假設(shè)值完全一致。

對于第2種損傷情況，若考慮3%的測量誤差時(測量誤差的添加方法為：在每個精確的位移值的基礎(chǔ)上，加上誤差水平(如3%)乘以一個[-1,1]范圍內(nèi)的隨機數(shù))，我們可以求出矩陣ΔU的三個奇異值為{1.073 6,0.028 0,0.001 3}，其中第3個奇異值很小可以近似認為等于0，仍然可以認為ΔU的秩應(yīng)為2，因此取前兩個位移差向量來計算柔度擾動量的限定秩解，然后應(yīng)用柔度擾動新方法計算出各單元損傷參數(shù)見圖4中。另外，為了討論誤差大小對本文方法的影響，添加10%的誤差時的損傷識別結(jié)果也同時列于圖4中。由圖4可見，3%噪聲水平下本文方法識別結(jié)果和假設(shè)值很接近，10%噪聲水平下識別精度有所降低，但仍然可以清楚判斷出損傷單元為10和15。

對于第3種損傷情況，不考慮測量誤差時ΔU的奇異值為{4.158 9,0.103 4,0.000 0}；考慮3%的測量誤差時，ΔU的奇異值為[4.084 7,0.103 4,0.016 0}，考慮10%的測量誤差時，ΔU的奇異值為{3.655 0,0.154 1,0.066 2}。以上三種情況下，均可以認為ΔU的秩應(yīng)為2，因此采用兩個位移差向量來計算柔度擾動的限定秩解，最后的損傷識別計算結(jié)果都列于圖5中。由圖5可見，無誤差時損傷程度計算值為α10=0.3，α15=0.4，和假設(shè)值完全一致。3%和10%噪聲水平下也清楚的表明單元10和15發(fā)生損傷。因此，由于計算過程中對秩的限制，所提方法對測量噪聲有著較好的魯棒性。

圖5 損傷識別結(jié)果(損傷情況3)

3 結(jié) 論

本文研究了結(jié)構(gòu)損傷與系統(tǒng)矩陣秩變化之間的對應(yīng)關(guān)系，并聯(lián)合利用最小秩和柔度擾動新方法提出了一種很精確的損傷識別方法。所提方法的優(yōu)勢在于：由于方法中對系統(tǒng)矩陣秩的限制，因此計算結(jié)果對測量誤差具有很好的魯棒性。數(shù)值算例結(jié)果表明：所提方法是合理可行的。

[1]楊秋偉, 劉濟科. 工程結(jié)構(gòu)損傷識別的柔度方法研究進展[J]. 振動與沖擊, 2011, 30(12): 147-153.

YANG Qiu-wei, LIU Ji-ke. Structural damage identification by flexibility change: a review[J]. Journal of Vibration and Shock, 2011, 30(12): 147-153.

[2]程遠勝, 王真. 基于提高頻率靈敏度的結(jié)構(gòu)損傷統(tǒng)計識別方法[J]. 計算力學學報, 2008, 25(6): 844-849.

CHENG Yuan-sheng, WANG Zhen. Statistical damage detection based on frequencies of sensitivity-enhanced structures[J]. Chinese Journal of Computational Mechanics, 2008, 25(6): 844-849.

[3]安永輝, 歐進萍. 大跨度鋼桁橋模型的精細化損傷定位模擬和試驗研究[J]. 計算力學學報, 2011, 28(5): 730-736.

AN Yong-hui, OU Jin-ping. Experimental and numerical studies on precise damage localization of a long-span steel truss bridge model[J]. Chinese Journal of Computational Mechanics, 2011, 28(5): 730-736.

[4]Farrar C R, Lieven N A J. Damage prognosis: the future of structural health monitoring[J]. Philosophical Transactions of the Royal Society A, 2007, 365: 623-632.

[5]Kammer D C. Optimum approximation for residual stiffness in linear system identification. AIAA J, 1988, 26(1): 104-112.

[6]Smith S W, Beattie C A. Secant-method adjustment for structural models[J]. AIAA J, 1991, 29(1): 119-126.

[7]Yu L, Cheng L, Yam L H, et al. Application of eigenvalue perturbation theory for detecting small structural damage using dynamic responses[J]. Composite Structures, 2007, 78: 402-409.

[8]Shi Z Y, Law S S, Zhang L M. Damage localization by directly using incomplete mode shapes[J]. Journal of Engineering Mechanics, 2000, 126(6): 656-660.

[9]Yang Q W. A new damage identification method based on structural flexibility disassembly. Journal of Vibration and Control, 2011, 17(7): 1000-1008.

[10]Zimmerman D C, Kaouk M. Eigenstructure assignment approach for structural damage detection[J]. AIAA Journal, 1992, 30(7): 1848-1855.

[11]Kiddy J, Pines D. Eigenstructure assignment technique for damage detection in rotating structures[J]. AIAA Journal, 1998, 36(9): 1680-1685.

[12]Zimmerman D C, Kaouk M.Structural damage detection using a minimum rank update theory[J]. Journal of Vibration and Acoustics, 1994, 116(2): 222-231.

[13]Zimmerman D C, Smith S W, Kim H M, et al. An experimental study of structural damage detection using incomplete measurements[J]. Journal of Vibration and Acoustics, 1996, 118(4): 543-550.

[14]Doebling S W. Minimum-Rank optimal update of elemental stiffness parameters for structural damage identification[J]. AIAA Journal, 1996, 34(12): 2615-2621.

[15]Souza K D, Epureanu B I. Damage detection in nonlinear systems using system augmentation and generalized minimum rank perturbation theory[J]. Smart Materials and Structures, 2005, 14: 989-1000.

[16]Zimmerman D C. Statistical confidence using minimum rank perturbation theory[J]. Mechanical Systems and Signal Processing, 2006, 20: 1155-1172.

[17]Wang X, Hu N, Fukunaga H, et al. Structural damage identification using static test data and changes in frequencies[J]. Engineering Structures, 2001, 23: 610-621.