基于OMA試驗模態參數的砌體結構有限元建模及修正

管曉明，傅洪賢，王夢恕，崔堃鵬，林凡濤

(北京交通大學 土木建筑工程學院，北京 100044)

在城市地鐵建設中較多為下穿隧道工程。用鉆爆法開挖隧道時爆破產生的地面振動極易導致臨近建筑物損傷甚至破壞。砌體結構房屋建筑因整體性差、應力分布不均勻，更易在隧道施工爆破振動作用下損傷致結構開裂，繼而威脅居民的正常工作及生活。因此，該種結構的損傷評價及振動控制技術研究頗受關注，而建立能準確、全面反映砌體真實動力學行為的結構模型為評價結構動力損傷的重要基礎。建筑結構模型大多采用有限元模型，但因建模時各種簡化及幾何、材料、邊界等參數取值不確定性，使有限元模型與實際結構有一定差距。目前大多采用模型修正技術建立準確的結構有限元模型。通過對建筑結構進行模態試驗獲取模態參數(固有頻率、振型和阻尼)，據結構固有頻率參數靈敏度分析選擇結構模型中靈敏度較大材料參數修正結構模型，從而使有限元模型分析所得模態參數與試驗模態參數趨于一致。該修正方法參數意義明確，結果準確可靠，工程實用性好，已廣泛用于機械、橋梁、航天等領域[1-3]；但對建筑結構尤其砌體結構模型修正研究較少。

本文以成渝客運專線新紅巖隧道工程為背景進行相關研究。該隧道出入口淺埋段下穿重慶市沙坪壩區建設新村與新民坡村，埋深僅15～30 m，隧道周圍2、3層砌體房屋密集，且房屋老舊，隧道施工爆破時其安全會直接受到威脅。在隧道上方選一典型2層磚房，通過對磚房進行隧道爆破振動激勵的OMA(Operational Modal Analysis)試驗獲取結構模態參數，采用等效體積單元法建立砌體結構有限元模型，以試驗所得模態參數為基準，采用基于結構固有頻率的參數靈敏度分析修正方法借助ANSYS軟件修正砌體結構有限元模型，為評價砌體結構在隧道施工爆破作用下的動力損傷提供準確、可靠的結構模型。

1 砌體結構有限元模型建立

1.1 砌體結構等效體積單元方法

砌體由磚、砂漿規則砌筑而成，需將砌體均質化為連續性介質，實現采用一種單元建立砌體結構有限元模型。目前大多采用等效體積單元(Representative Volume Element，RVE)法，其包含砌體所有集合與組成信息結構，適用于大規模砌體結構力學行為機理研究。砌體結構等效體積單元建模滿足條件[4]為：①包括組成砌體的所有材料，如磚塊、砂漿；②能按周期性、連續分布規律組成完整結構；③滿足上兩條件的最小單元。據此，砌體采用等效體積單元的均質化建模過程見圖1。等效體積單元應力應變值用單元中各組成部分的應力應變平均值。等效體積單元用正交各向異性材料模型時的材料參數[4]，即將等效體積單元中磚塊、砂漿在彈性階段用各向同性模型，在塑性階段用Drucker-Prager模型，分別對等效體積單元進行三向(X、Y、Z)單軸抗壓及純剪切試驗，獲得相應方向彈性模量、泊松比及剪切模量。采用有限元數值分析或室內試驗。據試驗結果，砌體結構等效體積單元正交各向異性模型參數計算公式[4]為

Ex=σx/εx，μxy=εy/εx，μxz=εz/εx

(1)

Ey=σy/εy，μyx=εx/εy，μyz=εz/εy

(2)

Ez=σz/εz，μzx=εx/εz，μzy=εy/εz

(3)

Gxy=τxy/γxy，Gyz=τyz/γyz，Gzx=τzx/γzx

(4)

圖1 砌體等效體積單元均質化過程[5]

對復雜條件下應力狀態，砌體采用正交各向異性模型時本構模型為

(5)

式中：εx，εy，εz為正應變；γxy，γyz，γzx為剪應變；σx，σy，σz為正應力；τxy，τyz，τzx為剪應力；Ex，Ey，Ez為彈性模量；μxy，μxz，μyx，μyz，μzx，μzy為泊松比；Gxy，Gyz，Gzx為剪切模量；x，y為水平坐標；z為垂直坐標。

