胡文哲 崔 闖 王 昊 張清華

(西南交通大學(xué)橋梁工程系,成都 610031)

隨著我國交通強(qiáng)國和海洋強(qiáng)國等國家戰(zhàn)略的推進(jìn),大跨度橋梁結(jié)構(gòu)越來越廣泛應(yīng)用于跨江、跨海等重要交通要道中。同時(shí),由于極端環(huán)境情況頻發(fā),大跨或超大跨橋梁力學(xué)性能分析的需要和要求逐步提升。但實(shí)際中,數(shù)值仿真模型與實(shí)際橋梁難免存在誤差,導(dǎo)致數(shù)值模型不能準(zhǔn)確地表征實(shí)際橋梁的行為[1]。通常情況下,模型計(jì)算響應(yīng)與實(shí)際結(jié)構(gòu)響應(yīng)之間的誤差可以分為以下兩種：1)在建模過程中,為簡化計(jì)算過程引入理想化假定造成的誤差;2)由于結(jié)構(gòu)在長期使用過程中受到眾多不利因素的影響,使得模型計(jì)算得到的結(jié)構(gòu)響應(yīng)與實(shí)際響應(yīng)間存在誤差,此誤差來源于模型結(jié)構(gòu)誤差、參數(shù)誤差以及階次誤差[2-3]。當(dāng)上述誤差超過閾值時(shí),認(rèn)為此模型無法正確反映真實(shí)結(jié)構(gòu)特性。有限元模型修正可以通過不斷調(diào)整設(shè)計(jì)參數(shù),使基于有限元模型計(jì)算的結(jié)構(gòu)特性值和實(shí)際結(jié)構(gòu)測量的真實(shí)觀測特性值差異盡可能減小。

目前常用的修正方法主要有矩陣型修正和參數(shù)型修正[4]。其中,前者修正過程拋棄了各類參數(shù)原本的物理含義,僅對數(shù)值進(jìn)行優(yōu)化,在實(shí)際場景分析中實(shí)用性較低;后者修正過程以模型某些物理參數(shù)作為修正目標(biāo),在工程分析中應(yīng)用廣泛[5]。傳統(tǒng)參數(shù)型方法計(jì)算量大且求解過程中易出現(xiàn)病態(tài)解,常見的代理模型算法包括神經(jīng)網(wǎng)絡(luò)、多項(xiàng)式響應(yīng)面和Kriging模型等[6-7]。馬印平等結(jié)合實(shí)測弦桿應(yīng)力數(shù)據(jù)并使用響應(yīng)面法對一座鋼管混凝土組合桁梁橋的全橋多尺度模型進(jìn)行了材料參數(shù)與車輛荷載的修正[4];唐煜等使用人工蜂群算法,以某斜拉橋的自振模態(tài)與靜載響應(yīng)轉(zhuǎn)角為目標(biāo)函數(shù),對塔梁連接處局部剛度進(jìn)行了修正[8];HA-SANCEBI等采用神經(jīng)網(wǎng)絡(luò)方法基于實(shí)測數(shù)據(jù)對某老化T形梁橋有限元模型進(jìn)行了修正[9];王蕾等提出一種基于徑向基函數(shù)神經(jīng)網(wǎng)絡(luò)的有限元模型修正方法,對預(yù)應(yīng)力混凝土大跨剛構(gòu)-連續(xù)梁橋的設(shè)計(jì)參數(shù)進(jìn)行了修正[10]。對于多尺度有限元模型,響應(yīng)面法無法全面考慮全局模型與局部模型之間的相互影響。針對此問題,文章提出一種使用Kriging元模型,基于R2指標(biāo)的差分高維多目標(biāo)進(jìn)化(Multi-Objective Evolutionary Algorithm,MOEA)算法[11-13],可以在考慮變量沖突關(guān)系的同時(shí)完成模型修正,并保證修正精度。

此處以某大跨度鋼斜拉橋?yàn)檠芯繉ο?通過使用殼/實(shí)體單元對結(jié)構(gòu)關(guān)鍵構(gòu)件進(jìn)行詳細(xì)建模生成局部模型,而使用簡單梁/桁架單元模擬其他構(gòu)件生成全局模型,并在二者邊界處滿足兼容耦合條件,形成結(jié)構(gòu)的多尺度模型,結(jié)合Kriging元模型提出了可反映荷載與結(jié)構(gòu)響應(yīng)映射關(guān)系的代理模型,并采用R2-MOEA多目標(biāo)優(yōu)化算法對有限元模型進(jìn)行了修正,以期為結(jié)構(gòu)運(yùn)維安全提供理論支持。

