夏慶福，余弘婧，朱 銳，章少輝，郭新蕾

（1.流域水循環(huán)模擬與調(diào)控國家重點實驗室 中國水利水電科學研究院，北京 100038；2.南水北調(diào)工程建設(shè)監(jiān)管中心，北京 100038）

1 研究背景

Saint-Venant方程組由雷諾平均的Navior-Stokes方程組沿垂向積分平均獲得，是描述地表淺水運動的基本方程［1-5］。在地表淺水流的數(shù)值模擬中，為獲得優(yōu)良的物理（質(zhì)量與動量）守恒性和有效捕捉淺水波（間斷波與稀疏波），在空間離散格式方面往往基于Saint-Venant方程組的雙曲守恒形式建立通量差分分裂（Flux Difference Splitting，F(xiàn)DS）類型的有限體積法［5-6］，例如Roe格式和HLL/HLLC格式等，其具有物理意義明確清晰的優(yōu)點，但結(jié)構(gòu)比較復雜，耗散項具有顯著的向量特征，致使效率較低［7-8］。與之相應(yīng)，時間格式方面多采用顯式，但嚴格的穩(wěn)定性條件導致了低效率。雖然已有學者建立了FDS類型有限體積法的全隱時間格式［9］，但空間離散格式的結(jié)構(gòu)復雜和較低效率，使其無法發(fā)揮隱時間格式高效率的優(yōu)點。Liou等指出［10］，在數(shù)值求解雙曲類方程時，隱時間格式適宜與結(jié)構(gòu)簡潔且高效率（耗散項具有典型的標量特征）的矢通量分裂（Flux Vector Splitting，F(xiàn)VS）類型的有限體積法相結(jié)合，從而達到高精度和高效率的統(tǒng)一。

本文基于Saint-Venant方程組的守恒形式，針對對流通量項建立了FVS類型的標量耗散有限體積法。通過增加代數(shù)項，使得地表水位相對高程梯度項的中心差分格式在整個計算區(qū)域內(nèi)均能嚴格地滿足靜水條件下其值為零的物理條件。借助雙時間步法，實現(xiàn)了Saint-Venant方程組空間離散的無條件穩(wěn)定全隱式解法。最后借助2個典型算例，驗證解法的模擬效果。

2 控制方程

Saint-Venant方程組包括質(zhì)量與動量守恒方程，其守恒形式表達如下［11］：

式中：A為過流斷面面積（m2）；Q為通過任意斷面的流量（m3/s）；z=zb+h地表水位相對高程，且zb為地表相對高程，h為地表水深（m），t為時間坐標（s）；x為空間坐標（m）；g為重力加速度（m/s2）；R為水力半徑（m）。

數(shù)值求解時往往采用式（1）和（2）的向量表達式：

式中：U為因變量向量；F為對流通量；Sζ和Sf，分別為地表水位相對高程梯度向量和地表糙率向量。各向量的表達式如下：

3 數(shù)值解法

3.1空間離散在任意空間單元格i上對式（3）進行空間積分如下：

式中的下標（i+1/2）和（i-1/2）分別用于標記單元格i的左右兩側(cè)邊界。

3.1.1 物理變量的二階精度重構(gòu) 對地表水位變量z在任意單元格i邊界（i+1/2）兩側(cè)的狀態(tài)值zi+1/2，L和zi+1/2，R進行空間重構(gòu)如下：

采用 min mod限制器［12］計算式（6）中的Δζi和Δζi+1。

對地表水流速度u、水深h等其它物理變量在任意單元格i邊界（i+1/2）、以及其他邊界兩側(cè)的相應(yīng)狀態(tài)值進行空間重構(gòu)的過程與地表水位相同。

單元格i邊界（i+1/2）兩側(cè)的地形相對高程值計算如下［12］：

若單元格無水、或有水單元格緊鄰無水區(qū)域，則其邊界變量值采用單元格中心值。

3.1.2 物理變量的黎曼狀態(tài)重構(gòu) 基于各淺水運動物理變量的二階精度重構(gòu)值，構(gòu)造黎曼狀態(tài)的目的，是正確計算通量值并捕捉間斷波（包括干濕邊界）［12］。地形相對高程、地表水深、水位和單寬流量的黎曼狀態(tài)分別重構(gòu)如下：

3.1.3 對流通量梯度項的數(shù)值離散 利用散度定理于式（5）中的對流通量梯度積分項，可獲得如下表達式：

式中：Fi+1/2和Fi-1/2為通過單元格邊界（i+1/2）和（i-1/2）的因變量對流數(shù)值通量。

Fi+1/2標量耗散有限體積法表達如下：

