河灘濕地對(duì)河道洪水影響的二維數(shù)值模擬

葉 群，計(jì) 勇，吳培軍，王 寅，張 潔，王 凱

(1.南昌工程學(xué)院水利與生態(tài)工程學(xué)院，南昌 330099；2.江西省鄱陽(yáng)湖水利樞紐建設(shè)辦公室，南昌 330046)

圖1 樟里河及濕地位置圖Fig.1 The location of Zhangli River and wetland

濕地作為一種重要的生態(tài)系統(tǒng)，它具有調(diào)節(jié)生態(tài)環(huán)境的特殊功能[1]，以及涵養(yǎng)水源[2]、生物多樣性[3]、凈化水質(zhì)[4]、地下水補(bǔ)給[5]和蓄滯洪水[6]等作用。國(guó)內(nèi)外主要采用水文分析與數(shù)學(xué)模型相結(jié)合的方法對(duì)濕地飽和或淹水時(shí)機(jī)、時(shí)長(zhǎng)和水深等水文屬性展開(kāi)研究[7,8]，以指導(dǎo)濕地保護(hù)與運(yùn)行。我國(guó)各地在河道治理中逐漸引入“濕地河道”的設(shè)計(jì)方案[9,10]，如在河道兩側(cè)原河灘地建立濕地公園，作為一種工程與非工程相結(jié)合的措施，以恢復(fù)河道走廊生態(tài)功能。樟里河位于江西省西部，萍鄉(xiāng)市域西南部，為淥水河支流，如圖1所示。樟里河主河道全長(zhǎng)9.36 km，控制流域面積26.6 km2，作為典型的山溪性河流，季節(jié)性強(qiáng)且發(fā)生頻率大。樟里河濕地位于主河道中游，濕地面積達(dá)49 800 m2。為探明該區(qū)域洪水規(guī)律，本文建立了平面二維水流數(shù)學(xué)模型對(duì)河道及周邊濕地水流特征進(jìn)行了模擬研究與分析，有利于明晰該區(qū)域濕地的存在對(duì)洪水的滯蓄作用，為該區(qū)域的防洪規(guī)劃設(shè)計(jì)提供參考。

1 水動(dòng)力數(shù)學(xué)模型

1.1 一維模型的率定

1.1.1 模型邊界條件推求

一維水動(dòng)力模型控制方程采用一維非恒定流Saint-Venant方程組，其中連續(xù)性方程反映了水流質(zhì)量守恒定律，動(dòng)量方程反映了動(dòng)量守恒定律。

(1)

式中：Q為流量，m/s；q為側(cè)向入流，m3/s；A為過(guò)水面積，m2；h為水位，m；R為水力半徑，m；C為謝才系數(shù)；α為動(dòng)量修正系數(shù)。

1.1.2 模型邊界條件推求

對(duì)樟里河河網(wǎng)進(jìn)行概化，如圖2所示。概化后的樟里河河網(wǎng)由樟里河主河道及樟里河支流甘泉河組成，其中樟里河主河道段3 650 m，支流段2 369 m。樟里河的斷面由實(shí)測(cè)資料進(jìn)行概化，斷面分布間隔大約為50 m，部分間隔100 m,模型中具體斷面數(shù)為干流67處斷面，支流49處斷面。由于濕地研究區(qū)位置位于干流處，故選取圖1中16處斷面進(jìn)行模擬水位計(jì)算，具體斷面樁號(hào)如表1所示。其中1-1斷面對(duì)應(yīng)樟里河匯入淥水處，即樁號(hào)K3+650處，16-16斷面對(duì)應(yīng)茶嶺下村X128縣道公路橋，即樁號(hào)K0+0處。河道水面曲線的計(jì)算采用一維恒定流數(shù)學(xué)模型來(lái)求解。模型上邊界設(shè)定為流量邊界，設(shè)計(jì)洪水的推求采用推理公式法，下邊界為下游出口斷面的水位流量關(guān)系，由曼寧公式Q=1/n×A×R2/3×J1/2進(jìn)行推求計(jì)算，如表2所示(水位為黃海高程基準(zhǔn)，下同) 。

圖2 MIKE11河網(wǎng)概化(單位：m)Fig.2 The generalization of MIKE 11 river network

表1 樟里河模擬水位與斷面設(shè)計(jì)水位成果分析(P=10%)Tab.1 Analysis the results of simulation and section design flood(P=10%)

表2 樟里河河口水位～流量關(guān)系Tab.2 The relationship between water level and discharge in Zhangli River estuary

1.1.3 一維模型水位率定結(jié)果

