三維泥沙動力數(shù)值模型的高效應(yīng)用
——準(zhǔn)三維模型“輕裝”效應(yīng)

楊 靜，黑鵬飛，假冬冬，尚毅梓

(1.中央民族大學(xué) 生命與環(huán)境科學(xué)學(xué)院， 北京 100081；2.南京水利科學(xué)研究院 水文水資源與水利工程科學(xué)國家重點實驗室， 江蘇 南京 210029;3.中國水科院 流域水循環(huán)模擬與調(diào)控國家重點實驗室, 北京 100038)

三維泥沙動力數(shù)值模型的高效應(yīng)用
——準(zhǔn)三維模型“輕裝”效應(yīng)

楊 靜1，黑鵬飛1，假冬冬2，尚毅梓3

(1.中央民族大學(xué) 生命與環(huán)境科學(xué)學(xué)院， 北京 100081；2.南京水利科學(xué)研究院 水文水資源與水利工程科學(xué)國家重點實驗室， 江蘇 南京 210029;3.中國水科院 流域水循環(huán)模擬與調(diào)控國家重點實驗室, 北京 100038)

近年來泥沙動力和河床沖淤三維計算模型(SB3D)取得了快速的發(fā)展，然而由于SB3D基于全三維水動力模型，計算量過大將會降低SB3D的工程實用性。從模型結(jié)構(gòu)、模型求解與程序代碼編譯三個角度，探討準(zhǔn)三維水動力模型和SB3D聯(lián)用的合理性和實用性，并采用實驗室模擬和工程應(yīng)用予以證明。結(jié)果表明：準(zhǔn)三維模型多數(shù)條件下可以為SB3D提供合理的三維水流速度；全三維水動力模型和程序模塊可以方便的由準(zhǔn)三維模型和模塊替換；準(zhǔn)三維水動力模型極大的提高了模型計算效率，具有理想的“輕裝”效應(yīng)，可以提高SB3D的工程實用性。

準(zhǔn)三維模型；計算河流動力學(xué)；泥沙

河流動力模型包含水動力模型和泥沙動力模型。近年來復(fù)雜形態(tài)河道演變數(shù)值模擬研究，促進了三維泥沙動力模型(SB3D)的迅速發(fā)展[1-7]。但是，這些模型都基于全三維水動力模型[3-6]，而當(dāng)前計算機速度尚無法滿足全三維模型工程計算的要求，限制了SB3D最新成果在工程中的應(yīng)用。準(zhǔn)三維水動力模型基于靜壓假定，極大的提高了三維模型的計算效率，增加了三維水動力模型的應(yīng)用可行性[2]。當(dāng)前準(zhǔn)三維模型越來越多的應(yīng)用于河流水庫的水質(zhì)模擬中[8-9]，遺憾的是，基于準(zhǔn)三維水動力模型的泥沙動力模型研究仍相對較少，準(zhǔn)三維河流動力模型的應(yīng)用優(yōu)勢尚未展現(xiàn)。工程應(yīng)用中，當(dāng)全三維模型因計算量過大而不具有可行性時，則直接選取二維[10-13]或一維泥沙動力模型[14-15]，而未將SB3D模塊應(yīng)用于準(zhǔn)三維模型，由此導(dǎo)致準(zhǔn)三維河流動力模型應(yīng)用的滯后性。

本文在對全三維和準(zhǔn)三維水動力模型進行理論分析和試驗對比(時均流速)的基礎(chǔ)上，將全三維模型的SB3D應(yīng)用于準(zhǔn)三維模型，探索其在準(zhǔn)三維河流動力學(xué)模型中的應(yīng)用可行性和優(yōu)越性。

1 準(zhǔn)三維水動力模型分析

1.1 準(zhǔn)三維水動力模型

河流全三維水動力模型求解主要難點源于自由表面和壓強項。在水位起伏不大的條件下，常由連續(xù)方程直接垂向積分獲得水位方程。自由表面獲得后，壓強并無獨立的控制方程，對此需要采用經(jīng)典流體力學(xué)求解方法，常用方法為SIMPLE法。每一次壓強修正迭代，都需要求解一個七對角方程。若采用GMRES方法求解，Im、Jm、Km分別表示x、y、z方向的網(wǎng)格數(shù)，則計算量為O(Im×Jm×Km)3，如每一時間步中壓力修正n步，則總計算量為O(n×(Im×Jm×Km)3)。

