基于HLLC近似Riemann求解器的天然河道水流運(yùn)動模擬

孫萬光，楊海滔，楊斌斌，房 巍，范寶山

（1.中水東北勘測設(shè)計(jì)研究有限責(zé)任公司，長春130061；2.水利部寒區(qū)工程技術(shù)研究中心，長春130061；3.遼寧省河庫管理服務(wù)中心（遼寧省水文局），沈陽110003）

0 引言

天然河道一維非恒定流數(shù)值模擬在水利規(guī)劃設(shè)計(jì)以及防洪減災(zāi)等領(lǐng)域應(yīng)用十分廣泛。天然河道特別是山區(qū)性河道斷面幾何形狀快速變化、河底比降變化大，水流流態(tài)可能出現(xiàn)急流、緩流交替現(xiàn)象；另外，高緯度嚴(yán)寒地區(qū)江河冰塞及冰壩現(xiàn)象頻繁發(fā)生，冰壩潰決后壩下將出現(xiàn)激波間斷，這些都屬于淺水間斷流動。因間斷處水力要素突變，具有較強(qiáng)非線性，傳統(tǒng)數(shù)值計(jì)算方法（如有限差分）往往失效［1］。

對于間斷問題，控制方程上，Lax［2］證明了守恒數(shù)值格式的解可收斂至守恒控制方程的弱解。Hou［3］證明了采用非守恒型控制方程求解間斷問題將會得到錯(cuò)誤結(jié)果，這些方程僅對光滑流動有效，對于間斷問題，根據(jù)跳躍條件計(jì)算出的激波波速是錯(cuò)誤的。數(shù)值計(jì)算方法上，Godunov［4］格式因具備較強(qiáng)的不連續(xù)問題處理能力而廣受歡迎［5-8］。在Godunov 格式下，將齊次淺水方程組的Riemann 問題定義為具有分段變量恒定的初值問題，采用精確Riemann 求解器［9］可求得單元界面處的通量值，精度高，但需要對非線性代數(shù)恒等式進(jìn)行迭代求解，增加了計(jì)算資源的消耗［10］。為了簡化計(jì)算，多種近似Riemann 求解器得到快速發(fā)展，如Roe 法［11］和HLL 法［12］等。Roe 線性化Riemann 求解器缺點(diǎn)如下［10］：①在跨臨界流或激波模擬中，為避免非物理解必須進(jìn)行熵修正；②在強(qiáng)稀疏波作用下（水流非常淺的區(qū)域），線性化Riemann 求解器計(jì)算水深出現(xiàn)負(fù)值；③在強(qiáng)波相互作用下，線性化Riemann 解一般缺乏魯棒性。HLL 近似Riemann 求解器將波族簡化為2 個(gè)，只適用于具有2 個(gè)方程的一維系統(tǒng)［10］，對接觸間斷和剪切波模擬精度較差。基于2 波族的HLL 近似Riemann 求解器，Einfeldt 提出了波速的計(jì)算方法，這被稱作HLLE［13］求解器，之后又進(jìn)一步將HLL 求解器單個(gè)中間狀態(tài)修改成線性分布狀態(tài)，稱作HLLM［14］求解器，特定參數(shù)條件下，HLLM 求解器降為線性化Roe 求解器。為了恢復(fù)HLL 求解器缺失的波族，Toro［15］提出了HLLC求解器，為3波族模型，可準(zhǔn)確求解接觸間斷和剪切波問題［16］。Alireza［17］采用HLLC 求解器對潰壩水流運(yùn)動進(jìn)行模擬；Zhou［18］建立了基于HLLC 求解器的二維水動力和水質(zhì)耦合模型。對比HLL 求解器，HLLC 求解器可擴(kuò)充濃度組分方程，為天然河道環(huán)境水力學(xué)數(shù)值模擬創(chuàng)造了有利條件，推廣應(yīng)用前景十分廣闊。以往的研究中，多采用HLLC 求解器對淺水方程進(jìn)行求解，而天然河道一維水動力通常以圣維南方程作為控制方程，采用HLLC 求解器對圣維南方程進(jìn)行求解至今少見報(bào)道。與淺水方程相比，守恒型圣維南方程需考慮河寬及其沿程變化產(chǎn)生的影響，實(shí)際應(yīng)用中，需重新推導(dǎo)基于守恒型圣維南方程的HLLC求解器中通量計(jì)算公式。

