伊丹河輸水河道封凍期冰情演變數(shù)值模擬

樊 霖,茅澤育,齊文彪,劉樹峰

(1.清華大學(xué)水利水電工程系,北京 100084; 2.吉林省水利水電勘測(cè)設(shè)計(jì)研究院,吉林 長(zhǎng)春 130021)

伊丹河輸水河道封凍期冰情演變數(shù)值模擬

樊 霖1,茅澤育1,齊文彪2,劉樹峰2

(1.清華大學(xué)水利水電工程系,北京 100084; 2.吉林省水利水電勘測(cè)設(shè)計(jì)研究院,吉林 長(zhǎng)春 130021)

根據(jù)水動(dòng)力學(xué)、熱力學(xué)、河冰水力學(xué)及固體力學(xué)等基本理論,針對(duì)吉林省伊丹河輸水河段具體特征,建立了伊丹河的河冰數(shù)值模型,并應(yīng)用1958—1988共30年的水文、氣象資料對(duì)該輸水河段進(jìn)行了數(shù)值模擬及分析。結(jié)果表明,伊丹河輸水河段冰期封凍形式為平封,可形成穩(wěn)定冰蓋從而實(shí)現(xiàn)冰期冰蓋下輸水,且滿足冰期正常輸水的要求。

冰期輸水;調(diào)水工程;冰情演變;數(shù)值模擬; 伊丹河

河冰是寒冷地區(qū)尤其是高緯度區(qū)域冬季河道中普遍存在的自然現(xiàn)象。北方河道冰期輸水時(shí),嚴(yán)寒的氣候往往出現(xiàn)冰害,導(dǎo)致運(yùn)行和管理上的許多困難。因此研究冰期輸水所涉及的相關(guān)冰害問題對(duì)于河渠道冬季安全運(yùn)行具有重要意義。

河冰生消演變?yōu)樗畡?dòng)力、機(jī)械力和熱力等相互耦合作用且極為復(fù)雜的物理過程。多年來,隨著河冰水力學(xué)理論的發(fā)展,國(guó)內(nèi)外學(xué)者們相繼提出了一些模擬河道冰情的數(shù)值模型。例如Lal等[1]根據(jù)熱力和兩層輸冰模型的理論分析框架,建立了冰情演變河冰數(shù)值模型—RICE模型。Shen等[2]提出了二維與一維結(jié)合的河冰動(dòng)力學(xué)模型DynaRICE,主要用于模擬冰壩形成的動(dòng)力學(xué)過程等。Zufelt等[3]建立了一維冰水耦合的冰塞動(dòng)力學(xué)模型,模擬封河期和開河期水位變化。茅澤育等[4-5]針對(duì)天然河道彎曲復(fù)雜的特點(diǎn),首次建立了適體坐標(biāo)下的二維河冰數(shù)值模型,對(duì)黃河河曲段、黃河萬(wàn)家寨水庫(kù)庫(kù)區(qū)河段、松花江上游白山河段、新疆額爾齊斯河段等進(jìn)行了數(shù)值模擬研究,取得了較好的效果。Fu等[6]在RICE和RICEN模型的基礎(chǔ)上,開發(fā)了用于模擬冰期洪水的YRIDM模型。靳國(guó)厚等[7]采用一維非恒定流水力學(xué)及熱力學(xué)相耦合的模型,預(yù)報(bào)了京密引水渠、引黃濟(jì)青等輸水工程的冰期輸水過程。卞雪軍等[8]以多元線性回歸理論與人工神經(jīng)網(wǎng)絡(luò)理論為基礎(chǔ),設(shè)計(jì)開發(fā)了黃河內(nèi)蒙古河段冰情預(yù)報(bào)系統(tǒng)。

如上所述,河流冰情演變受水力、氣象、熱力、河勢(shì)及河道特征等眾多因素的綜合影響,這些影響因素的不同組合必將導(dǎo)致發(fā)生不同的冰情。所以,對(duì)于任一具體河道,需根據(jù)具體的氣候、地形、水力等條件建立模型、率定參數(shù)并進(jìn)行數(shù)值計(jì)算。