為驗(yàn)證模擬結(jié)果的合理性，將模型計(jì)算結(jié)果與湘東區(qū)淥水樟里河綜合治理可行性研究報(bào)告設(shè)計(jì)結(jié)果進(jìn)行對(duì)比分析，各斷面模擬水位值與設(shè)計(jì)水位值對(duì)比結(jié)果見(jiàn)表1。由于淥水河上游無(wú)水文觀測(cè)站，設(shè)計(jì)水位值利用推理公式法查算暴雨手冊(cè)后，采用能量法自下而上逐段進(jìn)行水面線推算而來(lái)。經(jīng)率定計(jì)算值與設(shè)計(jì)值吻合較好，模型糙率取n=0.025～0.046，絕對(duì)誤差與相對(duì)誤差值均在可接受的范圍內(nèi)。

1.2 二維控制方程組

視水流為不可壓流，水流運(yùn)動(dòng)規(guī)律采用Navier-Stokes方程組進(jìn)行描述，符合布辛涅斯克Boussinesq假設(shè)和靜水壓力假設(shè)。模型的平面采用非結(jié)構(gòu)化網(wǎng)格，應(yīng)用有限體積法作為數(shù)值解法。連續(xù)性方程與動(dòng)量方程如下[11]：

(2)

其中c=(1/n)h1/6。

式中：x、y和t分別為空間和時(shí)間坐標(biāo)；h為水深；U、V分別為x、y方向流速；z為水位；g為重力加速度；u、v分別為垂線平均流速在x、y方向的分量；c為謝才系數(shù)；vt為紊動(dòng)黏性系數(shù)；n為曼寧糙率系數(shù)。

1.3 二維模型求解

1.3.1 邊界條件

二維模型的上下邊界取在濕地的進(jìn)出口處，沿河道模型總長(zhǎng)504 m。模型的左右邊界采用均為天然地形構(gòu)成，通過(guò)地形等高線按一定的高程取值來(lái)劃定[12]。上游模型進(jìn)口邊界條件：取p=10%設(shè)計(jì)頻率下的洪峰流量；下游模型出口邊界條件：采用曼寧公式推求出口斷面處的水位～流量關(guān)系，見(jiàn)表3；固定邊界：法向流速為 0，水流邊界的切向流速不為 0，采用壁面函數(shù)處理；初始條件：給出計(jì)算區(qū)域內(nèi)各點(diǎn)初始時(shí)刻的水深、流速等物理量的初值。

1.3.2 參數(shù)設(shè)置

糙率分區(qū)設(shè)置，其中河道糙率采用一維斷面糙率率定結(jié)果進(jìn)行設(shè)置，濕地糙率根據(jù)區(qū)域內(nèi)地形及濕地植被情況，參考《水力學(xué)》中灘地粗糙系數(shù)n值表[13],取值0.05，插值生成糙率場(chǎng)。水平渦黏系數(shù)取默認(rèn)值0.28，干水深0.005 m，淹沒(méi)水深0.05 m，濕水深度0.1 m。

表3 濕地出口斷面水位流量關(guān)系Tab.3 The relationship between water level and discharge in the export section of wetland

1.4 二維模型計(jì)算網(wǎng)格生成

根據(jù)研究區(qū)邊界與地形特征，采用非結(jié)構(gòu)三角網(wǎng)格進(jìn)行剖分。工況一(河道)網(wǎng)格數(shù)1 358，網(wǎng)格最大步長(zhǎng)為3m；工況二(河道加濕地)網(wǎng)格數(shù)10 815，對(duì)河道進(jìn)行網(wǎng)格加密，濕地網(wǎng)格最大步長(zhǎng)6 m，河道最大步長(zhǎng)2 m；工況三(河道加濕地加涵洞)網(wǎng)格數(shù)11 939，最大步長(zhǎng)同工況二，其中涵洞的設(shè)置的尺寸為長(zhǎng)5 m、寬2 m，涵洞的位置設(shè)置為由北到南，間隔約24 m設(shè)置一個(gè)涵洞，共設(shè)置了3個(gè)涵洞。各工況模型網(wǎng)格剖分如圖3。

圖3 網(wǎng)格剖分圖(單位：m)Fig.3 The graph of mesh generation

2 模型計(jì)算結(jié)果與分析

2.1 水位過(guò)程分析

如圖4所示，十年一遇洪峰流量情況下河段[圖4(a)]最高水位為82.73 m，最低水位為81.75 m，水位差達(dá)0.98 m。現(xiàn)河岸高程不能滿足防洪要求，需要加高對(duì)應(yīng)岸堤高程。增加河灘濕地[圖4(b)]后，最高水位降為82.3 m，降低了0.43 m，整體區(qū)域水位降低。濕地受公路的阻隔，被分為上下游兩塊區(qū)域，原本主河道行洪區(qū)域形成束窄，過(guò)水?dāng)嗝媸湛s，上下游濕地形成較明顯的水位落差約0.35 m，存在上游濕地?cái)嗝嫠慧崭攥F(xiàn)象，影響范圍為上游濕地全域，河道洪水進(jìn)入下游濕地，水流得到較大的擴(kuò)散，水位快速下降，水位基本維持81.95 m，較工況一下游水位下降最多處達(dá)0.43 m。在濕地的基礎(chǔ)上添加涵洞后[圖4(c)]，最高水位出現(xiàn)在上游濕地區(qū)域?yàn)?2.16 m，較工況二降低0.14 m，有效緩解上游水位壅高現(xiàn)象，回水水域減小，河道束窄區(qū)域由于涵洞的分流作用，上下游濕地水位落差降為0.16 m，有利于公路地基的穩(wěn)定。涵洞添加處出現(xiàn)局部水位壅高現(xiàn)象，但影響范圍較小。較未加涵洞情況，束窄處水位落差下降0.19 m，有利于濕地行洪，有利于束窄區(qū)域河岸安全。