Godunov 格式下，為了獲得高精度數(shù)值解，一般需要進(jìn)行變量空間重構(gòu)。以MUSCL［19］（Monotonic Upstream-centered Scheme for Conservation Laws）為代表的變量空間重構(gòu)方法應(yīng)用比較廣泛，但該方法僅適用于斷面漸變條件下的變量空間重構(gòu)。天然河道斷面幾何形狀復(fù)雜，河寬和水深沿程快速變化，直接采用MUSCL法進(jìn)行變量空間重構(gòu)會產(chǎn)生較大誤差。

本文采用守恒型圣維南方程作為天然河道一維非恒定流控制方程，基于Godunov格式，提出了天然河道斷面幾何形狀快速變化條件下的變量空間重構(gòu)方法，推導(dǎo)了基于守恒型圣維南方程的HLLC 求解器通量計(jì)算公式，為天然河道復(fù)雜水動力數(shù)值模擬提供一種高精度、簡便的方法。

1 控制方程

守恒型圣維南方程組表達(dá)式如下［20］：

式中：U為變量；F為通量；S為源項(xiàng)；t為時(shí)間；x為空間坐標(biāo)；A=A（x，t）為過流斷面面積；Q為流量；g為重力加速度；S0為床面比降；Sf為摩阻比降；I1和I2分別為靜力矩和側(cè)壓力，表達(dá)式如下：

式中：h為水深；b(x，η)為斷面寬度。

天然河道中，床面比降S0變化較大，在源項(xiàng)處理中通常避開底坡源項(xiàng)。Cunge［21］給出了對I1求導(dǎo)的表達(dá)式：

因此，式（2）中源項(xiàng)可改寫為［22］：

該源項(xiàng)處理方式是將側(cè)壓力項(xiàng)I2和床面比降S0轉(zhuǎn)化為對靜力矩I1的偏導(dǎo)，精度高且便于計(jì)算。

2 數(shù)值計(jì)算方法

2.1 離散方法

采用Godunov 格式的有限體積法對控制方程離散（見圖1），表達(dá)式如下：

圖1 中心格式有限體積法單元離散Fig.1 Finite volume element discretization with central scheme

式中：i為計(jì)算單元序號；n為時(shí)刻；Δxi為計(jì)算單元i的長度；Δt為計(jì)算時(shí)間步長，為n時(shí)刻最大波速；Fi+1/2和Fi-1/2分別為i+1/2和i-1/2界面處的通量。

2.2 HLLC近似Riemann求解器

HLLC近似Riemann求解器波的結(jié)構(gòu)見圖2。

圖2 HLLC近似Riemann求解器波結(jié)構(gòu)Fig.2 Wave structure of HLLC approximate Riemann solver

HLLC 求解器是在HLL 求解器2 波族的基礎(chǔ)上增加了中間波，中間波的波速為S*。HLLC通量表達(dá)式如下［9，16］：

式中：UL和UR分別為界面左側(cè)和右側(cè)變量；FL和FR分別為界面左側(cè)和右側(cè)通量；U*L和U*R分別為中間波左側(cè)和右側(cè)變量，為待求變量；F*L和F*R分別為中間波左側(cè)和右側(cè)通量；SL和SR分別為界面左側(cè)和右側(cè)波速。

在計(jì)算通量時(shí)，波速估計(jì)方法的選擇至關(guān)重要，因?yàn)樗鼈儠绊憯?shù)值結(jié)果的質(zhì)量［23］。文獻(xiàn)［9］建議采用如下方法估算波速：

式中：uL和uR分別為界面左側(cè)和右側(cè)流速；aL和aR分別為界面左側(cè)和右側(cè)重力波波速；qK（K=L，R）的表達(dá)式如下：

式中：aTR是基于界面左、右兩側(cè)均為稀疏波的假定對重力波波速a的估計(jì)值。

為了計(jì)算界面通量，還需U*L和U*R的值。這里引入如下假設(shè)［9，16］：

事實(shí)上，上述假設(shè)對于精確Riemann求解器也是正確的［16］。

本文推導(dǎo)得出基于守恒型圣維南方程的中間波變量U*=[A*，Q*]T以及中間波波速S*的表達(dá)式。

對跨越波速SL、SR和S*的狀態(tài)分別運(yùn)用Rankine-Hugoniot條件［24］，可得：

將式（2）中第一個(gè)分量（即質(zhì)量守恒方程）帶入式（12）中前兩項(xiàng)，并將式（11）帶入，可得：

上式是關(guān)于Q*和A*的方程組，經(jīng)求解可得：