式中：ci+1/2為地表淺水的重力波速（m/s）。

借助Liou［10］在空氣動力學中構(gòu)造的AUSM格式，將式（13）中的ci+1/2和Fri+1/2分別定義如下：

基于式（13），式（12）中的對流數(shù)值通量差值表達如下：

式（18）中的系數(shù)表達式如下：

3.1.4 水位梯度項的數(shù)值離散 通過改進借助Liang等［12］成果，地表水位相對高程梯度空間離散格式如下：

式（23）中右側(cè)第二和第三項分別表達如下：

在無水區(qū)域，式（23）中的右側(cè)所有項恰好等于零，這正確描述了無水區(qū)域內(nèi)沒有水位梯度驅(qū)動力的事實。

3.1.5 地表糙率向量的空間離散 在任意空間單元格i上地表糙率向量表達如下：

3.2時間離散采用二階精度時間離散格式對式（5）的空間離散表達式處理如下：

為無條件穩(wěn)定地數(shù)值求解式（29）的目的，采用雙時間步法對其處理如下：

式中：上標“p”為虛擬時間迭代步；Δτ為虛擬時間步長。

經(jīng)過合并同類項整理后，式（30）表達如下：

3.3初始和邊界條件設(shè)置無論真實時間步長Δt取何值，選取適當?shù)奶摂M時間步長Δτ，使得式（31）的系數(shù)矩陣始終保持對角占優(yōu)，則該式即能始終保持收斂，達到無條件穩(wěn)定。然而，較大時間步長意味著較大的數(shù)值耗散，故真實時間步長應(yīng)根據(jù)物理問題的特征進行選取。

初始條件分為給定水深與流速、和無水深與流速兩種情景。由于無水深h=0是糙率項的奇點，故假設(shè)初始水深 h0=10-10m［12］。

邊界條件包括入流、出流和無流三種條件。三種邊界條件均借助鏡像單元格方法給定。無流邊界條件下計算域邊界單元格的外邊界流速和水位梯度均設(shè)置為零。入流和出流邊界條件，則依據(jù)流態(tài)是否為亞臨界或超臨界流態(tài)，分別給定相應(yīng)的物理和數(shù)值邊界條件，可詳見文獻［12］。由此設(shè)置的邊界條件滿足二階精度［12］。

4 數(shù)值模擬驗證

選取文獻中2個標準算例，驗證上述建立的數(shù)值解法的模擬效果。任意算例均采用時間步長呈倍數(shù)遞減的三組模擬結(jié)果與解析或?qū)崪y值進行對比，以達到驗證模擬結(jié)果的時間收斂性［13］。

4.1矩形斷面下河床有水的潰壩過程若不考慮矩形河床糙率，描述河床有水潰壩過程的Saint-Ve?nant方程存在解析解［14］。假設(shè)河床長度2 000 m，寬度1 m，在距離上游1 000 m的位置有高度1m的水壩，壩后水深1 m，下游河床水深0.1 m。該算例用于驗證數(shù)值解法捕捉激波和稀疏波、以及模擬水流由亞臨界態(tài)向超臨界態(tài)演變的能力。

圖1 矩形斷面河床有水時潰壩水位解析解與模擬結(jié)果比較（t=200s）

采用尺寸Δx=1 m的單元格離散計算區(qū)域，圖1給出了模擬結(jié)果與解析值之間的比較?？梢钥闯觯煌瑔卧癯叽缦?，模擬結(jié)果與解析解之間均具有良好的擬合度；隨著時間步長以倍數(shù)遞減，模擬值迅速收斂至解析解，表明數(shù)值解法具有優(yōu)良的模擬效果和收斂性。

4.2矩形水槽斷面急劇收縮的潰壩過程矩形水槽的坡度為0.005，長度為122 m，寬度為1.22m［14］。糙率系數(shù)為0.009 s/m1/3。沿水槽縱向中點處有一個高0.305 m的水壩，壩前水位與壩高持平。在距離水槽上游30.5、45.75、68.625和106.75 m的位置設(shè)置測點，測點名稱分別為STA100、STA150、STA225和STA350。在這4個測點處觀測水位變化過程，并在STA225和STA350觀測淺水流的流速變化過程。

水壩突然在中間潰開一個寬0.732 m、高0.183 m的缺口，故本算例屬于典型的橫斷面急劇收縮的復雜物理問題，對數(shù)值解法的模擬效果提出了更嚴格驗證。

單元格尺寸Δx=0.1、0.2和0.4 m時，圖2給出了模擬結(jié)果與實測值之間的比較。在結(jié)果比較時，若單元格中心與測點坐標不重合，則采用線性插值的方式獲得觀測點處的水位和流速值。

