基于靜載試驗的PPC斜拉橋有限元模型修正分析

陳羽中，張家濱，顏東煌

（長沙理工大學 土木工程學院，湖南 長沙 410114）

有限元模型修正理論已有五十多年的發展，許多學者對有限元模型修正技術提出了改進和創新，提高了有限元模型的可靠度和精確性，能更好反映結構真實情況和狀態，常運用于橋梁工程的實際應用中[1-2]。斜拉橋通常依據設計參數建立初始有限元模型。由于有限元模型中理想化的邊界條件與實際工程結構的真實邊界條件存在差異，導致有限元模型不能反映實際結構的真實狀態。因此，需要采用有限元模型修正技術對初始有限元模型進行修正，將誤差控制在工程允許范圍內。基于靜力試驗有限元模型修正，通常使用靜力試驗中測量出的撓度值或靜態應變與仿真結果的差值作為目標函數的設定依據。相對于動力試驗中測量出的固有振型等數據，靜載試驗的數據更準確、更可靠[3-4]。夏品奇等人[5]以新加坡的Safti橋為研究背景，建立了有限元模型，并進行了修正，修正中考慮斜拉索垂度效應，基于靈敏度修正后得到Safti橋可靠的有限元模型。宗周紅等人[6]通過對比其他學者的修正方法，提出了聯合模態柔度和靜力位移的有限元模型修正方法。李義強等人[7]利用最小二乘法構建了目標函數，以混凝土梁為研究對象，提出了基于參數識別的靜力模型修正方法。謝海龍[8]對一座先簡支后連續預應力混凝土箱梁橋進行荷載試驗，搜尋目標函數最優解，達到有限元模型修正目的，修正過程中舍棄了GN法，引入了阻尼因子LM 法。劉天怡等人[9]利用模態頻率、模態振型和模態柔度組合指標作為神經網絡修正的輸入參數，并以實例進行驗證，比較不同輸入參數的模型修正精度。模態柔度作為一種損傷識別新指標，被學者們廣泛運用在實際工程中[10-12]。楊秋偉等人[13]研究發現利用廣義柔度矩陣靈敏度進行損傷識別時，并非階次越高識別結果越準確，因而通常采用一次或二次廣義柔度矩陣。本研究擬利用振型、固有頻率及模態柔度關于索及支座剛度之間的靈敏度系數，推導基于柔度的靈敏度矩陣，基于靜載試驗的有限元模型修正方法對模型進行修正，并以多種實際工況驗證修正后的有限元模型仿真實際結構。

1 有限元模型修正方法

1.1 模態柔度矩陣

多自由度體系振動特征方程通常可表示為：

式中：K為結構的總剛矩陣；M為結構的總質量矩陣；ωr為結構在第r階的圓頻率；φr為結構在第r階的特征振型，r=1,2,···,n。

對振型質量歸一化，且式（1）滿足正交條件，即：

式中：Φ= （φ1,φ2,···,φn）。

柔度矩陣D通常是剛度矩陣K的廣義逆矩陣，文獻[14]定義了柔度矩陣D和剛度矩陣K。運用模態參數的表達式為：

式中：φr是進行歸一化處理后的第r階振型向量。

對比式（3）和式（4）可知，柔度矩陣與頻率的平方成反比，隨著模態階數越高，模態對柔度的影響越小。實際工程中，有限元對高階模態仿真結果較差，因此，利用高階模態對柔度矩陣影響微弱的特性，在計算過程中只取有限元模型中低階模態參數構建柔度矩陣D進行模型校驗。

假設真實柔度矩陣為Da，初始有限元模型模態參數構建的柔度矩陣為Do。通過試驗得出模型撓度為Ua，假設試驗值為真實值，由有限元模型計算響應撓度為Uo。根據柔度矩陣與位移的關系，可得：

式中：f為荷載向量；Ua、Uo分別為位移向量。

將式（5）、（6）相減，即為試驗結果與仿真計算結果間的誤差。

對應的柔度矩陣間差值ΔD為：

通過對柔度矩陣與待修正參數間的靈敏度分析，得到撓度差值與待修正參數的關系表達式，實現有限元模型修正。

1.2 靈敏度分析

基于靜載試驗的有限元模型修正，對結構待修正參數進行靈敏度分析，柔度矩陣的靈敏度為有限元模型修正提供了快速搜索的方向。通過模態柔度法，推導模型關于某一單元參數Kj的靈敏度，Kj通常為結構邊界支撐的剛度、模型的設計參數等。

求解結構體系的固有頻率就是求解多自由度振動特征方程的特征解λ，λ=ω2，則為結構在第r階的特征值，r=1,2,···,n，將其帶入式（3），即：

將式（9）對單元參數Kj求偏導，得：

本研究修正目標為斜拉索的剛度和邊界支撐剛度，質量矩陣與斜拉索和邊界支承剛度無關，所以式（11）中，整理式（11）得：

根據差分法，柔度矩陣對單元參數Kj的導數可改寫為：