根據(jù)式（11）和式（12）的第三項(xiàng)可知，當(dāng)HLLC 法應(yīng)用于一維水動力方程數(shù)值求解時(shí)，U*L=U*R，并且F*L=F*R，通量計(jì)算表達(dá)式還可以寫成：

上式與HLL 通量計(jì)算表達(dá)式完全一致。由此可見，HLLC法引入的中間波僅作用于擴(kuò)充的濃度組分方程，一維水動力計(jì)算中，HLLC法和HLL法計(jì)算結(jié)果完全一致。

2.3 變量空間重構(gòu)

在變量空間重構(gòu)中，MUSCL方法將單元平均值線性外推到單元界面，表達(dá)式如下：

為了避免虛假振蕩，外推的斜率由斜率限制器函數(shù)限制。坡度限制器采用minmod 限制器，限制器中有2 個(gè)變量，當(dāng)其符號相同時(shí)取其最小值，否則取0［25］，因此單元坡度表達(dá)式如下：

為了保證重構(gòu)后的變量均為非負(fù)值，同時(shí)保證單元左右兩側(cè)重構(gòu)值的平均值與單元平均值相等［26］，表達(dá)式如下：

天然河道斷面幾何形狀不規(guī)則，河寬和水深沿程快速變化，若直接按上述方法進(jìn)行變量空間重構(gòu)，可導(dǎo)致通量計(jì)算結(jié)果嚴(yán)重失真，無法保證格式守恒。本文首先依據(jù)過流斷面面積和靜力矩等效原則將河道斷面概化成矩形，通過線性插值構(gòu)造單元界面處斷面幾何形狀，采用基于MUSCL法的水位空間二階精度重構(gòu)結(jié)果計(jì)算單元界面兩側(cè)過流斷面面積和靜力矩的重構(gòu)值，保證計(jì)算格式守恒。具體步驟如下：

（1）天然河道斷面概化方法。給定河道水位，首先計(jì)算過流斷面面積A和靜力矩I1（式（3）積分獲得），依據(jù)過流斷面面積和靜力矩等效原則將天然河道斷面概化成矩形，即：I1，其中分別為等效過流斷面面積和等效靜力矩。由此可推導(dǎo)出過流斷面等效水深h′、等效河寬B′和等效河底高程，表達(dá)式如下：

水深是波速計(jì)算的關(guān)鍵參數(shù)，天然河道通常采用h=A/B來計(jì)算斷面平均水深，其中B為水面寬度。然而對于復(fù)式斷面（見圖3），B和h均存在突變（見圖4），對變量空間重構(gòu)產(chǎn)生明顯影響。而等效水深和等效河寬則不存在突變。此概化方法保證單元水體受力條件不變，物理概念清晰，等效河寬和等效水深的概化結(jié)果唯一。

圖3 復(fù)式斷面示意圖（單位：m）Fig.3 Schematic diagram of compound section

圖4 斷面平均水深和等效水深對比Fig.4 Comparison of section average water depth and equivalent water depth

（2）單元界面處斷面幾何形狀構(gòu)造。根據(jù)概化后相鄰河道斷面的寬度和底高程，通過線性插值構(gòu)造單元界面處斷面幾何形狀。以i+1/2界面處斷面幾何形狀構(gòu)造為例，表達(dá)式如下：

（3）變量空間重構(gòu)。需空間重構(gòu)的變量包括A、Q和I1。先采用式（16）構(gòu)造界面處（以i+1/2 界面為例）的水位和流量，然后基于構(gòu)造界面處過流斷面面積和靜力矩，表達(dá)式如下：

2.4 源項(xiàng)及邊界條件

采用二階龍格—庫塔離散的時(shí)間分裂方法來處理源項(xiàng)［27］：

源項(xiàng)包括斷面靜力矩梯度和摩阻項(xiàng)。摩阻項(xiàng)可采用顯格式處理，離散表達(dá)式如下：

邊界條件有兩種：①固壁邊界，引入反射邊界條件，邊界處的水位與鄰近單元一致，而流速與鄰近單元相反；②允許波穿過的邊界條件，邊界處的水位和流速與鄰近單元一致。

3 算 例

3.1 天然河道洪水演進(jìn)算例

某河流經(jīng)盆地和峽谷，河道平面形態(tài)變化較大，常水位下河寬變化范圍為200～1 000 m，河道平面示意圖見圖5，此河段長度59.44 km，布置實(shí)測大斷面18條（0～17號）。

圖5 河道平面示意圖及斷面布置Fig.5 Plan and section layout of river