1.2 砌體與混凝土物理力學參數取值

砌體材料參數[6-9]為，密度1 600～2 000 kg/m3，彈性模量1.062～3.646 GPa，泊松比0.15～0.16。據文獻[10]燒結普通磚選MU20～MU10，砂漿選M2.5～M5，計算得砌體彈性模量為1.807～3.392 GPa，確定彈性模量為1.800～3.600 GPa，砌體剪切模量按砌體彈性模量的0.4倍選用。混凝土材料參數為，強度等級選C15～C20，彈性模量取22.0～25.5 GPa[11]，剪切變形模量可按彈性模量的0.40倍選用，泊松比 0.18～0.22；密度取 2 200～2 400 kg/m3。砌體、混凝土材料參數取值范圍及平均值見表1。

表1 砌體、混凝土材料參數取值表

1.3 砌體結構有限元模型建立

經現場調研、實測知，砌體結構有限元模型建模幾何參數為，樓層總高7 m，層高3 m，女兒墻高0.7 m，樓板厚0.15 m，外縱墻及承重橫墻厚0.24 m，部分隔墻厚0.12 m。房屋一層10.14×7.19 m，二層10.14×8.31 m，挑梁總長3.22 m，二層挑出長度1.12 m。砌體為橫向承重結構，二層地面與頂層為混凝土預制板，未設置圈梁及構造柱。二層磚房正、側面見圖2。

圖2 二層磚房正面、側面圖

圖3 砌體結構初始有限元模型

用ANSYS14.0建立砌體結構有限元模型，建模過程中需適當簡化，不考慮樓梯、窗戶及門。砌體與混凝土構件均采用solid185六面體單元，砌體結構有限元模型節點66887個，單元40940個，見圖3。

2 砌體結構OMA試驗與有限元模態分析

2.1 隧道爆破振動激勵砌體結構OMA試驗

結構模態試驗通用方法為試驗模態(EMA)法及運行模態(OMA)法。OMA法無需測量激勵力，只需測量結構在外界振動激勵的響應數據。本文用OMA法進行砌體隧道施工爆破振動激勵的模態試驗。

OMA模態試驗采用同步測試法，布置測點12個，其中一、二樓傳感器布置方式一致，縱向測點均4個(Z1，Z2，Z3，Z4)，橫向測點均2個(H1，H2)，見圖4。傳感器安裝時需保持與測量方向一致，底部用石膏固結，需防止振動或人為干擾。試驗用TST126型低頻水平速度傳感器，通過伺服放大器將采集信號傳送至數據處理系統進行轉換、存儲。模態試驗選INV3060S型智能數據采集處理分析儀與DASP軟件及模態分析軟件。試驗中測點典型振動波形見圖5，測點頻譜分析見圖6。

隧道爆破振動激勵較自然環境振動激勵雖持時較短(約1 s，圖5)，但能量大、頻帶寬，高頻成分豐富(圖6)，在測點布置足夠多條件下，除能激勵出結構低階振型外，亦能激勵出結構較高階振型，此對研究結構在隧道施工爆破作用的動力損傷極為重要。由于爆破地震波主頻較高，會致結構發生高階共振效應，造成局部應力過大，導致結構開裂。故隧道爆破振動激勵效果好于自然環境脈動激勵；但爆破激勵時由于振動能量豐富，易摻雜巖石振動能量或周圍房屋振動能量，結構模態參數識別時需謹慎，以確保獲取結構固有頻率而非周圍環境干擾振動頻率。

圖4 OMA測試傳感器一層布置圖

2.2 隨機子空間法模態參數識別

砌體結構OMA試驗后采用隨機子空間法(SSI)進行砌體結構模態參數識別。SSI法尤其適用于OMA試驗的結構固有頻率、阻尼及振型模態參數識別，識別精確度高。SSI法基本模型-離散時間隨機狀態空間模型[12]為

xk+1=Axk+wk,yk=Cxk+vk

(6)

式中：A為系統離散狀態矩陣；C為輸出矩陣；xk為離散時間狀態向量，yk為輸出狀態向量；wk，vk為環境激勵、測試過程誤差，通常設均值為零、互不相關的白噪聲協方差矩陣為