從圖2可以看出，4個測點的水位模擬值與實測值之間具有良好的擬合精度和收斂性。與之相比，流速擬合稍差，但仍能模擬出流速的變化特征。這表明，建立的模型能比較準確地模擬出橫斷面急劇變化下潰壩產(chǎn)生的激波變化和地表水的逐漸消退過程。

圖2 矩形水槽斷面急劇收縮下潰壩過程中不同測點模擬結(jié)果與實測值比較

5 結(jié)論

基于Saint-Venant方程組的守恒形式，重構(gòu)了各物理變量在單元格邊界的黎曼狀態(tài)值，實現(xiàn)了各變量在計算區(qū)域內(nèi)的二階精度分布?；谝蜃兞繉α魍渴歉等甑聰?shù)函數(shù)的事實，構(gòu)造了因變量對流通量的標量耗散有限體積法，具有結(jié)構(gòu)簡潔、計算效率和模擬精度高的優(yōu)點。通過增加構(gòu)造的代數(shù)項，使地表水位相對高程梯度項的中心差分離散格式能適用于有水、無水和干濕邊界等整個計算區(qū)域。采用雙時間步法，實現(xiàn)了Saint-Venant方程組空間離散式的全隱時間格式離散，實現(xiàn)了無條件穩(wěn)定求解。

選取文獻中2個標準算例，系統(tǒng)驗證了上述建立的數(shù)值解法的模擬效果。任意算例均數(shù)量呈倍數(shù)遞減的三個時間步長進行數(shù)值模擬，模擬結(jié)果與解析解及觀測值的對比表明，數(shù)值解法能以優(yōu)良的精度模擬出不同斷面形式下的潰壩過程，模擬結(jié)果表現(xiàn)出了良好的收斂性。

［1］HUBBARD M E，GARCIA-NAVARRO P.Flux difference splitting and the balancing of source terms and flux gradients［J］.Journal of Computational Physics，2000，165（1）：89-125.

［2］LIANG Q.Flood simulation using a well-balanced shallow flow model［J］.Journal of Hydraulic Engineering，2010，136（9）：669-675.

［3］WANG Y，LIANG Q，KESSERWANI G，et al.A 2D shallow flow model for practical dam-break simulations［J］.Journal of Hydraulic Research，2011，49（3）：307-316.

［4］LIANG D，LIN B，F(xiàn)ALCONER R A.A boundary-fitted numerical model for flood routing with shock-capturing capability［J］.Journal of Hydrology，2007，332（3）：477-486.

［5］HOU J，LIANG Q，SIMONS F，et al.A 2D well-balanced shallow flow model for unstructured grids with novel slope source term treatment［J］.Advances in Water Resources，2013，52：107-131.

［6］BENKHALDOUN F，ELMAHI I，SEAID M.Well-balanced finite volume schemes for pollutant transport by shal?low water equations on unstructured meshes［J］.Journal of Computational Physics，2007，226：180-203.

［7］BENKHALDOUN F，SAHMIM S，SEAID M.A two-dimensional finite volume morphodynamic model on unstruc?tured triangular grids［J］.International Journal of Numerical Methods Fluids，2010，63（11）：1296-1327.

［8］SIVAKUMAR P，HYAMS D，TAYLOR L，et al.A primitive-variable Riemann method for solution of the shal?low water equations with wetting and drying［J］.Journal of Computational Physics，2009，228（19）：7452-7472.

［9］YU H，HUANG G，WU C.Efficient finite-volume model for shallow-water flows using an implicit dual time-stepping method［J］.Journal of Hydraulic Engineering，2015，141（6）：04015004.

［10］LIOU M S，STEFFEN Jr C J.A new flux splitting scheme［J］.Journal of Computational Physics，1993，107（1）：23-39.

［11］XING Y，ZHANG X，SHU C-W.Positivity-preserving high order well-balanced discontinuous Galerkin meth?ods for the shallow water equations［J］.Advances in Water Resources，2010，33（12）：1476-1493.

［12］LIANG Q，MARCHE F.Numerical resolution of well-balanced shallow water equations with complex source terms［J］.Advances in Water Resources，2009，32（6）：873-884.

［13］SHAO S，LO E Y M.Incompressible SPH method for simulating Newtonian and non-Newtonian flows with a free surface［J］.Advances in Water Resources，2003，26（7）：787-800.

［14］YING X，KHAN A A，WANG S S Y.Upwind conservative scheme for the Saint Venant equations［J］.Journal of Hydraulic Engineering，2004，130（10）：977-987.