河道典型橫斷面見圖6。盆地河段為復(fù)式斷面，河寬較寬，而峽谷段為“V”字形斷面，河寬較小。河床縱斷面起伏劇烈，深泓最低點(diǎn)高程-48.67 m，深泓最高點(diǎn)高程為33.3 m，縱斷面見圖7。該河段1994年發(fā)生了頻率P=4%洪水，洪峰流量為44 400 m3/s（見圖8）。洪水過后對各斷面洪痕進(jìn)行了準(zhǔn)確測量，該場次洪水水面線數(shù)據(jù)完備、準(zhǔn)確，可以此作為河道糙率率定以及計(jì)算結(jié)果驗(yàn)證的依據(jù)。根據(jù)洪痕實(shí)測數(shù)據(jù)，對洪峰流量條件下的河道糙率進(jìn)行率定，糙率取值范圍為0.028～0.088。

圖7 河道縱斷面圖Fig.7 River profile

圖8 1994年典型洪水過程Fig.8 Typical flood process in 1994

將斷面計(jì)算最高水位與實(shí)測洪痕進(jìn)行對比，結(jié)果見圖9（本文方法簡稱“HLLC”）。為了便于對比，基于HEC-RAS 軟件，采用相同的斷面數(shù)據(jù)和邊界條件，建立一維非恒定流數(shù)值計(jì)算模型，對1994年典型洪水過程進(jìn)行計(jì)算（簡稱“HEC-RAS”）；HEC-RAS［28］軟件非恒定流控制方程為非守恒型，模型采用有限差分進(jìn)行求解。文獻(xiàn)［25］中非恒定流控制方程同為非守恒型，基于Godunov 格式，采用HLL 近似Riemann 求解器求解，采用該文方法同步進(jìn)行對比（簡稱“HLL”）。

從圖9 中可見，基于HLLC 算法的計(jì)算結(jié)果與實(shí)測值吻合良好，而HEC-RAS 和HLL法模擬結(jié)果與實(shí)測值有較大偏差，偏差的起始點(diǎn)位于9 號斷面，9～8 號斷面河道水面比降明顯高于實(shí)測值，9 號斷面上游水位計(jì)算結(jié)果均高于實(shí)測值，其中HECRAS 法偏差最大可達(dá)1.69 m（位于11 號斷面），HLL 法偏差最大可達(dá)1.42 m（位于10號斷面）。9～8號斷面河道由盆地向峽谷過渡（見圖6），河寬快速縮窄，斷面最高水位對應(yīng)的水面寬度分別為595 和269 m。以1994年典型洪水中最大洪峰起漲至消落過程為分析對象（6月16日-6月25日），將斷面間動量方程對應(yīng)的通量進(jìn)行比較，守恒型方程的通量為Q2/A+gI1，非守恒型方程的通量為Q2/A，結(jié)果見圖10。

圖6 河道典型實(shí)測斷面及概化斷面圖Fig.6 Typical measured section and generalized section of river

圖9 斷面最高水位計(jì)算值與實(shí)測值對比Fig.9 Comparison of calculated and measured maximum water level

從圖10 中可見，守恒型動量方程對應(yīng)的通量，低水位情況下9 號斷面小于8 號斷面，而高水位情況下恰恰相反；非守恒型動量方程對應(yīng)的通量，9 號斷面在整個(gè)分析過程中均小于8 號斷面。由式（6）可知，流量Q的大小取決于單元入口和出口通量差，以及源項(xiàng)大小。因斷面間流量過程差異較小（見圖11），在非守恒型動量方程下，單元入口通量明顯小于出口通量，只能通過加大源項(xiàng)中的水面比降來提高河段流量，這也是HECRAS 和HLL 法計(jì)算出9 號～8 號斷面河道水面比降明顯偏大的根本原因，而本文HLLC 法采用守恒型方程則不會出現(xiàn)此偏差。由此可見，對于斷面幾何形狀快速變化的河段，采用守恒型控制方程進(jìn)行數(shù)值計(jì)算是十分必要的。

圖10 斷面間動量守恒方程對應(yīng)通量比較Fig.10 Comparison of fluxes corresponding to momentum conservation equations between cross sections

圖11 流量計(jì)算結(jié)果Fig.11 Flow calculation results