靜壓假定后，動量方程中壓強p可由水位ζ簡單表示。ζ與動量方程常采用分離求解，先采用顯式求解ζn+1，再分別由動量方程和連續(xù)性方程先后求得x、z、y方向流速un+1、vn+1和wn+1，無需采用全三維模型常用的壓力修正法迭代求解[2]。其中最為簡單的是水位和動量方程全部顯式求解(如POM模型)?；蛘邇H對垂向擴散項隱式離散(如EFDC模型)，所得三對角方程也無需迭代求解，追趕法求解計算量僅為O(8Km×Im×Jm)。

水位方程

(1)

水平動量方程

(2)

(3)

垂向流速由下式

(4)

(5)

直接顯式離散求得。以上各式中t為時間，x、y、z為物理坐標(biāo)；σ=(z-zb)/H，zb為床面高程，H為水深；u、v分別為x、y方向的流速；ω、w分別為σ坐標(biāo)下和z坐標(biāo)下的垂向流速；f、g分別為科氏力和重力加速度；Fx和Fy分別是x和y方向的動量源項。

總之，雖然準(zhǔn)三維模型并不適用于垂向速度(w)變化較大情況，但準(zhǔn)三維和全三維水動力模型都可以為SB3D提供三維流速信息，且準(zhǔn)三維模型基于靜壓假定，降低了水動力模型求解難度，增加三維模型的求解效率。目前關(guān)于準(zhǔn)三維水動力模型開放型代碼較多，其中相對成熟且應(yīng)用較多的是EFDC。本文水動力學(xué)語言模塊(HYD)由Intel Visual Fortran 11.0.061專業(yè)版編譯器雙精度Fortran90語言編寫。

1.2 準(zhǔn)三維水動力精度分析

河道形態(tài)如圖1所示，彎道的底坡為0.0003。試驗流量為4.17 L/s、下游水深10.59 cm[16]。模型進口紊動動能k=0.02U2(U為垂向平均流速)，紊動黏性系數(shù)為νt=0.001 m2/s；模型橫向網(wǎng)格數(shù)10個，縱向網(wǎng)格數(shù)360個，垂向網(wǎng)格數(shù)10個。

理論分析已表明，由于垂向動量方程直接采用靜壓假定，因此準(zhǔn)三維模型無法準(zhǔn)確計算垂向流速(w)，不適合于垂向流速變化較大的情況，也就是說準(zhǔn)三維模型無法準(zhǔn)確計算涉及垂向流速的二次流。本文分別采用三維水動力和準(zhǔn)三維水動力方程，對圖1彎道水動力進行計算，驗證彎道中水流流速大小沿垂向的變化，結(jié)果表明，準(zhǔn)三維模型和三維模型流速大小相差甚小。圖2給出第三彎道45°和75°斷面流速沿垂向分布對比。由圖2可看出全三維、準(zhǔn)三維模型數(shù)值計算值與實測值基本符合。只是數(shù)值計算流速在近床面明顯降低，沿垂向呈對數(shù)分布，而實測流速在近床面略有波動，甚至出現(xiàn)近底層流速相對較大的現(xiàn)象，這可能由于彎道位置斷面螺旋流所致。模型驗證表明，雖然基于靜壓假定，但在動壓作用非主導(dǎo)條件下，準(zhǔn)三維模型流速計算結(jié)果與全三維模型基本一致。準(zhǔn)三維水動力模型可以為SB3D提供可靠的水流流速信息。

圖1 連續(xù)彎道模型平面布置圖

圖2 典型斷面流速沿水深分布驗證

2 三維泥沙動力模型及高效應(yīng)用

目前三維水動力模型和泥沙動力模型都采用分離求解。在每一時間步內(nèi)，先求解水動力學(xué)方程獲得un+1、vn+1和ωn+1，后求解泥沙輸運方程和河床變形方程獲得泥沙濃度(sn+1)、水位高程(zn+1)等(見圖3)。un+1、vn+1和ωn+1得出后，泥沙動力模塊自身封閉。若該泥沙動力模型的un+1、vn+1和ωn+1由準(zhǔn)三維水動力模型提供，則稱為準(zhǔn)三維河流動力學(xué)模型(見圖3)。

(1) 懸沙輸運模型。無論是全三維模型還是準(zhǔn)三維模型，懸沙輸運方程均采用單相多組分流體輸運方程[2-6]，只是在垂向?qū)α黜椫性黾恿四嗌炒瓜虺了佴豷：

(6)

