基于子結構的有限元模型修正方法

翁 順， 左 越， 朱宏平， 陳 波， 趙會賢， 田 煒， 顏永逸

(1.華中科技大學 土木工程與力學學院，控制結構湖北省重點實驗室，武漢 430074； 2.武漢理工大學 土木工程與建筑學院，武漢 430070)

基于子結構的有限元模型修正方法

翁 順1， 左 越1， 朱宏平1， 陳 波2， 趙會賢1， 田 煒1， 顏永逸1

(1.華中科技大學 土木工程與力學學院，控制結構湖北省重點實驗室，武漢 430074； 2.武漢理工大學 土木工程與建筑學院，武漢 430070)

提出了一種基于子結構的有限元模型修正方法。該方法將整體結構有限元模型劃分為多個子結構模型，求解獨立子結構的主模態特征解和特征靈敏度；通過位移協調條件和能量方程，約束相鄰獨立子結構，得到整體結構的特征解和特征靈敏度；并以整體結構模態和結構試驗模態的殘差為目標函數，通過調整子結構單元參數，完成有限元模型修正。當結構局部參數發生變化，通過分析某一個或幾個子結構即可求解整體結構特征解和特征靈敏度，而不需要分析其他未發生變化的子結構。由于子結構模型尺寸遠小于整體結構，該方法能夠極大地提高有限元模型修正方法的精度和效率。

有限元模型修正; 損傷識別; 子結構；特征解; 特征靈敏度

精確的有限元模型是結構健康評估、動態分析以及優化設計的基礎。建模過程中參數的不確定性、邊界條件的假定等，導致有限元模型的動態響應和實測的試驗數據存在誤差。設計規范規定, 有限元模型必須通過振動模態試驗或者地面共振試驗來檢驗。因此，近年來, 有限元模型修正技術得到了長足的發展[1-8]。根據修正對象的不同可將修正方法分為矩陣型方法和設計參數型方法。前者直接重建分析模型的剛度和質量矩陣。后者重復的修正有限元模型的物理參數來使模型的模態特性(頻率和振型)和現場試驗模態的殘差最小化。設計參數型方法物理意義明確，并且修正后的系統矩陣在對稱性、正定性和稀疏性等方面與修正前的矩陣保持一致。其中，基于靈敏度分析的設計參數型修正方法為模型修正過程提供有效地優化搜索方向，效率高，被廣泛的應用于實際工程結構中[9-13]。

大多數設計參數型修正方法采用優化技術，通過重復計算模型的特征解和特征靈敏度矩陣，尋找最接近實際工程結構的物理參數，實現模型修正過程。但是，土木工程結構體積龐大，為了有效地模擬實際工程結構，其有限元模型通常由大量的單元組成，并且包含許多需要被修正的結構設計參數，優化迭代過程通常需要花費很多時間，特別是從龐大的有限元模型的系統矩陣中重復計算特征解和特征靈敏度[14-16]。

子結構方法將龐大的整體結構分解為若干個較小的獨立子結構，通過位移協調條件和力的協調條件建立整體結構和子結構之間的關系。由于獨立子結構系統矩陣較小，未知參數少，在大型結構的模型修正及其相關應用中具有以下優勢[17]： ①子結構系統矩陣較小，分析精度和效率高;②各子結構相互獨立，可以對獨立子結構進行修正和重新分析，不需要分析整體結構;③針對當前的多核計算機硬件，進行并行計算，提高計算效率。

本文研究基于子結構的有限元模型修正方法。該方法通過組集子結構主模態的特征解和特征靈敏度，求解整體結構的特征解和特征靈敏度，用于有限元模型修正。在計算整體結構對結構單元參數的特征靈敏度時，只需要計算包含該單元參數的一個子結構的特征靈敏度，其他子結構的特征靈敏度為0。由于子結構的尺寸遠遠小于整體結構，該子結構方法極大的提高結構特征值和特征靈敏度的計算效率，從而有效地提高了模型修正過程的效率。基于正向子結構的有限元模型修正方法應用于一個實際工程橋梁的有限元模型修正中，驗證其精度和效率。

1 子結構方法求特征解

將N個自由度的整體結構分解為NS個子結構。以第j(j=1, 2, …,NS)個子結構為例，假如第j個子結構有n(j)個自由度，剛度矩陣為K(j)，質量矩陣為M(j)，其n(j)對特征值和特征向量為

[Φ(j)]TK(j)Φ(j)=Λ(j)

[Φ(j)]TM(j)Φ(j)=I(j)

