沖沙河道演變過程的數(shù)值研究

申 俊

(塔里木河干流管理局，新疆 庫(kù)爾勒 841000)

本文采用三維數(shù)值模型，模擬了新疆地區(qū)塔里木河流域希尼爾水庫(kù)2020年沖洗作業(yè)過程中河床隨表面速度場(chǎng)的演變趨勢(shì)，并將模擬結(jié)果與現(xiàn)場(chǎng)測(cè)量結(jié)果進(jìn)行比較。此外，還定量評(píng)估了沖出泥沙量在不同沖洗作業(yè)階段的變化情況。

1 材料與方法

1.1 案例說明

本文所研究的大壩水庫(kù)總庫(kù)容和有效庫(kù)容分別為9.01×106m3和1.66×106m3。圖1顯示了2020年6月進(jìn)行的水庫(kù)沖洗作業(yè)之前測(cè)量的河床水位，以及用于進(jìn)一步定量評(píng)估河床演變模式的橫截面位置(A-F)。將近2 km長(zhǎng)的研究案例分為三個(gè)區(qū)域，即區(qū)域Ⅰ、Ⅱ和Ⅲ，根據(jù)每個(gè)區(qū)域存在的河床物質(zhì)類型分析研究。Ⅰ區(qū)河床物質(zhì)類型為粗粒型，Ⅲ區(qū)和Ⅱ區(qū)后半部分河床物質(zhì)類型則變?yōu)榧?xì)粒型。研究區(qū)侵蝕沉積物的總體積為408.7×103m3。在該水庫(kù)中，中間區(qū)域的左岸(即區(qū)域Ⅱ)沉積的細(xì)泥沙(圖1中矩形框所示)無法有效清除。表1顯示了沖洗作業(yè)前不同橫截面的平均泥沙粒徑分布。選取0.37 ～ 316 mm的7種沉積物尺寸作為代表性粒徑。圖2顯示了2020年沖洗作業(yè)期間的進(jìn)入流量和水位變化，其中，主要的泥沙流入量是沖洗負(fù)荷。開始作業(yè)后，水位初步下降時(shí)間為8～24 h，自由流動(dòng)狀態(tài)為24～38 h。

圖1 水庫(kù)沖水作業(yè)前的實(shí)測(cè)河床地形以及用于評(píng)估河床演變的橫截面A-F位置

圖2 沖洗作業(yè)期間的水位和進(jìn)水流量變化

1.2 數(shù)值模型

數(shù)值模型求解質(zhì)量守恒方程和雷諾平均N-S方程來計(jì)算湍流的水運(yùn)動(dòng)，如式(1)、式(2)：

(1)

(2)

本文采用有限體積法進(jìn)行離散化[1-5]，采用二階迎風(fēng)格式求解雷諾平均N-S方程中的對(duì)流項(xiàng)。湍流由標(biāo)準(zhǔn)k-ε模型建模，網(wǎng)格是非正交的、非結(jié)構(gòu)化的和自適應(yīng)的，并隨著河床和自由水面的變化而垂直移動(dòng)。因此，僅對(duì)水體進(jìn)行建模。網(wǎng)格在每次時(shí)間步長(zhǎng)后重新生成，采用的潤(rùn)濕/干燥算法使模型能夠根據(jù)水和床位的變化具有不同數(shù)量的網(wǎng)格單元(例如在垂直方向上)。潤(rùn)濕/干燥算法消除了水深小于定義邊界的單元。

進(jìn)水采用第一類邊界條件(對(duì)數(shù)速度分布)，出水和泥沙采用零梯度邊界條件。對(duì)于河床和墻壁的邊界條件，在沒有水通量的情況下，采用壁面定律。河床形式，即沙丘效應(yīng)和波紋效應(yīng)，也通過經(jīng)驗(yàn)公式考慮。模擬形態(tài)變化的泥沙輸移計(jì)算分為懸移質(zhì)輸移和推移質(zhì)輸移。懸移質(zhì)通過求解瞬態(tài)對(duì)流-擴(kuò)散方程公式計(jì)算得出，推移質(zhì)可通過修正的MPM公式或Van Rijn公式來模擬。MPM公式更適用于陡峭的河流，這些河流主要輸送靠近河床的粗泥沙，如式(3)所示。

