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

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

0 引言

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

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

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

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

1 控制方程

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

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

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

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

因此，式（2）中源項可改寫為［22］：

該源項處理方式是將側壓力項I2和床面比降S0轉化為對靜力矩I1的偏導，精度高且便于計算。

2 數(shù)值計算方法

2.1 離散方法

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

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

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

2.2 HLLC近似Riemann求解器

HLLC近似Riemann求解器波的結構見圖2。

圖2 HLLC近似Riemann求解器波結構Fig.2 Wave structure of HLLC approximate Riemann solver

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

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

在計算通量時，波速估計方法的選擇至關重要，因為它們會影響數(shù)值結果的質量［23］。文獻［9］建議采用如下方法估算波速：

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

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

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

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

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

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

將式（2）中第一個分量（即質量守恒方程）帶入式（12）中前兩項，并將式（11）帶入，可得：

上式是關于Q*和A*的方程組，經求解可得：

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

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

2.3 變量空間重構

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

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

為了保證重構后的變量均為非負值，同時保證單元左右兩側重構值的平均值與單元平均值相等［26］，表達式如下：

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

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

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

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

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

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

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

2.4 源項及邊界條件

采用二階龍格—庫塔離散的時間分裂方法來處理源項［27］：

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

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

3 算 例

3.1 天然河道洪水演進算例

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

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

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

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

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

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

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

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

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

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

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

圖11 流量計算結果Fig.11 Flow calculation results

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

本算例HLLC法計算結果明顯優(yōu)于HLL法，主要原因為：①采用的控制方程不同，即守恒型圣維南方程優(yōu)于非守恒型方程；②變量空間重構方法不同，即本文變量空間重構方法優(yōu)于所有變量均直接采用MUSCL的變量空間重構方法。

3.2 喉口混合流算例

由于缺乏天然河道跨臨界流的實測數(shù)據(jù)，本文采用喉口混合流這一經典算例檢驗模型模擬跨臨界流以及淺水間斷等復雜明渠水流運動的能力。明渠長度為3 m，渠寬、底高程均沿程變化，表達式如下：

式中：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。計算空間步長Δx為0.01 m，不計河道摩阻項。

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

圖14 混合流計算結果Fig.14 Calculation results of mixed flow

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

4 結論

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

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

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

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

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