伊丹河全流域面積為483 km2,其河源至新立城水庫(kù)入庫(kù)斷面全長(zhǎng)47.71 km(圖1)。伊丹河為伊通河的最大支流,沿途流經(jīng)伊通縣二道鎮(zhèn)、伊通縣伊丹鎮(zhèn)、長(zhǎng)春市南關(guān)區(qū)新湖鎮(zhèn)等,在長(zhǎng)春市南關(guān)區(qū)新湖鎮(zhèn)榆樹屯注入伊通河。

圖1 伊丹河輸水河段示意圖

伊丹河冰期輸水工程為吉林省中部城市群引松供水長(zhǎng)春干線工程的重要組成部分。長(zhǎng)春干線全長(zhǎng)14.1 km,設(shè)計(jì)流量13.2 m3/s,由豐滿水庫(kù)經(jīng)110 km有壓隧洞(隧洞取水口位于水庫(kù)正常蓄水位以下21 m),經(jīng)馮家?guī)X分水樞紐進(jìn)入長(zhǎng)春干線,再在C09斷面處泄入伊丹河,并自流入新立城水庫(kù)。匯入斷面C09至新立城水庫(kù)入庫(kù)斷面河段長(zhǎng)度為7.3 m。由于工程地處嚴(yán)寒地區(qū),河道輸水一旦形成冰害,將抬高長(zhǎng)春干線末端即出口斷面C09的水位,減小長(zhǎng)春干線的輸水作用水頭,影響長(zhǎng)春干線設(shè)計(jì)輸水量。由于長(zhǎng)春干線與馮家?guī)X分水樞紐連通,還將影響四平泵站、遼源泵站設(shè)計(jì)參數(shù),并對(duì)馮家?guī)X調(diào)壓井結(jié)構(gòu)產(chǎn)生不利影響,進(jìn)而影響整個(gè)總干線的設(shè)計(jì)參數(shù)。此外,長(zhǎng)春干線伊丹河下游末端為水庫(kù),極易形成冰害,一旦形成冰害,對(duì)河岸、堤防產(chǎn)生影響,影響輸水穩(wěn)定性。因此,實(shí)現(xiàn)長(zhǎng)春干線伊丹河段冰期安全輸水對(duì)于長(zhǎng)春干線、馮家?guī)X分水樞紐乃至整個(gè)中部供水工程正常運(yùn)行至關(guān)重要。本文根據(jù)伊丹河河道的地形、氣象及研究河段實(shí)際情況,建立河冰數(shù)值模型,并應(yīng)用伊丹河1958—1988年的水文、氣象資料對(duì)伊丹河冰期輸水進(jìn)行了冰情數(shù)值模擬。

1 河冰數(shù)學(xué)模型及計(jì)算方法

封凍期河冰演變的模型主要包括河道水力學(xué)模型、輸冰模型、熱力學(xué)模型及冰凍模型等。

封凍河道的水流一維流動(dòng)近似由以下連續(xù)性方程和動(dòng)量方程描述[4]:

(1)

式中:Q為流量;A為過流斷面面積;ql為單寬側(cè)向入流量;x為距離;t為時(shí)間;g為重力加速度;z為水位;zb為河床高程;ti為冰蓋厚度;ρ、ρi分別為水和冰的密度;h為冰蓋下水深;nc為綜合糙率系數(shù);P為濕周。

根據(jù)熱量守恒原理,沿河水溫η的時(shí)空分布近似由以下對(duì)流-擴(kuò)散方程表述[4]:

(3)

其中η=ρCpATw

式中:v為斷面平均流速;B0為水面寬;∑S為水體與周圍環(huán)境的單位面積熱交換量(增熱或失熱);Cp為水的比熱容,取4 185 J/(kg·℃);Tw為斷面平均水溫。

