楊金波，段文剛，黃國兵，李訂芳

（1．長江科學(xué)院水力學(xué)研究所，武漢430010；2．武漢大學(xué)數(shù)學(xué)與統(tǒng)計學(xué)院，武漢430072）

一種Godunov算法及其在非平底流污染物傳輸中的應(yīng)用

楊金波1，段文剛1，黃國兵1，李訂芳2

（1．長江科學(xué)院水力學(xué)研究所，武漢430010；2．武漢大學(xué)數(shù)學(xué)與統(tǒng)計學(xué)院，武漢430072）

采用和諧的加權(quán)平均通量（Weighted Average Flux，簡稱WAF）算法，研究了淺水波方程的間斷解及污染物傳輸問題。該算法采用WAF格式和HLLC Riemann求解器近似單元邊界數(shù)值通量，中心差分格式離散地形源項，然后理論上證明了該算法是和諧的。最后利用WAF算法對非平底地形上淺水波間斷解及污染物傳輸問題進(jìn)行數(shù)值計算，精確地捕捉到了間斷解和污染物運(yùn)動過程，結(jié)果表明該算法滿足守恒性，具有高分辨率、無振蕩及捕捉污染物運(yùn)動邊界的能力。

淺水波方程；WAF格式；潰壩；污染物傳輸

1 研究背景

河流、湖泊和近海岸的污染問題越來越嚴(yán)重，一方面破壞了生態(tài)環(huán)境；另一方面對人類生命安全和經(jīng)濟(jì)發(fā)展帶來致命性的危害，因此有必要研究污染物傳輸方式的數(shù)學(xué)模型。

污染物傳輸?shù)臄?shù)值算法包括淺水波和污染物傳輸擴(kuò)散的計算。在淺水波方程計算方面，近年來空氣動力學(xué)中的Godunov［1］型有限體積算法越來越普遍地應(yīng)用到淺水波方程中，WAF方法［2］是Toro基于Godunov思想提出的顯示、二次精度的守恒有限體積格式。非平底河床是淺水方程計算的一個關(guān)鍵問題，為了構(gòu)造和諧的數(shù)值格式，Roe［3］提出迎風(fēng)離散雙曲型方程中源項的思想；Bermudez和Vazquez Cendon［4］用迎風(fēng)離散源項的方法計算了淺水波方程；Hubbard和Garcia Navarro［5］用迎風(fēng)離散源項的方法構(gòu)造了平衡源項和通量梯度的二階格式；Zhou［6］等人對二維淺水波方程提出了水位梯度法。對于污染物傳輸部分，考慮對流占優(yōu)的問題，文章采用垂直方向平均的傳輸方程，忽略了擴(kuò)散項。數(shù)值求解時，污染物傳輸方程作為整個系統(tǒng)的一個波系，同淺水方程耦合求解。

本文首先構(gòu)造了淺水波方程和污染物擴(kuò)散方程的二階精度和諧WAF格式，然后用該算法研究了非平底河床上淺水流的間斷解及污染物傳輸問題，將計算結(jié)果同精確解或其他文獻(xiàn)比對，證明該方法的計算結(jié)果精度高、分辨率高，能計算真實(shí)的潰壩和污染物傳輸問題。

2 數(shù)值格式

2．1 淺水波及污染物傳輸方程

控制方程由淺水波方程和污染物運(yùn)輸方程組成，其中，污染物部分主要研究污染物傳輸過程，忽略污染物擴(kuò)散項，如

式中：H＝h＋b（x）為水位高程；h為水深，b（x）為河床高程；q＝hu為單寬流量；u為流速，傳輸標(biāo)量r＝hc；c為污染物濃度；g為重力加速度。

2．2 WAF有限體積格式

對式（1）關(guān)于區(qū)間［xi－1／2，xi＋1／2］×［tn，tn＋1］積分，并應(yīng)用格林（Green）公式，則得有限體

項，這里采用TVD WAF格式給出數(shù)值通量，中心差分格式離散源項。如圖1所示，假設(shè)局部Riemann問題自相似解僅有激波結(jié)構(gòu)，TVD WAF格式［7］的數(shù)值通量為

圖1 局部Riemann問題自相似解的波結(jié)構(gòu)Fig．1 W ave structure of self sim ilar solution for local Riemann problem