(j=1, 2, …,NS)

(1)

KRON[18]的子結構方法基于虛功原理和幾何相容性原則，通過在子結構相鄰界面處施加約束，得到整體結構的特征方程為

(2)

其中，

Γ=[CΦp]T

Λp=Diag[Λ(1),Λ(2),…,Λ(NS)]

Φp=Diag[Φ(1),Φ(2),…,Φ(NS)]

(3)

從能量守恒的觀點出發，全部子結構的所有階模態都會對整體結構的特征模態有貢獻，也就是說，要組裝得到Λp和Φp需要計算全部子結構的所有階特征解。實際工程結構往往只計算部分低階模態，求解子結構的全部模態效率很低。為了克服這個困難，本文通過引入模態截斷方法來提高KRON子結構方法的效率。選取每個子結構的低階模態作為“主模態”，剩下的高階模態作為“從模態”。在組裝整體結構的特征方程時只需要計算主模態，從模態的能量貢獻用剩余柔度來補充。

(j=1,2,…,NS)

(4)

整體結構特征方程式(2)可以根據主模態和從模態重新分解為

(5)

根據式(5)的第二行，從模態的模態參與系數可以表示為

(6)

將式(6)代入式(5)中得到

(7)

(8)

(9)

式(9)的第二行將τ表示為zm，并且代入式(9)的第一行，可得到

(10)

(11)

(12)

簡化后的特征方程式(10)的尺寸等于NPm×NPm，遠小于原始特征方程式(2)的尺寸NP×NP，因此將極大地提高子結構方法求解特征解的效率。

2 子結構方法求解特征靈敏度

本節以第R個子結構的第r單元參數為例，推導整體結構第i階模態對參數r的特征靈敏度矩陣。式(10)用第i階模態表示為

(13)

(14)

(15)

其中,

(16)

(18)

將式(18)對參數r求偏導，可以得到第i階模態的特征向量靈敏度為

(19)

(20)

式中,ci為第i階模態{zi}的參與系數。

將式(20)代入式(14)得

(21)

求解式(21)，得到常向量{vi}。

求解式(13)的特征解，特征向量{zi}滿足正交條件

{zi}T{zi}=1

(22)

式(22)對r求偏導后得

(23)

將式(20)代入式(23)中，得到第i階模態的參與系數ci為

(24)

求解到常向量{vi}和第i階模態的參與系數ci后即可得到

(25)

3 基于靈敏度分析的有限元模型修正方法

在基于靈敏度分析的模型修正過程中，將有限元模型的計算模態同結構試驗模態的殘差作為目標函數

(26)

基于靈敏度分析的設計參數型修正方法，求解目標函數關于結構設計參數的靈敏度矩陣，提供最優的優化搜索方向。目標函數對結構設計參數的靈敏度矩陣為

(27)

其中特征值和特征向量對于參數r的靈敏度矩陣表示為

(28)

本研究采用子結構方法求解所有單元參數的一階偏導。在有限元模型修正過程中，基于子結構方法通過式(26)構建目標函數，通過式(29)建立靈敏度矩陣為優化算法提供搜索方向，采用常用的基于Trust-region的優化算法，完成有限元模型修正過程。

4 工程應用

為了驗證基于正向子結構的有限元模型修正方法在實際結構中的可行性和計算效率，將其應用于西澳大利亞的巴拉巴拉河大橋的有限元模型修正。根據設計圖紙可以建立如圖1所示有限元模型。該有限用模型的橋包含907個單元、947個節點。每個節點有6個自由度，共5 420個自由度。在現場的振動測試中，加速度計放置在7個縱梁對應的7列上，每列有19個測量點，共133個測量點。通過現場模態試驗提取前10階自振頻率和振型(見圖2)。

圖1 Balla Balla Bridge有限元模型和子結構劃分方法Fig.1 FE model of Balla Balla Bridge and the substructures

分別采用傳統的基于整體結構的有限元模型修正方法和本文所提出的基于子結構的有限元模型修正方法對有限元模型進行修正。一共選取了包括膈、梁、板的楊氏模量(E)，剪力連接件的拉壓剛度(EA)和抗彎剛度(EIxx,EIxy)在內的1 289個物理參數作為修正參數。目標函數根據式(26)定義為有限元模型的頻率振型與試驗模態的頻率振型之間的殘差。由于振型測量的誤差大于頻率，將振型的權重設定為0.1而頻率的權重設定為1.0。

