用于抽水蓄能電站大壩不穩(wěn)定溫度場仿真的水平分區(qū)異步長解法

孫 偉，何蘊龍，王 南，熊 堃

(1. 武漢大學(xué)水資源與水電工程科學(xué)國家重點實驗室，武漢 430072；2. 長江勘測規(guī)劃設(shè)計研究有限責任公司，武漢 430010)

用于抽水蓄能電站大壩不穩(wěn)定溫度場仿真的水平分區(qū)異步長解法

孫 偉1，何蘊龍1，王 南1，熊 堃2

(1. 武漢大學(xué)水資源與水電工程科學(xué)國家重點實驗室，武漢 430072；2. 長江勘測規(guī)劃設(shè)計研究有限責任公司，武漢 430010)

對于抽水蓄能電站(PSPS)的攔河大壩上游迎水面，由于水庫水位周期性升降變化，會出現(xiàn)與水、氣等不同介質(zhì)接觸的交替變化現(xiàn)象. 將具有這種現(xiàn)象的迎水壩面邊界稱之為特殊邊界. 這種交替變化的邊界條件對特殊邊界表面及其影響范圍內(nèi)的溫度場分布將產(chǎn)生較大的影響，在溫度仿真分析中不能忽略. 提出了考慮特殊邊界影響的抽水蓄能電站大壩不穩(wěn)定溫度場水平分區(qū)異步長解法，在分析特殊邊界特征的基礎(chǔ)上，對水平分區(qū)異步長解法的基本原理進行了推導(dǎo)，并提出了等效放熱系數(shù)的概念，從平均意義上近似考慮表面混凝土與不同介質(zhì)接觸時放熱系數(shù)處于變化的問題，給出了水平分區(qū)界限的位置范圍. 將該解法用于呼和浩特PSPS下水庫攔河壩的不穩(wěn)定溫度場的仿真分析中. 結(jié)果表明，提出的方法是有效和省時的，并為計算更為復(fù)雜耗時的溫度應(yīng)力場提供了可能. 最后，對嚴寒地區(qū)特殊邊界處可能出現(xiàn)的結(jié)冰與解凍問題進行了初步探討.

抽水蓄能電站；水平分區(qū)；異步長；不穩(wěn)定溫度場；等效放熱系數(shù)

抽水蓄能電站(pumped storage power station，PSPS)是電力系統(tǒng)發(fā)展到一定階段的產(chǎn)物．電力系統(tǒng)日負荷圖在每日上、下午各有一個高峰，午夜則有一個低谷．PSPS可利用夜間低谷負荷時火電站提供的剩余電能，從高程低的下水庫抽水到高程高的上水庫中，通過水體這一能量載體將電能轉(zhuǎn)化為水的位能．在日間出現(xiàn)高峰負荷時，再從上水庫放水發(fā)電，擔任負荷圖中的峰荷部分[1]．這意味著PSPS具有特殊的運行方式，在較短的周期內(nèi)承擔調(diào)峰填谷的任務(wù)，水庫水位變幅較大．由于存在水庫水位的周期性升降變化，電站的攔河大壩上游迎水面會出現(xiàn)與水、氣等不同介質(zhì)接觸的交替變化現(xiàn)象．本文將具有這種現(xiàn)象的迎水壩面邊界稱為特殊邊界．這種特殊邊界的邊界條件在與水接觸的第1類邊界條件、與大氣接觸的第3類邊界條件、寒冷季節(jié)甚至可能與冰接觸的第4類邊界條件之間不斷變化．特殊邊界的某一點屬于其中哪一類邊界條件取決于該點相對水庫水位的空間位置、PSPS泵站運行工況和電站運行工況的分配過程、仿真計算的時間位置以及混凝土表面水分傳輸機制等．由于水庫水位變幅大、周期短，特殊邊界所處的邊界條件一直處于變動中，所引起的特殊邊界表面及其影響范圍內(nèi)的溫度變化不能忽略．

在具有特殊邊界的PSPS大壩不穩(wěn)定溫度場的仿真分析中，由于運行周期的制約，常規(guī)水庫大壩運行期較長的計算步長不再適用．一方面，具有特殊邊界的這種情況與短期遭遇寒潮、晝夜溫差不同，電站在工作期間，水庫水位常年處于這種周期性升降的狀態(tài)．計算時常年都采用短步長，大壩不穩(wěn)定溫度場需要求解的計算量將是巨大的，再進一步顧及到計算量更大的溫度應(yīng)力場，實際上是難以實現(xiàn)的；另一方面，這種特殊邊界條件對壩體不穩(wěn)定溫度場的影響范圍有限，對那些大部分溫度變化不受特殊邊界影響的地方采用短步長計算是不經(jīng)濟的．

常規(guī)水庫大壩的溫度場及溫度應(yīng)力場的研究[2-5]表明，由于氣溫和水溫的年變化，運行期的溫度應(yīng)力仍可能引起混凝土的開裂．對于具有特殊邊界的PSPS攔河大壩，氣溫和水溫的變化周期減小，表面溫度的交替變化可能會帶來運行期溫度應(yīng)力分布特征的改變，例如在冬季，蓄水一般相當于保溫，運行后伴隨著水位的升降，壩面就如同短時間內(nèi)在保溫和沒有保溫之間不斷循環(huán)，這自然會引起特殊的溫度及溫度應(yīng)力分布特征．

