王 建，盧兆輝，路振剛

(1.河海大學(xué)水資源高效利用與工程安全國家工程研究中心，江蘇南京 210098;2.河海大學(xué)水文水資源與水利工程科學(xué)國家重點(diǎn)實(shí)驗(yàn)室，江蘇南京 210098;3.國網(wǎng)新源控股有限公司，北京 100761)

溫度對重力壩、拱壩等大體積混凝土的應(yīng)力具有重要的影響，因此溫度場的計算是大體積混凝土結(jié)構(gòu)設(shè)計的重要內(nèi)容之一。對于溫度場和溫度應(yīng)力，國內(nèi)外有關(guān)學(xué)者［1-2］開展了大量的研究，取得了廣泛而有價值的成果。準(zhǔn)確計算溫度場的前提之一是獲取可靠的混凝土熱學(xué)參數(shù)，然而，這些參數(shù)通常是由室內(nèi)試驗(yàn)甚至經(jīng)驗(yàn)得出的，有時難以真實(shí)反映實(shí)際情況。因此，如能根據(jù)現(xiàn)場監(jiān)測數(shù)據(jù)進(jìn)行反演分析，對于準(zhǔn)確把握熱學(xué)參數(shù)、指導(dǎo)類似工程具有重要的意義，國內(nèi)外已經(jīng)有學(xué)者對此進(jìn)行了較為深入的研究［3-11］。但現(xiàn)場監(jiān)測數(shù)據(jù)一般是基于不同時刻、不同空間位置的溫度測值，對于某些工程，如早期建設(shè)的老壩，常常缺乏壩體溫度監(jiān)測資料，如何獲取這些工程的熱學(xué)參數(shù)，如導(dǎo)溫系數(shù)等，目前仍缺乏有效的手段。

變形監(jiān)測是大壩最常規(guī)的監(jiān)測項(xiàng)目，對于中等規(guī)模以上的混凝土壩，通常有完善的變形監(jiān)測數(shù)據(jù)。如果能夠從這些變形監(jiān)測數(shù)據(jù)中挖掘出混凝土導(dǎo)溫系數(shù)的某些特征，對深入掌握大壩的工作性態(tài)是有幫助的。鑒于此，本文嘗試提出一種基于變形的壩體導(dǎo)溫系數(shù)反演方法。

1 壩體導(dǎo)溫系數(shù)反演方法

熱傳導(dǎo)基本方程為

式中:T——溫度;τ——時間;a——混凝土的導(dǎo)溫系數(shù);x，y，z——空間坐標(biāo);θ——混凝土絕熱溫升，對于老壩，由于水化熱已經(jīng)消散，溫度場計算時可不考慮θ的影響。

實(shí)踐表明:受氣溫周期性變化的影響，壩體位移也表現(xiàn)出明顯的周期性特征，且其相位滯后于氣溫，如圖1所示;而且，滯后時間Δ與混凝土導(dǎo)溫系數(shù)a具有很強(qiáng)的相關(guān)性。因此，如果能準(zhǔn)確計算出Δ，那么，結(jié)合溫度場及溫度變形的有限元計算，則可以反演得到壩體的綜合導(dǎo)溫系數(shù)。

1.1 溫度分量的提取

在大壩安全監(jiān)測資料分析中，常通過建立監(jiān)控模型，如統(tǒng)計模型，將觀測效應(yīng)量(如壩體位移)分解為水壓分量、溫度分量和時效分量［12］。其中，溫度分量用于刻畫效應(yīng)量受溫度變化的影響程度。對于運(yùn)行多年的老壩，壩體溫度變化進(jìn)入準(zhǔn)穩(wěn)定狀態(tài)，壩體變形的溫度分量呈周期性變化，且相位特征基本不受壩體材料性質(zhì)(如各種材料的分布及其彈性模量、熱膨脹系數(shù)、泊松比等)以及地基巖體地質(zhì)條件(結(jié)構(gòu)面產(chǎn)狀、結(jié)構(gòu)面充填物的力學(xué)參數(shù)、基巖變形模量及其分布等因素)等非周期性因素的影響，這種情況下，采用統(tǒng)計模型分離效應(yīng)量的溫度分量精度總體上較高。

當(dāng)對壩體變形監(jiān)測數(shù)據(jù)采用統(tǒng)計監(jiān)控模型進(jìn)行分析時，溫度分量可采用余弦函數(shù)表示，即