首先應用傳統的整體結構有限元模型修正方法完成有限元模型修正。采用Lanczos方法從整體結構系統矩陣中計算結構特征解構建目標函數，用NELSON方法計算特征靈敏度用于提供優化搜索方向。由于試驗模態只測到部分模態，在每次優化迭代中，計算有限元模型前30階模態并與實測的前10階模態進行匹配。通過模態置信準則MAC(Modal Assurance Creterion)

(29)

從有限元模型計算的30階模態中找出與試驗模態匹配的10階模態。

模型修正過程迭代69次后收斂到預先設定的閾值，在普通個人電腦上每一步迭代大約花費1.26 h，模型修正過程共花費86.16 h，其收斂過程如圖3所示。

然后，應用所提出的基于子結構的方法有限元模型修正方法完成有限元模型修正，并采用和上述整體結構方法相同的優化算法、修正參數、收斂準則等前提條件，并且在相同的個人電腦上運行程序。用本文所提出的子結構方法計算整體結構的特征值和特征靈敏度。將整體結構沿縱向分成11個子結構，如表1所示。每一步優化迭代過程中，通過組集獨立子結構的特征值求解整體結構的特征值。計算整體結構對某一單元參數的特征靈敏度，只需要求解包含這個單元參數的一個子結構的主模態特征靈敏度，而其他的子結構的特征靈敏度設定為0。

圖2 Balla Balla Bridge的測量頻率和振型Fig.2 Measured frequencies and mode shapes

圖3 有限元模型修正過程Fig.3 The finite element model udating process

比較有限元模型修正前后結構的頻率和振型，如表2所示。該基于子結構的有限元模型修正方法可以得到與整體結構非常相近的結果。修正后模型頻率和實驗頻率之間的平均差異<1%。修正后振型的相關系數MAC從0.85提高到0.93。對于這個中型結構，在達到相同計算精度的條件下，基于子結構的有限元模型修正過程為傳統方法的大約一半，極大地提高了有限元模型修正的效率。子結構方法應用于大型結構對效率的提高將更明顯。

由于子結構中保留的主模態的數量會影響求解特征值和特征靈敏度的精度和效率，本文在模型修正過程中動態選取主模態。在模型修正初始階段，選取每個子結構的前40個模態作為主模態來計算整體結構的前30個特征值和特征靈敏度。然后，隨著參數逐漸接近最優解，逐步增加子結構的主模態的數量。直到在最后幾步中，在每個子結構中保留90個主模態來提高特征值和特征向量計算的精度。該基于子結構的有限元模型修正過程迭代76次收斂到相同的閾值，如圖3所示?；谧咏Y構的有限元模型修正方法所花費的計算時間如表3所示，共花費48.07 h，所需時間大約是傳統的整體結構模型修正方法的56%。

表1 整體結構劃分為11個子結構后各子結構的信息

表2 修正前后橋的頻率和振型

表3 基于整體結構有限元模型修正方法和基于子結構的有限元模型修正方法計算時間對比

5 結 論

本文提出了一種基于子結構的有限元模型修正方法。將整體結構有限元模型劃分為多個獨立子結構有限元模型，通過求解一個或幾個發生變化的獨立子結構特征解，即可求解整體結構特征解，并用于有限元模型修正和損傷識別。并且，采用子結構方法只需要計算某一個子結構特征解靈敏度矩陣即可完成對整體結構特征靈敏度的求解。當結構局部發生損傷，只需要重復某一個或幾個子結構模型，避免對整體結構模型重復分析，從而有效地提高大型結構有限元模型修正的精度和效率?；谧咏Y構的有限元模型修正方法可以有效地用于對大型結構的動力分析和健康監測中。

[ 1 ] MOTTERSHEAD J E, FRISWELL M I. Model updating in structural dynamics: a survey[J]. Journal of Sound and Vibration,1993,167(2): 347-375.

[ 2 ] 朱宏平, 黃民水. 基于環境激勵的橋梁結構動力有限元模型修正研究[J]. 華中科技大學學報, 2009,26(1):1-11. ZHU Hongping, HUANG Minshui. Study on dynamic finite element model updating of bridge structures based on ambient excitation[J]. Journal of Huazhong University of Science and Technology, 2009, 26(1):1-11.

[ 3 ] 方圣恩, 基于有限元模型修正的結構損傷識別方法研究[D]. 長沙：中南大學,2010.

[ 4 ] 姜東, 費慶國, 吳邵慶. 基于攝動法的不確定性有限元模型修正方法研究[J].計算力學學報, 2014, 31(4): 431-437. JIANG Dong,FEI Qingguo, WU Shaoqing. A study on stochastic finite element model updating based on perturbation approach[J]. Chinese Journal of Computational Mechanics, 2014,31(4): 431-437.

