陳羽中，張家濱，顏東煌

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

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

1 有限元模型修正方法

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

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

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

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

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

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

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

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

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

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

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

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

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

1.2 靈敏度分析

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

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

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

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

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

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

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

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

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

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

當(dāng)ΔU= 0 時(shí)，表明：理論撓度值等于試驗(yàn)撓度值，且理論柔度矩陣等于真實(shí)柔度矩陣。

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

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

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

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

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

2 模型與試驗(yàn)布置

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

圖1 試驗(yàn)?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 加載方案

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

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

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

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

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

表1 各個(gè)測點(diǎn)撓度值變化Table 1 The change of deflection of measuring points

3 算例

3.1 初始有限元模型

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

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

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

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

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

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

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

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

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

3.3 修正前的撓度偏差

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

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

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

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

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

工況一和工況二作用下，撓度優(yōu)化結(jié)果取5個(gè)關(guān)鍵測點(diǎn)和試驗(yàn)值作對(duì)比，其詳細(xì)撓度數(shù)值見表7～8。修正后的撓曲線和擬合后的試驗(yàn)值曲線進(jìn)行對(duì)比，如圖9～10所示。

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

表8 工況二作用下修正后撓度對(duì)比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)過修正后，理論撓度值與試驗(yàn)實(shí)測出的撓度值基本一致，最大撓度誤差0.04 mm，遠(yuǎn)小于未優(yōu)化前的有限元模型計(jì)算的撓度理論值與試驗(yàn)值的差值。表明：基于靜載試驗(yàn)的有限元模型修正方法是可行性的，可以滿足實(shí)際工程對(duì)精度的要求。

4 結(jié)論

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

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

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