劉 洋

（江西省水利規(guī)劃設(shè)計(jì)院 江西 南昌 330029）

1 概述

對(duì)于河流的泥沙問(wèn)題，目前采用較多的研究方法主要有物理模型試驗(yàn)和數(shù)值模擬兩類(lèi)。物理模型具有直觀、可重現(xiàn)歷史狀況、可進(jìn)行未來(lái)問(wèn)題預(yù)測(cè)等優(yōu)點(diǎn)，但其實(shí)驗(yàn)周期長(zhǎng)、費(fèi)用高，且受比尺及實(shí)驗(yàn)環(huán)境、場(chǎng)地的影響大。自二十世紀(jì)五十年代以來(lái)，隨著計(jì)算機(jī)的出現(xiàn)和現(xiàn)代高速計(jì)算機(jī)技術(shù)的快速發(fā)展，數(shù)學(xué)模型在水利工程的許多領(lǐng)域已得到廣泛應(yīng)用，而計(jì)算水力學(xué)理論的成熟又促進(jìn)了求解流體力學(xué)方程方法的發(fā)展，使流動(dòng)現(xiàn)象的數(shù)值模擬達(dá)到了更高的水平。數(shù)學(xué)模型以計(jì)算水力學(xué)、河流動(dòng)力學(xué)與環(huán)境水力學(xué)等學(xué)科為基礎(chǔ)，將現(xiàn)代計(jì)算機(jī)與數(shù)值解技術(shù)相結(jié)合對(duì)河流進(jìn)行模擬研究，與物理模型試驗(yàn)相比，具有研究費(fèi)用較低、周期短、運(yùn)用靈活、不存在比尺效應(yīng)和試驗(yàn)環(huán)境影響等顯著優(yōu)點(diǎn)，應(yīng)用數(shù)學(xué)模型研究河流泥沙問(wèn)題無(wú)疑是該領(lǐng)域研究的發(fā)展方向。

水沙數(shù)學(xué)模型按所模擬的水沙運(yùn)動(dòng)在空間上的變化情況可分為一維模型、二維模型和三維模型。到目前為止，一、二維水沙數(shù)學(xué)模型的研究已相對(duì)成熟，三維水沙數(shù)學(xué)模型的研究尚處于探索階段。一維水沙數(shù)學(xué)模型只能求解河道斷面各水沙要素沿河道長(zhǎng)度方向的平均值，無(wú)法回答河段內(nèi)部沿河寬方向沖淤變化的分布問(wèn)題，對(duì)于很多實(shí)際工程問(wèn)題比如水庫(kù)壩前泥沙的運(yùn)動(dòng)、河床變形分析和過(guò)機(jī)泥沙預(yù)測(cè)等是解決不好的，而平面二維水沙數(shù)學(xué)模型卻能較好地解決這些問(wèn)題，其在水利工程上的應(yīng)用非常廣泛。

2 數(shù)學(xué)模型的構(gòu)建

2.1 模型簡(jiǎn)介

丹麥水力研究所（DanishHydraulic Institute，簡(jiǎn)稱(chēng)DHI公司）研發(fā)的DHI系列模型軟件，是目前世界上領(lǐng)先的并經(jīng)實(shí)際工程驗(yàn)證最多的，被水資源研究人員廣泛認(rèn)同和應(yīng)用的優(yōu)秀模型軟件之一。其涉及范圍廣泛，從一維模型到三維模型，從單一的水動(dòng)力模擬到水環(huán)境與生態(tài)系統(tǒng)相結(jié)合的復(fù)雜模型。其具有界面友好、模擬計(jì)算快速準(zhǔn)確、后處理方便等優(yōu)點(diǎn)。

MIKE21是DHI水動(dòng)力學(xué)模型系統(tǒng)之一，是一種通用的二維數(shù)值模擬系統(tǒng)，它可以用來(lái)研究河流、湖泊、河口、海灣、海岸及海洋等各類(lèi)環(huán)境下的水流、波浪和泥沙問(wèn)題。它既可以模擬簡(jiǎn)單的二維非恒定流，也可以對(duì)密度變化、水下地形變化、潮汐變化和降水蒸發(fā)等復(fù)雜條件下的流動(dòng)進(jìn)行分析計(jì)算。本文基于MIKE21中的水動(dòng)力模型和輸泥模型構(gòu)建一個(gè)平面二維水沙數(shù)學(xué)模型。

2.2 基本方程

