李明軍,王均星
(武漢大學 水資源與水電工程科學國家重點實驗室,湖北 武漢 430072)
大壩安全監測監控模型主要應用的有3種:統計模型、確定性模型、混合模型。確定性模型相對于其他兩種模型有更明確的物理概念,可以更好的聯系大壩和壩基的結構性態,具有更好的預測效果[1]。溫度位移分量的確定一直是建立位移確定性監控模型的難點。混凝土拱壩的整體性較好,大壩承載能力較強但壩體結構對溫度變化十分敏感,溫度邊界條件復雜且溫度位移分量在總位移中占有比重較大,在建立混凝土拱壩的位移確定性監控模型時一般要求有能全面反映壩體溫度變化的溫度實測資料。而目前很多大壩在溫度監測方面只有環境溫度資料,尤其是對于長期運行以后(一般50年及以上)的大壩,由于建設時壩體內溫度計埋設不合理,或是運行后溫度計失效,無法獲得完整的壩體溫度實測資料,限制了確定性模型的應用。
目前,國內外對大壩位移確定性模型的研究主要適用于在壩體內部溫度場確定的情況,當缺乏壩體溫度實測資料時,利用邊界溫度推求溫度位移分量值是建立大壩位移確定性模型的唯一方法,非此就不能建立確定性模型[2,10]。但是由于拱壩溫度初始條件和邊界條件復雜,其溫度荷載很難準確模擬,為確定溫度位移分量帶來了極大的困難。李步娟[3]提出用“單位邊界溫度模型”計算測點的位移溫度分量,但由于測點位移變化滯后于邊界溫度變化,需要深入分析不同測點的滯后規律,其應用具有明顯局限性。李端有[4]采用朱伯芳法以氣溫估算庫水溫度,再通過氣溫和水溫確定拱壩邊界條件,以單個的壩體測點實測溫度值檢驗大壩計算溫度場,在此基礎上使用有限元數值分析法建立拱壩位移確定性模型,然而此種方法要求在壩體內布置有效溫度計以獲取壩體實測溫度資料,在壩體內未埋設溫度計或者有效溫度計數量極少且測點分布不均勻的情況下并不適用。
對于運行多年的混凝土拱壩,其運行期任意時刻位移溫度分量由壩體內部溫度場決定。壩體內部溫度場受壩體外部溫度影響且滯后于壩體外部溫度。壩體外部的環境溫度是由氣溫、大壩上游水庫水溫所組成的[5][12]。當利用邊界溫度推求溫度分量值時,不可避免的提出一些假定,這些假定的提出讓溫度場計算變得簡單,卻使得由邊界溫度推求的拱壩溫度場與實際溫度場存在一定的差異,當缺乏實際壩溫監測資料來對拱壩計算溫度場進行修正時,如何校驗拱壩溫度場,提高確定性模型的精度成為一個有待解決的問題。而對運行多年的混凝土拱壩實測位移分析表明,拱壩位移的變化(尤其是徑向位移)受環境溫度變化的影響非常明顯,位移隨溫度做一定的周期性變化。在正常情況下,對于長期運行的混凝土拱壩而言,壩體自重保持不變,位移的改變主要是由于水位、溫度以及時效變形引起的,而運行期時效變形趨于穩定。通常情況下,在一年中,時效位移小于1%的總位移(例如,一個130 m的高度拱壩,熱位移振幅約為20 mm,而其不可逆的趨勢是每年約0.1 mm)[6]。因此可以認為在任意時段Δt內,當Δt比較小時(Δτ<1 a)運行期混凝土拱壩位移的改變主要是由于水位和溫度變化引起的,即在 時段內,由溫度和水位變化所引起的位移變化值與該時間段內實測位移變化值應該保持一致,而水位變化所引起的位移改變是確定的。以此對拱壩實際溫度場進行檢驗。
針對運行期拱壩這一變形特點,提出以環境溫度(主要是氣溫與水溫)建立大壩溫度邊界條件,采用有限元法分析大壩實時溫度場,選取短時間實測位移序列校驗大壩溫度場,建立運行期拱壩位移確定性模型。這為我國大部分缺乏壩體溫度監測資料的運行期混凝土拱壩的安全監控提供了一個新思路。
在區域Ω內,準穩定溫度場T(x,y,z,τ)應滿足熱傳導方程:
(1)
初始條件:
T(x,y,z,τ0)=T(x,y,z)
(2)
設區域Ω的邊界為Γ,則拱壩溫度場的邊界條件[7]:
(1)第一類邊界條件:
(3)
(2)第三類邊界條件:
T(x,y,z,τ)|(x,y,z),Γ1=g(x,y,z,τ)
(4)
式中:α為導溫系數,m2/h;T0(x,y,z)為初始溫度,℃;g(x,y,z,τ)為邊界溫度,℃。
根據變分原理,任意時刻τ溫度場為:

(5)
當溫度T在τ=τ0時,T=T0,邊界Γ1上的邊界溫度為g(x,y,z,τ)。當式(5)所示泛函I(T)取極小值,溫度場T在區域Ω內滿足熱傳導方程(1),且在邊界Γ2上滿足公式(3),準穩定溫度場T則確定了。
將區域Ω用有限元離散,則泛函I(T)變成在各單位內的積分,即:
I=∑Ie
(6)
其中,Ie為單元e內的積分,即:
(6)
式中:ΔΩ為單元e所包含的子域。
對式(7)在積分號內求導,則有:
(8)
泛函取極值,有:
(9)
設單元e的結點溫度為:
{T(τ)}e=[Ti(τ),Tj(τ),Tm(τ),LL]T
(10)
則單元內任一點的溫度:
T(x,y,z,τ)=[Ni(x,y,z),Nj(x,y,z),
Nm(xy,z),LL]{T(τ)}e=[N]{T(τ)}e
(11)
將公式式(11)代入公式(8),通過式(9)可得到方程組:
[K]{T}+[C]{T}=0
(12)
式中:[K]為熱傳導矩陣;[C]為熱傳導矩陣;[T]為結點溫度向量;{T}為結點溫度向量隨時間的變化率。
求解方程組(12),就可以得到所有結點的溫度值。
當缺乏詳細的溫度資料來描述大壩實時溫度場時,一般采用氣溫來估算壩前水溫、擬定壩體溫度場的邊界條件,再通過有限元法得到大壩實時溫度場,進而建立大壩位移確定性模型。
采用氣溫估算壩前水溫的常用方法有華東勘院法、朱伯芳法、數學模型法等。目前混凝土拱壩主要采用朱伯芳法,其計算公式簡單,通過月平均溫度和壩前水深,綜合水庫地理位置,估算壩前垂直方向的水溫分布,具有良好的工程應用效果[7,11,14]。具體計算方法如下:
壩前任意水深處y處的水溫變化為
T(y,τ)=Tm(y)+A(y)cosω(τ-τ0-ε)
(13)
式中:Tm(y)為y處的年平均溫度,Tm(y)=c+(b-c)e-ay;A(y)為y處的溫度年變幅,A(y)=A0e-βy;ω=2τ/p,p為溫度變化周期,這里p=12月;τ0=6.5月(因τ=τ0時氣溫最高);ε為水溫與氣溫的變化相位差,ε=e-fe-γ;α、β、γ、d、f等都為相關的參數。
朱伯芳法假設壩前水位是正常蓄水位,其計算公式對水庫壩前水深小于100 m的工程應用效果較好。當水位變幅較大且壩前水深超過100 m時,采用朱伯芳法計算會產生較大誤差。當運行水位高于或等于正常蓄水位時,朱伯芳法的估算精度能滿足工程需求,否則需要對朱伯芳庫水溫度估算公式進行改進[7]。
由于庫水溫度與水深 有關,當壩前水位為實際運行水位時的水溫用 表示,相應的氣溫用 表示, 為運行水位,當壩前正常水深H>100 m時, 用分段函數表示:
(14)
改進的朱伯芳法估算的庫水溫度能更好地符合實際溫度邊界條件,使有限元計算溫度場能更好地符合拱壩實際溫度場,為提高位移確定性模型精度提供更好的基礎。
按照前文所述方法所擬定的混凝土拱壩溫度邊界條件與實際情況并不相同,由此所得到的混凝土拱壩溫度場不可避免與實際溫度場存在一定的差異。不同混凝土拱壩對于這種差異的響應并不相同,主要表現為位移監控模型中的精度。
按照前文所述方法所擬定的混凝土拱壩溫度邊界條件與實際情況并不相同,由此所得到的混凝土拱壩溫度場不可避免與實際溫度場存在一定的差異。不同混凝土拱壩對于這種差異的響應并不相同,主要表現為位移監控模型中的精度差異[8,9]。為了驗證混凝土拱壩溫度邊界條件的合理性,本文針對缺乏壩體溫度實測資料的運行期(一般運行10年及以上)混凝土拱壩的變形特點,提出一種校驗混凝土拱壩擬定的溫度邊界條件合理性的新方法。
在正常運行條件下,對于長期運行的混凝土拱壩而言,壩體自重保持不變,位移的改變主要是由于水位、溫度以及時效變形引起的,而運行期時效變形趨于穩定。本文采用有限元法計算任意時間段Δt(Δt≤1 a)內拱壩的變化溫度場,由熱-結構的耦合得到在 時段內溫度變化所引起的位移變化值,由于水位變化引起的位移是確定的,其與該時段內水位變化所引起的位移變化值疊加與同時段內實測位移變化值應該保持一致,即:
(15)

