基于數(shù)據(jù)同化參數(shù)反演的華光潭拱壩結(jié)構(gòu)分析

林咸志,羅增浤,許 正,閔皆昇,曹傳波,趙 權(quán)

(1.浙江浙能華光潭水力發(fā)電有限公司,浙江 杭州 311322;2.浙江遠(yuǎn)算科技有限公司,浙江 杭州 310012)

拱壩穩(wěn)定性主要是靠兩岸壩肩的反力來維持,是一種經(jīng)濟(jì)性和安全性都非常好的壩型。在長期服役下,拱壩受復(fù)雜地質(zhì)條件和自身材料侵蝕老化等因素影響,存在一定的結(jié)構(gòu)安全問題,對(duì)國民經(jīng)濟(jì)和人民財(cái)產(chǎn)安全等具有重大影響。

隨著計(jì)算機(jī)技術(shù)的發(fā)展,有限元法成為大壩結(jié)構(gòu)安全狀態(tài)評(píng)估的重要手段之一[1-4],但是,該方法對(duì)大壩混凝土和基巖材料參數(shù)的準(zhǔn)確性要求較高[5-7]。通常的做法是通過參數(shù)反演來獲得較準(zhǔn)確的大壩材料參數(shù),當(dāng)前主流的參數(shù)反演方法包括Levenberg-Marquardt算法[8]、高斯-牛頓法[9]等傳統(tǒng)方法,以及神經(jīng)網(wǎng)絡(luò)算法、遺傳算法等智能反演方法[10]。傳統(tǒng)反演方法一般存在速度慢、精度差等缺陷,計(jì)算結(jié)果十分依賴于初值的選取;智能反演方法則存在神經(jīng)網(wǎng)絡(luò)全局搜索能力差,遺傳算法容易早熟收斂和結(jié)果依賴于初值選取等問題[11]。

本文旨在研究華光潭拱壩在長期服役過程中結(jié)構(gòu)的運(yùn)行狀態(tài),在已有監(jiān)測(cè)數(shù)據(jù)基礎(chǔ)上,對(duì)拱壩進(jìn)行基于數(shù)據(jù)同化參數(shù)反演的有限元計(jì)算。首先,建立分壩段的拱壩有限元模型;其次,在混凝土和基巖材料參數(shù)敏感性分析的基礎(chǔ)上,引入數(shù)據(jù)同化算法對(duì)其進(jìn)行分類反演;最后,根據(jù)反演結(jié)果校正有限元計(jì)算的材料參數(shù),并對(duì)有限元計(jì)算結(jié)果進(jìn)行殘差分析。本方法在確定觀測(cè)值和背景值后,通過Python腳本運(yùn)行,可實(shí)現(xiàn)自動(dòng)迭代優(yōu)化,在保證精度的前提下,大幅度提高參數(shù)反演的效率,對(duì)拱壩結(jié)構(gòu)安全評(píng)估具有重要意義。

1 工程背景介紹

華光潭一級(jí)水電站是一座以發(fā)電為主,兼有防洪作用的中型工程,水庫總庫容8 257萬m3。攔河壩為混凝土雙曲薄拱壩,壩頂高程449.85 m,最大壩高103.85 m,壩頂寬5.8 m,分13個(gè)壩段。該拱壩于2005年5月下閘蓄水,壩頂設(shè)有13個(gè)水平位移監(jiān)測(cè)點(diǎn),分別編號(hào)為TP1—13。華光潭拱壩在長期運(yùn)行中,由于多種原因,監(jiān)測(cè)設(shè)備發(fā)生故障,存在觀測(cè)數(shù)據(jù)不完整、混凝土和基巖材料力學(xué)性質(zhì)難以監(jiān)測(cè)等問題。

2 有限元分析