根據(jù)分層輸冰理論,河流中的流冰由面冰和懸浮冰兩部分組成,其質(zhì)量守恒方程可參見文獻(xiàn)[9]。根據(jù)冰蓋前緣斷面的水動(dòng)力條件,以及上游來冰情況,冰蓋的發(fā)展一般有并置推進(jìn)、水力增厚推進(jìn)、機(jī)械增厚推進(jìn)[1,10]。冰蓋下冰的輸移和積聚采用文獻(xiàn)[11]中的輸冰能力公式計(jì)算。當(dāng)積聚冰蓋體向上游發(fā)展過程中,一方面,冰蓋體中的孔隙水凍結(jié)而形成固狀冰體外殼,增大冰體強(qiáng)度;另一方面,由于冰蓋表面和底面發(fā)生熱交換,冰體厚度將發(fā)生熱力增厚(或消融)。冰蓋體的熱力增厚或消融過程由熱力平衡原理得到,即:

(4)

式中:ti為水溫度；Ts為冰層表面溫度;Ta為氣溫;Tm為冰層熱融度,即0℃;hia為大氣與冰層表面的熱交換系數(shù),取19.71;hwi為水體與冰層底面的熱交換系數(shù)。

模型中式(1)(2)(3)采用顯式有限差分MacCormack步進(jìn)格式求解[12],其在時(shí)間和空間上具有二階精度o(Δx2,Δt2)。MacCormack方法屬兩步預(yù)測(cè)-修正型,其預(yù)測(cè)和修正步驟,可以交替地向前和向后差分。各子模型組成的聯(lián)合方程組采用同步求解。為準(zhǔn)確模擬伊丹河輸水河段的冰情,筆者曾應(yīng)用白山河段完整詳盡的氣象、冰情原型觀測(cè)資料,采用工程類比驗(yàn)證方法對(duì)上述數(shù)學(xué)模型進(jìn)行了模型參數(shù)率定及驗(yàn)證[13]。

自豐滿水庫(kù)取水口至長(zhǎng)春干線出水口斷面C09之間的埋地輸水管道水溫,采用以下一維水溫對(duì)流-擴(kuò)散方程計(jì)算[14],得到長(zhǎng)春干線出水口斷面水溫隨時(shí)間的變化過程:

(5)

2 伊丹河輸水河段數(shù)值模擬

2.1 河段基本概況

對(duì)伊丹河(全長(zhǎng)47.7 km)上游河源、下游新立城水庫(kù)及豐滿水庫(kù)有壓隧洞取水口作全系統(tǒng)研究。長(zhǎng)春干線全長(zhǎng)14.1 km(PCCP管道,平均埋深6 m),匯入斷面C09;豐滿水庫(kù)至長(zhǎng)春干線的有壓隧洞全長(zhǎng)110 km。計(jì)算域共布置斷面154個(gè)，平均斷面間距約300 m。斷面1布置在伊丹河河源斷面(C09斷面),編號(hào)為84，距河源40.40 km；新立城水庫(kù)入庫(kù)斷面編號(hào)為154，距河源47.7 km；斷面127位于河源下游44.93 km。河段平均坡降3.2×10-4。

2.2 氣象水文資料選取與整理

采用新立城壩上水文站、伊通氣象站1958—1988年水文、氣象多年平均觀測(cè)數(shù)據(jù)。伊通氣象站冬季(11月至次年3月)多年日均氣溫如圖2所示(橫坐標(biāo)時(shí)間以11月1日為起點(diǎn)(下同))。計(jì)算總時(shí)長(zhǎng)151 d,時(shí)間步長(zhǎng)60 s,對(duì)于實(shí)測(cè)斷面間距過大河段進(jìn)行內(nèi)插,最大計(jì)算斷面間距取50 m。各實(shí)測(cè)橫斷面的過水面積與濕周通過微斷面法[15]進(jìn)行計(jì)算。

圖2 伊通氣象站冬季多年平均日均氣溫