圖1 壩頂水平位移滯后于氣溫示意圖Fig.1 Schematic diagram of phase lag of dam＇s horizontal displacement to air temperature

1.2 滯后時間的確定

類似的，氣溫也可以采用余弦函數(shù)來表示。設(shè)氣溫的形式為

式中:Ta——?dú)鉁?Tam——年平均氣溫;Aa——?dú)鉁卣穹?φ——?dú)鉁爻跏枷辔?，rad。

效應(yīng)量的溫度分量與氣溫之間的相位差則代表了滯后時間。如圖1所示，氣溫及效應(yīng)量的溫度分量均為余弦函數(shù)，當(dāng)規(guī)定壩體位移向下游為正時，溫度下降時位移增加，因此兩者呈負(fù)相關(guān)關(guān)系，且壩頂位移曲線略滯后于氣溫。

對比式(2)，(3)可知，溫度分量滯后于氣溫的相位差為

將相位角轉(zhuǎn)化為天數(shù)，則滯后時間為

1.3 導(dǎo)溫系數(shù)與滯后時間的關(guān)系

僅考慮溫度荷載作用，計算出壩體內(nèi)的準(zhǔn)穩(wěn)定溫度場以及由此產(chǎn)生的溫度變形，通過相位的比較，可以得到測點(diǎn)處壩體位移的峰值相對于氣溫的滯后時間Δ。假定不同的壩體材料導(dǎo)溫系數(shù)ai，則根據(jù)有限元模擬，可以得到多個反映導(dǎo)溫系數(shù)ai與滯后時間Δi關(guān)系的樣本點(diǎn)。計算表明:壩頂位移滯后于氣溫的時間通常隨著導(dǎo)溫系數(shù)的增加而降低，曲線如圖2所示。根據(jù)曲線的特征，可以采用多項(xiàng)式擬合壩體導(dǎo)溫系數(shù)與滯后時間的關(guān)系，一般采用二次多項(xiàng)式可以很好地擬合:

式中，c0，c1，c2為系數(shù)，可采用最小二乘法等優(yōu)化方法確定。

這里應(yīng)說明，計算壩頂位移時，所涉及的邊界溫度除了氣溫以外，還包括水溫。對于有水溫監(jiān)測的大壩，可以直接采用監(jiān)測值;若無水溫監(jiān)測資料，則可按經(jīng)驗(yàn)公式確定，如文獻(xiàn)［13］中的方法。

1.4 導(dǎo)溫系數(shù)反演

當(dāng)根據(jù)式(5)獲得壩頂位移相對于環(huán)境溫度的滯后時間后，則壩體混凝土的導(dǎo)溫系數(shù)為

式中，f-1為函數(shù)f的反函數(shù)。

由于f為二次多項(xiàng)式形式，即式(6)，因此，a可通過直接求解二次方程得到。

圖2 混凝土導(dǎo)溫系數(shù)與滯后時間的關(guān)系Fig.2 Relationship between phase lag and thermal diffusivity

2 計算實(shí)例

豐滿水電站是1949年中華人民共和國成立前我國修建的第一座大型水力發(fā)電工程。樞紐工程由混凝土重力壩、溢流壩、壩后式廠房及變電站、左岸泄洪洞(后改建成三期引水洞)及三期發(fā)電廠房組成，最大壩高91.7 m，壩頂長度1080 m，壩頂高程267.70m，共有60個壩段。這里以38～40號擋水壩段2001—2010年間壩頂順河向的水平位移為例，對豐滿老壩壩體材料的導(dǎo)溫系數(shù)進(jìn)行反演分析。

2.1 滯后時間確定

由于豐滿大壩為重力壩且規(guī)定位移向下游為正，溫降時壩頂位移增加，因此，最大位移出現(xiàn)時間與最低氣溫值出現(xiàn)時間之差即為滯后時間Δ。為了避免復(fù)雜邊界以及引水或泄水等不確定因素對壩體溫度場計算的影響，本文選擇接近河床中央且遠(yuǎn)離溢流壩段和廠房的38～40號3個擋水壩段的壩頂水平位移建立監(jiān)控模型:

式中:δ——壩頂位移;H——上游水位;a0，ai，b1，b2，ci——回歸系數(shù)。括號內(nèi)為溫度分量形式，通過和差化積可轉(zhuǎn)化為式(2)的形式。

