基于FLOW 3D的山丘區(qū)小流域山洪演進(jìn)三維數(shù)值模擬研究

易建州，楊培生，許小華

（江西省水利科學(xué)院，江西 南昌，330029）

0 前言

山洪災(zāi)害防治是我國防汛工作的短板，也是各學(xué)者高度關(guān)注的研究領(lǐng)域[1]。隨著丘陵地區(qū)人口的不斷增長以及氣候變化引起的局部強(qiáng)降雨頻率和強(qiáng)度的增加，相關(guān)研究越來越受到國際社會的高度重視。各學(xué)者對洪水的模擬均開展了一定的研究，劉昌軍[2]針對小流域暴雨山洪精細(xì)模擬問題，提出了小流域時空變源混合產(chǎn)流模型，以小流域為單元構(gòu)建了暴雨山洪分布式模擬模型。朱呈浩[3]比較了不同泄流計算方法的優(yōu)缺點(diǎn)及適用范圍，改進(jìn)泄流計算并構(gòu)建了城市雨洪模型。熊俊楠[4]利用重慶市1950～2015年歷史山洪災(zāi)害數(shù)據(jù)，采用平均中心、標(biāo)準(zhǔn)差橢圓、核密度分析和M-K突變檢測等方法分析了重慶市歷史山洪災(zāi)害時空分布規(guī)律，并在此基礎(chǔ)上分析了山洪與各影響因素的相關(guān)性。

近十年來，我國進(jìn)一步加大了山洪災(zāi)害防治力度，初步建成了適合我國國情的山洪災(zāi)害防御體系[5]，山洪災(zāi)害防治技術(shù)體系實現(xiàn)了從無到有的突破，在防災(zāi)減災(zāi)中發(fā)揮了重要作用。但由于山洪問題涉及氣象學(xué)、洪水水文學(xué)、沖刷與泥沙輸運(yùn)力學(xué)、水力學(xué)及動力學(xué)等多個學(xué)科，目前在機(jī)理研究方面對其認(rèn)識不夠全面和完善[6]。為了彌補(bǔ)資料獲取困難、物理模型試驗受模型和實驗條件限制等問題，數(shù)值模擬是研究山洪的有效手段[7]。由于山洪暴發(fā)涉及因素較多，成因較復(fù)雜，目前該方面的模擬軟件相對較少，主要有MIKE、FLOW 3D等。本文采用FLOW 3D數(shù)值軟件，基于動量方程、能量方程，進(jìn)行Navier-Stokes方程求解，RNG k-ε湍流模型對方程進(jìn)行封閉[8]，采用VOF方法處理自由面，通過有限體積法（FVM）對方程進(jìn)行離散，對水流進(jìn)行分析。

1 研究區(qū)概況

毛公洞小流域位于江西省宜春市靖安縣中源鄉(xiāng)合港村（中心坐標(biāo) E114.952559，N28.848591），流域面積約5.53km2，在兩山之間形成一條山洪溝，溝道長約1.84km，主河道比降較大，平均坡降i=9.6%，小流域影像如圖1所示，山洪溝現(xiàn)場斷面如圖2所示。本流域是暴雨區(qū)，短歷時暴雨強(qiáng)度較大。強(qiáng)降雨常發(fā)生于每年5、6月份，24h最大暴雨一般在200～400mm，一次降雨一般在300～500mm，最大超過800mm。

圖1 研究區(qū)影像圖

圖2 山洪溝斷面

根據(jù)山洪災(zāi)害調(diào)查評價成果，沿河兩側(cè)村莊防洪年限均低于20年一遇，受影響人口達(dá)450余人[9]，同時20世紀(jì)70年代靖安縣九嶺欣榮鎢礦落戶于此。由于該區(qū)域降雨相對集中，每年汛期都發(fā)生不同程度的山洪災(zāi)害；降雨引起的山洪多年沖刷以及人類活動，在溝道內(nèi)形成一條山洪溝，對沿河村民人身安全和鎢礦作業(yè)區(qū)造成嚴(yán)重威脅[10]。

1998年6月該處發(fā)生較大山洪災(zāi)害，11日開始連降暴雨，持續(xù)到25日結(jié)束，馬腦背站最大3h降雨118mm，最大6h降雨202mm，洪水量大、峰高。根據(jù)現(xiàn)場調(diào)研及走訪，本流域自20世紀(jì)50年代以來，特大洪水年份有1968年、1982年、1998年、2002年等，洪水多由峰面雨造成。

根據(jù)項目所處區(qū)域，查詢《江西省暴雨洪水查算手冊》（2010年版），獲取其降雨參數(shù)如表1所示，參照水力計算手冊[11]計算其20年一遇的設(shè)計洪水如圖3所示，山洪溝水位流量關(guān)系如表2所示。