(7)

式中：E為數學期望；δpq為Kronrcker delta(p=q時δpq=1，p≠q時δpq=0)；p，q，k為離散時間點；Q，S，R為wk及vk協方差矩陣分塊矩陣。

隨機子空間方法有協方差驅動隨機子空間(Cov-SSI)及數據驅動隨機子空間(Data-SSI)兩種算法。應用較多的數據驅動隨機子空間算法識別步驟[13]為

(1) 利用系統輸出響應構建Hankel矩陣：

(8)

式中：Y0|2i-1為Hankel矩陣，Y0|2i-1∈R2ij，每行代表塊行，即由l個輸出響應組成；下標p，f表示“過去”與“將來”，為Hankel矩陣劃分塊行方式。

(2) 計算“將來”輸入行空間在“過去”輸入行空間投影，并經QR分解在保持系統原信息情況下縮減數據：

(9)

式中：(·)+為·的廣義逆。

通過對Y0|2i-1的Hankel矩陣進行QR分解，Pi可表示為(n為Pi的秩)

Pi=RQT∈Rnij

(10)

投影計算為隨機子空間算法核心，可利用“過去”行空間信息預測“將來”；

(3) 對投影進行奇異值分解，并結合卡爾曼濾波理論計算系統狀態矩陣A及輸出矩陣C；

(4) 對A進行特征值分解，獲得特征值、特征向量求解系統模態參數。

采用DASP模態分析軟件中SSI方法計算獲得穩定圖見圖7。穩定圖含義[14]為，將不同階數模型模態參數繪于同一幅圖中，在相應于某階模態軸上將高一階模型識別模態參數與低一階參數比較，若特征頻率、阻尼比及模態振型差異均小于預設的限定值，則該點稱穩定點，組成的軸稱穩定軸，相應模態即為系統模態。限定值可據實際工程及經驗確定，一般設特征頻率為1%，阻尼比 5%，模態振型 2%。由圖7看出，穩定圖由兩部分組成，第一部分互譜，在譜峰對應的豎向位置會出現一排特征頻率。第二部分用字母“s”“d”“v”“f”“o”分別對應不同階數計算模型所得特征頻率特性。“s”表示頻率、阻尼及振型均穩定，“d”表示頻率、阻尼比穩定，“v”表示頻率、振型穩定，“f”表示頻率穩定，“o”表示新頻率。若穩定圖中譜峰對應的豎向位置出現由低到高的“s”，即對應結構一階模態。模態參數識別過程中在不明的顯譜峰處出現穩定圖可能由于計算誤差及干擾信號產生的虛假模態引起，需剔除該模態。由圖7識別的結構真實模態參數見表2。

圖7 SSI法計算的穩定圖

表2 砌體結構OMA模態試驗分析結果

2.3 有限元模態分析及實測結果對比

為分析砌體采用不同材料模型計算砌體結構模態參數的差別，分別采用各向同性、橫觀各向同性、正交各向異性材料模型進行模態分析。混凝土用各向同性模型，兩種材料物理力學參數取平均值(表1)。砌體材料參數取值為砌體用各向同性模型，選砌體材料參數均值；砌體材料用正交各向異性模型，包含Ex，Ey，Ez，μxy，μyz，μzx，Gxy，Gyz，Gzx9個獨立參數。令Ex=Ey=Ez=2.7 GPa，μxy=μyz=μzx=0.15，Gxy=Gyz=Gzx=1.08 GPa。砌體材料采用橫觀各向同性模型，考慮x，y向同性，包含Ex(或Ey)，Ez，μxy，μyz(或μzx)，Gyz(或Gzx)5個獨立參數。令Ex(Ey)=2.7 GPa，Ez=2.7 GPa，μxy=0.15，μyz(μzx)=0.15，Gyz(Gzx)=1.08 GPa。

模態分析時固定一樓底部節點6個自由度，用Block Lanczos法計算砌體結構1～4階自振頻率及振型，并對比頻率計算值與實測值相對誤差及總相對誤差，總相對誤差為前4階相對誤差平方和。初始有限元模型自振頻率計算值及實測值相對誤差見表3，振型見圖8。為分析有限元模型計算振型與試驗振型的相關性，采用模態置信準則MAC值(表3)進行驗證，MAC值計算公式為

