基于靜力殘余力向量的結構損傷評估方法

楊秋偉， 周 聰， 李翠紅， 羅 帥

(1.紹興文理學院 土木工程系，紹興 312000； 2.寧波工程學院 建筑與交通工程學院，寧波 315211；3.寧波工程學院 浙江省土木工程工業化建造工程技術研究中心,寧波 315211)

1 引 言

有限單元法已廣泛用于土木、機械、航空航天、汽車和船舶等眾多工程設計領域中，成為解決復雜工程分析計算問題的有效途徑。然而，由于大型工程結構的復雜性，用有限元軟件建立的結構模型與真實結構總是存在一定的差異，具體體現在，由有限元模型計算所得的結構響應參數(如位移、頻率和振型等)總是與儀器測試所得的結構響應參數存在差異。當這種差異較大時，則必須對結構的有限元模型進行修正，確保計算所得與測試所得的響應參數盡可能接近。只有修正后的有限元模型方可用于力學分析與計算。進一步，如果結構中存在損傷，通過修正完好結構的有限元模型來與測試所得的損傷結構的響應參數相匹配，根據模型的修正量就可以判斷出結構的損傷單元及損傷程度，這就是基于模型修正的結構損傷識別方法基本原理。近年來，模型修正法和以它為基礎的損傷識別方法已成為眾多工程技術領域的一項重要課題[1-4]。

殘余力向量法[5-16]是一類常用的損傷識別算法，由殘余力向量可以判定損傷位置，求解殘余力方程可以得到損傷程度。求解方法有最小范數法[5,6]、最小秩法[7,8]或者其他的迭代優化算法[14-18]。目前，已有的殘余力向量法都是基于動力測試的模態參數，因此可以統稱為動力殘余力向量法。和動力測試數據相比，靜力測試數據往往精度更高，且無需模態分析等復雜操作，易于實現。因此，本文基于靜力測試數據(位移參數)，結合有限元剛度矩陣，定義了靜力殘余力向量，根據其中不為零的元素來判斷損傷自由度，再根據自由度和單元之間的對應關系來確定損傷位置，在此基礎上，提出一種求解損傷參數的代數解法。另外，鑒于梁、剛架或實體結構中節點角位移難以測量的情況，進一步提出了一種靜力縮聚殘余力向量法，只需要應用節點線位移數據，從而有效拓寬了所提方法的應用范圍。數值算例結果說明所提方法合理有效。

2 靜力殘余力向量法

一般而言，結構的損傷主要導致結構剛度的變化，質量的變化基本可以忽略不計。不失一般性，假設K0和K分別為結構損傷前后的剛度矩陣(n×n維)，由于損傷所導致的剛度損失量為ΔK，即

(1)

式中Ki為第i個單元剛度矩陣，αi為第i個單元的損傷參數(αi∈[0,1])，αi=0表示第i個單元未損傷，αi=1表示第i個單元完全損傷，N為單元體總數目。結構損傷前后，在荷載向量l作用下的位移向量分別為

K0d0=l, (K0-ΔK)d=l

(2,3)

式中d0和d分別為完好結構和損傷結構的位移向量，其可以通過靜力測試得到。由方程(3)變形可得

ΔK·d={δ}， {δ}=K0d-l

(4，5)

方程(4)即為靜力殘余力方程，其中向量{δ}稱為靜力殘余力向量，由方程(5)可知，只要測試獲得損傷結構在已知荷載向量l作用下的位移d，即可由方程(5)計算得出其靜力殘余力向量。根據方程(4)，由靜力殘余力向量{δ}來進行損傷定位的原理在于，結構損傷一般只導致ΔK的少數行向量不為零，即和損傷對應的自由度所在的行向量不為零，因此，靜力殘余力向量{δ}中不為零的元素即對應著損傷自由度，由損傷自由度并根據幾何關系即可確定發生損傷的具體單元。具體而言，可將方程(4)改寫為以下方程組，

(6)

當損傷位置由上述方法判定出來以后，結合方程(1，6)，可以進一步求解損傷程度，即求解損傷參數αi。分兩種情況討論。

首先討論單個損傷情況，如果由靜力殘余力向量{δ}判定出來只有一個單元體發生損傷，則只有一個未知的損傷參數αi需要求解，此時方程(1)簡化為

ΔK=αiKi

(7)

相應地，方程(6)簡化為

(8)

(9)

方程(9)中δj取{δ}的任意一個非零元素即可。

(10)

方程(10)代入式(6)整理可得

Π{αi}′={δ}′

(11)

式中Π為整理所得的系數矩陣，{αi}′為q個未知損傷參數組成的列向量，{δ}′為q個非零的δj組成的列向量。由方程(11)即可求得損傷程度為

{αi}′=Π-1{δ}′

(12)

對于梁、剛架和實體結構，其節點位移既有線位移也有角位移，實踐中由于角位移難以測量，上述的靜力殘余力向量法需要進一步改良方可使用。為此，本文采用模型縮聚來對這些角位移進行縮減處理，從而使得縮聚后的相關公式都只包含各節點線位移，拓寬了靜力殘余力向量法的應用范圍。縮聚方法具體過程如下。

首先,將結構的全部自由度分為平動自由度(保留自由度)和轉動自由度(縮減自由度)兩類，相應地將方程(2)改寫為

(13)

式中上標m和s分別表示平動自由度和轉動自由度，m和s同時也表示平動自由度和轉動自由度的數目。若所有的轉動自由度上均未施加外荷載，即ls=0，則由方程(13)的第二行可得

(14)

由方程(14)可得

(15)

(16，17)

式中矩陣Θ為自由度轉換矩陣。將方程(16)代入式(13)，并兩邊左乘ΘT可得