表1 降雨參數(shù)

表2 山洪溝水位流量關(guān)系表

圖3 山洪20年一遇設(shè)計洪水（推理公式法）

2 數(shù)值模擬

2.1FLOW 3D模型簡介

FLOW 3D模型是由Flow Science公司開發(fā)的三維流體力學(xué)軟件。與其他的流體動力學(xué)軟件不同，F(xiàn)LOW 3D軟件有其獨(dú)特的自由流體表面跟蹤算法（VOF），能夠追蹤液-液或液-固交界面并結(jié)合有限差分法求解三維N-S方程，可采用多網(wǎng)格體的方法進(jìn)行建模[12]。FLOW 3D可以將STL文件直接導(dǎo)入，便于建立模型，軟件功能強(qiáng)大；能夠模擬復(fù)雜的物理現(xiàn)象，包括對空氣模型、沖刷模型、自由表面流、各種流、凝固和融化、固體顆粒運(yùn)動等復(fù)雜物理現(xiàn)象的模擬；還可將整體結(jié)果直接導(dǎo)入支持的繪圖軟件中，實現(xiàn)實時連續(xù)結(jié)果。該模型在理論上能較為全面和精細(xì)地反映河道三維流場，運(yùn)行結(jié)果更接近于實際，因而廣泛應(yīng)用于水利、水工結(jié)構(gòu)、水環(huán)境等的三維流體仿真領(lǐng)域[13]。

2.2 控制方程

FLOW 3D模擬采用笛卡爾坐標(biāo)系，其主要控制方程包括連續(xù)性方程，動量方程，紊動能K方程和耗散率ε方程[14]，具體如下：

連續(xù)性方程：

動量方程：

紊動能K方程：

紊動能耗散率ε方程：

2.3 模型構(gòu)建

FLOW 3D建模的三維模型分為兩部分。首先，以毛公洞小流域的地形數(shù)據(jù)建立數(shù)值模型，地形范圍為5km×4km，網(wǎng)格數(shù)量采用760萬個，毛公洞小流域高程圖如圖4所示；再將DEM數(shù)據(jù)經(jīng)過預(yù)處理后，轉(zhuǎn)換為FLOW 3D模型可識別的STL格式，生成下墊面的三維模型，模型結(jié)構(gòu)分兩層，最上層為疏松層，厚度約5m，下層為基巖層（將一個多孔介質(zhì)上層疊加在不滲水的基層巖上）；最后將其導(dǎo)入 FLOW 3D后，再精確放置在實際的位置即可。

圖4 毛公洞小流域高程圖

過程中，用到的兩個主要模型是質(zhì)量源模型和多孔介質(zhì)模型。本文定義這個多孔介質(zhì)層為上層表面液體的來源（即通過土體中均勻涌出水來模擬降雨）。較低一層為基巖層，流體不僅滲透到上層和可滲透層，也出現(xiàn)對流現(xiàn)象/納維斯托克斯行為。一旦可滲透層飽和之后，表層水在地表運(yùn)動，表層水的翻滾導(dǎo)致了峽谷中海拔較低處兇猛的洪水。

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

由于山洪暴發(fā)涉及因素較多，國內(nèi)外多考慮降雨為主，因此，本章假設(shè)其他參數(shù)不變，把降雨量作為山洪災(zāi)害防治非工程措施的預(yù)警指標(biāo)，分析不同降雨條件下的山洪發(fā)生過程，進(jìn)而探討其適用性。過程中，采用多孔介質(zhì)模塊模擬降雨，將降雨強(qiáng)度轉(zhuǎn)化為單位面積的流量，以雨水落地后的徑流過程模擬山洪暴發(fā)的過程[15]。

2.4.1 邊界條件設(shè)置

本數(shù)值模擬中設(shè)置了固壁不滑移（壁面邊界條件）、自由表面和出流三種邊界條件，具體邊界條件如下：

（1）自由表面邊界條件。本文中的液體本身沒有任何運(yùn)動，自由液面是空間與水的交界面，在邊界條件中，空氣對流體的作用也要考慮。因此邊界條件可取為：P=Pa（即為大氣壓力），并且要求在該表面上速度分量沿法向的梯度為0，表面切應(yīng)力為0。本文數(shù)值模擬的邊界條件設(shè)置中，上部表面采用自由表面邊界條件。

（2）固壁不滑移邊界條件。一般在在實際情況中，流體與模型底部直接接觸，不會穿過這些材料向下滲透或者運(yùn)動。在FLOW 3D軟件中，提供的固壁不滑移（wall）邊界條件恰好滿足了要求。在此邊界上流體法向速度為0，本文中數(shù)值模擬底部邊界條件的設(shè)置中，將模型底部均設(shè)為固壁邊界。