國內(nèi)外有關(guān)PSPS的這種周期短、交替變化快的特殊邊界的處理，鮮見報導(dǎo)．文獻[6]提出了高壩溫度仿真分析的分區(qū)異步長解法，該方法用于解決高壩仿真分析中由于新老混凝土溫度敏感性不同而建立的減小計算量的有效方法，在大壩空間域看來是一個由于澆筑時間、邊界條件、內(nèi)源熱量變化劇烈程度等不同而建立的豎直分區(qū)異步長解法．本文借鑒分區(qū)異步長解法的思想，將其應(yīng)用到考慮PSPS特殊邊界影響的運行期不穩(wěn)定溫度場計算中．在大壩空間域上看來，這是由于受特殊邊界影響程度不同而建立的水平分區(qū)異步長解法．與豎直分區(qū)異步長解法不同的是，水平分區(qū)異步長解法特殊邊界處不僅接觸的溫度在變化，而且接觸介質(zhì)不同，表面的放熱系數(shù)也在變化，因此本文在不穩(wěn)定溫度場分解過程中提出等效放熱系數(shù)的概念，從平均意義上近似考慮表面混凝土與不同介質(zhì)接觸時放熱系數(shù)也處于變化的問題．通過某PSPS的工程實例分析，將這種水平分區(qū)異步長解法的計算結(jié)果與整個壩體統(tǒng)一短步長解法的計算結(jié)果進行比較，證明本文解法的合理性和有效性．最后針對嚴寒地區(qū)PSPS寒冷季節(jié)迎水面水位升降變化過程中存在的更為特殊的結(jié)冰與解凍問題進行了初步探討．

1 特殊邊界特征分析

本節(jié)所述及的特殊邊界特征的分析均針對普通常溫地區(qū)PSPS而言，特殊邊界只可能是水邊界或氣邊界，而不涉及大壩上游迎水面的結(jié)冰與解凍的問題.有關(guān)冰層覆蓋與融化的問題將在第4節(jié)進行探討.

PSPS攔河大壩特殊邊界特征分析如圖1所示.設(shè)點A1為大壩壩面特殊邊界上任一點，點A1距離死水位高度為H，正常蓄水位與死水位高度差為H′. 下面從兩個方面研究點A的邊界特征．

1)點A1為水邊界或氣邊界的條件

水庫水位高程和干濕交替變化條件下混凝土(或混凝土表面附著的永久保溫材料)表層水分的傳輸機制決定了點A1的邊界類型．

圖1 PSPS攔河大壩特殊邊界特征分析Fig.1 Analysis of the characteristics of special boundary of PSPS dam

水庫水位高程的變化與PSPS的泵站運行工況和電站運行工況的分配過程有關(guān)．不失一般性，將PSPS運行周期(1,d)分為4個時段：①下水庫抽水時段(零時至凌晨T1時)；②下水庫低水位時段(凌晨T1時至下午T2時)；③下水庫蓄水時段(下午T2時至晚間T3時)；④下水庫高水位時段(晚間T3時至24時) ．水庫的水位-時間曲線如圖2所示．

圖2 水庫的水位-時間曲線Fig.2 Curve of water level-time of reservoir

干濕交替條件下混凝土表層水分的傳輸機制屬于混凝土濕度場的研究范疇．Bazant等[7]在1972年以混凝土內(nèi)部相對濕度為基本變量，用擴散方程的形式描述混凝土在干燥過程中內(nèi)部相對濕度的變化，即

式中：H為混凝土相對濕度；t為時間；DH(H)為隨H變化的擴散系數(shù)．式(1)為一個非線性方程，需迭代求解．利用該擴散方程，在給定的初始濕度和邊界濕度條件下，可求出混凝土多孔介質(zhì)的濕度場分布．然而，一方面DH(H)與溫度、混凝土組成材料、配合比(特別是水膠比)、開裂情況、應(yīng)力狀況等均有關(guān)系，加上壩體施工期永久保溫材料(蓄水后運行期可能處于部分失效狀態(tài))的影響，故難以準確確定；另一方面，濕度場、溫度場相互影響，兩者的耦合計算更加大了特殊邊界條件下的非穩(wěn)定溫度場的計算難度，不便應(yīng)用于實際工程．

文獻[8]對混凝土近表面層一維水分傳輸問題進行了實驗研究．結(jié)果表明：升溫干燥時，水分傳輸?shù)慕砻鎸雍穸群苄?，近表面水蒸氣含量梯度很快消失；等溫吸濕時，近表面不存在明顯的濕度梯度．基于此，為了便于工程應(yīng)用，本文采用以下簡單原則確定點A1的水氣邊界條件：①忽略水庫水位上升過程中吸濕過程的濕度梯度，認為點A1在水位高程以下，則視為水邊界；②忽略水位下降中干燥過程的水蒸氣含量梯度，認為點A1在水位高程以上，則視為氣邊界.