其中：

式（4）中變量p可以為任意一個在波k處變化的物理量，這里取p（k）＝H，k＝1，3…；p（2）＝r，k＝2。

為了計算FWAFi＋1／2，采用HLLC解法器［6］來近似Riemann問題自相似解的中間狀態(tài)，

下標(biāo)（k）表示通量向量的第k項，其中HLL求解器定義為

左右波速度S1和S3通過文獻(xiàn)［7］給出，即

式中

當(dāng)出現(xiàn)干濕邊界時S1，S3分別調(diào)整：

左干河床

右干河床

此時，式（8）中的S3會自動和估計的干河床式（9）或式（10）重合。

為保證數(shù)值格式的穩(wěn)定性，時間步長

其中CFL（0＜CFL≤1）是Courant數(shù)。

地形源項采用中心差分格式進(jìn)行離散

2．3 和諧性質(zhì)證明

為了證明WAF格式的和諧性，只需證明靜態(tài)條件q＝0 m，H＝常數(shù)及r＝常數(shù)時，下式成立

由TVD WAF格式的數(shù)值通量定義式（3）和靜水條件，單元Ii右邊界xi＋1／2處的數(shù)值通量為

將式（14），式（15）和式（16）分別代入式（13）很清楚地證明格式滿足和諧性。

3 算 例

這里利用WAF格式對非平底間斷淺水波方程中污染物傳輸問題進(jìn)行了數(shù)值計算，研究了潰壩波傳播過程中污染物的傳輸過程和機(jī)理，試驗結(jié)果表明該算法可以精確地捕捉激波和污染物傳輸邊界，無振蕩。

3．1 梯形河床上的污染物傳輸問題

本算例［8］梯形河床定義為

計算步長Δx＝0．08 m，CFL＝0．9。T＝1 s時計算結(jié)果如圖2所示。WAF算法計算的水位和HLLC格式計算的水位吻合較好，同時由圖3可知，WAF格式對污染物運(yùn)動的計算比一階精度的HLLC格式分辨率要高，WAF算法的計算數(shù)值耗散低。

圖2 水面線Fig．2 W ater surface profile

圖3 污染物分布Fig．3 Pollutant distribution

3．2 潰壩波中的污染物傳輸問題

本算例研究了潰壩過程中污染物傳輸規(guī)律，計算區(qū)域為2 000 m長的平底無摩擦河道，潰壩初始條件如下：

采用100個計算單元，CFL＝0．9。T＝250 s時計算水位和單寬流量如圖4所示，可以看出WAF格式比HLLC格式計算的更逼近解析解。在污染物計算方面，如圖5，WAF計算的污染物區(qū)域比HLLC格式窄，濃度更高，說明WAF比HLLC在數(shù)值計算方面的耗散低，精度高。

圖4 水位和單寬流量Fig．4 W ater level and discharge per unit w idth

圖5 污染物中心分布Fig．5 Distribution of pollutant center

3．3 帶干濕河床的污染物傳輸問題

該算例利用WAF算法研究干濕河床上的污染物傳輸問題［9］，初始條件為：

其中地形定義為

空間步長Δx＝0．1 m，CFL＝0．8。t＝0．25 s時，計算結(jié)果如圖6所示。WAF格式計算的水位能很好的和參考值吻合（該算例沒有解析解，可采用HLLC格式在2 000個網(wǎng)格下的計算結(jié)果作為參考值）。并且精確地捕捉到阻礙物左邊x＝25／3 m處的激波和x＝6．5 m處的水面涌起。同時，分析t＝0．18 s時污染物運(yùn)動狀態(tài)，如圖7?？梢钥闯鯳AF比HLLC格式精度高。進(jìn)一步，由初始條件產(chǎn)生的向左和向右運(yùn)動的2個稀疏波中間形成干河床，改進(jìn)后的WAF算法可以很好處理干濕邊界問題。該算例說明WAF算法具有高精度、高分辨率，及能處理復(fù)雜干河床問題，具有真實(shí)問題的處理能力。

圖6 T＝0．25 s時刻水面線Fig．6 W ater surface profile at T＝0．25 s

圖7 T＝0．18 s時污染物分布Fig．7 Pollutant distribution at T＝0．18 s