天然河道斷面幾何形狀快速變化，斷面間水深和過流斷面面積關(guān)系差別較大。以8號斷面和9號斷面為例（見圖6），最高水位條件下8 號斷面過流斷面面積為9 595 m2，等效水深為47.13 m；而9 號斷面過流斷面面積為17 133 m2，為8 號斷面的1.79 倍，等效水深為33.9 m，為8 號斷面的0.72 倍。若直接進(jìn)行變量空間重構(gòu)則會產(chǎn)生較大誤差，通量計(jì)算結(jié)果嚴(yán)重失真，同時(shí)也會產(chǎn)生格式不守恒問題。人為構(gòu)造單元界面處斷面幾何形狀，采用界面處基于MUSCL法的水位重構(gòu)值計(jì)算與之對應(yīng)的過流斷面面積和靜力矩，以此作為變量空間重構(gòu)結(jié)果，成功避開了斷面幾何形狀快速變化對變量空間重構(gòu)的影響，同時(shí)保證了計(jì)算格式的守恒。

本算例HLLC法計(jì)算結(jié)果明顯優(yōu)于HLL法，主要原因?yàn)椋孩俨捎玫目刂品匠滩煌?，即守恒型圣維南方程優(yōu)于非守恒型方程；②變量空間重構(gòu)方法不同，即本文變量空間重構(gòu)方法優(yōu)于所有變量均直接采用MUSCL的變量空間重構(gòu)方法。

3.2 喉口混合流算例

由于缺乏天然河道跨臨界流的實(shí)測數(shù)據(jù)，本文采用喉口混合流這一經(jīng)典算例檢驗(yàn)?zāi)Ｐ湍M跨臨界流以及淺水間斷等復(fù)雜明渠水流運(yùn)動的能力。明渠長度為3 m，渠寬、底高程均沿程變化，表達(dá)式如下：

式中：B（x）為x處河寬；Zb（x）為x處河道底高程。

圖12 河道底高程沿程變化Fig.12 Variation of river bottom elevation along the river

圖13 河寬沿程變化Fig.13 River width variation along the river

渠道初始水位1.0 m，初始流量為0，入口流量Q=1.878 m3/s，下游水位Z=1.0 m。計(jì)算空間步長Δx為0.01 m，不計(jì)河道摩阻項(xiàng)。

計(jì)算結(jié)果表明（見圖14），水位的計(jì)算值與解析解吻合較好：當(dāng)0≤x≤1 時(shí)，水位計(jì)算值為1.094 1 m，而解析解為1.094 4 m；當(dāng)x=1.92 m時(shí)，計(jì)算水位最低，其值為0.500 2 m，而解析解中當(dāng)x=1.95 m 時(shí)，水位最低，其值為0.499 7 m。流量方面，計(jì)算值與解析解完全一致。

圖14 混合流計(jì)算結(jié)果Fig.14 Calculation results of mixed flow

該算例中既有緩流又有急流，還存在跨臨界流動，即水跌和水躍，是典型的淺水間斷流動。本文算法對激波進(jìn)行精確捕捉，計(jì)算精度高，適應(yīng)性強(qiáng)。

4 結(jié)論

本文采用守恒型圣維南方程作為控制方程，基于Godunov格式，采用HLLC 近似Riemann 求解器對模型進(jìn)行求解，主要結(jié)論如下。

（1）推導(dǎo)得出了基于守恒型圣維南方程的HLLC 近似Riemann 求解器通量計(jì)算表達(dá)式，將該求解器推廣應(yīng)用至天然河道一維非恒定流計(jì)算中。

（2）提出了針對天然河道復(fù)雜斷面幾何形狀下的變量空間重構(gòu)方法：提出了基于過流斷面面積和靜力矩等效原則的天然河道斷面概化方法，通過線性插值構(gòu)造單元界面處斷面幾何形狀，基于水位的空間二階精度重構(gòu)結(jié)果計(jì)算界面兩側(cè)過流斷面面積和靜力矩的重構(gòu)值，避免變量空間重構(gòu)產(chǎn)生較大誤差，同時(shí)保證計(jì)算格式守恒。

（3）天然河道斷面幾何形狀快速變化條件下，本文方法計(jì)算結(jié)果與實(shí)測值吻合良好，明顯優(yōu)于HEC-RAS 法和HLL 法（非守恒型圣維南方程，MUSCL變量空間重構(gòu)方法），同時(shí)對于混合水流運(yùn)動以及淺水間斷具有較高的模擬精度。

（4）本文將HLLC 近似Riemann 求解器由淺水方程拓展至守恒型圣維南方程，較HLL 求解器增加了中間波，可擴(kuò)充濃度組分方程，為天然河道一維水動力及環(huán)境水力學(xué)高精度數(shù)值模擬提供了新的方法?！?/p>