綜上所述，可得如下結(jié)論：①點A1在水位高程以下，則視為水邊界，邊界溫度等于相應(yīng)時刻的水溫；②點A1在水位高程以上，則視為氣邊界，邊界溫度等于相應(yīng)時刻氣溫．

2)在一個運行周期里，點A1為水邊界或氣邊界的時段分配

在PSPS一個運行周期1,d里，點A1(距離死水位高度為H)為水邊界的時間為1t，為氣邊界的時間為．以下討論t1、t2是如何分配的．假設(shè)在下水庫抽水及蓄水過程中，水位均勻變動．水邊界和氣邊界時段分配如圖3所示．

圖3中，1M、2M為以點A1處的水位H做平行于時間軸的直線與水庫水位-時間曲線的兩個交點的時間坐標，即點A1在0～1M和2M～24為水邊界，在

圖3 水邊界和氣邊界時段分配Fig.3 Allocation of the period of water and gas boundary

2 水平分區(qū)異步長解法基本原理

2.1 公式推導(dǎo)

水平分區(qū)異步長解法如圖4所示．特殊邊界影響區(qū)記為R1，非影響區(qū)記為R2，特殊邊界記為S，其余常規(guī)水邊界、氣邊界分別記為S1～S4，影響區(qū)和非影響區(qū)接觸面記為B1．周圍介質(zhì)溫度分別為：死水位以下S1邊界水溫w1T，S2、S3邊界氣溫gT，尾水位以下S4邊界水溫w2T，特殊邊界S溫度f(t)．蓄水后運行期混凝土水化熱溫升的影響已不必考慮．圖4中，點A～E(影響區(qū)上游點A、影響區(qū)內(nèi)部點B、分區(qū)界限位置點C、非影響區(qū)內(nèi)部點D、非影響區(qū)下游點E)是為了便于分析而取出的典型點．

2.1.1 總溫度場——T場

T場應(yīng)滿足下列條件．

(1) 在影響區(qū)域R1內(nèi)熱傳導(dǎo)方程為

式中：T為總溫度場；a為導(dǎo)溫系數(shù)．

式(4)的初始條件為

式(4)的邊界條件為

式中：λ為導(dǎo)熱系數(shù)；wβ為大壩混凝土與水接觸時的表面放熱系數(shù)；gβ為大壩混凝土與空氣接觸時的表面放熱系數(shù)．

式(4)的特殊邊界條件為

式中：β為特殊邊界處大壩混凝土與介質(zhì)接觸時的表面放熱系數(shù)．當與不同介質(zhì)接觸時，β值不同．

(2) 在接觸面上溫度滿足

式中1T、2T為接觸面兩側(cè)的溫度．

(3) 在非影響區(qū)域R2內(nèi)熱傳導(dǎo)方程為

式(11)的初始條件為

由于待求解的熱傳導(dǎo)方程(4)、(11)是線性的，故可將總溫度場T進行分解，即

2.1.2 分解溫度場—U場

分解的溫度場—U場應(yīng)滿足下列條件．

(1) 在影響區(qū)域R1內(nèi)熱傳導(dǎo)方程為

式(16)的初始條件為

式(16)的邊界條件為

式(16)的特殊邊界條件為

式中(x, y, z)∈ΓS．

(2) 在接觸面上溫度滿足

(3) 在非影響區(qū)域R2內(nèi)熱傳導(dǎo)方程為

式(23)的初始條件為

式(23)的邊界條件為

2.1.3 分解溫度場—V場

分解的溫度場—V場應(yīng)滿足下列條件．

(1) 在影響區(qū)域R1內(nèi)熱傳導(dǎo)方程為

式(27)的初始條件為

式(27)的邊界條件為

式(27)的特殊邊界條件為

式中：(x, y, z)∈ΓS；βe為等效放熱系數(shù)，具體分析見第2.2節(jié)．

在接觸面上滿足

在非影響區(qū)域R2內(nèi)熱傳導(dǎo)方程為

式(34)的初始條件為

式(34)的邊界條件為

由于待求解的熱傳導(dǎo)方程都是線性的，上述分解是嚴格的．

在U場中，溫度場的劇烈變化主要局限于影響區(qū)域R1，在影響區(qū)域R1與非影響區(qū)域R2的接觸面B1上，其影響已趨于0，即U1=U2=0，x∈ΓB1．非影響區(qū)域R2滿足下列條件：

熱傳導(dǎo)方程為

初始條件為

邊界條件為

顯然，滿足式(38)～式(42)的解為：在非影響區(qū)域R2中0U=．

因此，U場只需要在影響區(qū)域R1中求解．由于U場中特殊邊界條件隨時間而急劇變化，故計算中需要采用較小的時間步長1τΔ．V場在全區(qū)域R1＋R2內(nèi)求解，因其變化平緩，計算中采用較大的時間步長2τΔ．總溫度場的分解如表1所示．

需要說明的是，在取時間步長時需要注意，1τΔ和2τΔ是有關(guān)聯(lián)的，2τΔ必須為n個1τΔ的代數(shù)和．本文取1=1hτΔ，2=1dτΔ．

在接觸面上有