(11)

式中：φai，φei分別為第i階振型計算及實測振型向量。

表3 砌體不同材料模型下的結構固有頻率和MAC值表

圖8 砌體結構有限元計算振型

由表3看出：① 無論采用何種材料模型，前4階自振頻率相對誤差不一致，第1階頻率相對誤差均較小，第2階頻率相對誤差均較大，各向同性模型相對誤差最大，第3、4階相對誤差為中間值；② 橫觀各向同性及正交各向異性模型前4階自振頻率值相差不大，且總相對誤差小于各向同性模型；③ 不同材料模型的MAC值差別不大，第2階振型相關性最高，第1、3階相關程度均大于77%，說明前3階振型相關性良好，而第4階實測振型不明顯，誤差較大，故MAC值較小。

綜上分析，砌體采用橫觀各向同性及正交各向異性模型模態分析結果與實測值更接近，說明各向異性模型能更好反映砌體結構的動力特性，但有限元模型的固有頻率與實測值存較大差距，需采用基于結構固有頻率的參數靈敏度分析修正。

3 砌體結構有限元模型修正

3.1 結構有限元模型修正方法

砌體結構模型修正為基于OMA模態試驗所得模態參數對結構有限元模型材料的物理、力學參數修正。通過對結構有限元模型固有頻率參數靈敏度分析，選靈敏度大的材料參數作為修正參數，用ANSYS軟件的優化方法進行不斷修正，直到有限元計算值與試驗值誤差最小為止。修正結構模型時先用隨機搜索法初步確定最優設計序列，并以此為初始值用零階或一階算法進行迭代修正。滿足收斂條件時優化迭代結束，所得全局最優設計序列即為材料修正參數的最佳值。

3.2 參數靈敏度分析及修正參數選擇

參數靈敏度分析[15]指將結構模態參數表示為模型物理、力學參數的函數。設F(pm)為關于pm=(1,2,…,n)的多元函數，l階微分靈敏度及差分靈敏度統稱為F對pm的l階靈敏度：

(12)

(13)

式中：F為實測模態參數；pm為材料物理、力學參數，也可為結構幾何尺寸。

本文砌體采用三種不同材料模型分別進行結構固有頻率的參數靈敏度分析及模型修正，以求獲得最優砌體材料模型及參數修正值。砌體各向同性模型參數包括密度(MD)、彈性模量(ME)、泊松比(MP)3個獨立參數；砌體橫觀各向同性模型材料參數包括密度(MD)、彈性模量(MEX、MEZ)、剪切模量(MGXZ)、泊松比(MPXY、MPXZ)6個獨立參數；砌體正交各向異性模型材料參數包括密度(MD)、彈性模量(MEX、MEY、MEZ)、剪切模量(MGXY、MGYZ、MGXZ)、泊松比(MPXY、MPYZ、MPXZ)10個獨立參數；混凝土用各向同模型材料參數包括密度(CD)、彈性模量(CE)及泊松比(CP)3個獨立參數。

本文采用差分靈敏度計算公式，每次計算只改變一個材料參數值，使其增大20%，其它參數不變，計算材料參數對結構前4階固有頻率(FREQ1～FREQ4)靈敏度，以百分比表示。砌體結構固有頻率參數靈敏度分析結果見圖9～圖11。由三圖看出，① 砌體密度、混凝土密度及泊松比的增大使結構頻率降低，其它參數增大使結構頻率增大；②砌體采用各向同性模型時對結構頻率影響的參數靈敏度大小依次為ME、MD、CD，其它參數影響較小；橫觀各向同性模型時對結構頻率影響的參數靈敏度大小依次為MD、MGXZ、MEZ、CD、CE、MEX，其它參數影響較小；砌體采用正交各向異性模型時對結構頻率影響的參數靈敏度大小依次為MD、MGXZ、 MEZ、CD、CE、MEX、MGYZ，其它參數影響較小。有限元模型修正參數應優先選擇靈敏度較大的材料參數。

圖9 各向同性模型參數靈敏度分析

3.3 設計變量、狀態變量及目標函數

