邵錦焯

(廣州市水務(wù)規(guī)劃勘測設(shè)計研究院有限公司，廣東 廣州 510640)

0 引言

增江及分洪道位于廣州市增城區(qū)石灘鎮(zhèn)，本研究分別采用TUFLOW與MIKE 21建立增江及分洪道的二維水動力數(shù)學(xué)模型，對河道水流流場、水位等進行數(shù)值模擬研究，分析河道水動力，揭示不同模型軟件的特點。

1 數(shù)學(xué)模型簡介

1.1 TUFLOW模型

TUFLOW由某公司聯(lián)合昆士蘭大學(xué)在1989—1990年共同開發(fā)，在隨后的多年持續(xù)改進、完善和擴充功能，已在全球推廣應(yīng)用。應(yīng)用TUFLOW軟件可以構(gòu)建精細化的分布式、直接降雨法、水文水動力耦合、一維二維耦合暴雨洪水模型[1]，并采用GPU并行計算技術(shù)提高計算速度，采用Quadtree網(wǎng)格四杈剖分技術(shù)、SGS(Sub-Grid-Sampling)地形數(shù)據(jù)亞網(wǎng)格提取技術(shù)，兩者技術(shù)結(jié)合使得對河道可以構(gòu)建內(nèi)外連成一體的二維洪水模型，能夠有效提高構(gòu)建洪水模型的效率，獲得高精度的洪水模擬成果。

針對二維淺水方程的求解，TUFLOW使用的是基于Stelling網(wǎng)格的交替方向隱格式有限差分法(ADI)。采用追趕法求解，首先利用連續(xù)方程和x方向動量方程求解水位ζ和流速u：將動量方程代入到連續(xù)方程并消去流速u，獲得1個三對角矩陣，從而求得水位ζ，接著將其反代入到動量方程中求出流速u，更換求解方向即可求得v。

1.2 MIKE21模型

MIKE 21是專業(yè)的二維水動力計算模塊，用于模擬河流、湖泊、河口等水流、波浪狀態(tài)，高級圖形用戶界面與高效的計算引擎結(jié)合，為工程應(yīng)用、規(guī)劃設(shè)計提供完備、有效的設(shè)計環(huán)境[2-5]。MIKE 21具有用戶界面友好、強大的前后處理功能、多種計算網(wǎng)格模塊、干濕節(jié)點單元設(shè)置、集成式多種控制性結(jié)構(gòu)設(shè)置、可視化建模、集成展示等特點。

MIKE 21的水流連續(xù)方程如下：

(1)

水流運動方程如下：

(2)

(3)

式中，x、y、t—分別為空間、時間坐標；z—水位；h—水深；u、v—分別為垂線平均流速在x、y方向上的分量；M、N—分別為單寬流量在x、y方向上的分量，M=hu，N=hv；n—曼寧糙率系數(shù)；c—謝才系數(shù)；vt—紊動粘性系數(shù)；g—重力加速度。

2 增江及分洪道二維水動力模型

2.1 工程區(qū)域基本情況

2.1.1流域概況

研究工程區(qū)域位于增江流域下游。增江發(fā)源于廣東省新豐縣七星嶺，全長203km，流域面積3160km2，增城區(qū)內(nèi)長66km，河寬90～220m，區(qū)間集雨面積971km2，占全區(qū)面積的53%。增江上游建有天堂山水庫及梅州、白沙河、七星墩水庫，區(qū)內(nèi)有百一座中型水庫(花林水庫)，59座小型水庫，主要支流96條，總長416.62km，石灘鎮(zhèn)內(nèi)主要支流有縣江河、水門頭涌等。

2.1.2氣象條件

研究工程區(qū)域位于廣州珠江三角洲地區(qū)，屬亞熱帶氣候，區(qū)域受東南亞季風(fēng)影響很大，且處于低緯度地區(qū)，太陽輻射強，日照時數(shù)多，平均氣溫高，氣候炎熱多雨，夏季綿長。

區(qū)域為廣州市的暴雨區(qū)，由于地形多山，南部臨近海洋，氣流的抬升經(jīng)常帶來強度大、歷時長的暴雨。多年平均降雨量為1820mm，但年內(nèi)分配不均，4—6月多季風(fēng)雨，占全年降雨量的46.7%，7—9月多臺風(fēng)雨，占全年雨量的36.27%，其余10月—次年3月降雨量只占全年的17.03%。據(jù)新家埔站實測，最大24h降雨量為476mm(1981年6月29日)。多年平均蒸發(fā)量為1232mm。流域內(nèi)水汽充沛，濕度較大，平均相對濕度達84%，極端最大相對濕度99%。

2.1.3工程區(qū)域現(xiàn)狀問題