(3)

式中：qb，i為單位寬度推移質(zhì)第i部分的泥沙輸送速率，kg/s；g為重力加速度，m/s2；ρs為泥沙密度，kg/m3；ρw為水的密度，kg/m3；d50為特征泥沙粒徑，mm；r為水力半徑，mm；I為能量線的斜率。

Van Rijn公式如式(4)為：

(4)

式中：di為第i部分的直徑，m；τ為剪切應(yīng)力，N/m2；τc，i為di的臨界剪切應(yīng)力，N/m2；v為運(yùn)動(dòng)黏度，Pa·s。模型中引入了在單個(gè)時(shí)間步長(zhǎng)內(nèi)可被侵蝕的上部活性層厚度以及為活性層提供足夠沉積物的下部非活性層厚度。

該模型利用初始運(yùn)動(dòng)臨界剪應(yīng)力的折減函數(shù)和沙滑算法來考慮邊坡效應(yīng)。當(dāng)河床坡度超過河床物質(zhì)的休止角時(shí)，沙滑算法將修正河床坡度。如果沒有內(nèi)聚力影響，泥沙顆粒可以在邊坡上相互獨(dú)立移動(dòng)，并可以適當(dāng)模擬河岸侵蝕。

2 結(jié)果與討論

2.1 模型校準(zhǔn)

根據(jù)圖1所示沖洗前測(cè)量的河床水位，制作計(jì)算網(wǎng)格。水流方向和橫向的網(wǎng)格單元尺寸分別為10～20 m和5～10 m。假定河床物質(zhì)為2650 kg/m3。圖2顯示了模擬中用作水動(dòng)力邊界條件的水位和進(jìn)水流量波動(dòng)。此外，還利用7種具有代表性的沉積物尺寸將水庫(kù)內(nèi)的非均勻河床物質(zhì)尺寸分布引入模型。根據(jù)可用河床物質(zhì)樣品的尺寸插值獲得庫(kù)區(qū)內(nèi)的河床物質(zhì)尺寸分布。由于沖洗負(fù)荷是假定在水庫(kù)中無沉積的情況下被輸送，因此，在計(jì)算中忽略了沖洗負(fù)荷對(duì)河床演變的影響。

為了揭示以經(jīng)驗(yàn)為主參數(shù)的影響趨勢(shì)(如活性層厚度、河床物質(zhì)的含水量以及沙滑算法的臨界休止角)，首先，建立了一個(gè)參考案例，假設(shè)上述參數(shù)的一般值。 然后，改變假設(shè)值，計(jì)算沖出泥沙的總體積(以下稱為TVFS)并與測(cè)量值(即4008.7×103m3)進(jìn)行比較，以校準(zhǔn)模型。表2顯示了與參考案例相比，TVFS對(duì)選定的主要經(jīng)驗(yàn)參數(shù)的敏感性分析。

表2 參考案例對(duì)所選經(jīng)驗(yàn)參數(shù)的敏感性分析 m3