采用ANSYS軟件優化設計時需確定設計變量、狀態變量及目標函數。設計變量據參數靈敏度分析結果確定：砌體采用各向同性模型時設計變量為ME、MD、CD，砌體采用橫觀各向同性模型時設計變量為MD、MGXZ、MEX、MEZ、CD、CE，砌體采用正交各向異性模型時設計變量為MD、MGXZ、MGYZ、 MEZ、MEX、CD、CE，設計變量約束條件按表1中材料參數取值范圍。狀態變量為1～4階計算頻率(FREQ1～FREQ4)，分別增大、減小1～4階實測頻率值的10%作為狀態變量約束條件，MAC值作為驗證計算、實測振型的相關性參考值，滿足MACi=1～3≥75%，MAC1i=4≥60%。目標函數用于評價模型修正效果，用前4階有限元計算頻率與實測頻率相對誤差的平方和作為目標函數，因其修正效率更高、修正效果更明顯[16]，即

(14)

式中：MBHS為目標函數；fei為試驗所得頻率；fai為有限元分析所得頻率。

3.4 模型修正結果分析

用ANSYS程序優化算法進行模型修正時，優化迭代收斂條件為：當前設計序列與前一設計序列的目標函數值小于目標函數容差ε：

|MBHS(j)-MBHS(j-1)|≤ε

(15)

最佳合理設計序列與當前設計序列目標函數值小于目標函數容差ε：

|MBHS(j)-MBHS(b)|≤ε

(16)

式中：j為迭代次數；ε為較小數：b為最佳設計序列時迭代次數。

砌體結構模型修正后，三種不同模型設計變量及目標函數計算結果及與初始值的差值見表4，計算頻率與實測頻率相對誤差及MAC值見表5。

表4 有限元模型設計變量及目標函數修正結果

表5 模型修正后有限元計算與試驗模態參數比較表

由表4、表5看出，① 結構有限元模型修正采用ANSYS優化設計模塊多次迭代分析完成，模型修正以目標函數收斂于最小值作為修正目標，故模型修正結果使有限元模型計算所得模態參數總體與實測值趨于一致，不易實現計算與實測每階模態參數均保持一致或計算值與實測值完全一致。對三種不同材料模型有限元模型修正后目標函數值及計算得結構頻率與實測頻率相對誤差均有減小，且以正交各向異性模型目標函數值為最小。② 有限元模型中需修正的參數通常為難以準確確定的參數，如砌體密度、彈性模量、剪切模量。由修正后材料參數改變量看出，砌體密度、彈性模量及剪切模量改變較大。而在各向異性模型中，正交各向異性模型砌體彈模及剪模的改變量總體小于橫觀各向同性模型，原因為正交各向異性性模型考慮砌體3個主軸向不同力學參數值，相比橫觀各向同性模型及各向同性模型而言其修正效果更好。③有限元模型修正后MAC值略有提高，但提高量并不明顯，因此模型修正時MAC值驗證計算振型與實測振型相關性可起一定參考作用。④ 有限元模型修正后，雖計算頻率與實測頻率相對誤差均有所減小，但部分階次相對誤差仍較大，如第2階頻率，其原因可能為有限元模型未考慮房屋內家具、物品等，測試中存在噪音、干擾信號等，造成有限元模型與實際結構邊界條件存在一定誤差。

4 結 論

本文通過對砌體結構有限元模型建模及基于OMA模態參數進行結構模型修正研究，結論如下：

(1)采用隨機子空間法(SSI)獲得二層砌體結構前4階固有頻率位于9～25 Hz。

(2)砌體采用3種材料模型建模時，橫觀各向同性及正交各向異性模型前4階固有頻率計算值較各向同性模型與實測值更接近，振型相關性更好。

(3)對三種不同材料有限元模型修正后，目標函數值均有減小，計算所得結構頻率與實測頻率相對誤差均有減小，且正交各向異性模型目標函數值最小。MAC值可驗證計算振型與實測振型的相關性，有一定參考作用。

(4)有限元模型修正中由于正交各向異性模型考慮砌體3個主軸向不同力學參數值，較橫觀各向同性模型及各向同性模型修正效果更好。