1 模型修正方法

基于Kriging元模型的R2-MOEA修正方法可分為兩個(gè)主要步驟：1)使用靈敏度分析確定出對目標(biāo)函數(shù)影響最顯著的材料參數(shù),使用拉丁超立方采樣方法建立參數(shù)樣本集,并代入初始多尺度模型計(jì)算得到對應(yīng)目標(biāo)響應(yīng)值,以此形成有效數(shù)據(jù)集建立Kriging代理模型,滿足精確度條件后進(jìn)行后續(xù)優(yōu)化計(jì)算;2)使用多目標(biāo)進(jìn)化與演化控制算法相結(jié)合,迭代計(jì)算出與實(shí)測數(shù)據(jù)更加接近的有限元模型最優(yōu)參數(shù),完成模型修正。全流程見圖1。

圖1 基于Kriging元模型的R2-MOEA多目標(biāo)優(yōu)化全流程Fig.1 The whole process of multi-objective optimization based on Kriging meta model R2-MOEA

1.1 基于Kriging元模型的結(jié)構(gòu)響應(yīng)代理模型建立

Kriging元模型是一種非線性插值模型[14],在單目標(biāo)問題中,其主要形式如式(1),由一個(gè)基本多項(xiàng)式擬合函數(shù)fT(x)β和一個(gè)隨機(jī)誤差函數(shù)z(x)組成：

y(x)=fT(x)β+z(x)

(1)

式中：β為多項(xiàng)式函數(shù)的系數(shù)矩陣;誤差函數(shù)z(x)服從正態(tài)分布N(0,σ2),并且滿足以下關(guān)系：

cov[z(xi),z(xj)]=σ2R(xi,xj)

(2)

式中：R(xi,xj)表示樣本點(diǎn)xi與xj之間的相關(guān)度,其相關(guān)函數(shù)R由二階多項(xiàng)式函數(shù)及高斯函數(shù)確定：

(3)

式中：n為隨機(jī)變量個(gè)數(shù);θk需要根據(jù)樣本值由最大似然估計(jì)確定。假定樣本集S=[x1,x2,…,xn]和其對應(yīng)的有限元模型計(jì)算值Y=[y1,y2,…,yn],θk可由式(4)計(jì)算：

(4)

Kriging元模型的預(yù)測值根據(jù)最佳線性無偏估計(jì)方法確定：

(5)

(6)

研究中采用Latin-Hypercube取樣方法建立樣本集[15],以滿足樣本的隨機(jī)性與樣本規(guī)模的可控性。Kriging元模型的精度采用均方誤差MSE與R2確定：

(7a)

(7b)

式中：Var表示數(shù)組yi的方差值。R2取值范圍是[0,1],當(dāng)模型誤差可以忽略時(shí),R2=1。

結(jié)合建立的有限元模型,采用Kriging元模型建立荷載與結(jié)構(gòu)響應(yīng)之間的映射關(guān)系,以建立基于Kriging元模型的結(jié)構(gòu)響應(yīng)代理模型。

1.2 R2-MOEA多目標(biāo)優(yōu)化算法

參照橋梁監(jiān)測系統(tǒng)實(shí)際測量數(shù)據(jù),選取橋梁結(jié)構(gòu)宏觀響應(yīng)和局部響應(yīng)所對應(yīng)的結(jié)構(gòu)物理參量為修正對象。以某種類型響應(yīng)指標(biāo)i為例(應(yīng)力指標(biāo)、位移指標(biāo)等),模型修正旨在使得模型響應(yīng)與實(shí)際響應(yīng)更加接近,即響應(yīng)i的模型計(jì)算值FEi與實(shí)際測量值TEi之間的殘差ei最小：

(8)

式中：Fi為單種類型目標(biāo)響應(yīng)i的殘差和值。

模型修正問題可表達(dá)為含有n類響應(yīng)目標(biāo)函數(shù)的多目標(biāo)優(yōu)化問題：

F(x)={F1(x),F2(x),F3(x),…,Fn(x)}

(9)

其中x∈(xmin,xmax)

在多目標(biāo)問題中,若在有限范圍內(nèi)的兩組參數(shù)值x,y滿足式(10)關(guān)系,則稱x支配y,也即x更接近最優(yōu)解。

(10)

由此關(guān)系在樣本范圍內(nèi)循環(huán)計(jì)算可得出最優(yōu)支配解,該解稱為Pareto最優(yōu)解x*,其對應(yīng)的目標(biāo)函數(shù)值稱為Pareto前沿[16],根據(jù)x*的定義,應(yīng)存在以下推論：