根據式（5）、（6），撓度U對單元參數Kj的靈敏度分析，得到：

在式（15）中，力向量f與Kj無關，因此，整理得到單元參數K與撓度U的關系j式為：

式（17）中ΔU與ΔKj的具體表達式為：

式（20）給出撓度殘差與單元參數修正值的關系表達式。每次修正后，迭代計算理論撓度向量，并再次計算新的撓度殘差，循環迭代。

修正有限元模型，撓度殘差ΔU需無限趨近于0，則代表有限元模型與真實結構相一致，聯立式（5）～（8），整理可得：

當ΔU= 0 時，表明：理論撓度值等于試驗撓度值，且理論柔度矩陣等于真實柔度矩陣。

為實現有限元模型能仿真實際結構，使理論撓度值與試驗撓度值差別最小。ΔU為一列向量，實質上是模型n個節點上的殘差位移，可表示為：

利用最小二范數原理，構造有限元模型修正的優化目標函數：

撓度修正實質上是對柔度矩陣修正，撓度殘差最小化實質上是模態柔度殘差最小化。撓度殘差在單元參數ΔKj附近進行攝動，對其過程進行一階泰勒級數線性展開，可得：

式（25）中，ΔUk為第k次迭代時的位移殘差；為迭代時靈敏度矩陣，每次迭代時都會發生變化；ΔKj為迭代時的單元參數增量。

聯立式（24）、（25），可得到每次迭代入式（24），即：

2 模型與試驗布置

本試驗以廣州南二環李家沙大橋為工程背景，采用（38+72+220+72+38）m 的橋跨布置，邊中跨比為0.5，橋梁橫向寬度50 m。該橋設置等效的節段縮尺模型，砼強度等級為C50，混凝土彈性模量為 3.5×104 MPa，混凝土密度為 2.55×103 kg·m3，泊松比為 0.3，K1 剛度為 3.39×103 kN/m，K2 剛度為6.71×103 kN/m，斜拉索彈性模量為1.95×105 MPa，斜拉索截面面積為1.39×10-4m2，10#索的長度為12.58 m，11#索的長度為13.49 m，12#索的長度為14.42 m，13#索的長度為15.34 m，14#索的長度為16.28 m。模型上張拉的5 根斜拉索編號從左至右分別為10#～14#。節段模型試驗現場如圖1 所示。模型主梁兩端支撐在彈簧上為彈性支撐，左側彈簧支撐為K1，右側彈簧支撐為K2。梁的左端設置軸承與止推裝置，并借助鋼橫截梁和混凝土支墩直接頂在左側剪力墻上，以限制縱橋向的位移。主梁兩側焊接2 個H 型鋼，以限制橫向位移。模型具體情況參考文獻[16]，截面示意圖如圖2 所示，模型結構示意圖如圖3所示。

圖1 試驗模型的布置Fig.1 The layout of the test model

圖2 模型橫截面示意（單位：mm）Fig.2 The cross-section（unit:mm）

圖3 模型結構示意Fig.3 The model structure diagram

2.1 加載方案

本試驗以撓度殘差作為目標函數，對初始有限元模型進行修正。因此，在部分預應力混凝土（partially prestressed concret， 簡稱為PPC）斜拉橋縮尺模型上進行靜力荷載試驗，收集記錄各測點的撓度值。在模型跨中270 m 處，以2 種集中荷載方式進行加載：

工況一，選取16 kN 集中荷載作用；工況二，選取12 kN集中荷載作用。

2.2 撓度測點布置及測試結果

試驗中，以梁左端支座為原點，在0、1.35、2.70、4.40、6.00 m 處，分別布置5 個動態位移數據采集點，其余8 個靜態位移計按照位置均勻布置，位移計精度為0.01 mm，布置情況如圖4所示。以2種工況的數據為實例進行修正，試驗中實測撓度值見表1。撓度方向向下為正值。

圖4 位移計布置情況Fig.4 The layout of the displacement gauge

表1 各個測點撓度值變化Table 1 The change of deflection of measuring points

3 算例

3.1 初始有限元模型

根據模型的設計參數和竣工圖紙，采用MATLAB 軟件，按照結構的實際尺寸建立PPC 斜拉橋有限元模型。模型中的梁段部分，單元剛度矩陣采用平面彎曲梁單元，單元質量矩陣采用協調質量矩陣。掛籃體系中，單元質量矩陣采用集中質量矩陣。斜拉索由軸向桿單元構建。模型中，梁段部分、支撐邊界和斜拉索聯接處共用節點實現共同協調受力。

斜拉索的設計參數見表2，彈性模量E=1.95×105MPa，橫截面積A=1.39×10-4m2，根據斜拉索的剛度計算公式K=EA/L，可計算出斜拉索的初始剛度。5根斜拉索施加在梁單元上的初始剛度見表2。

表2 斜拉索設計參數Table 2 The design parameters of stay cables

