Saint-Venant方程組全隱式標(biāo)量耗散有限體積法

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

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

1 研究背景

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

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

2 控制方程

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

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

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

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

3 數(shù)值解法

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

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

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

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

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

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

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

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

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

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

Fi+1/2標(biāo)量耗散有限體積法表達(dá)如下：

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

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

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

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

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

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

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

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

3.2時(shí)間離散采用二階精度時(shí)間離散格式對(duì)式（5）的空間離散表達(dá)式處理如下：

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

式中：上標(biāo)“p”為虛擬時(shí)間迭代步；Δτ為虛擬時(shí)間步長(zhǎng)。

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

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

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

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

4 數(shù)值模擬驗(yàn)證

選取文獻(xiàn)中2個(gè)標(biāo)準(zhǔn)算例，驗(yàn)證上述建立的數(shù)值解法的模擬效果。任意算例均采用時(shí)間步長(zhǎng)呈倍數(shù)遞減的三組模擬結(jié)果與解析或?qū)崪y(cè)值進(jìn)行對(duì)比，以達(dá)到驗(yàn)證模擬結(jié)果的時(shí)間收斂性［13］。

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

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

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

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

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

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

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

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

5 結(jié)論

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

選取文獻(xiàn)中2個(gè)標(biāo)準(zhǔn)算例，系統(tǒng)驗(yàn)證了上述建立的數(shù)值解法的模擬效果。任意算例均數(shù)量呈倍數(shù)遞減的三個(gè)時(shí)間步長(zhǎng)進(jìn)行數(shù)值模擬，模擬結(jié)果與解析解及觀測(cè)值的對(duì)比表明，數(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.