表1 總溫度場的分解Tab.1 Decomposition of total temperature field

水平分區(qū)異步長解法的求解步驟為：①求解U場，求解域為影響區(qū)域R1，時間步長為1,h；②求解V場，求解域為影響區(qū)域R1和非影響區(qū)域R2，時間步長為1,d；③將U場24,h末時刻的計算結(jié)果與V場1,d末時刻的計算結(jié)果相疊加，即得到最終T場1,d的分布．在實際工程中，疊加過程可以通過外部程序調(diào)用有限元軟件自動連續(xù)地加以實現(xiàn)．

2.2 等效放熱系數(shù)

水平分區(qū)異步長解法與文獻[6]中豎直向分區(qū)異步長解法不同之處在于，這種特殊邊界處不僅接觸的溫度在變化，而且由于接觸介質(zhì)不同，其表面放熱系數(shù)也在變化．在進行U場小步長計算時，可以隨著接觸介質(zhì)的變化設(shè)定相應(yīng)的表面放熱系數(shù)，而在V場長步長計算時，一個長步長中已經(jīng)包含了與不同介質(zhì)接觸的變化過程，因此需要設(shè)定一個合理的等效放熱系數(shù)，從平均意義上近似考慮表面混凝土與不同介質(zhì)接觸時的放熱系數(shù)處于變化的問題．

與第1節(jié)研究方法相同，仍假設(shè)點A為水邊界的時間為t1，為氣邊界的時間為t2，t1+t2=24 h ．

V場的特殊邊界處混凝土表面?zhèn)鬟f的總熱量為

式中：wQ為混凝土與水接觸時，混凝土表面?zhèn)鬟f的熱量；gQ為混凝土與空氣接觸時，混凝土表面?zhèn)鬟f的熱量；aV為混凝土與水接觸時的表面溫度；bV為混凝土與空氣接觸時的表面溫度．

此時邊界溫度無論是水溫還是氣溫，均保持不變．在氣邊界時，相當于空氣溫度變化非常緩慢，周期趨于無限大，虛擬厚度的影響可以忽略(實際上，變化周期為年的時候可以忽略[2])，可按第1類邊界條件近似處理，即假設(shè)Vb≈Va(如此假設(shè)是由于βw遠大于βg(βw≈320βg)，而且特殊邊界附近交替變化的溫度已經(jīng)在U場得到反映，特殊邊界附近U場占T場的比例很大，所以V場的近似處理所帶來的影響也會較小，因此取Vb≈Va影響也將不大) ．此時有

式中ΔV為V場中壩面混凝土表面溫度與外界水溫之差，ΔV=Va．

將總熱量Qtotal表示成熱流量與時間的乘積，式(46)可寫成

由式(47)有 e

β即為具有平均意義的特殊邊界等效放熱系數(shù).

特殊邊界范圍內(nèi)，位于不同高程的點，與水、氣接觸的時間長短不同(具體接觸時間的分析見第1節(jié))，因而其等效放熱系數(shù)也不同．為了提高精度，需將特殊邊界按澆筑層細分，根據(jù)每一層的水、氣接觸時間計算式(2)～(3)和等效放熱系數(shù)計算式(49)，求得每一層的等效放熱系數(shù)值．

2.3 分區(qū)位置的合理確定

在離特殊邊界較遠處，特殊邊界處交替的溫度變化對不穩(wěn)定溫度場的影響很?。∧硞€很小的溫度影響值ε0作為判別標準，對于溫度影響值ε＜ε0的部分，特殊邊界對其影響可忽略不計，以此可以確定特殊邊界的影響深度．

在以下的分析中，取溫度影響值ε0=0.01℃．影響深度d的確定可通過理論計算確定．

熱傳導(dǎo)方程為

初始條件為

邊界條件為

式中：P為溫度變化周期；B為表面溫度變幅．

拉普拉斯變換得到滿足上述條件的解為

對于氣邊界，考慮到虛擬厚度，g(,)T x t=因此氣溫邊界的影響衰減得更快，影響深度比水溫邊界小．

對于這種周期短、變幅大的溫度變化，應(yīng)取較小的周期P．取溫度變化周期為1,d，t以1,h為單位步長，計算得到水溫邊界的影響深度如表2所示．

表2 水溫邊界影響深度Tab.2 Influence depth of the water temperature boundary

取B＝20,℃，影響深度d＝1.36,m處溫度影響值滿足ε＜ε0．結(jié)合網(wǎng)格劃分及分區(qū)方便等因素，選定分區(qū)界限位置范圍為特殊邊界表面附近2～3,m．

3 工程實例分析

以呼和浩特PSPS為例，將考慮特殊邊界影響的水平分區(qū)異步長解法應(yīng)用于其下水庫攔河壩的不穩(wěn)定溫度場的仿真分析中．該電站為國內(nèi)第1例位于嚴寒地區(qū)的PSPS，位于內(nèi)蒙古自治區(qū)呼和浩特市東北的大青山脈中麓哈拉沁溝下游，距呼和浩特市中心公路里程約20,km．下水庫庫區(qū)所在流域?qū)儆谥袦貛Ъ撅L(fēng)亞干旱氣候區(qū)，具有冬長夏短、寒暑變化急劇的特征．冬季長達5個月，漫長而嚴寒．

