基于數據同化參數反演的華光潭拱壩結構分析

林咸志,羅增浤,許 正,閔皆昇,曹傳波,趙 權

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

拱壩穩定性主要是靠兩岸壩肩的反力來維持,是一種經濟性和安全性都非常好的壩型。在長期服役下,拱壩受復雜地質條件和自身材料侵蝕老化等因素影響,存在一定的結構安全問題,對國民經濟和人民財產安全等具有重大影響。

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

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

1 工程背景介紹

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

2 有限元分析

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

圖1 拱壩有限元計算網格

圖2 基于內聚力算法的連接單元

表1 有限元計算材料參數設計值

2.1 拱壩封拱狀態仿真

拱壩建造初期的應力狀態對運行期的結構計算有較大的影響,因此需要對封拱應力狀態進行計算。不同壩段在重力作用下發生變形,導致橫縫面產生復雜的應力狀態,通過奇數壩段和偶數壩段在重力作用下的變形疊加計算,可以較好地模擬拱壩的封拱建造過程,從而獲得封拱狀態的初始應力場[13]。

2.2 拱壩運行狀態仿真

拱壩運行狀態仿真本質上是熱力耦合計算,以封拱狀態仿真獲得的應力場作為初始物理場,將熱應力場和靜水壓力等外載荷產生的應力場進行耦合計算,獲得最終的拱壩結構狀態。拱壩運行狀態仿真中施加的載荷有靜水壓力、泥沙淤積壓力、氣溫和水溫,拱壩左右岸、上下游和底面各邊界均設置為法向約束。

3 材料參數反演

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

(y-H[x])

(1)

式中xb——先驗值;B——背景誤差協方差矩陣;y——觀測值;H——觀測算子;R——觀測誤差協方差矩陣。

3.1 材料參數敏感性分析

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

表2 參數敏感性分析的變量范圍

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

圖3 熱導率對徑向位移的影響

圖4 熱膨脹系數對徑向位移的影響

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

圖6 對流換熱系數對徑向位移的影響

圖7 壩基泊松比對徑向位移的影響

圖8 壩體泊松比對徑向位移的影響

3.2 數據同化參數反演

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

δ=δH+δS+δT

(2)

(3)

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

(4)

(5)

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

數據同化材料參數反演根據參數類別分為彈性結構材料參數反演和熱學材料參數反演2個步驟(圖9)。

注：W—水位;T—氣溫和水溫,其中下標表示不同日期對應的值。

a)彈性結構材料參數反演。考慮位移參考初值的影響,將2個不同日期下的壩頂水平徑向位移水位分量差值作為數據同化算法的觀測值。以4個彈性結構材料參數作為背景場,背景場初值為對應參數的設計值。以僅靜水壓力和泥沙淤積壓力作用下的水位分量仿真作為代價函數,進行迭代計算,獲得一組最優彈性結構材料參數:EL0、ER0、E0、P0。

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

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

表3 壩頂水平徑向位移各分量及數據同化參數反演觀測值

4 結果與討論

4.1 數據同化參數反演結果

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

圖10 數據同化參數反演迭代次數

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

數據同化參數反演后,對2012年6月15日工況下水位分量仿真和季節分量仿真的結果進行分析。

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

注：帶網格云圖為拱壩變形后的位置。

注：帶網格云圖為拱壩變形前的位置。

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

圖13 拱冠梁截面溫度

4.2 位移殘差分析

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

注：帶網格云圖為拱壩變形后的位置。

a)參數反演前

計算同一日期下TP8測點水平徑向位移監測值和仿真值的差值,即殘差。反演前后殘差呈對稱分布,但是反演前的殘差均值約為-12,總體標準差為5.81,而反演后的殘差均值約為0,并符合正態分布,總體標準差為1.09(圖16)。經過數據同化參數反演后殘差標準差下降了81.2%,極大地提高了計算結果的準確性。

a)反演前殘差分布

4.3 拱壩應力分析

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

圖17 上游壩面第一主應力

圖18 下游壩面第一主應力

5 結論

為了準確地研究華光潭拱壩結構的運行狀態,本文采用基于數據同化參數反演的有限元結構仿真方法,對華光潭拱壩結構位移和應力狀態進行了分析,得到如下結論。

a)使用數據同化參數反演獲得的混凝土和基巖材料參數,能夠較大地提高有限元仿真的準確性,為華光潭拱壩結構的安全評估提供一個重要依據。

b)基于數據同化參數反演的華光潭拱壩有限元仿真,可以根據實時采集的水位、外部溫度等數據,同步計算獲得整個拱壩位移和應力等物理量在時空上的分布,在一定程度上彌補監測點損壞和數據采集頻率低的不足,為華光潭拱壩管理人員提供實時且全面的結構狀態信息。

c)拱壩在外載荷綜合作用下發生結構變形,具有一定規律性:拱冠梁中心附近區域向下游方向位移較大,以水位的影響為主導,向上游方向位移較大,以外界溫度導致的熱膨脹為主導;拱壩最大拉壓應力產生位置多為壩踵和壩趾位置,并且拱壩向下游變形時,壩踵多為拉應力,壩趾多為壓應力。