式中，ωs為σ坐標(biāo)下的泥沙沉速；Qs為泥沙源項。

圖3 全三維河流動力學(xué)模型和準(zhǔn)三維河流動力學(xué)模型的轉(zhuǎn)換

由于單相多組分控制方程應(yīng)用的廣泛性，多數(shù)準(zhǔn)三維模型程序語言都包含相關(guān)模塊，只需要在程序代碼中添加垂向沉速，并修改床面邊界條件以及相應(yīng)參數(shù)即可用于懸沙計算。床面邊界可采用

(7)

式中，νt是紊動擴散系數(shù)；σs為泥沙的Schmidt數(shù)；sb和sb*分別為近底層泥沙濃度和飽和濃度。

(2) 河床變形模型。河床縱向變形模型的一般形式

(8)

式中

(9)

其中角標(biāo)l表示粒級級配分組；lM為泥沙級配組數(shù)；p′為泥沙的空隙比；δb為推移質(zhì)泥沙運動層的厚度；qb,l,x和qb,l,y分別在x和y方向泥沙推移通量。

不同文獻中上式各項具體形式不同，但具體計算都是采用顯式求解，因此全三維模型研究成果，都可直接用于準(zhǔn)三維模型，且程序代碼的實現(xiàn)也可統(tǒng)一。此外，河床橫向演變模型也都采用顯式求解[1,3-6]，因此同樣可直接應(yīng)用于準(zhǔn)三維模型。

3 模型工程應(yīng)用

3.1 工程概況

北本水電站位于湄公河上游河段，為湄公河水電開發(fā)規(guī)劃的第一級電站，電站位于老撾北部烏多姆賽省北本縣境內(nèi)，壩址在北本縣城上游約14 km處。電站采用徑流式開發(fā)，樞紐工程由混凝土重力壩、泄洪沖沙閘、沖沙底孔、船閘和魚道等建筑物組成(見圖4)。水庫為日調(diào)節(jié)水庫，水庫正常蓄水位340 m，死水位334 m，總庫容約7.8×108m3，電站裝機容量912 MW。

圖4 電站及下游航道分布平面示意圖

根據(jù)規(guī)劃，當(dāng)入庫流量小于3 a一遇洪水流量13 200 m3/s時，水庫水位盡量維持在340 m運行，通過泄洪沖沙閘(12孔)、沖沙底孔來維持；當(dāng)入庫流量大于3 a一遇洪水流量13 200 m3/s、小于5 a一遇洪水流量14 900 m3/s時，船閘停止使用，水庫水位盡量維持在340 m運行，通過泄洪沖沙閘(12孔)、沖沙底孔和機組過流、航道沖沙閘打開泄水(2孔)來調(diào)節(jié)。電站推移質(zhì)泥沙大都被上游電站攔截，懸沙中值粒徑為0.007 7 mm，電站下游床沙中值粒徑0.025 mm。

3.2 離散網(wǎng)格和邊界條件

計算區(qū)域上游選自大壩位置，出口取至引航道以下500 m位置。模型采用了曲線正交網(wǎng)格，縱向網(wǎng)格數(shù)為129個，橫向網(wǎng)格數(shù)為60個，垂向分為5層。網(wǎng)格劃分主要考慮主流流向以及水流、泥沙參數(shù)的變化。網(wǎng)格與引航道邊界擬合很好，主流區(qū)網(wǎng)格走向與河道深泓線一致，網(wǎng)格間距為10 m～20 m，引航道區(qū)域網(wǎng)格局部加密。區(qū)域網(wǎng)格夾角為80°～90°。時間步長為10 s。計算區(qū)域上游邊界給定水庫下泄流量和沙量，下游給定水位(見圖5)。

圖5 下游出口水位變化過程曲線

基于所建立準(zhǔn)三維模型，計算2003年—2009年水庫修建之前的河道沖淤過程，參數(shù)率定表明，模型飽和系數(shù)采用常數(shù)，沖刷時取2.3，淤積時取0.55，床面粗糙高度ks=0.01 m。

3.3 引航道表層橫向流速和通航安全分析

