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

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

(1.清華大學水利水電工程系,北京 100084; 2.吉林省水利水電勘測設計研究院,吉林 長春 130021)

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

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

(1.清華大學水利水電工程系,北京 100084; 2.吉林省水利水電勘測設計研究院,吉林 長春 130021)

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

冰期輸水;調水工程;冰情演變;數值模擬; 伊丹河

河冰是寒冷地區尤其是高緯度區域冬季河道中普遍存在的自然現象。北方河道冰期輸水時,嚴寒的氣候往往出現冰害,導致運行和管理上的許多困難。因此研究冰期輸水所涉及的相關冰害問題對于河渠道冬季安全運行具有重要意義。

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

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

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

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

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

1 河冰數學模型及計算方法

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

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

(1)

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

根據熱量守恒原理,沿河水溫η的時空分布近似由以下對流-擴散方程表述[4]:

(3)

其中η=ρCpATw

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

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

(4)

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

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

自豐滿水庫取水口至長春干線出水口斷面C09之間的埋地輸水管道水溫,采用以下一維水溫對流-擴散方程計算[14],得到長春干線出水口斷面水溫隨時間的變化過程:

(5)

2 伊丹河輸水河段數值模擬

2.1 河段基本概況

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

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

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

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

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

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

圖3 新立城入庫斷面冬季流量水位關系線

圖4 豐滿水庫有壓隧洞取水口處實測水溫隨時間變化

2.4 計算結果及分析

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

圖5 不同時間沿程水位變化

圖6 C09斷面冬季水位過程

圖7 入庫斷面冬季流量過程

圖8 入庫斷面水溫變化過程

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

圖10 冰蓋前緣位置隨時間變化

圖11 入庫斷面冰厚隨時間變化

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

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

3 結 論

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

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

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

[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] 茅澤育,吳劍疆,張磊,等.天然河道冰塞演變發展的數值模擬[J].水科學進展,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] 茅澤育,許昕,王愛民,等.基于適體坐標變換的二維河冰模型[J].水科學進展,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] 靳國厚,高霈生,呂斌秀.明渠冰情預報的數學模型[J].水利學報,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] 卞雪軍,冀鴻蘭,姜新華,等.黃河內蒙古段冰情預報系統的開發[J].水利水電科技進展,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] 王永填.河冰數值模型及河渠冬季輸水冰情特性研究[D].北京:清華大學,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] 樊霖,茅澤育,吳劍疆，等.松花江白山河段冰情數值計算及分析[J].水科學進展,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] 崔慧,吳長春.熱油管道非穩態工況傳熱與流動的耦合計算模型[J].中國石油大學學報(自然科學版),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

國家重點研發計劃(2016YFC0402504);清華大學博士生短期出國訪學基金(2013118)

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

茅澤育(1962—),男,教授,博士,主要從事水力學及河流動力學研究。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 編輯：鄭孝宇)