1) 解集x*不被任何其他解集支配;

2) 解集x*中的任意兩個(gè)解之間不相互支配。

在實(shí)際計(jì)算中,準(zhǔn)確的Pareto最優(yōu)解是很難得到的,常用其近似集AF代替Pareto解集[17-19]。為了評估及控制算法精度,設(shè)置目標(biāo)向量理想值為零向量組f,優(yōu)化算法的收斂性則為近似解集AF的誤差R2。為表明各目標(biāo)函數(shù)對于優(yōu)化問題所占權(quán)重,引入權(quán)重系數(shù)矩陣w,采用標(biāo)準(zhǔn)加權(quán)切比雪夫函數(shù)定義[20]：

R2(AF)=minx∈A{max{wi|fi-Fi(x)|,

i=1,2,…,n}}

(11)

假定AF集維度為m?1,根據(jù)比較目標(biāo)值之間的支配性,可得到最優(yōu)解xopt：

(12a)

(12b)

優(yōu)化過程表現(xiàn)為反復(fù)迭代取舍的過程,所需計(jì)算次數(shù)極其龐大,需采用建立的代理Kriging模型反復(fù)迭代。由于誤差限的設(shè)定及考慮到Kriging模型建立過程的反復(fù)迭代,導(dǎo)致Kriging元模型與有限元模型間的誤差不可避免,使用演化控制可減小此誤差。當(dāng)樣本組預(yù)測值的均方根誤差RMSE較大時(shí),使用真實(shí)值代替Kriging預(yù)測值以重建代理模型提高精度[21]：

(13)

Rmin(xopt)=smin

(14)

該指標(biāo)s即可用于表示每組預(yù)測值的不確定性,設(shè)定此不確定性失效范圍為計(jì)算得到近似解集的前20%,認(rèn)為此范圍內(nèi)的預(yù)測值失去可靠性。為控制此循環(huán)過程的終止,根據(jù)式(13)與式(14)計(jì)算出最優(yōu)解的RMSE值Rmin,以此作為優(yōu)化過程評定指標(biāo),設(shè)置演化控制次數(shù)及Rmin界限值φ=1.0,當(dāng)0

2 基于實(shí)測數(shù)據(jù)的多尺度有限元模型修正

2.1 多尺度有限元模型建立

以某大跨度鋼斜拉橋作為研究對象,建立橋梁多尺度模型,橋梁立面及斷面如圖2所示。大橋主橋?yàn)殡p塔雙索面鋼箱梁斜拉橋。主橋跨徑布置為：48+204+460+204+48 m,主橋總長為964 m,橋面總寬度為38.8 m。鋼箱梁標(biāo)準(zhǔn)梁段長度為12 m,單邊由18根不同規(guī)格的斜拉索組成。斜拉橋索塔的塔頂標(biāo)高為171.000 m(黃海高程,下同),承臺頂面標(biāo)高為7.500 m,總高度為163.5 m。根據(jù)設(shè)計(jì)參數(shù),采用ANSYS有限元計(jì)算軟件建立的多尺度有限元模型如圖3所示。主梁采用正交異性鋼橋面板,其中頂板厚度為18 mm,U肋厚度為8 mm。

圖2 大跨斜拉橋斷面 mmFig.2 Section of the long-span cable-stayed bridge

圖3 大跨斜拉橋多尺度有限元模型Fig.3 Multi-scale finite element model of the long-span cable-stayed bridge

全橋共劃分43 486個(gè)結(jié)點(diǎn),48 792個(gè)單元。其中,跨中6.8 m節(jié)段為局部關(guān)注區(qū)域,通過板殼-梁單元的混合模型進(jìn)行模擬形成整體模型。局部節(jié)段模型與整體雙塔斜拉橋通過使用多點(diǎn)約束方法(MPC)連接形成多尺度模型。橋塔與橫梁、橫梁與輔助墩之間采用彈簧單元Combin40連接;橋墩與橫梁之間使用CP耦合,可確保局部鋼橋面板細(xì)觀模型和整體宏觀模型之間剛度連續(xù)和變形協(xié)調(diào)。