本文采用的數(shù)學(xué)模型的基本方程主要包括水動(dòng)力模型和泥沙模型控制方程兩部分。

水動(dòng)力模型的控制方程為：質(zhì)量守恒方程:

x方向動(dòng)量方程:

y方向動(dòng)量方程:

其中：h為水深，m；d為水位隨時(shí)間的變化量，m；ζ為表面水位,m；p和q分別為x與y方向的通量密度,m3/sm；u和v分別為沿水深的平均流速,m/s；C為謝才系數(shù),m1/2/s；f為風(fēng)摩擦系數(shù)；V、Vx和Vy為風(fēng)速及其在x與y方向的分量,m/s；Ω為科氏系數(shù)，s-1；Pa為大氣壓，kg/m/s2；ρw為水的密度，kg/m3；t為時(shí)間，s；τxx、τxy、τyy為有效切應(yīng)力分量，N/m2。

該方程的數(shù)值解法主要采用交替方向隱式技術(shù)（ADI）分別對(duì)質(zhì)量守恒與動(dòng)量方程進(jìn)行時(shí)空上的積分，根據(jù)雙掃描法求解每個(gè)方向及每個(gè)單獨(dú)的網(wǎng)格產(chǎn)生的方程矩陣。

泥沙模型的計(jì)算是基于水動(dòng)力計(jì)算的結(jié)果，其控制方程為:

其中：c為垂向平均含沙量，kg/m3；x與y為空間坐標(biāo)，m；u、v分別為x與y方向的垂向平均流速，m/s；h為水深，m；Dx、Dy分別為x與y方向的泥沙紊動(dòng)擴(kuò)散系數(shù)，m2/s；QL為水平單位面積源項(xiàng)流量，m3/s/m2；CL為源含沙量，kg/m3；S為沖淤項(xiàng)，kg/m3/s。

該方程為一對(duì)流擴(kuò)散方程，采用以QUICKEST格式為基礎(chǔ)的三階有限差分顯式格式（即ULTIMATE格式）求解。該求解格式已經(jīng)在許多涉及對(duì)流擴(kuò)散方程的紊流模型、環(huán)境模型以及其他類(lèi)似問(wèn)題中得到了很好的應(yīng)用。與其他數(shù)值解法相比，它具有諸多優(yōu)越性，尤其避免了與對(duì)流項(xiàng)中心差分相關(guān)的“蠕動(dòng)”的不穩(wěn)定性，同時(shí)也在一定程度上減小了一階迎風(fēng)格式的數(shù)值耗散。

圖1 河道淤積厚度分布

圖2 各斷面淤積厚度歷時(shí)變化曲線

3 數(shù)學(xué)模型在灞河橡膠壩庫(kù)區(qū)的應(yīng)用

3.1 模擬區(qū)域概況

灞河屬黃河流域，渭河水系，是渭河的一級(jí)支流，發(fā)源于秦嶺北麓藍(lán)田、渭南、華縣交界處的藍(lán)田縣灞源鄉(xiāng)箭峪嶺南九道溝，由南向北流，經(jīng)灞源后西行，到馮家灣出峪口。河流全長(zhǎng)104km，流域面積2581km2（含浐河760km2），河床平均比降6.0‰。

灞河流域地形南高北低，屬秦嶺土石山區(qū)，巖石裸露，土層較薄，植被良好，平時(shí)水質(zhì)清澈，洪水時(shí)挾帶少量泥沙。上游（藍(lán)田故京）平均比降9.0‰，洪水猛漲暴落，水流湍急；中游（故京至浐灞交匯口）河道為平原彎曲型河道，河床斷面形態(tài)為寬淺式，比降為2.35‰；下游（浐灞交匯口至灞河入渭口）屬于平原河道，河床寬淺，比降較緩，為1.58‰。

目前，在灞河中下游修建有大量橡膠壩工程。為了城市防洪和維持水面景觀，橡膠壩建成以后長(zhǎng)期立壩運(yùn)行，加之該處河道比降較小、水面寬淺，使得河道水體流動(dòng)性變差。水土流失問(wèn)題與上述問(wèn)題的綜合作用使得該處庫(kù)區(qū)泥沙淤積問(wèn)題非常突出，采用數(shù)值模擬的方法對(duì)該河段泥沙淤積問(wèn)題進(jìn)行研究不失為一種有效的途徑。

3.2 模擬計(jì)算結(jié)果分析