監(jiān)控模型的復(fù)相關(guān)系數(shù)、擬合標(biāo)準(zhǔn)差、溫度分量等主要成果見表1。從表1可以看出，各測點(diǎn)溫度分量的初始相位比較接近，平均值約為0.4348。

表1 監(jiān)控模型主要成果Table 1 General results of regression equation for safety monitoring

氣溫初始相位根據(jù)豐滿壩址區(qū)氣溫實(shí)測資料按余弦函數(shù)擬合后得到，擬合函數(shù)為

本文始測日取2010年1月1日。

根據(jù)各測點(diǎn)的平均相位以及式(9)中的相位，由式(5)可以得到Δ=11.8 d。

2.2 溫度場計算

分別假設(shè)壩體導(dǎo)溫系數(shù)為 0.050 m2/d，0.075 m2/d，0.100 m2/d，0.125 m2/d，0.150 m2/d，進(jìn)行溫度場計算。各工況下，壩體初始溫度均取多年平均氣溫。為了消除初始溫度對準(zhǔn)穩(wěn)定溫度場的影響，溫度場計算時間較長，共30 a;而計算壩頂變形，為了減小計算量，僅取最后1年(2010年)進(jìn)行計算。計算軟件采用ABAQUS。

溫度場計算時的邊界條件包括氣溫和水溫。其中，壩體上游面在2001—2010年間的平均庫水位以下采用水溫邊界，其余采用氣溫邊界。氣溫邊界按式(9)計算，水溫邊界利用豐滿水庫溫度實(shí)測資料的擬合公式計算［13］:

式中為y水深，m。

圖3為導(dǎo)溫系數(shù)取0.075m2/d，氣溫最低時的壩體溫度場分布云圖。為盡量提高壩體顯示的比例尺，增加云圖的清晰度，壩基部分未在圖中給出。

2.3 導(dǎo)溫系數(shù)的反演

根據(jù)不同導(dǎo)溫系數(shù)下溫度場的成果，假設(shè)彈性模量為25GPa，線膨脹系數(shù)為1.0×10-5，采用ABAQUS有限元軟件分別對2010年的壩體溫度變形進(jìn)行了計算，并輸出各個時刻測點(diǎn)位置的變形量，由此可確定最大位移出現(xiàn)時間。將其與氣溫最低時間進(jìn)行對比，即可得到不同導(dǎo)溫系數(shù)下壩頂位移滯后于氣溫的時間。這里要說明的是，盡管位移計算時彈性模量與線膨脹系數(shù)是假設(shè)的，但對壩體位移周期性的相位特征，如最大位移出現(xiàn)時間是沒有影響的，因此不會給滯后時間的確定帶來誤差。

分別取導(dǎo)溫系數(shù)為0.050 m2/d，0.075 m2/d，0.100 m2/d，0.125 m2/d，0.150 m2/d，得到對應(yīng)的滯后時間為14.1 d，12.8 d，11.7 d，10.8 d，10.0 d，圖2為對應(yīng)曲線。對上述結(jié)果按式(6)進(jìn)行擬合，得到回歸方程(11)，其復(fù)相關(guān)系數(shù)R=0.999，接近1.0，擬合精度很高。

將Δ=11.80代入式(11)，求解二次方程并舍去不合理的解，得到壩體的導(dǎo)溫系數(shù)為0.0978 m2/d。

圖3 低溫時壩體溫度場分布云圖Fig.3 Contour of dam temperature when air temperature is low

3 結(jié) 論

a.對于運(yùn)行多年的老壩，壩體溫度變化進(jìn)入準(zhǔn)穩(wěn)定狀態(tài)，壩體變形的溫度分量呈周期性變化，且相位特征基本不受壩體材料性質(zhì)(如各種材料的分布及其彈性模量、熱膨脹系數(shù)、泊松比等)以及地基巖體地質(zhì)條件(結(jié)構(gòu)面產(chǎn)狀、結(jié)構(gòu)面充填物的力學(xué)參數(shù)、基巖變形模量及其分布等因素)等非周期性因素的影響，這種情況下，采用統(tǒng)計模型分離溫度分量的精度總體上較高。

b.導(dǎo)溫系數(shù)a的反演主要是利用壩體位移的溫度分量滯后于氣溫的時間Δ，因此，Δ的精度將直接影響a的反演精度。本文通過建立監(jiān)控模型分離出余弦形式的壩體溫度分量，同時對氣溫也采用余弦函數(shù)擬合，利用兩者間的相位差來確定壩體位移滯后于氣溫的時間Δ，而不是以目測的方法，保證了Δ的準(zhǔn)確性。