呼和浩特PSPS一般從凌晨零時開始到凌晨7時從下水庫抽水至上水庫，調(diào)峰發(fā)電泄水主要是在16：00—20：00，下水庫庫盆內(nèi)水位在死水位(1,355.00,m)到正常蓄水位(1,400.00,m)之間變化，晝夜升降變幅為45,m，水位變幅大，周期短．

選取呼和浩特PSPS下水庫攔河壩為研究對象，驗證水平分區(qū)異步長算法有效性．大壩采用“金包銀”形式，即壩體上游側(cè)包有厚2～3,m的常態(tài)混凝土．為簡單起見，將影響區(qū)域R1和非影響區(qū)域R2的接觸面劃分在上游常態(tài)混凝土與碾壓混凝土交界面處，即B1取為上游常態(tài)混凝土與碾壓混凝土分界線，影響區(qū)域R1和非影響區(qū)域R2范圍見圖4．

設(shè)計算例1 計算時間取運行期夏季6月．以施工期及運行期蓄水開始至電站以抽水蓄能方式運行穩(wěn)定的溫度分布為初始條件．為反映特殊邊界處溫度變化是由于其與水溫和氣溫交替接觸所引起，而排除晝夜溫差等其他因素的影響，假設(shè)在此短時間(3,d)內(nèi)，水溫和氣溫不變．根據(jù)文獻[9]分析結(jié)果，水溫取為定值17.0,℃，氣溫取為定值19.5,℃．采用兩種解法對比分析．解法1精確模擬一天內(nèi)水庫實際變化過程，時間步長均勻為1,h，時長為72,h(以下簡稱整個壩體統(tǒng)一短步長解法，記為0T場)，計算結(jié)果為每天24,h末時刻溫度值．解法2采用水平分區(qū)異步長解法，即：①求解U場，求解域為影響區(qū)域R1，時間步長為1,h，時長為24,h；②求解V場，求解域為影響區(qū)域R1＋非影響區(qū)域R2，時間步長為1,d；③取U場24,h末時刻溫度計算結(jié)果與V場第1天末時刻的計算結(jié)果相疊加，得到最終第1天的T場．重復(fù)操作①、②、③，共計算3,d．計算結(jié)果為U場24,h末時刻與V場第1天末時刻計算結(jié)果疊加的溫度．

取1,377,m高程點A～E(見圖4)，算例1兩種解法典型點溫度值對比見表3．兩種解法壩體不穩(wěn)定溫度場分布對比見圖5．

表3 算例1兩種解法典型點溫度值對比Tab.3 Comparison of temperature value of typical points under two algorithms in example 1, ℃

圖5 算例1兩種解法壩體不穩(wěn)定溫度場分布對比(單位：℃)Fig.5Comparison of the unsteady temperature filed of the dam under two algorithms in example 1(unit：℃)

分析相應(yīng)的圖表，可得到以下結(jié)論．

(1) 從特殊邊界表面到分區(qū)界限處，U場越來越小，而V場越來越大，即特殊邊界附近U場占T場的比例較大，V場中采用等效放熱系數(shù)的近似處理引起的誤差較小，不超過0.15,℃．

(2) 總體上來說，采用水平分區(qū)異步長解法與整個壩體統(tǒng)一短步長解法在同一時刻大壩不穩(wěn)定溫度場分布十分接近，非影響區(qū)域R2溫度分布吻合較好，影響區(qū)域R1則在上游邊界和分區(qū)界限位置處由于V場的簡化處理存在誤差，但溫度差值不超過0.15,℃．

(3) 解法1時間步長均勻為1,h，時長為72,h，在此期間，大壩上游面不穩(wěn)定溫度場變化不大，變幅為0.55,℃，最大值約為18.30,℃，出現(xiàn)在表面特殊邊界處1,360,m高程附近．而解法2只得到了3,d中3個時刻(24 ht=，48 h，72 h)的計算結(jié)果，大壩上游面不穩(wěn)定溫度場變幅為0.35,℃，最大值約為18.05,℃，變幅和最大值均較解法1小，原因是解法2每天的溫度分布是由U場24,h末時刻溫度計算結(jié)果與V場該天末時刻的計算結(jié)果相疊加得到的，因此只能反映周期為1,d的溫度變幅和最大值．通常來說，運行期不穩(wěn)定溫度場仿真分析通常以1,d(甚至更大)為計算步長，關(guān)注周期為1,d的溫度變幅和最大值已經(jīng)足夠滿足工程分析的需要．因此，水平分區(qū)異步長解法在溫度變幅和最大值方面雖有誤差，但能夠滿足分析的需要．

(4) 水平分區(qū)異步長解法只需要在特殊邊界影響的較小范圍內(nèi)采用短步長計算劇烈變化的U場，因而比整個壩體統(tǒng)一短步長解法節(jié)約計算時間．同時，采用該算法，先假定影響區(qū)與非影響區(qū)的分界面暫時固定，各個區(qū)按不同的時間步長分別計算，然后再放松分界面的約束，可實現(xiàn)應(yīng)力場的簡化計算，即水平分區(qū)異步長解法為復(fù)雜耗時的PSPS攔河大壩蓄水后正常運行期溫度應(yīng)力場的計算提供了可能．