根據(jù)大跨斜拉橋的基本受力特性,為獲得結(jié)構(gòu)整體和局部響應(yīng),進(jìn)行現(xiàn)場結(jié)構(gòu)試驗(yàn),加載試驗(yàn)現(xiàn)場如圖5所示,荷載試驗(yàn)工況如表1所示。整體位移的數(shù)據(jù)取自布置在橋梁兩邊跨跨中及中跨1/4、中跨跨中、中跨3/4共5個(gè)位移測點(diǎn)(D1～D5);局部應(yīng)力數(shù)據(jù)取自全橋跨中截面頂、底板局部應(yīng)力測點(diǎn),各測點(diǎn)布置詳細(xì)位置如圖4所示。選取跨中斷面中線位置處測點(diǎn)L26、L30作為數(shù)據(jù)提取點(diǎn)。

表1 加載試驗(yàn)工況設(shè)置Table 1 Load test case settings

圖4 位移測點(diǎn)與跨中斷面局部應(yīng)力測點(diǎn)布置Fig.4 Arrangements of displacement measuring points and local stress measuring points across the interrupted surface

圖5 加載試驗(yàn)現(xiàn)場Fig.5 Load test site

對應(yīng)5種工況,各點(diǎn)位移及應(yīng)力數(shù)據(jù)可對應(yīng)設(shè)置目標(biāo)函數(shù)40個(gè),包括結(jié)構(gòu)體系前5階自振頻率、5種工況作用下頂板測點(diǎn)L26應(yīng)力、底板測點(diǎn)L30應(yīng)力、D1～D5測點(diǎn)位移。在有限元模型中計(jì)算時(shí),荷載時(shí)間步設(shè)置與試驗(yàn)加載方案相同。全局模型中施加節(jié)點(diǎn)荷載模擬工況一、二及工況四、五,對于工況三,荷載面加載在局部模型橋面板對應(yīng)位置模擬真實(shí)輪載。

2.2 參數(shù)選取與目標(biāo)函數(shù)分組

根據(jù)實(shí)際結(jié)構(gòu)受力特征,選取結(jié)構(gòu)頻率、位移響應(yīng)及應(yīng)力響應(yīng)作為優(yōu)化目標(biāo)函數(shù)。隨著橋梁服役期的增長,材料的固有屬性發(fā)生變化,對于整體結(jié)構(gòu)的模態(tài)分析以及局部結(jié)構(gòu)應(yīng)力均存在影響,選擇橋梁不同部位的材料參數(shù)作為更新參數(shù)。在多尺度模型中,由于對結(jié)構(gòu)構(gòu)件細(xì)部構(gòu)造的假定,例如對焊縫、連接件、加勁肋等部位的簡化,使得局部模型與實(shí)際細(xì)節(jié)存在差異,雖然對于全局結(jié)構(gòu)模態(tài)響應(yīng)以及位移響應(yīng)影響甚微,但會使局部構(gòu)造細(xì)節(jié)處的應(yīng)力、應(yīng)變偏離正常值。因此,在設(shè)定多尺度模型的更新參數(shù)時(shí),應(yīng)分別選取整體和局部更新參數(shù)。

初選取結(jié)構(gòu)不同位置處的彈性模量與密度作為更新參數(shù)。為方便說明,將橋面分為10 個(gè)節(jié)段,編號沿縱橋向,由于斜拉索布置呈雙向?qū)ΨQ,設(shè)置單個(gè)索面斜拉索由內(nèi)向外編號依次增大共16 根,為保證參數(shù)對該橋的代表性,全局模型中初選取參數(shù)共32 個(gè),分別為：主梁鋼材10 個(gè)節(jié)段、橋塔混凝土、橋墩混凝土處彈性模量與密度,8 根斜拉索彈性模量;局部模型中初選取參數(shù)共8 個(gè),分別為：頂板及其加勁肋、底板及其加勁肋的彈性模量與密度。采用靈敏度分析得到待定參數(shù)對目標(biāo)函數(shù)的影響程度如圖6,設(shè)定閾值選定對目標(biāo)函數(shù)影響最為顯著的19 個(gè)更新參數(shù)如表2所示。

表2 更新參數(shù)初始值匯總Table 2 Update parameter initial value summary

圖6 靈敏度分析Fig.6 Sensitivity analysis

根據(jù)結(jié)構(gòu)響應(yīng)類型將上述40 個(gè)目標(biāo)函數(shù)進(jìn)行分組,可以將目標(biāo)函數(shù)間的相互影響充分考慮,具體采用的兩種分組方案a)和b),其中Ff、Fd、Fs分別表示頻率、位移、應(yīng)力響應(yīng)對應(yīng)殘差和,方案b)中F2～F6分別表示工況一至工況五對應(yīng)的位移響應(yīng)殘差和：

a.F1=Ff,F2=Fd,F3=Fs;