圖4 水位等值線圖(單位：m)Fig.4 The contour map of water level

2.2 流場(chǎng)過(guò)程與分析

如圖5所示，河道[圖5(a)]最大流速發(fā)生在河道匯流束窄處約為2.48 m/s，分析表明，河道平均流速為1.04 m/s。增加河灘濕地[圖5(b)]模擬結(jié)果顯示，由于受束窄區(qū)域影響，最大流速同樣發(fā)生在河道束窄處，流速為2.6 m/s，局部流速的增大，不利于河岸的穩(wěn)定性，該處區(qū)域易受到洪水沖刷侵蝕，同時(shí)束窄處水流情況復(fù)雜，出現(xiàn)上下邊界繞流、回流、斷面水位雍積現(xiàn)象。濕地的存在，對(duì)水流起到擴(kuò)散作用，主河道流速明顯下降。經(jīng)分析計(jì)算，區(qū)域的平均流速為0.3 m/s，極大地降低了流速。在濕地的基礎(chǔ)上增加涵洞[圖5(c)]，模擬結(jié)果顯示，主流河道上下游的流速均得到較明顯的降低，經(jīng)分析計(jì)算，平均流速降為0.26 m/s，較濕地情況，流速得到進(jìn)一步降低。同時(shí)涵洞的存在連通了上下游濕地，進(jìn)一步進(jìn)行分流，雖流速極值較濕地沒(méi)有太大變化，但束窄流速得到較大降低，流速降為1.65 m/s，有利于防止沖刷，防止河岸邊坡不穩(wěn)定等情況發(fā)生，同時(shí)整體流速差進(jìn)一步降低。

圖5 流速等值線圖(長(zhǎng)度單位：m；流速：m/s)Fig.5 The contour map of velocity

2.3 控制點(diǎn)對(duì)比分析

沿主河道沿程選取5個(gè)控制點(diǎn)，控制點(diǎn)位置見(jiàn)圖3(d)。從水位分析，如圖6(a)所示，濕地顯著降低河道的水位，涵洞的添加，進(jìn)一步降低上游濕地的水位，同時(shí)減少上下游的水位差，由原先的河道上下游水位差由0.6 m降至0.15 m，使得水流更為平緩。由于控制點(diǎn)2處水面束窄，使得濕地工況下存在水位頂托效應(yīng)，涵洞的添加緩解了頂托效應(yīng)，降低了上游水位0.2 m。工況2、3水面最大比降發(fā)生在控制點(diǎn)2到控制3這段，這是由于下游濕地的存在，水流得到充分的擴(kuò)散，削弱了洪峰，降低了沿程的洪水位，充分體現(xiàn)出濕地的防洪減洪，削峰作用。

圖6 控制點(diǎn)分析 Fig.6 Analysis of control point

從流速場(chǎng)分析，如圖6(b)所示，濕地的運(yùn)用對(duì)洪水起了分流作用，各控制點(diǎn)的流速下降作用明顯，特征點(diǎn)4處下降最大，達(dá)1.5 m/s。在濕地的基礎(chǔ)上添加涵洞，主要效果體現(xiàn)在束窄區(qū)域流速得到下降。濕地、涵洞的運(yùn)用有效地緩和急流，充分體現(xiàn)出濕地對(duì)于河流主流減緩流速的作用。

3 結(jié) 語(yǔ)

(1)建立的水動(dòng)力數(shù)學(xué)模型，可以較好地模擬復(fù)雜的邊界，經(jīng)定性分析，數(shù)學(xué)模型計(jì)算成果合理，基本能模擬該河段復(fù)雜的水流條件。

(2)計(jì)算區(qū)域水流形態(tài)主要取決于計(jì)算區(qū)域的邊界條件，邊界條件變化較小的上下游，水流流態(tài)平穩(wěn)，而在束窄區(qū)域，水流流動(dòng)復(fù)雜，上下邊界繞流、回流、斷面水位壅積等現(xiàn)象出現(xiàn)較多。

(3) 從水位分析，增加人工濕地可以顯著降低上、下游水位高度。從流速場(chǎng)分析，由于濕地的擴(kuò)散作用，主河道內(nèi)的平均流速顯著降低。增加涵洞可以改變局部流場(chǎng)形態(tài)，減少回水水域，降低束窄區(qū)域的流速。

□