近年來，受氣候變化和城市下墊面影響，城市熱島效應(yīng)嚴重，增城區(qū)降雨范圍、量級存在不確定性[6]。近年來極端天氣頻發(fā)，特別是2020年發(fā)生的“5.22”和“6.8”兩場特大暴雨，具有強度大、范圍廣、面雨量大的特點，增江堤岸出現(xiàn)漫堤，水深超1m，造成嚴重的經(jīng)濟社會損失。增江分洪道進口位于初溪水利樞紐下游約6km，分洪道進口至出口段的增江長6.3km、分洪道長7.6km，出口至增江口約2.3km。分洪道進口前增江河道寬約400m，增江石灘鎮(zhèn)中心段河道最窄處僅有99m，過洪能力不足，且河道內(nèi)存在鐵路橋、公路橋等橋梁存在壅水現(xiàn)象。分洪道河灘地高程較高并有農(nóng)田、魚塘、果林等，經(jīng)多年運行，分洪道存在分洪比例不協(xié)調(diào)的情況。如圖1—2所示。

圖1 增江及分洪道位置示意圖

圖2 增江及分洪道堤圍分布圖

從歷史水系分析，分洪道自20世紀80年代已形成，河道走勢較為穩(wěn)定，至今保留較為天然的河道形態(tài)[7]。根據(jù)《增江中下游干流設(shè)計洪潮水面線復(fù)核》[8]，2008年6月27日增江出現(xiàn)十年一遇洪水，省水文局廣州分局進行了現(xiàn)場測驗，左汊(分洪道)過流440m3/s，右汊(增江)過流2220m3/s，分流比例16.7%。

2.2 建模范圍及網(wǎng)格劃分

本研究的建模范圍為增江及分洪道，下游為新家埔水位站，上游為甩洲水位站，增江河道建模范圍為0+710～9+200，總長8.49km，分洪道建模范圍為0+110～8+210，總長8.10km，模擬水域面積4.69km2，并考慮莞高速橋、增江大橋、鐵路橋、犁耙渡橋等橋梁。

TUFLOW采用2D矩形網(wǎng)格嵌套的模式，先通過堤防確定計算范圍，對河道測量數(shù)據(jù)進行水下地形格柵化，設(shè)置堤內(nèi)地面高程14m，并生成asc格式文件用于模型處理，通過tgc控制文件設(shè)置模型區(qū)域x=9000，y=3000，網(wǎng)格大小Cell Size=5，通過tcf文件設(shè)置啟用SGS技術(shù)。如圖3—6所示。

圖3 建模范圍

圖4 水位站分布圖

圖5 TUFLOW網(wǎng)格地形

圖6 橋墩處局部示意圖

MIKE 21采用2D三角網(wǎng)格模式，先通過堤線確定計算范圍，對河道測量數(shù)據(jù)進行水下地形矢量化，對堤內(nèi)地面和橋墩位置不成生網(wǎng)格，設(shè)置網(wǎng)格大小不大于50m2，河道局部變化及橋墩處加密網(wǎng)格，生成網(wǎng)格總數(shù)148975個，平均網(wǎng)格大小31m2，最小網(wǎng)格大小0.08m2，生成mesh格式文件用于模型處理。如圖7—8所示。

圖7 MIKE 21網(wǎng)格地形

圖8 鐵路橋橋墩處局部示意圖

2.3 模型率定

采用2020年的6.8增江洪水進行率定。收集到甩洲、石灘圍、新家埔3個水位站當次洪水前后三日的水位實測數(shù)據(jù)，新家埔最高水位為3.06m(珠基，下同)，石灘圍為5.66m，甩洲為6.27m(超20年一遇)。通過甩洲站水位及增江河道大斷面，計算相應(yīng)洪峰流量為3700m3/s。

由于增江及分洪道水位高時，周邊地面主要靠排澇泵站排水，泵排流量較小，相比于增江、分洪道流量可忽略。本次研究下游采用新家埔的水位邊界，上游采用甩洲處流量作為邊界條件率定。如圖9—11所示。通過兩種模型軟件分別計算，計算結(jié)果見表1。

表1 6.8洪水兩種模型軟件計算成果對比表 單位：m

圖9 甩洲、石灘圍、新家埔水位站三日水位過程線

圖10 6.8洪水TUFLOW計算結(jié)果

圖11 6.8洪水MIKE 21計算結(jié)果

由上述對比圖可見，在地形基礎(chǔ)數(shù)據(jù)、邊界條件相同，網(wǎng)格大小相似的情況下，采用TUFLOW及MIKE 21兩種不同的模型軟件計算的水位分布相近，計算水位與實測水位吻合，TUFLOW計算結(jié)果略比MIKE 21的高。由于對橋墩的處理方式不同，以及結(jié)果輸出展示軟件的不同，水位分布效果略有差異。

2.4 模型應(yīng)用

經(jīng)過率定后的模型，可應(yīng)用于汛期洪水水位預(yù)報。增江若出現(xiàn)5年一遇左右洪水，麒麟咀站出現(xiàn)2350m3/s的最大流量，東江博羅站出現(xiàn)近5年一遇洪水。根據(jù)上述邊界條件，折合至本模擬范圍，下游邊界水位采用3.85m，上游甩洲流量取2380m3/s。通過兩種模型軟件分別計算，計算結(jié)果如下。如圖12—13所示。