若能滿足此條件,則證明溫度邊界條件擬定的合理,若不能滿足此條件則根據實際資料對混凝土拱壩溫度邊界條件進行修正。
大壩位移確定性模型通過物理理論計算成果來構造環境自變量與大壩變形聯系的確定性關系式,其建立的基本思路為假設壩體和地基材料的力學參數與時效的數學模型,采用有限元法計算大壩觀測控制點在荷載作用下的壩體和壩基的位移場,通過此計算值與實測值進行最小二乘法擬合,求得大壩變形影響因子的調整參數以及時效位移分量參數,從而建立大壩位移確定性模型[13,16]。影響混凝土拱壩變形主要有三個因素:水位,溫度以及時效,其相關表達式為:
(16)
(17)
δθ=C1θ+C2lnθ
(18)
δ=XδH+YδT+δθ
(19)
式中:δ為測點總位移;δH、δT及δθ分別為水壓、溫度及時效位移分量;H為壩前水深;T為一定時間段內的氣溫平均值;θ為時間;C1、C2為待定系數;X、Y為調整系數,其中系數在總位移模型時根據實測資料利用回歸方程求出。
位移確定性監控模型的難點是有效的確定溫度邊界條件,得到與實際情況符合的拱壩實時溫度場,從而提高確定性模型的精度。在利用1中所論述的方法,擬定拱壩的溫度邊界條件,計算拱壩初始溫度場,選取時間段Δt(Δt<1 a),計算此時間段內的溫度位移變化值,將其與 時段內實測位移變化值和水位變化引起的位移變化值之差進行比較,當誤差值小于10%,則認為此次溫度場邊界條件擬定合理,否則根據實測資料對溫度邊界條件進行修正,重新計算直到滿足條件。
東江大壩為混凝土雙曲拱壩,其左右岸基本對稱,壩頂高程294.0 m,壩底高程137.0 m,壩高157.0 m,設計正常蓄水位285 m。東江拱壩于1986年開始蓄水,已安全運行二十余年,荷載基本于左右對稱,其中11#,為中間壩段。由于拱冠15號壩段所在測點監測資料部分丟失,因此本文擬將靠近拱冠梁位置的11#壩段的前方交會的8號測點作為研究對象,認為此位置測點在外部荷載作用下能反映整個拱壩的工作性態。
依據東江大壩的壩型及基礎地質資料,選擇壩體和一定范圍內的壩基為研究對象,確定壩體和壩基的彈性模量、線膨脹系數、熱傳導系數等物理力學參數,將測點作為單元節點來劃分計算單元網格,從而建立拱壩和壩基的三維有限元分析整體模型。模型的計算范圍為:壩軸線上、下游分別取350、250 m;壩底向基巖延伸150 m;左、右岸各取250 m。其中左岸方向為x軸正方向,下游方向為y軸正方向,z軸向上為正。整個模型可以分成兩個部分:壩體部分以及拱壩基礎部分。其中壩體部分除壩身外主要包括的細部結構有:左右岸溢洪道,壩間分縫,閘墩,壩后背管以及橫向廊道。拱壩基礎部分則包含重力墩以及巖體。由于本文主要研究重點為壩體,故而將左右溢洪道只模擬其與壩體相關聯部分,忽略廊道作用,而拱壩的分縫則使用薄層來模擬。由東江拱壩多年運行經驗可知東江拱壩基礎中存在的巖層軟弱帶與巖層破碎帶范圍小,對大壩運行影響可以忽略不計,因此,基礎除重力墩以外部分視作整體,不考慮巖層破碎帶的影響。東江拱壩三維有限元數值模型如圖1,上、下游方向取y向法向約束,左右岸方向取x向法向約束,底面取z向法向約束。其單元數71 061個,節點數96 471個,其中壩體單元43 318個,壩基單元27 743個。