2.3 初始條件及邊界條件確定

初始條件為計(jì)算開始時(shí)刻計(jì)算域內(nèi)流量、水位、水溫、浮冰濃度、水內(nèi)冰濃度和冰蓋體厚度。伊丹河河道初始流量及水位采用上游邊界及下游邊界對(duì)應(yīng)的流量和水位值,通過長(zhǎng)時(shí)間迭代運(yùn)算直至穩(wěn)定,最終得到沿程流量和水位值;初始水體溫度采用6℃。根據(jù)水文、氣象資料,初始時(shí)刻沿程浮冰濃度、水內(nèi)冰濃度分布及冰蓋體厚度均設(shè)定為零。流場(chǎng)計(jì)算邊界條件：上游邊界為伊丹河河源斷面處的流量過程以及中途匯入的長(zhǎng)春干線設(shè)計(jì)引水流量(13.2 m3/s),下游邊界為入庫(kù)斷面(斷面154)流量水位關(guān)系曲線(圖3)。水流溫度場(chǎng)計(jì)算的邊界條件為伊丹河河源斷面處以及豐滿水庫(kù)取水口處實(shí)測(cè)水溫值(圖4),初始條件通過模型迭代計(jì)算,直至初始溫度場(chǎng)基本達(dá)到熱交換的準(zhǔn)平衡狀態(tài)得到。

圖3 新立城入庫(kù)斷面冬季流量水位關(guān)系線

圖4 豐滿水庫(kù)有壓隧洞取水口處實(shí)測(cè)水溫隨時(shí)間變化

2.4 計(jì)算結(jié)果及分析

圖5為封河過程中不同時(shí)刻輸水河段沿程水位計(jì)算結(jié)果，圖中橫坐標(biāo)以C09斷面為起點(diǎn)(下同)。隨著氣溫降低,冰蓋前緣逐漸向上游推進(jìn),前緣斷面附近水位隨之壅高,同時(shí)冰體厚度增長(zhǎng)。1月29日水位壅高達(dá)到最大;3月下旬隨著氣溫回暖、水溫升高,冰體開始消融,沿程水位降低。圖6給出的是C09斷面封凍期(11月至次年3月)水位變化過程。由圖6可知,12月末1月初,由于下游河段冰蓋的生成,C09斷面水位開始壅高,最高壅水達(dá)0.16 m。隨著冰蓋進(jìn)一步固結(jié),冰蓋糙率系數(shù)降低,計(jì)算壅水位逐漸降低。入庫(kù)斷面封凍期流量變化過程如圖7所示,12月初,入庫(kù)斷面流量有所減小,這是由于河道槽蓄水量增加所致,隨后入庫(kù)斷面流量穩(wěn)定在13.4 m3/s。入庫(kù)斷面水溫隨時(shí)間的變化過程如圖8所示。由于受長(zhǎng)春干線輸水水體水溫的影響,伊丹河輸水河段的水溫變化滯后于日均氣溫變化。從圖8可知,12月15日受大氣溫度下降及河道水體與大氣熱交換的影響,水溫持續(xù)降至0℃以下,水體進(jìn)入過冷卻狀態(tài),隨著冰層生成,水溫又逐漸穩(wěn)定在0℃附近。

圖5 不同時(shí)間沿程水位變化

圖6 C09斷面冬季水位過程

圖7 入庫(kù)斷面冬季流量過程

圖8 入庫(kù)斷面水溫變化過程

圖9 1月29日河道弗勞德數(shù)沿程分布

圖10 冰蓋前緣位置隨時(shí)間變化

圖11 入庫(kù)斷面冰厚隨時(shí)間變化

圖12 1月29日冰厚沿程分布