(5) 由以上分析可知，本文提出的考慮特殊邊界影響的PSPS大壩不穩(wěn)定溫度場水平分區(qū)異步長解法是合理有效．

4 嚴寒地區(qū)特殊問題初探

如果外界氣溫過低，殘余水分長時間暴露在空氣中會結(jié)冰，水位回升后，冰層浸泡在水中又會解凍．比如第3節(jié)中位于高寒地區(qū)的呼和浩特PSPS，冬季長達5個月，11月、12月、1月、2月的月平均氣溫均低于0,℃，1月則達到-12.2,℃，因此其特殊邊界迎水面存在結(jié)冰與解凍的問題．

由于國內(nèi)同類型工程較少，故目前對迎水面結(jié)冰與解凍問題認識不足．位于俄羅斯莫斯科的扎戈爾PSPS同樣地處高寒地區(qū)，它的研究結(jié)論可為嚴寒地區(qū)PSPS混凝土表面結(jié)冰與解凍的問題提供參考．

扎戈爾PSPS[10-11]位于俄羅斯莫斯科東北角的謝爾捷耶夫-波薩德斯克地區(qū)的庫尼亞河上，位于高寒地區(qū)．扎戈爾抽水蓄能電站上水庫護坡混凝土板工作溫度狀態(tài)[12]如表4所示．當水庫水位下降時，脫離水面而暴露在空氣中的護坡混凝土，在－10～0,℃的空氣中暴露時間超過3.5,h就會結(jié)冰，而要將該冰層融化則需要在水中浸泡4,h，循環(huán)時間不小于7.5,h．

類比分析作為本文工程實例的呼和浩特PSPS，考慮其壩址處氣溫分布規(guī)律及水庫運行特性，可得出以下類比結(jié)論：①呼和浩特PSPS下水庫在11月、12月、1月、2月4個月份運行時，迎水面在空氣中暴露不少于3.5,h就會結(jié)冰；②一個周期內(nèi)，水、氣交替區(qū)至少浸水4,h，迎水面冰層會融化．

表4 迎水面混凝土結(jié)冰與解凍規(guī)律Tab.4Rule of freezing and thawing of upstream boundary concrete

另外，扎戈爾PSPS上水庫護坡混凝土板溫度狀態(tài)的研究表明，冰下的混凝土溫度在-2～0,℃之間，如圖6所示．結(jié)冰對混凝土壩面是有利的，因為冰層起到了一定的保溫效果．為了簡單起見，結(jié)冰后的邊界溫度視為0,℃．

圖6 結(jié)冰后混凝土表面溫度分布(12月)Fig.6Temperature distribution of concrete surface after freezing (December)

現(xiàn)在來回答與第1節(jié)相同的兩個問題．

1) 點A為水邊界或氣邊界的條件

對于嚴寒地區(qū)冬季月平均氣溫小于0,℃的月份，按照以下原則近似確定邊界條件：①若點A在下水庫水位以上，在0～3.5,h內(nèi)，附著有殘余水分，仍視為水邊界；在3.5,h以上，點A殘余水分已結(jié)冰，則點A邊界溫度按0,℃處理．②若點A在下水庫水位以下，在0～4,h內(nèi)，附著的殘余水分結(jié)冰沒有解凍，則點A邊界溫度仍按0,℃處理；在4,h以上，結(jié)冰已解凍，則視為水邊界．

2) 在一個運行周期里，點A為水邊界或氣邊界的時段的分配

在抽水蓄能電站一個運行周期1,d里，令點A為水邊界的時段為1t，與冰接觸(稱之為冰邊界)時段為．假設(shè)在下水庫抽水及蓄水過程中，水位均勻變動．水邊界和冰邊界時段分布如圖7所示，則有

寒冷季節(jié)冰邊界，因為考慮到冰層的保護作用大壩混凝土表面溫度近似取0,℃，相當于忽略混凝土與冰層之間的熱阻，與冰接觸的放熱系數(shù)可取代入式(49)得到等效放熱系數(shù)

仍以呼和浩特PSPS下水庫攔河壩為例，驗證水平分區(qū)異步長算法對于寒冷地區(qū)結(jié)冰與解凍特殊問題的適用性．

設(shè)計算例2 計算時間取為運行期冬季12月，其他條件與算例1相同．算例2兩種解法壩體不穩(wěn)定溫度場分布對比如圖8所示．

分析相關(guān)計算成果，可得到以下結(jié)論．

(1) 總體上來說，對此嚴寒地區(qū)寒冷季節(jié)的大壩采用水平分區(qū)異步長解法，與整個壩體統(tǒng)一短步長解法在同一時刻大壩不穩(wěn)定溫度場分布十分接近，非影響區(qū)域R2溫度分布吻合較好，影響區(qū)域R1則在上游邊界和分區(qū)界限位置處由于V場的簡化處理存在誤差，但溫度差值不超過0.05,℃．