為了獲得華光潭拱壩的結(jié)構(gòu)狀態(tài),本文使用開源的結(jié)構(gòu)有限元分析軟件code_aster進(jìn)行拱壩有限元模型的建立與計(jì)算。拱壩有限元模型在上下游和左右岸方向范圍均為1.5倍壩高,豎直方向壩基為2倍壩高。選用3D六面體線性實(shí)體單元進(jìn)行結(jié)構(gòu)化網(wǎng)格劃分,共148 032個(gè)節(jié)點(diǎn),134 100個(gè)單元(圖1)。在壩體和壩基交界面位置設(shè)置基于內(nèi)聚力算法的連接單元,模擬混凝土和巖石材料中的微小孔隙和裂縫,避免出現(xiàn)較大的應(yīng)力集中問題(圖2)。有限元計(jì)算采用以下假定:①不考慮結(jié)構(gòu)配筋的影響,把結(jié)構(gòu)視為素混凝土;②混凝土的彈性模量假定為常值;③線性徐變理論[12]。有限元計(jì)算材料參數(shù)設(shè)計(jì)值見表1。

圖1 拱壩有限元計(jì)算網(wǎng)格

圖2 基于內(nèi)聚力算法的連接單元

表1 有限元計(jì)算材料參數(shù)設(shè)計(jì)值

2.1 拱壩封拱狀態(tài)仿真

拱壩建造初期的應(yīng)力狀態(tài)對(duì)運(yùn)行期的結(jié)構(gòu)計(jì)算有較大的影響,因此需要對(duì)封拱應(yīng)力狀態(tài)進(jìn)行計(jì)算。不同壩段在重力作用下發(fā)生變形,導(dǎo)致橫縫面產(chǎn)生復(fù)雜的應(yīng)力狀態(tài),通過奇數(shù)壩段和偶數(shù)壩段在重力作用下的變形疊加計(jì)算,可以較好地模擬拱壩的封拱建造過程,從而獲得封拱狀態(tài)的初始應(yīng)力場(chǎng)[13]。

2.2 拱壩運(yùn)行狀態(tài)仿真

拱壩運(yùn)行狀態(tài)仿真本質(zhì)上是熱力耦合計(jì)算,以封拱狀態(tài)仿真獲得的應(yīng)力場(chǎng)作為初始物理場(chǎng),將熱應(yīng)力場(chǎng)和靜水壓力等外載荷產(chǎn)生的應(yīng)力場(chǎng)進(jìn)行耦合計(jì)算,獲得最終的拱壩結(jié)構(gòu)狀態(tài)。拱壩運(yùn)行狀態(tài)仿真中施加的載荷有靜水壓力、泥沙淤積壓力、氣溫和水溫,拱壩左右岸、上下游和底面各邊界均設(shè)置為法向約束。

3 材料參數(shù)反演

由于拱壩長期運(yùn)行過程中,混凝土與基巖材料參數(shù)發(fā)生了變化,為了獲得準(zhǔn)確的材料參數(shù),就需要進(jìn)行材料參數(shù)反演。本文的數(shù)據(jù)同化參數(shù)反演是指從拱壩的位移監(jiān)測(cè)值和數(shù)據(jù)同化算法出發(fā),確定拱壩混凝土與基巖的材料參數(shù)。數(shù)據(jù)同化算法是一種考慮數(shù)據(jù)的時(shí)空分布,量化模型和觀測(cè)方法誤差,將觀測(cè)信息動(dòng)態(tài)融合到模型中的分析方法。三維變分法(3DVAR)是應(yīng)用最廣泛的數(shù)據(jù)同化算法之一,可滿足材料參數(shù)校準(zhǔn)的高精度要求,并且計(jì)算量相對(duì)較小,能夠較好地應(yīng)用于實(shí)際的工程[14-19],其代價(jià)函數(shù)的形式見式(1)。

(y-H[x])

(1)

式中xb——先驗(yàn)值;B——背景誤差協(xié)方差矩陣;y——觀測(cè)值;H——觀測(cè)算子;R——觀測(cè)誤差協(xié)方差矩陣。

3.1 材料參數(shù)敏感性分析

由于混凝土和基巖材料參數(shù)眾多,在參數(shù)反演之前,需要通過敏感性分析,評(píng)估各個(gè)材料參數(shù)對(duì)拱壩位移的影響大小,從而提高材料參數(shù)反演的效率。拱壩混凝土和基巖的彈性模量對(duì)其結(jié)構(gòu)狀態(tài)的影響是比較大的,因此,不對(duì)其進(jìn)行敏感性分析而直接作為背景值之一進(jìn)行參數(shù)反演。根據(jù)工程經(jīng)驗(yàn)初步篩選,擬對(duì)以下參數(shù)進(jìn)行敏感性分析:熱導(dǎo)率、熱膨脹系數(shù)、熱膨脹參考溫度、壩體泊松比、壩基泊松比和對(duì)流換熱系數(shù)(表2)。