4 結(jié) 語

本文利用二階精度和諧WAF格式研究了淺水波方程和污染物擴(kuò)散問題，研究了污染物在潰壩、非平底河床情況下的傳播機(jī)理計算結(jié)果同精確解或其他文獻(xiàn)比對，證明該方法是高精度、高分辨率，能計算真實(shí)的潰壩和污染物傳輸問題。

［1］ GODUNOV SK．Finite Difference Methods for the Com putation of Discontinuous Solutions of the Equations of Fluid Dynamics［J］．MatematicheskiǐSbornik，1959，47：271－306．

［2］ TORO E F．AWeighted Average Flux Method for Hyper bolic Conservation Laws［J］．Proceedings of the Royal Society A：Mathematical，Physical and Engineering Sci ences，1989，423（1865）：401－418．

［3］ ROE L P．Upwind Differencing Schemes for Hyperbolic Conservation Laws with Source Terms［C］∥Springer Verlag．Proceedings of an Advanced Research Workshop：Nonlinear Hyperbolic Problems：Lecture Notes in Mathe matics，Volume 1270．St．Etienne，F(xiàn)rance January 13－17，1986：41－51．

［4］BERMUDEZA，VAZQUEZ CENDON EM．Upwind Meth ods for Hyperbolic Conservation Laws with Source Terms［J］．Computers and Fluids，1994，23（8）：1049－1071．

［5］ HUBBARD E M，GARCIA NAVARRO P．Flux Differ ence Splitting and the Balancing of Source Terms and Flux Gradients［J］．Journal of Computational Physics，2001，165（1）：89－125．

［6］ ZHOU G J，CAUSON M D，MINGHAN G C，et al．The Surface Gradient Method for the Treatment of Source Terms in the Shallow Water Equations［J］．Journal of Computational Physics，2001，168（1）：1－25．

［7］ GALLOU?T T，H?RARD JM，SEGUIN N．Some Ap proximate Godunov Schemes to Compute Shallow Water Equations with Topography［J］．Computers and Fluids，2003，32（4）：479－513．

［8］TORO E F．Shock Capturing Methods for Free Surface Shallow Flows［M］．Chichester：Wiley，2001．

［9］FERNANDEZ NIETO E D，NARBONA REINA G．Exten sion of WAF Type Methods to Non homogeneous Shallow Water Equations with Pollutant［J］．Journal of Scientific Computing，2008，36（2）：193－217．

（編輯：王 慰）

Application of a Godunov Type A lgorithm to the Pollutant Transport w ith Irregular Bedslope

YANG Jin bo1，DUANWen gang1，HUANG Guo bing1，LIDing fang2
（1．Hydraulics Department，Yangtze River Scientific Research Institute，Wuhan 430010，China；2．School of Mathematics and Statistics，Wuhan University，Wuhan 430072，China）

A well balanced weighted average flux（WAF）scheme is proposed for the nonhomogeneous shallow water equations with pollutant transport．Thewater surface elevation and water discharge are used as the conserved varia bles in the shallow water equations，while the conserved variables and the bed elevation are set on the staggered mesh．Here，theWAF scheme combined with the HLLC Riemann solver is directly performed for the nonhomoge neous shallow water equationswith pollutant transport，and the spatial second order central difference approximation is applied to source terms．Then，it’s proved that the extended WAF scheme is harmonious in theory．Finally，several numerical tests about the steady and nonsteady flows over the irregular bottom or with the wet／dry front are simulated．The results obviously suggest that the WAF scheme is a well balanced，robust，efficient，and practical scheme with high resolution．

shallow water equations；WAF scheme；dam break；pollutant transport

X52

A

1001－5485（2013）08－0097－05

10．3969／j．issn．1001－5485．2013．08．021

2013，30（08）：97－101

2013－04－16；

2013－06－19

國家自然科學(xué)基金資助項目（11102027）；中央級基本科研業(yè)務(wù)費(fèi)項目（CKSF2011014／SL，CKSF2013025／SL，CKSF2012010／SL）

楊金波（1982－），男，山東德州人，工程師，博士，主要從事計算流體力學(xué)，（電話）027－82823698（電子信箱）yangjb1225＠163．com。