試驗模型的設計方案中，左側的邊界情況較為復雜，除了豎向彈簧支撐K1，在梁的左端還設置軸承與止推裝置；右側的邊界條件主要是豎向彈簧支撐K2；梁的左右兩側焊接了H 型鋼，限制梁的橫向位移，如圖5所示。

圖5 支座布置Fig.5 The arrangement of supports

根據結構的特點、邊界支撐和斜拉索的布置情況，考慮到集中荷載在加載中的位置和有限元模型的修正目標，將梁段模型結構劃分為32 個有限元模型單元，33 個節點，單元長度并不完全一致，單元劃分情況如圖3所示。

3.2 待修正參數的確定

在進行試驗時，兩側彈簧支座放置較久，兩側彈簧支撐的剛度可能發生改變，且邊界情況較為復雜。本試驗將邊界條件等效為兩側，僅有彈性支撐K1、K2 和縱橋向的水平約束。斜拉索受到重力作用出現垂度效應，導致設計參數與真實值有一定誤差。因此，模型修正的對象是5根斜拉索的剛度和兩側豎向彈性支撐K1、K2 的剛度。其中，為保證斜拉索剛度的修正符合實際情況，具有物理意義，對索的剛度修正主要是對索的設計參數EA進行修正。

表3 待修正參數初值Table 3 The initial value of parameters

3.3 修正前的撓度偏差

在建立的初始有限元模型中，分別輸入相應工況后，獲得試驗結構任意位置的撓度值，再將理論值和試驗值進行對比。工況一和工況二其結果各取5個關鍵測點與試驗測試值進行對比，分別見表4～5。由于測點數量較少，小于其計算響應的位移數。因此，根據試驗數據擬合撓度曲線，擬合后的試驗位移數與計算響應的位移數一致，在撓度變化曲線中更清晰呈現出理論值與試驗值的撓度偏差，擬合后位移數為31個。2種工況的結果如圖6～7所示。

表4 工況一作用下的撓度對比Table 4 The comparison of deflection under working condition one mm

表5 工況二作用下的撓度對比Table 5 The comparison of deflection under working condition two mm

圖6 工況一的撓度Fig.6 The deflection curve diagram under the action of working condition one

圖7 工況二的撓度Fig.7 The deflection curve diagram under the action of working condition two

3.4 修正的過程及結果

根據目標函數，利用MATLAB 軟件編寫設計參數修正程序，對初始有限元模型進行修正，7 個參數βi（i= 1,2,···,7）為待修正參數，進行迭代計算，得出最優解。迭代過程中，由于目標函數達到最優解，迭代結束，共迭代113次，計算響應出左側彈簧支撐K1、右側彈簧支撐K2 和5 根斜拉索剛度的修正結果，見表6。目標函數變化曲線如圖8所示。

表6 待修正參數的修正結果Table 6 The correction results of parameters 103 kN·m-1

圖8 目標函數迭代趨勢Fig.8 The iterative trend graph of objective function

工況一和工況二作用下，撓度優化結果取5個關鍵測點和試驗值作對比，其詳細撓度數值見表7～8。修正后的撓曲線和擬合后的試驗值曲線進行對比，如圖9～10所示。

表7 工況一作用下的修正后撓度對比Table 7 The comparison of modified deflection under the action of working condition one mm

表8 工況二作用下修正后撓度對比Table 8 The comparison of modified deflection under the action of working condition two mm

圖9 修正后工況一作用下的撓度Fig.9 The deflection curve diagram under the action of working condition one after modification

圖10 修正后工況二作用下的撓度Fig.10 The deflection curve diagram under the action of working condition two after modification

從圖9～10中可以看出，有限元模型經過修正后，理論撓度值與試驗實測出的撓度值基本一致，最大撓度誤差0.04 mm，遠小于未優化前的有限元模型計算的撓度理論值與試驗值的差值。表明：基于靜載試驗的有限元模型修正方法是可行性的，可以滿足實際工程對精度的要求。

4 結論

1）根據設計方案或者竣工圖紙中的設計參數，建立PPC 斜拉橋縮尺模型的初始有限元模型。在相應工況下，計算的理論值與靜載試驗的實測值有一定誤差。表明：初始有限元模型需要修正后才能反映試驗中的實際結構，滿足工程校驗的要求。

2）運用低階模態參數構建柔度矩陣，因柔度矩陣在推導過程中與待修正參數直接相關，所以可以直接建立撓度變形與參數之間的靈敏度關系。PPC 斜拉橋縮尺模型的待修正參數，主要有5 根斜拉索的剛度和左、右兩側彈簧支撐K1、K2 的剛度。

3）通過2 種不同工況下的試驗數據，驗證了該方法的可行性，修正后的PPC 斜拉橋有限元模型值與靜載試驗實測值更吻合。基于靜載試驗的有限元模型修正技術，利用模態柔度原理有效修正了初始有限元模型，建立了該結構基準有限元模型，具有效率高、擬合快、計算誤差較小等優點，可滿足工程應用要求。