1月29日輸水河段沿程弗勞德數(shù)如圖9所示。由圖9可知,輸水河段的弗勞德數(shù)均在第二臨界弗勞德數(shù)(0.09)以下,表明在氣象及水溫等條件下,河道大部分河段都可形成冰蓋。圖10給出的是冰蓋推進(jìn)前緣位置隨時(shí)間的變化，可以看出,12月30日下游入庫(kù)斷面首先形成冰蓋,6 d后冰蓋前緣到達(dá)127斷面,由于水溫的影響,冰蓋前緣到達(dá)127斷面之后便不在向前發(fā)展。圖11給出的是封凍期入庫(kù)斷面冰蓋平均厚度隨時(shí)間的變化，可見入庫(kù)斷面的初封時(shí)間為12月30日。隨后,冰蓋體在冰面與大氣及冰底面與水體熱交換的共同作用下,開始熱力增厚。2月20日入庫(kù)斷面冰厚達(dá)最大,為0.86 m,此后受氣溫回暖影響,冰厚開始減小。圖12為伊丹河輸水河段1月29日冰蓋體厚度沿程分布。

3 結(jié) 論

a. 伊丹河輸水河段冰期沿程弗勞德數(shù)均在第二臨界弗勞德數(shù)(0.09)以下,在氣象以及水溫條件等適宜條件下,河道大部分河段都可形成冰蓋,進(jìn)而實(shí)現(xiàn)河道冰期冰蓋下輸水,輸水河段沿程弗勞德數(shù)大多小于0.06,河段封凍形式為平封形式。

b. 河段封凍始于入庫(kù)斷面,然后逐漸上溯,6 d之后冰蓋前緣到達(dá)127斷面,冰期最大冰厚為0.86 m。

c. 在封凍過程中,冰蓋前緣推進(jìn)速度及冰體厚度分布受諸多因素影響,且隨冰蓋前緣逐步向上游推進(jìn),沿程水位也逐漸壅高,最高壅水達(dá)0.16 m。冰期入庫(kù)斷面流量值在冰蓋形成期間略有下降(河道槽蓄水量增加所致),隨后一直維持在13.4 m3/s附近,伊丹河輸水河段能夠滿足冰期設(shè)計(jì)流量下的正常輸水要求。

[1] LAL A M,SHEN H T.A mathematical model for river ice processes[J].Journal of Hydraulic Engineering,1993,117(7): 851-867.

[2] SHEN H T,SU J,LIU L.SPH simulation of river ice dynamics[J].Journal of Computational Physics,2000,165(2): 752-770.

[3] ZUFELT J E,ETTEMA R.Fully coupled model of ice-jam dynamics[J].Journal of Cold Regions Engineering,2000,14(1): 24-41.

[4] 茅澤育,吳劍疆,張磊,等.天然河道冰塞演變發(fā)展的數(shù)值模擬[J].水科學(xué)進(jìn)展,2003,14(6): 700-705.(MAO Zeyu,WU Jianjiang,ZHANG Lei,et al.Numerical simulation of river ice jam[J].Advance in Water Science,2003,14(6): 700-705.(in Chinese))

[5] 茅澤育,許昕,王愛民,等.基于適體坐標(biāo)變換的二維河冰模型[J].水科學(xué)進(jìn)展,2008,19(2):214-223.(MAO Zeyu,XU Xin,WANG Aimin,et al.2-D numerical model for river-ice processes based upon body-fitted coordinate[J].Advance in Water Science,2008,19(2):700-705.(in Chinese))

[6] FU C,POPESCU I,WANG C,et al.Challenges in modeling ice floods on the Ningxia-Inner Mongolia reach of the Yellow River,China[J].Hydrology and Earth System Sciences Discussions,2013,10(10): 12293-12329.

[7] 靳國(guó)厚,高霈生,呂斌秀.明渠冰情預(yù)報(bào)的數(shù)學(xué)模型[J].水利學(xué)報(bào),1997,28(10):1-9.(JIN Guohou,GAO Peisheng,LU Binxiu.Ice regime forecast model in open channel[J].Journal of Hydraulic Engineering,1997,28(10):1-9.(in Chinese))