可以看出，TVFS隨著活性層厚度和含水量的增加而增強(qiáng)，隨著臨界休止角的增加而降低。當(dāng)大量的河床物質(zhì)(即更厚的活性層)可以在一個(gè)時(shí)間步長(zhǎng)內(nèi)被侵蝕時(shí)，預(yù)計(jì)泥沙的侵蝕量會(huì)更大。假設(shè)泥沙沉積中較高的含水量(如50%)會(huì)降低河床物質(zhì)的淹沒密度，從而導(dǎo)致較高的夾帶泥沙。在每個(gè)時(shí)間步長(zhǎng)后，較高的臨界休止角會(huì)使沖沙河道的側(cè)岸更陡，導(dǎo)致水道進(jìn)一步加深，這并不是增加TVFS的有效方法。相反，較低的臨界休止角可以實(shí)現(xiàn)沖沙河道的橫向發(fā)展，有利于增加TVFS。除了TVFS，最終將模擬的河床地形模式與測(cè)量結(jié)果進(jìn)行了定性比較，以提供Ⅲ區(qū)下部周圍侵蝕的信息。之后，對(duì)參考案例進(jìn)行了修改并使用更新后的值進(jìn)行驗(yàn)證(例如活性層厚度、含水量和臨界休止角分別設(shè)置為1 m、40%和33°)。

2.2 沖沙河道演變模擬

2020年6月沖洗作業(yè)后數(shù)值模擬的河床水位與沖洗后的實(shí)測(cè)水位對(duì)比結(jié)果：在水位下降過程中，數(shù)值模型在自由流動(dòng)條件下觀察到獨(dú)特的沖沙河道，其特征與原型相似。

圖3顯示了采用MPM公式、Van Rijn公式(分別為圖3(a)和3(b))和實(shí)測(cè)結(jié)果(圖3(c))的沖洗后的模擬河床水位。可以看出，應(yīng)用Van Rijn公式可以降低Ⅲ區(qū)的侵蝕，而Ⅱ區(qū)的侵蝕被高估了。使用MPM公式計(jì)算的TVFS為313.14×103m3，而使用Van Rijn公式計(jì)算的TVFS為333.91 ×103m3。

圖3 模擬河床水位

當(dāng)采用MPM公式時(shí)，Ⅲ區(qū)被侵蝕的粗粒型物質(zhì)可能會(huì)再次沉積到Ⅰ區(qū)。Ⅰ區(qū)水位較深導(dǎo)致流速降低。此外，在Ⅰ區(qū)，無論所用推移質(zhì)輸沙公式類型如何，模擬結(jié)果都顯示出與觀測(cè)結(jié)果相比更窄的沖沙河道。這是由于在靠近大壩的大彎曲區(qū)域是由二次流和反向流形成的復(fù)雜流場(chǎng)。該區(qū)域形成的河床形式(即沙丘和波紋)將增加水流的復(fù)雜性，進(jìn)而增加泥沙輸移模式。此外，該區(qū)域的輸沙能力達(dá)到了最高水平，但仍低于原型中的原始值。因此，河床水位與實(shí)測(cè)水位不同。

由于水庫(kù)不同河段的推移質(zhì)類型和水位不同，每個(gè)推移質(zhì)輸運(yùn)公式(即MPM和 Van Rijn)都可以對(duì)特定段的河床演變模式提供更合理的預(yù)測(cè)。因此，將MPM和Van Rijn公式模擬的最終河床高程與圖4中橫截面A-A和E-E的測(cè)量值進(jìn)行了比較。選擇MPM公式對(duì)數(shù)值結(jié)果進(jìn)行進(jìn)一步的定量評(píng)估。圖5還給出了使用MPM公式在橫截面A-A、B-B、D-D和F-F中數(shù)值模擬的河床演變模式。

圖4 最終河床高程

圖5 各個(gè)橫截面的河床演變模式

