林咸志,羅增浤,許 正,閔皆昇,曹傳波,趙 權(quán)
(1.浙江浙能華光潭水力發(fā)電有限公司,浙江 杭州 311322;2.浙江遠算科技有限公司,浙江 杭州 310012)
拱壩穩(wěn)定性主要是靠兩岸壩肩的反力來維持,是一種經(jīng)濟性和安全性都非常好的壩型。在長期服役下,拱壩受復(fù)雜地質(zhì)條件和自身材料侵蝕老化等因素影響,存在一定的結(jié)構(gòu)安全問題,對國民經(jīng)濟和人民財產(chǎn)安全等具有重大影響。
隨著計算機技術(shù)的發(fā)展,有限元法成為大壩結(jié)構(gòu)安全狀態(tài)評估的重要手段之一[1-4],但是,該方法對大壩混凝土和基巖材料參數(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)反演方法一般存在速度慢、精度差等缺陷,計算結(jié)果十分依賴于初值的選取;智能反演方法則存在神經(jīng)網(wǎng)絡(luò)全局搜索能力差,遺傳算法容易早熟收斂和結(jié)果依賴于初值選取等問題[11]。
本文旨在研究華光潭拱壩在長期服役過程中結(jié)構(gòu)的運行狀態(tài),在已有監(jiān)測數(shù)據(jù)基礎(chǔ)上,對拱壩進行基于數(shù)據(jù)同化參數(shù)反演的有限元計算。首先,建立分壩段的拱壩有限元模型;其次,在混凝土和基巖材料參數(shù)敏感性分析的基礎(chǔ)上,引入數(shù)據(jù)同化算法對其進行分類反演;最后,根據(jù)反演結(jié)果校正有限元計算的材料參數(shù),并對有限元計算結(jié)果進行殘差分析。本方法在確定觀測值和背景值后,通過Python腳本運行,可實現(xiàn)自動迭代優(yōu)化,在保證精度的前提下,大幅度提高參數(shù)反演的效率,對拱壩結(jié)構(gòu)安全評估具有重要意義。
華光潭一級水電站是一座以發(fā)電為主,兼有防洪作用的中型工程,水庫總庫容8 257萬m3。攔河壩為混凝土雙曲薄拱壩,壩頂高程449.85 m,最大壩高103.85 m,壩頂寬5.8 m,分13個壩段。該拱壩于2005年5月下閘蓄水,壩頂設(shè)有13個水平位移監(jiān)測點,分別編號為TP1—13。華光潭拱壩在長期運行中,由于多種原因,監(jiān)測設(shè)備發(fā)生故障,存在觀測數(shù)據(jù)不完整、混凝土和基巖材料力學(xué)性質(zhì)難以監(jiān)測等問題。
為了獲得華光潭拱壩的結(jié)構(gòu)狀態(tài),本文使用開源的結(jié)構(gòu)有限元分析軟件code_aster進行拱壩有限元模型的建立與計算。拱壩有限元模型在上下游和左右岸方向范圍均為1.5倍壩高,豎直方向壩基為2倍壩高。選用3D六面體線性實體單元進行結(jié)構(gòu)化網(wǎng)格劃分,共148 032個節(jié)點,134 100個單元(圖1)。在壩體和壩基交界面位置設(shè)置基于內(nèi)聚力算法的連接單元,模擬混凝土和巖石材料中的微小孔隙和裂縫,避免出現(xiàn)較大的應(yīng)力集中問題(圖2)。有限元計算采用以下假定:①不考慮結(jié)構(gòu)配筋的影響,把結(jié)構(gòu)視為素混凝土;②混凝土的彈性模量假定為常值;③線性徐變理論[12]。有限元計算材料參數(shù)設(shè)計值見表1。
圖1 拱壩有限元計算網(wǎng)格
圖2 基于內(nèi)聚力算法的連接單元
表1 有限元計算材料參數(shù)設(shè)計值
拱壩建造初期的應(yīng)力狀態(tài)對運行期的結(jié)構(gòu)計算有較大的影響,因此需要對封拱應(yīng)力狀態(tài)進行計算。不同壩段在重力作用下發(fā)生變形,導(dǎo)致橫縫面產(chǎn)生復(fù)雜的應(yīng)力狀態(tài),通過奇數(shù)壩段和偶數(shù)壩段在重力作用下的變形疊加計算,可以較好地模擬拱壩的封拱建造過程,從而獲得封拱狀態(tài)的初始應(yīng)力場[13]。
拱壩運行狀態(tài)仿真本質(zhì)上是熱力耦合計算,以封拱狀態(tài)仿真獲得的應(yīng)力場作為初始物理場,將熱應(yīng)力場和靜水壓力等外載荷產(chǎn)生的應(yīng)力場進行耦合計算,獲得最終的拱壩結(jié)構(gòu)狀態(tài)。拱壩運行狀態(tài)仿真中施加的載荷有靜水壓力、泥沙淤積壓力、氣溫和水溫,拱壩左右岸、上下游和底面各邊界均設(shè)置為法向約束。
由于拱壩長期運行過程中,混凝土與基巖材料參數(shù)發(fā)生了變化,為了獲得準(zhǔn)確的材料參數(shù),就需要進行材料參數(shù)反演。本文的數(shù)據(jù)同化參數(shù)反演是指從拱壩的位移監(jiān)測值和數(shù)據(jù)同化算法出發(fā),確定拱壩混凝土與基巖的材料參數(shù)。數(shù)據(jù)同化算法是一種考慮數(shù)據(jù)的時空分布,量化模型和觀測方法誤差,將觀測信息動態(tài)融合到模型中的分析方法。三維變分法(3DVAR)是應(yīng)用最廣泛的數(shù)據(jù)同化算法之一,可滿足材料參數(shù)校準(zhǔn)的高精度要求,并且計算量相對較小,能夠較好地應(yīng)用于實際的工程[14-19],其代價函數(shù)的形式見式(1)。
(y-H[x])
(1)
式中xb——先驗值;B——背景誤差協(xié)方差矩陣;y——觀測值;H——觀測算子;R——觀測誤差協(xié)方差矩陣。
由于混凝土和基巖材料參數(shù)眾多,在參數(shù)反演之前,需要通過敏感性分析,評估各個材料參數(shù)對拱壩位移的影響大小,從而提高材料參數(shù)反演的效率。拱壩混凝土和基巖的彈性模量對其結(jié)構(gòu)狀態(tài)的影響是比較大的,因此,不對其進行敏感性分析而直接作為背景值之一進行參數(shù)反演。根據(jù)工程經(jīng)驗初步篩選,擬對以下參數(shù)進行敏感性分析:熱導(dǎo)率、熱膨脹系數(shù)、熱膨脹參考溫度、壩體泊松比、壩基泊松比和對流換熱系數(shù)(表2)。
表2 參數(shù)敏感性分析的變量范圍
拱壩位移監(jiān)測值中的水平徑向位移受溫度變化和水位調(diào)節(jié)影響較大,是拱壩安全監(jiān)測關(guān)注的重點。通過觀察壩頂水平徑向位移(DR)在拱壩橫河向位置坐標(biāo)(X)上的空間分布(圖3—8)可知,熱傳導(dǎo)系數(shù)、對流換熱系數(shù)和壩基泊松比對徑向位移的結(jié)果影響比較小,而熱膨脹系數(shù)、熱膨脹參考溫度和壩體泊松比對徑向位移結(jié)果的影響較大??紤]拱壩左右岸基巖物性差異,將其彈性模量根據(jù)拱壩左右岸進行區(qū)分。因此,最后確定數(shù)據(jù)同化參數(shù)反演的背景場有6個參數(shù):左岸彈性模量EL、右岸彈性模量ER、壩體彈性模量E、熱膨脹系數(shù)Alp、熱膨脹參考溫度Tref和壩體泊松比P。
圖3 熱導(dǎo)率對徑向位移的影響
圖4 熱膨脹系數(shù)對徑向位移的影響
圖5 熱膨脹參考溫度對徑向位移的影響
圖6 對流換熱系數(shù)對徑向位移的影響
圖7 壩基泊松比對徑向位移的影響
圖8 壩體泊松比對徑向位移的影響
為了獲得更加準(zhǔn)確的材料參數(shù),將材料參數(shù)敏感性分析獲得的6個混凝土和基巖材料參數(shù)分為2類:①彈性結(jié)構(gòu)材料參數(shù)(EL、ER、E、P);②熱學(xué)材料參數(shù)(Alp、Tref),對2類材料參數(shù)分別進行數(shù)據(jù)同化參數(shù)反演。根據(jù)不同載荷的影響,采用逐步回歸的統(tǒng)計學(xué)方法對水平徑向位移進行分解見式(2)—(5)[20]:
δ=δH+δS+δT
(2)
(3)
δS=b1sinD+b2cosD+b3sinD·cosD+b4sin2D+b5cos2D+…
(4)
(5)
式中δ——水平徑向位移;δH——徑向位移水位分量;δS——徑向位移季節(jié)分量;δT——徑向位移時效分量;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個步驟(圖9)。
注:W—水位;T—氣溫和水溫,其中下標(biāo)表示不同日期對應(yīng)的值。
a)彈性結(jié)構(gòu)材料參數(shù)反演。考慮位移參考初值的影響,將2個不同日期下的壩頂水平徑向位移水位分量差值作為數(shù)據(jù)同化算法的觀測值。以4個彈性結(jié)構(gòu)材料參數(shù)作為背景場,背景場初值為對應(yīng)參數(shù)的設(shè)計值。以僅靜水壓力和泥沙淤積壓力作用下的水位分量仿真作為代價函數(shù),進行迭代計算,獲得一組最優(yōu)彈性結(jié)構(gòu)材料參數(shù):EL0、ER0、E0、P0。
b)熱學(xué)材料參數(shù)反演。將2個不同日期下的壩頂水平徑向位移季節(jié)分量差值作為數(shù)據(jù)同化算法的觀測值。以2個熱學(xué)材料參數(shù)作為背景場,背景場初值為對應(yīng)參數(shù)的設(shè)計值。以僅在溫度作用下的季節(jié)分量仿真作為代價函數(shù),進行迭代計算,季節(jié)分量仿真本質(zhì)上是結(jié)構(gòu)熱應(yīng)力計算,需要用到①中的4個彈性結(jié)構(gòu)材料參數(shù):EL0、ER0、E0、P0。最終獲得一組最優(yōu)熱學(xué)材料參數(shù):Alp0、Tref0。
上述步驟中的觀測值數(shù)據(jù),是基于2012年壩頂水平位移測點TP3、TP5、TP8、TP9、TP11的徑向位移監(jiān)測值,通過逐步回歸后進一步處理得到的(表3)。在確定了背景場、觀測值和代價函數(shù)之后,通過Pyhon腳本運行數(shù)據(jù)同化算法,最終獲得6個混凝土和基巖材料參數(shù):EL0、ER0、E0、P0、Alp0、Tref0。
表3 壩頂水平徑向位移各分量及數(shù)據(jù)同化參數(shù)反演觀測值
在數(shù)據(jù)同化參數(shù)反演過程中,迭代次數(shù)達到50次后,各個參數(shù)值的變化基本趨于穩(wěn)定,在保證準(zhǔn)確性的同時具有較快的計算速度(圖10)。數(shù)據(jù)同化參數(shù)反演獲得的6個材料參數(shù)相較于設(shè)計值均有所變化(表4)。壩體彈性模量反演值比設(shè)計值有所提高,是因為隨著時間的增加,混凝土硬化,剛度逐漸增加,同時,混凝土材料設(shè)計參數(shù)與實際施工時的材料參數(shù)存在一定差異。而左右岸彈性模量減小是由于基巖風(fēng)化和隨時間軟化作用導(dǎo)致的。
圖10 數(shù)據(jù)同化參數(shù)反演迭代次數(shù)
表4 數(shù)據(jù)同化反演后的拱壩混凝土和基巖材料參數(shù)
數(shù)據(jù)同化參數(shù)反演后,對2012年6月15日工況下水位分量仿真和季節(jié)分量仿真的結(jié)果進行分析。
a)水位分量仿真結(jié)果。在水位分量仿真中,由于拱壩是一個薄殼結(jié)構(gòu),在殼結(jié)構(gòu)中部位置剛度最低,此處向下游方向位移最大為22 mm(圖11)。拱冠梁截面壩踵位置呈現(xiàn)向上抬升的趨勢,將其等效成底部固支的梁結(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 拱冠梁截面溫度
在華光潭拱壩的位移監(jiān)測中,每個測點均有一個位移參考初值,后續(xù)實測數(shù)據(jù)均為相對于該位移參考初值的位移。由于位移初始參考值數(shù)據(jù)缺失,選取2013年1月15日拱壩在外載荷作用下產(chǎn)生的位移值作為位移參考初值,實現(xiàn)仿真結(jié)果與監(jiān)測數(shù)據(jù)的可對比性。使用拱壩混凝土和基巖材料參數(shù)設(shè)計值和反演值,分別進行2013—2018年的拱壩運行狀態(tài)仿真,對比壩頂水平位移測點TP8的徑向位移數(shù)據(jù),從曲線圖可知,材料參數(shù)反演后仿真結(jié)果更接近于監(jiān)測值(圖15)。
注:帶網(wǎng)格云圖為拱壩變形后的位置。
a)參數(shù)反演前
計算同一日期下TP8測點水平徑向位移監(jiān)測值和仿真值的差值,即殘差。反演前后殘差呈對稱分布,但是反演前的殘差均值約為-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%,極大地提高了計算結(jié)果的準(zhǔn)確性。
a)反演前殘差分布
選取2018年6月15日拱壩的有限元仿真結(jié)果進行應(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,是相對危險區(qū)域,壩趾最大壓應(yīng)力3 MPa,拱壩的總體應(yīng)力水平均在合理工程范圍內(nèi)。
圖17 上游壩面第一主應(yīng)力
圖18 下游壩面第一主應(yīng)力
為了準(zhǔn)確地研究華光潭拱壩結(jié)構(gòu)的運行狀態(tài),本文采用基于數(shù)據(jù)同化參數(shù)反演的有限元結(jié)構(gòu)仿真方法,對華光潭拱壩結(jié)構(gòu)位移和應(yīng)力狀態(tài)進行了分析,得到如下結(jié)論。
a)使用數(shù)據(jù)同化參數(shù)反演獲得的混凝土和基巖材料參數(shù),能夠較大地提高有限元仿真的準(zhǔn)確性,為華光潭拱壩結(jié)構(gòu)的安全評估提供一個重要依據(jù)。
b)基于數(shù)據(jù)同化參數(shù)反演的華光潭拱壩有限元仿真,可以根據(jù)實時采集的水位、外部溫度等數(shù)據(jù),同步計算獲得整個拱壩位移和應(yīng)力等物理量在時空上的分布,在一定程度上彌補監(jiān)測點損壞和數(shù)據(jù)采集頻率低的不足,為華光潭拱壩管理人員提供實時且全面的結(jié)構(gòu)狀態(tài)信息。
c)拱壩在外載荷綜合作用下發(fā)生結(jié)構(gòu)變形,具有一定規(guī)律性:拱冠梁中心附近區(qū)域向下游方向位移較大,以水位的影響為主導(dǎo),向上游方向位移較大,以外界溫度導(dǎo)致的熱膨脹為主導(dǎo);拱壩最大拉壓應(yīng)力產(chǎn)生位置多為壩踵和壩趾位置,并且拱壩向下游變形時,壩踵多為拉應(yīng)力,壩趾多為壓應(yīng)力。