北本水電站通航建筑物為單線單級船閘，船閘級別為Ⅳ級，要求口門區(qū)垂直航線的橫向流速不能大于0.30 m/s，回流區(qū)不能大于0.40 m/s。本文選取通航最低流量、年平均流量、2 a一遇流量以及通航最大流量四種條件計算引航道的橫向流速。計算結(jié)果可以清楚的反映表層和底層的流速差異(見圖6)，表層流速可為航道通航安全提供有力依據(jù)。引航道回流區(qū)位于航道內(nèi)口門上游100 m～300 m處，如圖7中C#區(qū)。航道內(nèi)橫向流速隨著河道下泄流量增加而增大。在最低通航流量和年平均流量時，C#區(qū)橫向流速為0.10 m/s左右，對船只的安全航行無不利影響，回流流速小于0.40 m/s，尚能滿足通航要求，但已有礙航趨勢；當(dāng)流量達到2 a一遇流量11 600 m/s時，航道內(nèi)C#區(qū)最大橫向流速達0.50 m/s～0.60 m/s，超過口門區(qū)最大橫向流速限值0.40 m/s，無法滿足通航條件。上述計算結(jié)果表明，相較于二維模型，準(zhǔn)三維模型可以方便的提供引航道內(nèi)表層流速，進一步判斷不同水流條件下船只航行的安全性。

圖6 計算全域三維流場分布

3.4 電站下游泥沙沖淤計算

模型基于Dell/戴爾T420服務(wù)器(六核E5-2430V2/16G/2T)，對電站運行后下游泥沙沖淤進行計算。電站下游主河道及引航道沖淤示意圖如圖8所示。計算表明，電站運行后下游河道沿深泓線均有不同程度的沖刷，斷面沖刷由深泓線向兩岸遞減。最大沖深位于壩下游300 m處深槽位置，水庫運行1 a后，最大沖深達1.0 m左右，后沖刷的速度逐漸減緩，運行1 a～5 a內(nèi)，該點年平均沖刷深度0.5 m，至第5年時沖淤基本平衡，河道深泓線平均沖刷1.4 m～2.0 m左右，局部達3.0 m。當(dāng)河道泄洪沖沙閘關(guān)閉時，泄洪沖沙閘下游形成局部回流和淤積，淤積厚度約為0.7 m。在電站運行10 a內(nèi)，下引航道呈累積性淤積趨勢。淤積程度由口門向船閘呈遞減趨勢，5 a末船閘附近淤積厚度約為1.0 m左右，口門最大淤積厚度達2.6 m左右，將會出現(xiàn)礙航現(xiàn)象。以上計算表明，準(zhǔn)三維模型不僅可以區(qū)分垂向速度差異，而且可以成功模擬電站下游河道長期沖淤過程，對比全三維模型在計算效率方面具有更大的優(yōu)勢。計算機運行時間約為214 h，模型適用于中小型的工程的應(yīng)用。

圖7 引航道流速分布

圖8 電站下游主河道及引航道沖淤示意圖

4 結(jié) 論

近年來泥沙輸運和河床沖淤計算方法取得了較快的發(fā)展，但這些模型的水動力計算多是采用全三維水動力模型，計算量過大而限制了成果的工程實用性。本文研究表明：

(1) 準(zhǔn)三維水動力模型可以極大的降低模型的計算量，同時較為準(zhǔn)確的模擬非動水壓強主導(dǎo)的三維水流運動特征。

(2) 全三維河流動力模型中的水動力模型和模塊，可方便的被準(zhǔn)三維模型或模塊取代，變?yōu)闇?zhǔn)三維模型，增加其工程應(yīng)用可行性。

(3) 準(zhǔn)三維水動力模型可以成功應(yīng)用于中長尺度河流模擬，特別是對于部分既需要考慮三維垂向流速分布，計算時空尺度又較大的泥沙動力模擬。

[1] Olsen N R B, Kjellesvig H M. Three-dimensional numerical flow modeling for estimation of maximum local scour depth[J]. Journal of Hydraulic Research, 1998,36(4):579-590.

[2] 黑鵬飛,假冬冬,冶運濤,等.計算河流動力學(xué)理論體系框架探討[J].水科學(xué)進展,2016,27(1):152-164.

[3] 黃國鮮,周建軍,吳偉華.彎曲河道螺旋流作用下的物質(zhì)輸運三維模擬[J].清華大學(xué)學(xué)報(自然科學(xué)版),2008，48(6):977-982.

[4] 假冬冬,邵學(xué)軍,王 虹,等.考慮河岸變形的三維水沙數(shù)值模擬研究[J].水科學(xué)進展,2009,20(3):311-317．

[5] 假冬冬,黑鵬飛,邵學(xué)軍,等.分層岸灘側(cè)蝕坍塌過程及其水動力響應(yīng)模擬[J].水科學(xué)進展,2014,25(1):83-89.