表2 參數(shù)敏感性分析的變量范圍

拱壩位移監(jiān)測(cè)值中的水平徑向位移受溫度變化和水位調(diào)節(jié)影響較大,是拱壩安全監(jiān)測(cè)關(guān)注的重點(diǎn)。通過觀察壩頂水平徑向位移(DR)在拱壩橫河向位置坐標(biāo)(X)上的空間分布(圖3—8)可知,熱傳導(dǎo)系數(shù)、對(duì)流換熱系數(shù)和壩基泊松比對(duì)徑向位移的結(jié)果影響比較小,而熱膨脹系數(shù)、熱膨脹參考溫度和壩體泊松比對(duì)徑向位移結(jié)果的影響較大?？紤]拱壩左右岸基巖物性差異,將其彈性模量根據(jù)拱壩左右岸進(jìn)行區(qū)分。因此,最后確定數(shù)據(jù)同化參數(shù)反演的背景場(chǎng)有6個(gè)參數(shù):左岸彈性模量EL、右岸彈性模量ER、壩體彈性模量E、熱膨脹系數(shù)Alp、熱膨脹參考溫度Tref和壩體泊松比P。

圖3 熱導(dǎo)率對(duì)徑向位移的影響

圖4 熱膨脹系數(shù)對(duì)徑向位移的影響

圖5 熱膨脹參考溫度對(duì)徑向位移的影響

圖6 對(duì)流換熱系數(shù)對(duì)徑向位移的影響

圖7 壩基泊松比對(duì)徑向位移的影響

圖8 壩體泊松比對(duì)徑向位移的影響

3.2 數(shù)據(jù)同化參數(shù)反演

為了獲得更加準(zhǔn)確的材料參數(shù),將材料參數(shù)敏感性分析獲得的6個(gè)混凝土和基巖材料參數(shù)分為2類:①彈性結(jié)構(gòu)材料參數(shù)(EL、ER、E、P);②熱學(xué)材料參數(shù)(Alp、Tref),對(duì)2類材料參數(shù)分別進(jìn)行數(shù)據(jù)同化參數(shù)反演。根據(jù)不同載荷的影響,采用逐步回歸的統(tǒng)計(jì)學(xué)方法對(duì)水平徑向位移進(jìn)行分解見式(2)—(5)[20]:

δ=δH+δS+δT

(2)

(3)

δS=b1sinD+b2cosD+b3sinD·cosD+b4sin2D+b5cos2D+…

(4)

(5)

式中δ——水平徑向位移;δH——徑向位移水位分量;δS——徑向位移季節(jié)分量;δT——徑向位移時(shí)效分量;D=2π·t/365;t——天數(shù);H——水位;ai、bi、ci——常數(shù)(i=1,2,…n)。

數(shù)據(jù)同化材料參數(shù)反演根據(jù)參數(shù)類別分為彈性結(jié)構(gòu)材料參數(shù)反演和熱學(xué)材料參數(shù)反演2個(gè)步驟(圖9)。

注：W—水位;T—?dú)鉁睾退疁?其中下標(biāo)表示不同日期對(duì)應(yīng)的值。

a)彈性結(jié)構(gòu)材料參數(shù)反演?？紤]位移參考初值的影響,將2個(gè)不同日期下的壩頂水平徑向位移水位分量差值作為數(shù)據(jù)同化算法的觀測(cè)值。以4個(gè)彈性結(jié)構(gòu)材料參數(shù)作為背景場(chǎng),背景場(chǎng)初值為對(duì)應(yīng)參數(shù)的設(shè)計(jì)值。以僅靜水壓力和泥沙淤積壓力作用下的水位分量仿真作為代價(jià)函數(shù),進(jìn)行迭代計(jì)算,獲得一組最優(yōu)彈性結(jié)構(gòu)材料參數(shù):EL0、ER0、E0、P0。