(2) 呼和浩特PSPS下水庫庫容小，進出水口相對當?shù)貛斓赘叱虄H為5,m，電站裝機容量、出入庫流量大，因此認為電站運行方式引起的水體摻混作用足以使整個水域溫度分布均勻，不存在豎直高程上的梯度，因此12月整個水域均取為0,℃[9]，加上結(jié)冰后邊界溫度也視為0,℃，因此U場溫度分布梯度小，較為均勻，且U場占T場比例不大．對于其他庫容大、進出水口高程高的嚴寒地區(qū)PSPS水庫，水溫在庫表和庫底高程之間存在分布梯度，此時U場分布不再均勻，具有算例1中“從特殊邊界表面到分區(qū)界限處，U場越來越小，而V場越來越大，特殊邊界附近U場占T場的比例較大”的特點．

(3) 兩種解法計算結(jié)果對比表明，基于對結(jié)冰與解凍問題的這種近似處理原則，本文提出的考慮特殊邊界影響的PSPS大壩不穩(wěn)定溫度場水平分區(qū)異步長解法對于寒冷地區(qū)結(jié)冰與解凍特殊問題也適用．

圖8 算例2兩種解法壩體不穩(wěn)定溫度場分布對比(單位：℃)Fig.8Comparison of unsteady temperature filed of the dam under two algorithms in example 2(unit：℃)

5 結(jié) 論

(1) 本文提出了考慮特殊邊界影響的PSPS大壩不穩(wěn)定溫度場水平分區(qū)異步長解法，并在不穩(wěn)定溫度場分解過程中提出了等效放熱系數(shù)的概念，并通過理論計算，確定分區(qū)界限位置范圍大致為特殊邊界表面附近2～3,m．

(2) 對呼和浩特PSPS的工程實例進行了分析，將這種水平分區(qū)異步長解法的計算結(jié)果與整個壩體統(tǒng)一短步長解法計算結(jié)果進行了比較，二者的影響區(qū)域和非影響區(qū)域的不穩(wěn)定溫度場分布均十分接近，證明本文提出的算法是合理有效的．

(3) 水平分區(qū)異步長解法比整個壩體統(tǒng)一短步長解法節(jié)約了計算時間．同時，也為具有特殊邊界的PSPS攔河大壩蓄水后正常運行期溫度應(yīng)力場的計算提供了可能．具有特殊邊界的PSPS攔河大壩蓄水后正常運行期溫度應(yīng)力場的計算今后有待研究．

(4) 水平分區(qū)異步長解法在建模與求解過程中需要注意：分區(qū)界限位置既需要足夠大，以滿足U場和V場分開計算的要求，又希望盡量小，以節(jié)省計算資源，文中經(jīng)過理論推導(dǎo)給出了大致的范圍，可供參考；為了提高精度，需將特殊邊界按澆筑層細分，根據(jù)每一層的水、氣接觸時間，求得每一層的等效放熱系數(shù)值；總溫度場分解后得到的U場和V場，初始條件和邊界條件有很大的不同，求解過程需要加以注意；疊加過程可以通過外部程序調(diào)用有限元軟件自動連續(xù)的加以實現(xiàn)；初步計算時可與文中算例一樣，選用幾天，采用本文解法與整個壩體統(tǒng)一短步長解法進行比較驗證，以檢查或調(diào)整計算模型．

(5) 本文在處理特殊邊界與水、氣交替接觸時，對混凝土表面水分傳輸以及寒冷季節(jié)迎水面結(jié)冰與解凍問題進行了一定程度上的近似處理．筆者認為需要在以后有關(guān)實驗或者工程經(jīng)驗更加豐富時，對此近似處理進行驗證或改進(如考慮表層水分傳輸機制、進行濕度場和溫度場耦合計算等)，使得仿真條件更加接近真實情況．

［1］ 陸佑楣，潘家錚. 抽水蓄能電站[M]. 北京：水利電力出版社，1992. Lu Youmei，Pan Jiazheng. Pumped Storage Power Station[M]. Beijing：Water Resources and Electric Power Press，1992(in Chinese).

［2］ 朱伯芳. 大體積混凝土溫度應(yīng)力與溫度控制[M]. 北京：中國電力出版社，1999. Zhu Bofang. Thermal Stresses and Temperature Control of Mass Concrete[M]. Beijing：China Electric Power Press，1999(in Chinese).

［3］ 郭之章，傅 華. 水工建筑物的溫度控制[M]. 北京：水利電力出版社，1990. Guo Zhizhang，F(xiàn)u Hua. Temperature Control of Hydraulic Structure[M]. Beijing：Water Resources and Electric Power Press，1990(in Chinese).

［4］ 李松輝，張湘濤，張國新，等. 高混凝土重力壩關(guān)鍵部位溫控防裂研究[J]. 水力發(fā)電學(xué)報，2013，32(3)：182-186. Li Songhui，Zhang Xiangtao，Zhang Guoxin，et al. study on temperature control for key components of high concrete dam[J]. Journal of Hydroelectric Engineering，2013，32(3)：182-186(in Chinese).