[6] 胡德超,鐘德鈺,張紅武,等.三維懸沙模型及河岸邊界追蹤方法:Ⅱ—河岸邊界追蹤[J].水力發(fā)電學(xué)報,2010,29(6):106-113.

[7] Motta D, Abad J D, Langendoen E J, et al. A simplified 2-D model for meander migration with physically-based bank evolution[J]. Geomorphology, 2012,163/164:10-25.

[8] Tang C Y, Li Y P, Acharya K. Modeling the effects of external nutrient reductions on algal bloomsin hyper-eutrophic Lake Taihu, China[J]. Ecological Engineering, 2016，94:164-173.

[9] Chen L B, Yang Z F, Liu H F. Assessing the eutrophication risk of the Danjiangkou Reservoir based on the EFDC model. Ecol. Eng.[EB/OL]. http://dx.doi.org/10.1016/j.ecoleng.2016.02.021.

[10] 方紅衛(wèi),王光謙.平面二維全沙泥沙輸移數(shù)學(xué)模型及其應(yīng)用[J].應(yīng)用基礎(chǔ)與工程科學(xué)學(xué)報,2000,8(2):165-178.

[11] 夏軍強,宗全利,許全喜,等.下荊江二元結(jié)構(gòu)河岸土體特性及崩岸機理[J].水科學(xué)進展,2013,24(6):810-820.

[12] 周 剛,王 虹,邵學(xué)軍,等.河型轉(zhuǎn)化機理及其數(shù)值模擬—Ⅰ模型建立[J].水科學(xué)進展,2010,21(2):145-152.

[13] 鐘德鈺,張紅武.考慮環(huán)流橫向輸沙及河岸變形的平面二維擴展數(shù)學(xué)模型[J].水利學(xué)報,2004,35(7):14-20.

[14] 李義天,尚全民.一維不恒定流泥沙數(shù)學(xué)模型研究[J].泥沙研究,1998(1):81-87.

[15] 方紅衛(wèi),王光謙.一維全沙泥沙輸移數(shù)學(xué)模型及其應(yīng)用[J].應(yīng)用基礎(chǔ)與工程科學(xué)學(xué)報,2000,8(2):154-164.

[16] Wang B, Jia D D, Zhou G, et al. An experimental investigation on flow structure in channel with consecutive bends[C]//Proceedings of 16th IAHR-APD Congress and 3rd Symposium of IAH R-ISHS. Berlin: Springer, 2009:1811-1816.

3D River Fluvial Numerical Model Optimized by the Quisi-3D Hydrodynamic Numeric Model

YANG Jing1, HEI Pengfei1, JIA Dongdong2, SHANG Yizi3

(1.CollegeofLifeandEnvironmentalSciences,MinzuUniversityofChina,Beijing100081,China;2.StateKeyLaboratoryofHydrology-WaterResourcesandHydraulicEngineering,NanjingHydraulicResearchInstitute,Nanjing,Jiangsu210029,China;3.DepartmentofWaterResources,ChinaInstituteofWaterResourcesandHydropowerResearch,Beijing100038,China)

Quasi-3D numerical model has been widely used in the water flow modeling and water quality modeling, but its significance in the river fluvial computation has not be fully realized. In this paper, the quasi-3D hydrodynamics model and fluvial model used in the 3D model were combined. After the model calibration, the model were then used to simulate the hydrodynamic and fluvial process in downstream of the Beiben power station. Results manifested that quasi-3D numerical model for the river dynamics shows great advantage than the 3D numerical model in the computation efficiency while retaining the capability for simulating the vertical distribution of the flow and sediment dynamics. The quasi-3D numerical model plays the irreplaceable role when the vertical distribution of the flow and sediment are needed while 3D was not feasible due to the computation efficiency.

quasi-3D numerical model, compuational river dynamics; sediment

10.3969/j.issn.1672-1144.2017.02.002

2016-12-11

2017-02-01

中央民族大學(xué)校級自主科研項目(2016SHXY03)；國家自然科學(xué)基金項目(51209239)

楊 靜(1982—)，女，內(nèi)蒙古赤峰人，博士生，研究方向為水生態(tài)。E-mail:yangjing123456@126.com

黑鵬飛(1979—)，男，陜西榆林人，講師，主要從事地表水環(huán)境和水生態(tài)研究。E-mail:heipf06@mails.tsinghua.edu.cn

TV143

A

1672—1144(2017)02—0009—07