圖12 5年一遇洪水TUFLOW計算結(jié)果

圖13 5年一遇MIKE 21計算結(jié)果

增江石灘鎮(zhèn)中心區(qū)段路堤結(jié)合(樁號6+125～7+000)、堤頂高程為4.48～4.92m，2020年6.8期間曾出現(xiàn)1m左右的漫堤。根據(jù)《增江中下游干流設(shè)計洪潮水面線復(fù)核》[3]此段河道相應(yīng)5年一遇設(shè)計水位為4.65～4.89m，通過5年一遇洪水計算，TUFLOW模型預(yù)測水位為4.83～5.12m，MIKE 21模型預(yù)測水位為4.71～5.04m，均比原設(shè)計水位提高，分析原因為河道演變、分洪道多年運行過流減小、本次計算考慮最新實測橋墩等精細化計算。

根據(jù)上述計算結(jié)果，增江石灘鎮(zhèn)區(qū)段堤防在遭遇5年一遇洪水時，預(yù)計將出現(xiàn)約40cm的漫堤風(fēng)險，據(jù)此可在汛期準備臨時應(yīng)急沙包以應(yīng)對可能出現(xiàn)的增江高水位風(fēng)險。

3 TUFLOW與MIKE 21模型的對比

3.1 建模方式

TUFLOW的模型文件主要歸類于bc_dbase、check、model、results、runs等，分別存放邊界文件、檢查文件、模型數(shù)據(jù)文件、結(jié)果文件、運算文件等，依托GIS進行各類矢量數(shù)據(jù)處理、應(yīng)用電子表格對邊界文件和地類文件處理，通過Notepad++進行控制文件編寫，對文件管理嚴格，建模過程靈活。MIKE 21的界面友好，模型文件主要為mdf、mesh、m21fm等，分別用于數(shù)據(jù)處理、網(wǎng)格文件生成修改、模型總運行文件，建模過程需結(jié)合CAD、Dxf2xyz工具、電子表格等軟件處理，流水線式建模，文件較少、集成度高，可視化程度高。

3.2 計算結(jié)果

TUFLOW與MIKE 21在地形原數(shù)據(jù)、邊界條件等一致的情況下，計算結(jié)果略有差異，適用與二維模型計算各有特點，且在考慮河道彎道、水下地形變化、橋墩流態(tài)變化等都較為優(yōu)異，比一維水動力模型計算結(jié)果更符合實際。在橋墩處理上，本次處理方式不一樣，TUFLOW可采用單獨橋涵設(shè)置的方式，但本次研究采用抬高局部河底高程的方式，MIKE 21采用不生成網(wǎng)格的方式進行，因此在流態(tài)上略有差異。總體來說，兩種模型在本次計算實例中，都能得到很好的應(yīng)用。

3.3 計算速度

TUFLOW的GPU并行計算技術(shù)、SGS技術(shù)等有效提高運算效率，修改網(wǎng)格大小簡便，能根據(jù)實際情況提高計算速度。MIKE 21亦能啟用GPU計算，計算速度總體較慢。

本研究采用兩臺不同配置的臺式計算機，分別啟用中央處理器(CPU)和圖形處理器(GPU)組合，測試TUFLOW兩種網(wǎng)格大小與MIKE 21模型的運行時間。見表2。

表2 兩種模型軟件計算速度對比表

(1)網(wǎng)格大小接近的情況下，TUFLOW總體比MIKE 21快15～25倍；

(2)同樣的計算機，啟用GPU相比不啟用GPU，MIKE 21計算快3.2～7.8倍，TUFLOW計算快3.2～16.8倍。

(3)TUFLOW網(wǎng)格從5m增大至20m，網(wǎng)格的大小不同引起的計算速度不同，采用CPU+GPU運算差異是13～26倍，采用CPU運算差異是43～63倍。

(4)CPU、GPU是提高計算速度的關(guān)鍵，計算機的獨立顯卡宜更新至最新的驅(qū)動程序；二維水動力模型在計算過程中需讀取和生成大量數(shù)據(jù)，固態(tài)硬盤比機械硬盤讀寫數(shù)據(jù)快，配備固態(tài)硬盤的計算機在一定程度上提高運行速度。

(5)TUFLOW在運行大范圍、高精度的模型，計算速度優(yōu)勢更明顯。

4 結(jié)語

本文以增江及分洪道為例，分別采用TUFLOW和MIKE 21研究河道二維水動力數(shù)學(xué)模型，兩者在河道二維水動力模擬上都能較好的反映真實情況，在計算速度上，TUFLOW總體比MIKE 21較快。

作為兩款應(yīng)用廣泛的水動力計算商業(yè)軟件，在實際應(yīng)用中，建模方式、計算速度、結(jié)果展示等方面各有特點。在大范圍、高精度、高效預(yù)警預(yù)報的模擬應(yīng)用中，采用TUFLOW更合適；在小范圍、對計算速度要求不高、模擬結(jié)果展示效果方面，MIKE 21更有優(yōu)勢。