Krdm=lr,Kr=ΘTK0Θ,lr=ΘTl

(18～20)

利用方程(18)縮聚后的剛度矩陣和荷載向量，按照前述靜力殘余力向量同樣的推導過程，可推導得到縮聚后的殘余力向量計算公式如下，

{δ}r=Krdm-lr

(21)

式中Kr為完好結構縮聚后的剛度矩陣，lr為縮聚后的荷載向量，dm為損傷后結構各平動自由度處的線位移向量。相應地，也可以得到模型縮聚后類似于方程(9，12)的計算損傷程度公式為

(22)

(23)

3 數值算例

3.1 桁架結構

圖1 桁架結構及靜力加載

圖2 單元10損傷時的靜力殘余力向量(單位：kN)

表1 損傷結構位移和殘余力向量計算結果

接下來討論多個損傷情況，不失一般性，假設單元16和31剛度分別損傷15%和30%，這種情況下損傷結構的位移數據和殘余力向量分別列于表1的第4列和第5列，為了更直觀，將表1的第5列數據繪入圖3。可以看出，殘余力不為零的自由度編號為9/13，25/26/31/32，根據圖1，自由度9和13對應的節點為第5和第7號節點，顯然，這兩個節點之間的桿件剛好為單元16，而自由度25/26/31/32分別對應著13和16號節點，這兩個節點之間為單元31。因此，根據殘余力向量中不為零的元素最終可以判定單元16和31為發生損傷的單元。接下來計算損傷程度，由于欲求的損傷參數個數為2，需要列兩個方程(方程(11))，這里選擇第9個自由度(對應損傷單元16)和第32個自由度(對應損傷單元31)來列方程，方程(11)的具體計算結果為

圖3 單元16和31損傷時的靜力殘余力向量(單位：kN)

(24)

進一步由方程(24)計算可得損傷參數為α16=0.15和α31=0.3，剛好與假設值15%和30%相等。

需要說明的是，在上述模擬過程中，并未考慮有限元模型誤差與測試數據誤差，因此計算所得的損傷程度值和假設值完全相等，這也直觀地說明了本文所提的損傷程度代數解是一種理論上的精確解。如果考慮數據誤差，如對于單元10損傷情況，在損傷結構各個位移數據中添加誤差水平為3%的隨機數，則可得靜力殘余力向量如圖4所示，可以看出第7和第11自由度(對應單元10)處殘余力明顯偏大，故仍然可以判斷單元10為損傷單元，其損傷程度計算值為0.1989，與假設值0.2很接近。如果有些情況難以根據殘余力向量圖形做出準確的損傷定位，可以結合殘余力向量圖形和損傷程度計算結果來綜合進行判定，具體見接下來的梁結構算例。

圖4 利用含噪聲數據計算所得靜力殘余力向量(單元10損傷，單位：kN)

3.2 梁結構

圖5 梁結構

不失一般性，假設單元15剛度損傷20%，圖6給出了無噪聲時靜力縮聚殘余力向量的計算結果。可以看出，絕對值偏大的殘余力集中于節點13，14，15和16處。根據單元體和平動自由度編號之間的對應關系，位于這4個節點之間的是第14，15和16號單元。由此可見，采用靜力縮聚殘余力向量，由于模型縮聚所帶來的誤差，不能精確定位于真正發生損傷的第15號單元處，而是定位于以損傷單元為中心的局部小區域。為了進一步準確判斷，分別計算第14，15和16號單元的損傷程度值。首先，計算第14號單元，該單元對應的節點編號為13和14，這里采用節點13對應的殘余力來進行計算(由于節點14的殘余力可能由于第14號單元或第15號單元引起，故此處不宜采用節點14)，所得結果為α14=-0.3581，據此判斷第14號單元并非真正損傷單元。接下來，計算第16號單元，該單元對應的節點編號為15和16，這里采用節點16對應的殘余力進行計算，結果為α16=-0.3402，據此判斷第16號單元也并非真正損傷單元。最后，計算第15號單元，該單元對應的節點編號為14和15，為了可靠起見，分別采用這兩個節點處殘余力對應的方程來進行計算，所得結果分別為α15=0.2156(節點14)和α15=0.2159(節點15)，這兩個結果均說明第15號單元為真正損傷單元，損傷程度計算值可取為二者的平均值，顯然與假設值0.2很接近。考慮在位移數據中添加3%噪聲水平的隨機數，所得靜力縮聚殘余力向量如圖7所示，可見絕對值偏大的殘余力仍然集中于節點13，14，15和16處，利用和上述相同的損傷程度計算過程，所得結果分別為α14=-0.3005，α16=-0.3714，α15=0.2012(節點14)和α15=0.2235(節點15)，故仍然可以判定出第15號單元為損傷單元，說明所提的靜力縮聚殘余力向量法是合理可行的。

圖6 無噪聲時靜力縮聚殘余力向量(單元15損傷，單位：kN)

圖7 有噪聲時靜力縮聚殘余力向量(單元15損傷，單位：kN)

4 結 論

提出一種靜力殘余力向量法用于結構損傷評估，利用結構靜力測試數據和有限元模型剛度矩陣，根據靜力殘余力向量中不為零的元素來判斷損傷自由度；再根據自由度和單元之間的對應關系來確定發生損傷的位置，并提出一種求解損傷參數的代數解法；另外，針對實踐中角位移難以測量的情況，進一步提出了一種靜力縮聚殘余力向量法，有效拓寬了所提方法的應用范圍。數值算例結果表明，所提方法合理可行，對單個損傷和多個損傷情況均可適用，可為解決結構損傷識別問題提供參考。