圖1 東江拱壩有限元計算模型Fig.1 The finite element calculation model of Dongjiang arch dam


表1 292 m高程處8號測點計算溫度與實測溫度對比Tab.1 Comparison of the calculated temperature and measured temperature at 8 points of 292 m height

表2 292 m高程處8號測點計算位移與實測位移對比 mTab.2 Comparison of the calculated displacement and measured displacement at 8 points of 292 m height
由表1與表2中的溫度變化所引起位移變化的分析表明,采用有限元法分析計算拱壩壩體變化溫度場是可行的。當計算變化溫度場與實測溫度場較符合時,且 較小時,由此計算的由水位變化和溫度變化引起的位移變化值與實測位移變化值相差很小。當 增大的時候,誤差絕對值增大,但是其數值依然小于2 mm,這是由于運行期混凝土拱壩的時效變形比較穩定的緣故。因此,在壩體沒有埋設溫度計,或者溫度監測資料不完整的條件下,可以通過有限元法計算壩體溫度場,再采用本文所述方法,即公式(3)檢驗假定壩體溫度邊界條件,從而提高位移確定性模型精度。
通過有限元數值分析方法計算東江拱壩溫度場,由熱-結構耦合分析求得溫度荷載作用條件下的壩體位移。11號壩段前方交會8號測點在2003年2月-2012年12月時段內的溫度位移分量變化過程線如圖2所示。

圖2 8號測點溫度位移分量變化過程線(單位:mm)Fig.2 Temperature displacement component change process line of point 8
選取溫度位移計算日之前的某幾個特征時段內的氣溫平均溫度做因子。根據東江拱壩變形統計分析的經驗,參考不同測點的實測壩溫與氣溫過程線的滯后情況,選取觀測日前1、5、15、30、60、100天共6個特征時段的平均氣溫作為溫度預報因子,特征氣溫因子是相對于位移計算基準日對應特征氣溫的增量。根據公式(5)采用多元回歸方法與理論計算所得之溫度位移過程進行線性擬合,所得計算結果如下:

表3 溫度位移分量回歸系數Tab.3 Temperature displacement component regression coefficient
其中相關系數R2=0.924 7接近于1,統計量F=108.5487,與統計量F對應的概率P=0,小于顯著性水平 ,說明回歸在 水平上顯著,回歸方程有效。則溫度位移分量的回歸模型為:
δT=1.096 1-0.159 6T1+0.148 7T5+0.252 8T15+
0.555 6T30-1.265 8T60+1.994 9t100
(20)
式中:T1、T5、T15、T30、T60、T100分別表示觀測日前1、5、15、30、60、100 d 6個特征時段的平均氣溫。由表3中發現,B5的置信區間不包含零點,這說明對于東江拱壩而言,在考慮拱壩的氣溫滯后效應計算壩體溫度場時,計算日前100 d的氣溫平均值對于壩體溫度場有較大影響。同時也說明用氣溫變化近似代替環境溫度變化對大壩的影響是可行的。
通過有限元數值分析方法計算東江拱壩在靜水壓力荷載作用下壩體位移。11號壩段前方交會8號測點在2003年2月-2012年12月時段內的水壓位移分量變化過程線如圖3所示。

圖3 8號測點水壓位移分量變化過程線(單位:mm)Fig.3 Hydraulic displacement component change process line of point 8
由前文所述,水位位移分量是由于壩體或基巖在庫水位作用下產生的彈性變形所引起的,水荷載一經加上,變形便立即產生,因此不考慮位移對庫水位的滯后,而取位移觀測當天的日平均水位進行計算。根據東江大壩變形統計分析經驗,取壩體水平位移為水位的4次多項式函數。根據公式(4)采用線性回歸分析方法得8號測點對應部位的計算水壓水平位移與壩前水深的擬合結果如表4。