b)熱學(xué)材料參數(shù)反演。將2個(gè)不同日期下的壩頂水平徑向位移季節(jié)分量差值作為數(shù)據(jù)同化算法的觀測(cè)值。以2個(gè)熱學(xué)材料參數(shù)作為背景場(chǎng),背景場(chǎng)初值為對(duì)應(yīng)參數(shù)的設(shè)計(jì)值。以僅在溫度作用下的季節(jié)分量仿真作為代價(jià)函數(shù),進(jìn)行迭代計(jì)算,季節(jié)分量仿真本質(zhì)上是結(jié)構(gòu)熱應(yīng)力計(jì)算,需要用到①中的4個(gè)彈性結(jié)構(gòu)材料參數(shù):EL0、ER0、E0、P0。最終獲得一組最優(yōu)熱學(xué)材料參數(shù):Alp0、Tref0。

上述步驟中的觀測(cè)值數(shù)據(jù),是基于2012年壩頂水平位移測(cè)點(diǎn)TP3、TP5、TP8、TP9、TP11的徑向位移監(jiān)測(cè)值,通過逐步回歸后進(jìn)一步處理得到的(表3)。在確定了背景場(chǎng)、觀測(cè)值和代價(jià)函數(shù)之后,通過Pyhon腳本運(yùn)行數(shù)據(jù)同化算法,最終獲得6個(gè)混凝土和基巖材料參數(shù):EL0、ER0、E0、P0、Alp0、Tref0。

表3 壩頂水平徑向位移各分量及數(shù)據(jù)同化參數(shù)反演觀測(cè)值

4 結(jié)果與討論

4.1 數(shù)據(jù)同化參數(shù)反演結(jié)果

在數(shù)據(jù)同化參數(shù)反演過程中,迭代次數(shù)達(dá)到50次后,各個(gè)參數(shù)值的變化基本趨于穩(wěn)定,在保證準(zhǔn)確性的同時(shí)具有較快的計(jì)算速度(圖10)。數(shù)據(jù)同化參數(shù)反演獲得的6個(gè)材料參數(shù)相較于設(shè)計(jì)值均有所變化(表4)。壩體彈性模量反演值比設(shè)計(jì)值有所提高,是因?yàn)殡S著時(shí)間的增加,混凝土硬化,剛度逐漸增加,同時(shí),混凝土材料設(shè)計(jì)參數(shù)與實(shí)際施工時(shí)的材料參數(shù)存在一定差異。而左右岸彈性模量減小是由于基巖風(fēng)化和隨時(shí)間軟化作用導(dǎo)致的。

圖10 數(shù)據(jù)同化參數(shù)反演迭代次數(shù)

表4 數(shù)據(jù)同化反演后的拱壩混凝土和基巖材料參數(shù)

數(shù)據(jù)同化參數(shù)反演后,對(duì)2012年6月15日工況下水位分量仿真和季節(jié)分量仿真的結(jié)果進(jìn)行分析。

a)水位分量仿真結(jié)果。在水位分量仿真中,由于拱壩是一個(gè)薄殼結(jié)構(gòu),在殼結(jié)構(gòu)中部位置剛度最低,此處向下游方向位移最大為22 mm(圖11)。拱冠梁截面壩踵位置呈現(xiàn)向上抬升的趨勢(shì),將其等效成底部固支的梁結(jié)構(gòu),最大位移發(fā)生于梁的自由端,該位移是由靜水壓力和泥沙淤積壓力主導(dǎo)的(圖12)。

注：帶網(wǎng)格云圖為拱壩變形后的位置。

注：帶網(wǎng)格云圖為拱壩變形前的位置。

b)季節(jié)分量仿真結(jié)果。在季節(jié)分量仿真中,上游水位線以下部分溫度約為6～10℃,沿著拱壩厚度方向,溫度基本呈線性分布(圖13)。由于拱壩上部和下部的溫差較大約為20℃,引起拱壩不均勻的熱膨脹,壩頂向上游方向位移最大為21 mm,外界溫度的影響起到了主導(dǎo)作用(圖14)。

圖13 拱冠梁截面溫度

4.2 位移殘差分析