[ 5 ] 任偉新, 陳華斌. 基于響應面的橋梁有限元模型修正[J]. 土木工程學報, 2008,41(12): 73-78. REN Weixin, CHEN Huabin. Response-surface based on finite element model updating of bridge structures[J].China Civil Engineering Journal, 2008, 41(12): 73-78.

[ 6 ] 郭勤濤, 張令彌, 費慶國. 結構動力學有限元模型修正的發展-模型確認[J].力學進展, 2006, 36(1): 36-42. GUO Qintao, ZHANG Lingmi, FEI Qingguo. From FE model updating to model validation: advances in modeling of dynamic structures[J]. Advances in Mechanics, 2006, 36(1): 36-42.

[ 7 ] BROWNJOHN J M W, MOYO P, OMENZETTER P, et al. Assessment of highway bridge upgrading by dynamic testing and finite-element model updating[J]. Journal of Bridge Engineering,2003,8(3): 162-172.

[ 8 ] JAISHI B, REN W X. Structural finite element model updating using ambient vibration test results[J]. Journal of Structural Engineering,2005, 131(4): 617-628.

[ 9 ] ZHU H P, MAO L, WENG S. A sensitivity-based structural damage identification method with unknown input excitation using transmissibility concept[J]. Journal of Sound and Vibration, 2014, 333(26): 7135-7150.

[10] WENG S, ZHU A Z, ZHU H P, et al. Dynamic condensation approach to the calculation of eigensensitivity[J]. Computers and Structures, 2014, 132(1): 55-64.

[11] WENG S, ZHU H P, XIA Y, et al. Substructuring approach to the calculation of higher-order eigensensitivity[J]. Computers and Structures, 2013, 117(2): 23-33.

[12] WENG S, ZHU H P, XIA Y, et al. Damage detection using the eigenparameter decomposition of substructural flexibility matrix[J]. Mechanical Systems and Signal Processing, 2013, 34(1/2): 19-38.

[13] WENG S, XIA Y, XU Y L, et al. Improved substructuring method for eigensolutions of large-scale structures[J]. Journal of Sound and Vibration, 2009, 323(3/4/5): 718-736.

[14] CHOI K K, KIM N H. Structural sensitivity analysis and optimization 1: linear systems[M].New York: Springer Science and Business Media, 2005.

[15] FOX R L, KAPOOR M P. Rate of change of eigenvalues and eigenvectors[J]. AIAA Journal,1968, 6(12): 2426-2429.

[16] NELSON R B. Simplified calculation of eigenvector derivatives[J].AIAA Journal,1976,24(9):823-832.

[17] CRAIG R R. Coupling of substructures for dynamic analysis:an overview[C]//Atlanta: Proceedings of 41st AIAA/ASME/ASCE/ASC Structures, Structural Dynamics, and Materials Conference, 2000:171-179.

[18] KRON G. Diakoptics[M].London:Macdonald and Co., 1963: 46-87.

Model updating based on a substructuring method

WENGShun1，ZUOYue1，ZHUHongping1,CHENBo2，ZHAOHuixian1，TIANWei1，YANYongyi1

(1. School of Civil Engineering and Mechanics,Huazhong University of Science and Technology，Wuhan 430074, China； 2. School of Civil Engineering and Architecture,Wuhan University of Technology，Wuhan 430070, China)

This paper proposes a substructure-based model updating method. The global structure was divided into independent manageable substructures, which were analyzed to calculate the substructural eigensolutions and eigensensitivity. Afterwards, the independent substructures were constrained to calculate the eigensolutions and eigensensitivity of the global structure. Finally, the eigensolutions were used for the objective function and the eigensensitivity was used to indicate the searching direction to achieve the model updating. When the local area of a structure was changed, only one or more substructures were analyzed whereas the other substructures were untouched. Since the substructures were much smaller than the global structure, the proposed method significantly improved the accuracy and efficiency of the model updating process.

model updating; substructure; damage identification; substructure; eigensolution; eigensensitivity

國家自然科學基金(51328802;51108205);中央高校基本科研業務費專項資金(2014TS130;2015MS064);武漢城建委科研項目(201511)

2015-11-23 修改稿收到日期：2016-01-31

翁順 女，博士，副教授，1982年生

TH212；TH213.3

A

10.13465/j.cnki.jvs.2017.04.016