［5］ 李守義，趙基花. 碾壓混凝土重力壩溫度場與溫度徐變應(yīng)力仿真分析[J]. 西安理工大學(xué)學(xué)報，2004，20(1)：58-62. Li Shouyi，Zhao Jihua. Simulation analysis of temperature field and thermal stress of RCC gravity dam[J]. Journal of Xi’an University of Technology，2004，20(1)：58-62(in Chinese).

［6］ 朱伯芳. 不穩(wěn)定溫度場數(shù)值分析的分區(qū)異步長解法[J]. 水利學(xué)報，1995(8)：46-52. Zhu Bofang. A method using different time increments in different regions for solving unsteady temperature filed by numerical method[J]. Journal of Hydraulic Engineering，1995(8)：46-52(in Chinese).

［7］ Bazant Z P，Najjar L J. Nonlinear water diffusion in nonsaturated concrete[J]. Materials and Structures/Materiauxet Constructions，1972，5(25)：3-20.

［8］ 李春秋. 干濕交替下表層混凝土中水分與離子傳輸過程研究[D]. 北京：清華大學(xué)土木水利學(xué)院，2009. Li Chunqiu. Study on Water and Ionic Transport Processes in Cover Concrete Under Drying-Wetting Cycles[D]. Beijing：School of Civil Engineering，Tsinghua University，2009(in Chinese).

［9］ 孫 偉，何蘊龍，趙 鑫，等. 抽水蓄能電站水庫水溫簡化模型與應(yīng)用[J]. 武漢大學(xué)學(xué)報：工學(xué)版，2013，46(4)：449-457. Sun Wei，He Yunlong，Zhao Xin，et al. A simplified method of reservoir water temperature for pumped storage power station and its application[J]. Engineering Journal of Wuhan University，2013，46(4)：449-457(in Chinese).

［10］ 庫列紹夫A II. 扎戈爾抽水蓄能電站的問題及解決辦法[J]. 劉正啟，譯. 水利水電快報，2001，22(12)：14-15. Kuleshov A II. Problems and solutions in Zhageer pumped storage power station[J]. Liu Zhengqi，Trans. Water Conservancy and Hydropower Express，2001，22(12)：14-15(in Chinese).

［11］ 庫德林B II. 扎戈爾抽水蓄能電站專項工程施工[J].趙秋云，譯. 水利水電快報，2000，21(15)：10-14. Kudrin B II. The special engineering construction in Zhageer pumped storage power station[J]. Zhao Qiuyun，Trans. Water Conservancy and Hydropower Express，2000，21(15)：10-14(in Chinese).

［12］ 夏爾庫諾夫С В. 扎戈爾抽水蓄能電站上庫護坡混凝土板工作溫度狀態(tài)的研究[J]. 劉正啟，譯. 水利水電快報，2002，23(12)：17-19. Charles Kuno С В. Study on the working temperature state of concrete plate in the upper reservoir protection slope in Zhageer pumped storage power station[J]. Liu Zhengqi，Trans. Water Conservancy and Hydropower Express，2002，23(12)：17-19(in Chinese).

(責任編輯：樊素英)

A Method Using Different Time Steps in Different Horizontal Regions for Simulation of Unsteady Temperature Field of Pumped Storage Power Station Dam

Sun Wei1，He Yunlong1，Wang Nan1，Xiong Kun2
(1. State Key Laboratory of Water Resources and Hydropower Engineering Science，Wuhan University，Wuhan 430072，China；2. Changjiang Survey Planning Design and Research Limited Liability Company，Wuhan 430010，China)

Because of the cyclical fluctuation of the water level of the reservoir，the upstream boundary of the pumped storage power station(PSPS)dam would encounter alternate contact with water，gas and other different media. The upstream dam boundary with this kind of characteristics is called the special boundary. These alternating boundary conditions have a big impact on the distribution of the temperature field of special boundary surface and its impact range，and cannot be ignored in the temperature simulation analysis. This paper proposes a method of solving unsteady temperature field of pumped storage power station dam with different time steps in different horizontal regions，and analyzes in detail the characteristics of special boundary and gives the formula derivation of the basic principle of the method with different time steps in different horizontal regions. The concept of equivalent heat transfer coefficient is also put forward to averagely and approximately consider the variation of heat transfer coefficient when the surface concrete contacting with different media and the position of horizontal partitioning boundary is given as well. The method is used in the simulation analysis of unsteady temperature field of Hohhot PSPS reservoir dam. The results show that the proposed method is effective and provides the possibility of calculating the temperature stress field which is more complex and time-consuming. Finally，the paper tentatively discusses the problem of freezing andthawing at the special boundary in cold regions.

pumped storage power station；horizontal partition；different time steps；unsteady temperature field；equivalent heat transfer coefficient

TV697.2

A

0493-2137(2015)09-0817-10

10.11784/tdxbz201401032

2014-01-11；

2014-06-18.

中國博士后科學(xué)基金資助項目(2012M511594)；國家大壩安全工程技術(shù)研究中心開放基金資助項目(NDSKFJJ1103).

孫 偉（1990— ），男，博士研究生，sunweimoving@126.com.

何蘊龍，ylhe2002@aliyun.com.cn.