為方便說明,使用RM-3和RM-7分別表示分組方案a)和b)。

2.3 多尺度有限元模型修正結(jié)果

在模型參數(shù)修正過程中,使用Latin-hypercube方法生成樣本集維度為2 000?19,即一共生成2 000組初始樣本,控制優(yōu)化過程精度的終止指標(biāo)φ=1.0。在此基礎(chǔ)上,利用基于Kriging元模型的R2-MOEA多目標(biāo)優(yōu)化算法對多尺度模型進(jìn)行修正,得到兩種不同分組方案下各目標(biāo)響應(yīng)的優(yōu)化結(jié)果。在對建立的Kriging元模型進(jìn)行精度檢驗(yàn)時(shí),最終根據(jù)式(7)得出Kriging模型MSE和R2指標(biāo)列于表3。

表3 Kriging元模型精度Table 3 Kriging metamodel accuracy

根據(jù)式(8)計(jì)算可知,Kriging元模型精度滿足要求,表示此模型可以代替多尺度有限元模型進(jìn)行后續(xù)計(jì)算。根據(jù)式(14)計(jì)算得到RM-3和RM-7兩種分組方案的優(yōu)化指標(biāo)Rmin數(shù)值列于表4。

表4 不同分組方案下Rmin值Table 4 Rmin values under different grouping schemes

由前述可以得知,Rmin表征最優(yōu)解與實(shí)測結(jié)果的誤差大小,根據(jù)表4可以看出,在兩種分組方案下Rmin計(jì)算值均在有限范圍[0,1]內(nèi),優(yōu)化效果較好,且Rmin(RM-3)

表5 設(shè)計(jì)參數(shù)最優(yōu)解Table 5 Optimal solutions of design parameters

由表5可知,主梁彈性模量變化沿縱向關(guān)于跨中呈對稱關(guān)系,跨中主梁容重增大,橋塔混凝土容重減小,彈性模量增大;斜拉索彈性模量增減幅度相同;局部模型中,頂板與底板部位彈性模量減小,加勁肋部位彈性模量增大。修正后的設(shè)計(jì)參數(shù)在可接受范圍內(nèi)。根據(jù)最優(yōu)解改變初始有限元模型的設(shè)計(jì)參數(shù),計(jì)算得到修正后模型的結(jié)構(gòu)各類響應(yīng)值。初始模型、實(shí)測、修正模型前三階自振頻率及各測點(diǎn)最大位移、應(yīng)力響應(yīng)值對比如圖7～圖9所示,修正結(jié)果對比如表6所示。

表6 最大位移與應(yīng)力響應(yīng)誤差值Table 6 Maximum displacement and stress response error values

圖7 RM-7分組自振頻率修正結(jié)果Fig.7 Correction results of natural frequency values for group RM-7

圖8 RM-7分組位移最大響應(yīng)值修正結(jié)果Fig.8 Correction results of maximum displacement response values for group RM-7

圖9 RM-7分組頂板與底板最大應(yīng)力響應(yīng)值修正結(jié)果Fig.9 Correction results of maximum stress response values of top and bottom plates of group RM-7

研究結(jié)果表明,模型修正后的自振頻率、位移響應(yīng)值和應(yīng)力響應(yīng)值與實(shí)測值之間的誤差相較于修正前的初始值顯著降低。修正后的全局自振頻率平均相對誤差由6.59%降低至3.10%;修正后的全局模型位移最大相對誤差為8.07%,相較于初始模型的18.27%,降低了10.20%,全局位移平均相對誤差由8.00%降低至3.55%;修正后局部模型的應(yīng)力最大相對誤差為1.61%,相較于初始模型的6.77%,降低了5.16%,局部模型應(yīng)力平均相對誤差由5.50%降低至1.06%?？傮w而言,采用基于Kriging元模型和R2-MOEA優(yōu)化算法適用于大跨度斜拉橋多尺度有限元模型修正。

3 結(jié)束語

1)基于Kriging元模型和R2-MOEA多目標(biāo)優(yōu)化算法建立了大跨度斜拉橋多尺度有限元模型修正方法,并結(jié)合實(shí)際橋梁健康監(jiān)測的實(shí)測值,驗(yàn)證了所提方法的正確性和適用性。

2)所建立的多尺度有限元模型更新后的響應(yīng)值與實(shí)測值誤差大幅降低,其中更新后的自振頻率誤差為3.10%,結(jié)構(gòu)位移響應(yīng)誤差為3.55%,應(yīng)力響應(yīng)誤差為1.06%。