在華光潭拱壩的位移監(jiān)測(cè)中,每個(gè)測(cè)點(diǎn)均有一個(gè)位移參考初值,后續(xù)實(shí)測(cè)數(shù)據(jù)均為相對(duì)于該位移參考初值的位移。由于位移初始參考值數(shù)據(jù)缺失,選取2013年1月15日拱壩在外載荷作用下產(chǎn)生的位移值作為位移參考初值,實(shí)現(xiàn)仿真結(jié)果與監(jiān)測(cè)數(shù)據(jù)的可對(duì)比性。使用拱壩混凝土和基巖材料參數(shù)設(shè)計(jì)值和反演值,分別進(jìn)行2013—2018年的拱壩運(yùn)行狀態(tài)仿真,對(duì)比壩頂水平位移測(cè)點(diǎn)TP8的徑向位移數(shù)據(jù),從曲線圖可知,材料參數(shù)反演后仿真結(jié)果更接近于監(jiān)測(cè)值(圖15)。

注：帶網(wǎng)格云圖為拱壩變形后的位置。

a)參數(shù)反演前

計(jì)算同一日期下TP8測(cè)點(diǎn)水平徑向位移監(jiān)測(cè)值和仿真值的差值,即殘差。反演前后殘差呈對(duì)稱分布,但是反演前的殘差均值約為-12,總體標(biāo)準(zhǔn)差為5.81,而反演后的殘差均值約為0,并符合正態(tài)分布,總體標(biāo)準(zhǔn)差為1.09(圖16)。經(jīng)過數(shù)據(jù)同化參數(shù)反演后殘差標(biāo)準(zhǔn)差下降了81.2%,極大地提高了計(jì)算結(jié)果的準(zhǔn)確性。

a)反演前殘差分布

4.3 拱壩應(yīng)力分析

選取2018年6月15日拱壩的有限元仿真結(jié)果進(jìn)行應(yīng)力分析,上游壩面在中部區(qū)域最大壓應(yīng)力為3.1 MPa,壩肩區(qū)域產(chǎn)生較大的拉應(yīng)力,最大為0.91 MPa(圖17)。下游壩面的中部區(qū)域主要是拉應(yīng)力,最大值為0.36 MPa,壩肩區(qū)域主要是壓應(yīng)力,最大值為7 MPa(圖18)。壩踵主要承受拉應(yīng)力,最大為0.69 MPa,是相對(duì)危險(xiǎn)區(qū)域,壩趾最大壓應(yīng)力3 MPa,拱壩的總體應(yīng)力水平均在合理工程范圍內(nèi)。

圖17 上游壩面第一主應(yīng)力

圖18 下游壩面第一主應(yīng)力

5 結(jié)論

為了準(zhǔn)確地研究華光潭拱壩結(jié)構(gòu)的運(yùn)行狀態(tài),本文采用基于數(shù)據(jù)同化參數(shù)反演的有限元結(jié)構(gòu)仿真方法,對(duì)華光潭拱壩結(jié)構(gòu)位移和應(yīng)力狀態(tài)進(jìn)行了分析,得到如下結(jié)論。

a)使用數(shù)據(jù)同化參數(shù)反演獲得的混凝土和基巖材料參數(shù),能夠較大地提高有限元仿真的準(zhǔn)確性,為華光潭拱壩結(jié)構(gòu)的安全評(píng)估提供一個(gè)重要依據(jù)。

b)基于數(shù)據(jù)同化參數(shù)反演的華光潭拱壩有限元仿真,可以根據(jù)實(shí)時(shí)采集的水位、外部溫度等數(shù)據(jù),同步計(jì)算獲得整個(gè)拱壩位移和應(yīng)力等物理量在時(shí)空上的分布,在一定程度上彌補(bǔ)監(jiān)測(cè)點(diǎn)損壞和數(shù)據(jù)采集頻率低的不足,為華光潭拱壩管理人員提供實(shí)時(shí)且全面的結(jié)構(gòu)狀態(tài)信息。

c)拱壩在外載荷綜合作用下發(fā)生結(jié)構(gòu)變形,具有一定規(guī)律性:拱冠梁中心附近區(qū)域向下游方向位移較大,以水位的影響為主導(dǎo),向上游方向位移較大,以外界溫度導(dǎo)致的熱膨脹為主導(dǎo);拱壩最大拉壓應(yīng)力產(chǎn)生位置多為壩踵和壩趾位置,并且拱壩向下游變形時(shí),壩踵多為拉應(yīng)力,壩趾多為壓應(yīng)力。