[8] 卞雪軍,冀鴻蘭,姜新華,等.黃河內(nèi)蒙古段冰情預(yù)報(bào)系統(tǒng)的開發(fā)[J].水利水電科技進(jìn)展,2014,34(4):62-65.(BIAN Xuejun,JI Honglan,JIANG Xinhua,et al.Development and application analyze of forecasting system of ice conditions in Inner Mongolia reach of Yellow River[J].Advances in Science and Technology of Water Resources,2014,34(4):62-65.(in Chinese))

[9] 王永填.河冰數(shù)值模型及河渠冬季輸水冰情特性研究[D].北京:清華大學(xué),1999.

[10] PARISET E,HAUSSER R.Formation and evolution of ice covers on rivers[J].Transaction of Engineering Institute of Canada,1961,5(1): 41-49.

[11] SHEN H T,WANG D S,LAL A M.Numerical simulation of river ice processes[J].Journal of Cold Regions Engineering,1995,9(3): 107-118.

[12] MACCORMACK R W.A numerical method for solving the equations of compressible viscous flow[J].AIAA Journal,1982,20(9): 1275-1281.

[13] 樊霖,茅澤育,吳劍疆，等.松花江白山河段冰情數(shù)值計(jì)算及分析[J].水科學(xué)進(jìn)展,2016,27(6):890-897.(FAN Lin,MAO Zeyu,WU Jianjiang,et al.Study on the numerical model of river ice for the Baishan reach of Songhua River,China[J].Advance in Water Science,2016,27(6):890-897.(in Chinese))

[14] 崔慧,吳長(zhǎng)春.熱油管道非穩(wěn)態(tài)工況傳熱與流動(dòng)的耦合計(jì)算模型[J].中國(guó)石油大學(xué)學(xué)報(bào)(自然科學(xué)版),2005,29(3):101-105.(CUI Hui,WU Changchun.Heat transfer and flow coupling calculation model of transient scenario for hot oil pipeline[J].Journal of China University of Petroleum,2005,29(3):101-105.(in Chinese))

[15] FAN L,MAO Z Y,BAO J,et al.A river-ice model for Wanjiazhai reach of the Yellow River[C]//Proceedings of the 22nd IAHR International Symposium on Ice.Singapore:NUS,2014:869-876.

Numericalsimulationofriver-iceprocessesinawater-transferriverofYidanRiverduringfreezingperiod//

FAN Lin1, MAO Zeyu1, QI Wenbiao2, LIU Shufeng2

(1.DepartmentofHydraulicEngineering,TsinghuaUniversity,Beijing100084,China; 2.JilinInvestigationandDesignInstituteofWaterResourcesandHydropower,Changchun130021,China)

Based upon the theory of hydrodynamics, river mechanics, thermodynamics, river-ice hydraulics and mechanics of solids etc., a mathematic model of the river-ice processes for the Yidan reach is constructed. 30-year long-term series comprehensive field observation data, including geomorphy, meteorology, hydrology, thermo-dynamics and ice regime, are used to simulate the water transfer process of the Yidan River during freezing period. The results show that steady ice cover can be formed by juxtaposition on the Yidan River, which achieves the goal of water transfer beneath the ice cover. In addition, the design discharge of the Yidan River water transfer project can be guaranteed during freezing period.

water transfer during freezing period; water diversion project;river ice process; numerical simulation; Yidan River

國(guó)家重點(diǎn)研發(fā)計(jì)劃(2016YFC0402504);清華大學(xué)博士生短期出國(guó)訪學(xué)基金(2013118)

樊霖(1990—),男,博士研究生,主要從事河冰水力學(xué)研究。E-mail:l-fan11@mails.tsinghua.edu.cn

茅澤育(1962—),男,教授,博士,主要從事水力學(xué)及河流動(dòng)力學(xué)研究。E-mail:maozeyu@tsinghua.edu.cn

10.3880/j.issn.1006-7647.2017.06.003

TV143

A

1006-7647(2017)06-0014-05

2017-02-21 編輯：鄭孝宇)