孫 偉，何蘊(yùn)龍，王 南，熊 堃

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

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

孫 偉1，何蘊(yùn)龍1，王 南1，熊 堃2

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

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

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

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

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

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

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

1 特殊邊界特征分析

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

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

1)點(diǎn)A1為水邊界或氣邊界的條件

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

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

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

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

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

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

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

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

2)在一個(gè)運(yùn)行周期里，點(diǎn)A1為水邊界或氣邊界的時(shí)段分配

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

圖3中，1M、2M為以點(diǎn)A1處的水位H做平行于時(shí)間軸的直線與水庫水位-時(shí)間曲線的兩個(gè)交點(diǎn)的時(shí)間坐標(biāo)，即點(diǎn)A1在0～1M和2M～24為水邊界，在

圖3 水邊界和氣邊界時(shí)段分配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)．蓄水后運(yùn)行期混凝土水化熱溫升的影響已不必考慮．圖4中，點(diǎn)A～E(影響區(qū)上游點(diǎn)A、影響區(qū)內(nèi)部點(diǎn)B、分區(qū)界限位置點(diǎn)C、非影響區(qū)內(nèi)部點(diǎn)D、非影響區(qū)下游點(diǎn)E)是為了便于分析而取出的典型點(diǎn)．

2.1.1 總溫度場(chǎng)——T場(chǎng)

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

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

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

式(4)的初始條件為

式(4)的邊界條件為

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

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

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

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

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

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

式(11)的初始條件為

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

2.1.2 分解溫度場(chǎng)—U場(chǎng)

分解的溫度場(chǎng)—U場(chǎng)應(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 分解溫度場(chǎng)—V場(chǎng)

分解的溫度場(chǎng)—V場(chǎng)應(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)方程都是線性的，上述分解是嚴(yán)格的．

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

熱傳導(dǎo)方程為

初始條件為

邊界條件為

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

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

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

在接觸面上有

表1 總溫度場(chǎng)的分解Tab.1 Decomposition of total temperature field

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

2.2 等效放熱系數(shù)

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

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

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

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

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

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

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

由式(47)有 e

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

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

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

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

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

熱傳導(dǎo)方程為

初始條件為

邊界條件為

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

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

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

對(duì)于這種周期短、變幅大的溫度變化，應(yīng)取較小的周期P．取溫度變化周期為1,d，t以1,h為單位步長，計(jì)算得到水溫邊界的影響深度如表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 工程實(shí)例分析

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

另外，扎戈?duì)朠SPS上水庫護(hù)坡混凝土板溫度狀態(tài)的研究表明，冰下的混凝土溫度在-2～0,℃之間，如圖6所示．結(jié)冰對(duì)混凝土壩面是有利的，因?yàn)楸鶎悠鸬搅艘欢ǖ谋匦Ч疄榱撕?jiǎn)單起見，結(jié)冰后的邊界溫度視為0,℃．

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

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

1) 點(diǎn)A為水邊界或氣邊界的條件

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

2) 在一個(gè)運(yùn)行周期里，點(diǎn)A為水邊界或氣邊界的時(shí)段的分配

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

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

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

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

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

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

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

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

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

5 結(jié) 論

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

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

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

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

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

［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é)報(bào)，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］ 李守義，趙基花. 碾壓混凝土重力壩溫度場(chǎng)與溫度徐變應(yīng)力仿真分析[J]. 西安理工大學(xué)學(xué)報(bào)，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)定溫度場(chǎng)數(shù)值分析的分區(qū)異步長解法[J]. 水利學(xué)報(bào)，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ùn)龍，趙 鑫，等. 抽水蓄能電站水庫水溫簡(jiǎn)化模型與應(yīng)用[J]. 武漢大學(xué)學(xué)報(bào)：工學(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. 扎戈?duì)柍樗钅茈娬镜膯栴}及解決辦法[J]. 劉正啟，譯. 水利水電快報(bào)，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. 扎戈?duì)柍樗钅茈娬緦ｍ?xiàng)工程施工[J].趙秋云，譯. 水利水電快報(bào)，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］ 夏爾庫諾夫С В. 扎戈?duì)柍樗钅茈娬旧蠋熳o(hù)坡混凝土板工作溫度狀態(tài)的研究[J]. 劉正啟，譯. 水利水電快報(bào)，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).

(責(zé)任編輯：樊素英)

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é)基金資助項(xiàng)目(2012M511594)；國家大壩安全工程技術(shù)研究中心開放基金資助項(xiàng)目(NDSKFJJ1103).

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

何蘊(yùn)龍，ylhe2002@aliyun.com.cn.