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

陳羽中，張家濱，顏東煌

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

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

1 有限元模型修正方法

1.1 模態(tài)柔度矩陣

多自由度體系振動特征方程通?？杀硎緸椋?/p>

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

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

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

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

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

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

假設(shè)真實柔度矩陣為Da，初始有限元模型模態(tài)參數(shù)構(gòu)建的柔度矩陣為Do。通過試驗得出模型撓度為Ua，假設(shè)試驗值為真實值，由有限元模型計算響應(yīng)撓度為Uo。根據(jù)柔度矩陣與位移的關(guān)系，可得：

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

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

對應(yīng)的柔度矩陣間差值ΔD為：

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

1.2 靈敏度分析

基于靜載試驗的有限元模型修正，對結(jié)構(gòu)待修正參數(shù)進行靈敏度分析，柔度矩陣的靈敏度為有限元模型修正提供了快速搜索的方向。通過模態(tài)柔度法，推導模型關(guān)于某一單元參數(shù)Kj的靈敏度，Kj通常為結(jié)構(gòu)邊界支撐的剛度、模型的設(shè)計參數(shù)等。

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

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

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

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

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

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

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

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

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

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

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

利用最小二范數(shù)原理，構(gòu)造有限元模型修正的優(yōu)化目標函數(shù)：

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

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

聯(lián)立式（24）、（25），可得到每次迭代入式（24），即：

2 模型與試驗布置

本試驗以廣州南二環(huán)李家沙大橋為工程背景，采用（38+72+220+72+38）m 的橋跨布置，邊中跨比為0.5，橋梁橫向?qū)挾?0 m。該橋設(shè)置等效的節(jié)段縮尺模型，砼強度等級為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#。節(jié)段模型試驗現(xiàn)場如圖1 所示。模型主梁兩端支撐在彈簧上為彈性支撐，左側(cè)彈簧支撐為K1，右側(cè)彈簧支撐為K2。梁的左端設(shè)置軸承與止推裝置，并借助鋼橫截梁和混凝土支墩直接頂在左側(cè)剪力墻上，以限制縱橋向的位移。主梁兩側(cè)焊接2 個H 型鋼，以限制橫向位移。模型具體情況參考文獻[16]，截面示意圖如圖2 所示，模型結(jié)構(gòu)示意圖如圖3所示。

圖1 試驗?zāi)Ｐ偷牟贾肍ig.1 The layout of the test model

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

圖3 模型結(jié)構(gòu)示意Fig.3 The model structure diagram

2.1 加載方案

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

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

2.2 撓度測點布置及測試結(jié)果

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

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

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

3 算例

3.1 初始有限元模型

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

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

表2 斜拉索設(shè)計參數(shù)Table 2 The design parameters of stay cables

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

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

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

3.2 待修正參數(shù)的確定

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

表3 待修正參數(shù)初值Table 3 The initial value of parameters

3.3 修正前的撓度偏差

在建立的初始有限元模型中，分別輸入相應(yīng)工況后，獲得試驗結(jié)構(gòu)任意位置的撓度值，再將理論值和試驗值進行對比。工況一和工況二其結(jié)果各取5個關(guān)鍵測點與試驗測試值進行對比，分別見表4～5。由于測點數(shù)量較少，小于其計算響應(yīng)的位移數(shù)。因此，根據(jù)試驗數(shù)據(jù)擬合撓度曲線，擬合后的試驗位移數(shù)與計算響應(yīng)的位移數(shù)一致，在撓度變化曲線中更清晰呈現(xiàn)出理論值與試驗值的撓度偏差，擬合后位移數(shù)為31個。2種工況的結(jié)果如圖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 修正的過程及結(jié)果

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

表6 待修正參數(shù)的修正結(jié)果Table 6 The correction results of parameters 103 kN·m-1

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

工況一和工況二作用下，撓度優(yōu)化結(jié)果取5個關(guān)鍵測點和試驗值作對比，其詳細撓度數(shù)值見表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中可以看出，有限元模型經(jīng)過修正后，理論撓度值與試驗實測出的撓度值基本一致，最大撓度誤差0.04 mm，遠小于未優(yōu)化前的有限元模型計算的撓度理論值與試驗值的差值。表明：基于靜載試驗的有限元模型修正方法是可行性的，可以滿足實際工程對精度的要求。

4 結(jié)論

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

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

3）通過2 種不同工況下的試驗數(shù)據(jù)，驗證了該方法的可行性，修正后的PPC 斜拉橋有限元模型值與靜載試驗實測值更吻合?；陟o載試驗的有限元模型修正技術(shù)，利用模態(tài)柔度原理有效修正了初始有限元模型，建立了該結(jié)構(gòu)基準有限元模型，具有效率高、擬合快、計算誤差較小等優(yōu)點，可滿足工程應(yīng)用要求。