根據(jù)實(shí)測(cè)地形和水文資料，以灞河B號(hào)壩庫(kù)區(qū)為例，進(jìn)行數(shù)值模擬計(jì)算。模擬的河段長(zhǎng)約2.5km，平均比降為2.35‰，河床糙率為0.033。根據(jù)實(shí)測(cè)地形進(jìn)行非結(jié)構(gòu)網(wǎng)格劃分，共劃分1016個(gè)網(wǎng)格，該河段上游入口寬約356m，下游出口寬度采用的是B號(hào)壩的壩長(zhǎng)，為450m，立壩運(yùn)行，在該河段中心處有一河心島。

采用的水沙過(guò)程為枯水期資料，即1955年至1965年灞河逐日水沙過(guò)程。模擬總時(shí)長(zhǎng)為10年，時(shí)間步長(zhǎng)為50s，對(duì)模型參數(shù)進(jìn)行反復(fù)調(diào)整率定，進(jìn)行模擬計(jì)算后，最后達(dá)到?jīng)_淤平衡后的河床淤積厚度的分布如圖1所示。

從圖1中可以看出，隨著時(shí)間的推移，在該河段上游形成了一狹長(zhǎng)的主河槽，主河槽兩側(cè)的白色區(qū)域?yàn)橛俜e作用造成的露出水面的河灘地，水流主要沿著主河槽繞過(guò)中間河心島向下游流動(dòng)。主河槽處深色區(qū)域明顯大于兩側(cè)，說(shuō)明高濃度懸沙在枯水期橡膠壩庫(kù)區(qū)長(zhǎng)期立壩造成河道水體流動(dòng)性變差以后，更容易落淤，最大淤積厚度達(dá)到了3.5m以上，該河段在長(zhǎng)時(shí)間系列中的淤積作用還是相當(dāng)嚴(yán)重的。

提取模擬所得該河段各斷面各時(shí)間段淤積厚度數(shù)據(jù)，繪制淤積厚度歷時(shí)變化曲線，如圖2所示。從圖2中可以看出，河床的淤積厚度隨著時(shí)間的推移逐漸增大，最終于相應(yīng)的年限達(dá)到平衡狀態(tài)。達(dá)到?jīng)_淤平衡的時(shí)間表現(xiàn)為上游早下游晚，最終的淤積厚度主要分布在3m～3.5m之間。其中，灞4～灞27各斷面計(jì)算值達(dá)到平衡的年限分別為4年、4.5年、4.5年、4.5年、4.5年、5年、5年、6年，其平衡時(shí)的淤積厚度值分別為3.00m、3.11m、3.20m、3.32m、3.21m、3.03m、3.41m、3.42m，整個(gè)庫(kù)區(qū)的淤死年限為6年。

4 結(jié)語(yǔ)

本文首先對(duì)水沙數(shù)學(xué)模型的發(fā)展現(xiàn)狀進(jìn)行了簡(jiǎn)單介紹，然后闡述了所采用的數(shù)學(xué)模型的控制方程和其數(shù)值方法。應(yīng)用所構(gòu)建的數(shù)學(xué)模型在灞河B號(hào)橡膠壩庫(kù)區(qū)河段進(jìn)行了模擬計(jì)算，通過(guò)分析河道的沖淤演變規(guī)律發(fā)現(xiàn)，在枯水期橡膠壩庫(kù)區(qū)長(zhǎng)期立壩情況下，該河段淤積作用是相當(dāng)嚴(yán)重的。利用本文所構(gòu)建的平面二維水沙數(shù)學(xué)模型可以對(duì)河道的沖淤演變趨勢(shì)進(jìn)行預(yù)測(cè)分析，得出河道的演變規(guī)律。可見(jiàn)，采用數(shù)值模擬方法來(lái)研究并處理灞河橡膠壩庫(kù)區(qū)泥沙淤積問(wèn)題無(wú)疑是一種可供參考的有效途徑。

[1]李義天,鄧金運(yùn),孫昭華,等.河流水沙災(zāi)害及其防治[M].武漢:武漢大學(xué)出版社,2004.

[2]白玉川，楊建民，黃本勝.二維水沙數(shù)學(xué)模型在復(fù)雜河道治理中的應(yīng)用[J].水利學(xué)報(bào)，2003，（9）：25-30.

[3]孫小軍,黃志文,魏炳乾.灞河流域水文分析[J].黑龍江水專(zhuān)學(xué)報(bào)，2005，(4):40-42.

[4]邵學(xué)軍，王興奎.河流動(dòng)力學(xué)概論[M].北京:清華大學(xué)出版社,2005:49-148.