表4 水壓位移分量回歸系數Tab.4 Hydraulic displacement component regression coefficient
表4中SSE是擬合誤差的平方和,SSE=3.00說明本次擬合結果比較準確。R-square代表實測數據與擬合模型計算的數據之間相關系數的平方值,它越接近1,說明模型能更好地解釋變量間的比例關系。R-square=0.998 4說明本次擬合的兩組數據的相關性很好。RMSE是均方差,RMSE=0.223 8近于0說明本次擬合結果比較好,即回歸模型成立。其中相關系數R2=0.999 9接近于1,統計量F=5.147 9,與統計量F對應的概率P=0,小于顯著性水平 ,說明回歸在α水平上顯著,回歸方程有效。
水壓位移分量的回歸模型為:
δH=-0.272 5-0.968 3H-0.013H2-1.86e-04H3-
7.585e-06H4
(21)
式中:H為計算壩前水位。
由表4可知,回歸模型系數的置信區間都不包含零點,這說明該影響因子對于模型因變量的影響顯著。
由前文計算得到溫度位移分量以及水壓位移分量確定性模型中的待定系數,即Ai與Bj。然后根據東江大壩1/3拱壩處8號前方交會測點在2003年2月-2013年12月位移實測數據,采用數學回歸分析方法計算出公式(6)中的時效位移分量系數C1、C2以及公式(7)中的調整系數X、Y,得到東江大壩的位移確定性模型為:
δ=1.090 1δH+1.042 7δT+0.458 5lnθ-0.000 3θ-0.105 2
(22)
模型的復相關系數R2=0.916,統計量F=262.5,與統計量F對應的概率P接近于0,小于顯著性水平α=0.05,說明回歸計算值在α水平上顯著,回歸方程有效,此模型成立。
當不考慮時效位移時,測點總位移δ′為溫度位移分量(δH)以及水壓位移分量(δT)之和,即:
δ′=δH+δT
(23)
式(22)是利用東江大壩前方交會8號測點基于2003年2月-2012年12月的實測資料所建立的位移確定性監控模型。為了校驗此位移確定性模型精度,將東江大壩2003年2月-2013年12月的實測氣溫及水位代入模型,得到一組位移計算值,其為模型擬合預測值,其中2003年2月-2012年12月時段的位移計算值為模型擬合值,2012年12月-2013年12月時段的位移計算值為模型預測值。圖4為11#壩段8號測點位移實測值、模型擬合預測值與不考慮時效擬合值的對比圖。從圖中可以看出,確定性模型的擬合精度及預測精度較好,能滿足工程精度要求,證明本文所述建立確定性模型的方法可行。同時,圖4表明,對于長期運行的混凝土拱壩而言,不考慮時效的計算位移與位移實測值相差很小,證明時效位移在拱壩總位移中所占比例很小,幾乎可忽略不計。

圖4 8號測點位移實測值與模型擬合值、預測值對比圖Fig.4 Displacement contrast diagram of measured values and model fitting values and predicted values of point 8
(1)本文采用壩區氣溫資料構造溫度位移分量,提出以短時間序列實測資料來檢驗混凝土拱壩溫度場的新方法,能較好的消除估算混凝土拱壩溫度邊界條件對混凝土確定性模型的影響,這對于大壩監控是非常重要的。
(2)針對運行期混凝土拱壩的變形特點,采用三維有限元數值分析方法計算混凝土拱壩位移的水壓位移分量和溫度位移分量,建立混凝土拱壩位移的確定性模型。此位移確定性模型可以實現對水壓、溫度及時效位移分量的分離,其精度與預測監控性能較好,為我國大部分缺乏壩體溫度監測資料的運行期混凝土拱壩的安全監控提供了一個新的思路。
(3)對于運行多年的混凝土拱壩而言,時效位移分量相對于總位移所占比例很小,有時可以忽略不計,因此建立位移確定性模型時可不必包含時效分量,將時效分量和殘差放在一起與模型分開研究。