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

楊金波，段文剛，黃國(guó)兵，李訂芳

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

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

楊金波1，段文剛1，黃國(guó)兵1，李訂芳2

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

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

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

1 研究背景

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

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

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

2 數(shù)值格式

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

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

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

2．2 WAF有限體積格式

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

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

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

其中：

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

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

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

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

式中

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

左干河床

右干河床

此時(shí)，式（8）中的S3會(huì)自動(dòng)和估計(jì)的干河床式（9）或式（10）重合。

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

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

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

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

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

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

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

3 算 例

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

3．1 梯形河床上的污染物傳輸問(wèn)題

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

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

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

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

3．2 潰壩波中的污染物傳輸問(wèn)題

本算例研究了潰壩過(guò)程中污染物傳輸規(guī)律，計(jì)算區(qū)域?yàn)? 000 m長(zhǎng)的平底無(wú)摩擦河道，潰壩初始條件如下：

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

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

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

3．3 帶干濕河床的污染物傳輸問(wèn)題

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

其中地形定義為

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

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

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

4 結(jié) 語(yǔ)

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

［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

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

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