在河床物質(zhì)較粗的Ⅲ區(qū)，泥沙侵蝕被低估，模擬的沖沙河道寬度和深度小于測(cè)量值(圖5(a)和圖5(b))。由于粗糙河床物質(zhì)的復(fù)雜流速剖面和有效水動(dòng)力作用，特別是在低水頭的自由流動(dòng)狀態(tài)下，模型無法很好地模擬該區(qū)域的侵蝕模式。另一個(gè)原因是在Ⅲ區(qū)下部左岸的計(jì)算域內(nèi)消除了具有淺水頭(即干涸)的單元，其中已經(jīng)測(cè)量到明顯的河床退化(圖5(b))。在Ⅱ區(qū)，即水庫(kù)的寬中段，通過應(yīng)用MPM公式得到的沖沙河道形狀和位置與觀測(cè)結(jié)果基本一致。然而，在數(shù)值結(jié)果中，沿右路堤的河床退化被高估了(圖5(c))。Ⅱ區(qū)的沖沙河道特征對(duì)保持水庫(kù)蓄水能力和提高沖洗效率有很大貢獻(xiàn)。在自由流動(dòng)沖洗作業(yè)期間，狹窄的沖沙河道位置偏向Ⅱ區(qū)的右岸區(qū)域，因此，左岸區(qū)域中大量沉積物的侵蝕可能性較低。還應(yīng)考慮到研究水庫(kù)是在雨季進(jìn)行沖沙，并且在沖洗前后不進(jìn)行河床水位測(cè)量。因此，測(cè)量的河床水位可能存在不確定性。此外，沿非活性層也沒有高空間分辨率的泥沙粒徑分布測(cè)量，但這會(huì)顯著影響數(shù)值模型結(jié)果的質(zhì)量。

圖6顯示了在自由流動(dòng)條件下具有二次流矢量的水流速度。在沖沙河道中，流速升至4.5 m/s，多個(gè)區(qū)域出現(xiàn)超臨界流(圖6(a)和圖6(b))。三維模型還捕捉了河道彎道中河床發(fā)育的特征(即彎道外側(cè)的侵蝕和內(nèi)側(cè)的沉積)，并再現(xiàn)了寬度上的非對(duì)稱速度剖面，以及彎道頂點(diǎn)的側(cè)向水面傾斜(圖6(c))。

圖6 自由流動(dòng)條件下具有二次流矢量的水流速度

2.3 沖沙量變化

根據(jù)沖洗過程中每隔一小時(shí)沖洗掉的侵蝕河床物質(zhì)的計(jì)算通量可知，平均流速和湍流波動(dòng)的增加，即初流階段和自由流動(dòng)階段的最后部分，往往會(huì)增強(qiáng)泥沙的流動(dòng)性和隨后的泥沙夾帶量。 較小的沉積物尺寸往往比較大的沉積物尺寸更早被侵蝕和沖走，如直徑為3.7 mm的沉積物在t=5 h的洪峰后立即出現(xiàn)高濃度。直徑較小的沉積物對(duì)水位和流量波動(dòng)都很敏感。另一方面，直徑較大的沉積物對(duì)水位變化比對(duì)流量變化更敏感。

3 結(jié) 論

(1)三維數(shù)值模型能夠較好地模擬真實(shí)邊界條件(如非穩(wěn)態(tài)水動(dòng)力條件)下的沖沙河道演變模式。Van Rijn和MPM推移質(zhì)輸沙公式都可以在一定范圍內(nèi)代表沖洗過程中沖出的泥沙量。

(2)MPM公式能較好地預(yù)測(cè)由較粗物質(zhì)覆蓋的上游區(qū)域的河床水位變化。Van Rijn公式在靠近大壩的區(qū)域表現(xiàn)出更好的性能，其中，河床主要被細(xì)粒物質(zhì)覆蓋，而且水庫(kù)也更深。

(3)三維數(shù)值模型可以識(shí)別彎道中的螺旋流以及由此產(chǎn)生的侵蝕和沉積模式，比二維建模更具優(yōu)勢(shì)。

(4)就沖洗過程中沖出的泥沙量來看，較粗的河床物質(zhì)主要在初步水位下降結(jié)束時(shí)和自由流動(dòng)狀態(tài)下沖出。因此，在自由流動(dòng)條件下引入額外的流量和小幅波動(dòng)的水位降低可能會(huì)增加泥沙侵蝕的可能性。