[1] 夏益霖．基于試驗模態參數的結構有限元模型修正[J]．振動與沖擊,1993,12(1):61-65．

XIA Yi-lin.Finite element structural model updating based on experiment modal parameters[J].Journal of Vibration and Shock, 1993, 12(1): 61-65.

[2] 韓芳,鐘冬望,汪君．基于貝葉斯法的復雜有限元模型修正研究[J]．振動與沖擊,2012,31(1):39-43．

HAN Fang，ZHONG Dong-wang，WANG Jun.Complicated finite element model updating based on Bayesian method[J].Journal of Vibration and Shock, 2012, 31(1): 39-43.

[3] 殷海濤,姜金輝,張方,等．基于試驗模態參數及結構動力學優化設計的有限元建模[J]．國外電子測量技術,2012,31(9)：18-22．

YIN Hai-tao, JIANG Jin-hui, ZhANG Fang, et al.Finite element modeling based on experimental modal parameters and structural dynamics optimization[J].Foreign Electronic Measurement Technology, 2012, 31(9):18-22.

[4] Wu Cheng-qing, Hao Hong.Derivation of 3D masonry properties using numerical homogenization technique[J].International Journal For Numerical Methods In Engineering, 2006, 66(11): 1717-1737.

[5] Wei Xue-ying, Hao Hong.Numerical derivation of homogenized dynamic masonry material properties with strain rate effects[J].International Journal of Impact Engineering, 2009, 36(3): 522-536.

[6] 仲繼壽,楊舜臣,肖躍軍．礦山采動區磚砌體房屋三維有限元模型的建立與分析[J]．中國礦業大學學報，1993，22(2)：40-47．

ZHONG Ji-shou, YANG Shun-chen, XIAO Yue-jun.The establishment of a three dimensional finite element model of brick masonry buildings subjected to mining ground deformation[J].Journal of China University of Mining & Technology, 1993, 22(2): 40-47．

[7] 魏海霞,陳士海,張安康．基于動力有限元方法的典型砌體結構爆破振動安全標準的探討[J]．振動與沖擊,2011,30(5)：49-53．

WEI Hai-xia, CHEN Shi-hai, ZHANG An-kang.Safety standards discussion for blasting vibration of typical masonry buildings with dynamic finite element method[J].Journal of Vibration and Shock, 2011, 30(5): 49-53.

[8] 劉富君,朱玉華,馬曉輝．內構造柱加固砌體墻抗震性能有限元分析[J]．結構工程師,2011,27(3):62-66．

LIU Fu-jun, ZHU Yu-hua, MA Xiao-hui.Finite element analysis for seismic performances of unreinforced masonry walls strengthened with inside structural column[J].Structural Engineers, 2011, 27(3): 62-66.

[9] 譚曉晶,吳斌,辛文杰．林甸縣農村砌體房屋抗震性能調查與分析[J]．建筑科學與工程學報,2012,29(2):36-42．

TAN Xiao-jing, WU Bin, XIN Wen-jie.Investigation and analysis of seismic behavior of masonry buildings in rural areas located in Lindian country[J].Journal of Architecture and Civil Engineering, 2012, 29(2): 36-42．

[10] GB 50003-2001,砌體結構設計規范[S]．

[11] GB 50010-2010,混凝土結構設計規范[S]．

[12] 肖祥,任偉新．實時工作模態參數數據驅動隨機子空間識別[J]．振動與沖擊,2009,28(8):148-153．

XIAO Xiang, REN Wei-xin.Improved data-driven stochastic subspace identification of online operational modal parameters[J].Journal of Vibration and Shock, 2009, 28(8): 148-153.

[13] Van O P, Moor D B.Subspace identification for linear systems: theory implementation applications[M].Dordrecht: The Netherlands: Kluwer Academic Publishers, 1996.

[14] 禹丹江．土木工程結構模態參數識別-理論、實現與應用[D]．福州：福州大學，2006．

[15] 沃德·海倫，斯蒂芬·拉門茲，波爾·薩斯,著．白化同，郭繼忠,譯．模態分析理論與試驗[M]．北京:北京理工大學出版社,2001．

[16] 卞廣為．大跨斜拉橋有限元模型修正的數值模擬與試驗研究[D]．哈爾濱：哈爾濱工業大學,2008．