c.從豐滿大壩導(dǎo)溫系數(shù)的反演成果來看，a=0.0978 m2/d在合理范圍之內(nèi)，可以認(rèn)為反演的成果是可信的。因此，本文方法是可行的。不過，若將該數(shù)值與常規(guī)混凝土導(dǎo)溫系數(shù)總體均值［1］相比，略顯偏低，這是否與壩體混凝土老化有關(guān)目前尚無法定論。

d.由于壩體取不同導(dǎo)溫系數(shù)時算出的滯后時間的變化范圍較窄，使得導(dǎo)溫系數(shù)對滯后時間較為敏感，這是本文方法的不足之處。因此，當(dāng)壩體內(nèi)布置有溫度計時，宜優(yōu)先采用基于壩體溫度測值的反演方法。但對于缺乏壩體溫度測值的老壩，本文方法則具有明顯優(yōu)勢，可為大壩安全監(jiān)測的分析評估提供有益信息。

［1］朱伯芳.大體積混凝土溫度應(yīng)力與溫度控制［M］.北京:中國電力出版社，1998.

［2］CARSLAW H S，JAEGER J C.Conduction of heat in solids［M］.2th ed.Oxford:Oxford University Press，1959.

［3］KOHN R，VOGELINSM.Determining conductivity by boundary measurements［J］.Common Pure Applied Mathematics，1984，37(3):289-298.

［4］KUNIHITO M，NORIYUKI N，YOSHIRERU D，et al.Inverse method for estimation of thermal properties of mass concrete［J］.Transactions of the Japan Concrete Institute，1994，16:131-138.

［5］NAKAMURE M，TANAKA M，ADACHI O.Simulation of temperature control for transient heat conduction fields using boundary element inverse analysis［C］//Proceeding of the 1995 Joint ASME/JSME Pressure Vessels and Piping Conference.New York:ASME，1995:19-24.

［6］劉寧，張劍，趙新銘.大體積混凝土結(jié)構(gòu)熱學(xué)參數(shù)隨機(jī)反演方法初探［J］.工程力學(xué)，1996，20(5):114-120.(LIU Ning，ZHANG Jian，ZHAO Xinming.Stochastic back analysis of thermal parameters of mass concrete structures［J］.Engineering Mechanics，1996，20(5):114-120.(in Chinese))

［7］LESNIC D，ELLIOTT L，INGHAM D B.An inverse coefficient identification problem for the heat equation［C］//International Symposium on Inverse Problems in Engineering Mechanics.Nagano，Japan:［s.n.］，1998:11-16.

［8］NOWAK I，NOWAK A J，WROBEL L C.Boundary and geometry inverse thermal problems in continuous casting［C］//International Symposium on Inverse Problems in Engineering Mechanics Ⅳ.Nagano，Japan:［s.n.］，2003:21-31.

［9］張宇鑫，宋玉普，王登剛.基于遺傳算法的混凝土三維非穩(wěn)態(tài)溫度場反分析［J］.計算力學(xué)學(xué)報，2004，21(3):338-342.(ZHANG Yuxin，SONG Yupu，WANG Denggang.Inverse analysis of three dimensional unsteady temperature field with genetic algorithms［J］.Chinese Journal of Computational Mechanics，2004，21(3):338-342.(in Chinese))

［10］JEONG J H，ZOLLINGER D G.Finite-element modeling and calibration of temperature prediction of hydrating Portland cement concrete pavements［J］.Journal of Materials in Civil Engineering，2006，18(3):317-324.

［11］KUZMANOVIC V，SAVIC L，STEFANAKOS J.Long-term thermal two-and three-dimensional analysis of roller compacted concrete dams supported by monitoring verification［J］.Canadian Journal of Civil Engineering，2010，37(4):600-610.

［12］吳中如.水工建筑物安全監(jiān)控理論及其應(yīng)用［M］.北京:高等教育出版社，2003.

［13］朱伯芳.庫水溫度估算［J］.水利學(xué)報，1985(2):12-21.(ZHU Bofang.Prediction of water temperature in reservoirs［J］.Journal of Hydraulic Engineering，1985(2):12-21.(in Chinese))