（3）出流邊界條件。水流自由讓流體自由通過出流邊界，而不造成對計算區(qū)域內(nèi)流場的影響。由于本文模擬的是降雨引發(fā)的山洪暴發(fā)，模型四周需要水流自由進(jìn)出邊界恰好滿足了要求。本文中數(shù)值模擬四周邊界條件的設(shè)置中，將模型四周均設(shè)為出流邊界。

2.4.2 試驗工況

計算過程中，設(shè)置土的孔隙率參數(shù)為0.40，阻力系數(shù)為100，根據(jù)相關(guān)文獻(xiàn)[16]中河床糙率的推薦值：兩岸雜草叢生、河中有植被，河床糙率取0.035。本文主要分析不同降雨時的山洪現(xiàn)象，有3種工況，工況一：135mm/min；工況二：13.5mm/min；工況三：1.35mm/min（如圖 5）。

圖5 山洪演進(jìn)模擬過程圖

3 仿真結(jié)果分析

3.1 工況一計算結(jié)果

本次模擬時間120s，初始水位為37m。當(dāng)降雨量強(qiáng)度為135mm/min時，由計算結(jié)果可知，降雨隨著時間增加區(qū)域內(nèi)計算點(diǎn)水深總體呈上升趨勢，但是也存在特異性，在t=60s到t=90s，水深出現(xiàn)微弱回落趨勢，主要是由于水流的粘滯性，在此時間段，上游水流具有一定的流速，對計算點(diǎn)附近的水流沖擊、攜帶的作用使其流動，導(dǎo)致局部的水深出現(xiàn)回落現(xiàn)象。在t=100s左右時，水深明顯，主要是降雨隨著時間的持續(xù)，流域內(nèi)水流在流域出口處匯集（計算點(diǎn)處），使得水深明顯升高。在t=100s到120s左右，水流出現(xiàn)微弱的回落，主要是由于100s左右時，水流匯集，使局部水深急劇升高，但是隨著流速穩(wěn)定，水流的粘滯力及攜帶作用，使水深在短時間內(nèi)微輻回落。由此可知在降雨條件下山洪的形成過程是，山洪水深先上升，后出現(xiàn)微輻回落，再上升后回落，呈周期性上升的特點(diǎn)。

3.2 工況二計算結(jié)果

由計算結(jié)果可知，降雨強(qiáng)度減弱后，區(qū)域內(nèi)水深較降雨強(qiáng)度大時更低，但是水深隨著降雨的持續(xù)呈周期性增長的規(guī)律與前者相似。在此工況下，也具有其特異性，與降雨強(qiáng)度較大的工況相比，此工況沒有出現(xiàn)水位的急劇升高，水位上升較緩和，主要是由于降雨強(qiáng)度明顯減弱，水流流量不大，匯集后水位并沒有急劇上升。

3.3 工況三計算結(jié)果

從計算結(jié)果分析得出，降雨強(qiáng)度從13.5mm/min降至1.35mm/min，由于存在初始水深，隨著時間增加，區(qū)域內(nèi)水深都小于初始水深，且水深呈周期性變化。主要是由于水流具有粘滯性，上游水流的匯集對計算點(diǎn)范圍內(nèi)水流沖擊，攜帶流入下游，導(dǎo)致水深短暫的降低。在持續(xù)來水及降雨的情況下，水位變高，水深加深又回到初始水深左右，如此反復(fù)呈周期性的變化。但是隨著時間的推移水流逐步穩(wěn)定，水深最終穩(wěn)定在一個范圍內(nèi)。

4 結(jié)論

（1）采用FLOW 3D，利用多孔介質(zhì)模塊，模擬不同降雨條件的山洪演進(jìn)，得出演進(jìn)結(jié)果如下：①降雨強(qiáng)度較小時，隨著時間增加，上游水流的匯集對計算點(diǎn)范圍內(nèi)水流沖擊、攜帶后流入下游，使得水深降低。隨著時間增加，最終水深達(dá)到穩(wěn)定。②隨著降雨強(qiáng)度增加，水深隨著降雨的持續(xù)呈周期性變化，與降雨強(qiáng)度小時相比，監(jiān)測點(diǎn)位水深增加，山洪的流量增大。③當(dāng)降雨強(qiáng)度大時，流域的水位呈周期性上升，在上升期間也會有短暫的回落，當(dāng)水流匯集后，水深急劇升高。

（2）由于三維模型模擬范圍增大后網(wǎng)格數(shù)過多，計算時間太長，本研究僅涉及河道兩岸300m范圍內(nèi)的山洪演進(jìn)模擬。今后的研究應(yīng)結(jié)合其他二維水動力模型，開展更大范圍的、甚至整條流域的山洪演進(jìn)研究。