基于HLLC格式的異重流高分辨率立面二維數學模型

盧新華 秦超 張小峰

摘要：異重流是自然界中的常見現象，異重流運動過程中交界面附近可能存在物理量的間斷。為較好地捕捉這種間斷，建立了異重流高分辨率立面二維數學模型，該模型基于同位網格的Godunov型有限體積法求解σ坐標下的雷諾時均Navier-Stokes方程組。模型中水平方向界面數值通量采用HLLC近似黎曼求解器計算，湍流封閉采用非線性K-ε模型。選用3個經典的開閘式平坡和反坡異重流試驗對模型性能進行了檢驗。結果表明：該模型能較好地模擬異重流在平整或非平整床面上的運動過程，并具有較高的模擬精度。

關 鍵 詞：

異重流; HLLC; 反坡異重流; σ坐標; 非線性K-ε模型

中圖法分類號： TV131.2

文獻標志碼： A

DOI：10.16232/j.cnki.1001-4179.2021.06.021

異重流一般由兩種或多種流體因較小的密度差而引起。自然界中常出現因鹽度、溫度、泥沙濃度差異而引起的鹽水異重流、溫差異重流和泥沙異重流等。在異重流運動的數值模擬研究中常涉及物理量的間斷問題，較典型的如開閘式異重流，其在流體交界面附近存在較大梯度，數學模型若不在間斷處進行特殊處理較易產生數值震蕩[1-2]。目前，常用的異重流數學模型主要包括單層或深度平均數學模型，以及三維或立面二維數學模型[3-6]。對于單層或深度平均異重流數學模型而言，不同學者處理間斷方式不同，如Klemp等[7]及Ungarish等[8]在計算中采用顯式追蹤間斷位置，另有研究則采用激波捕捉數值格式以自動捕捉數值間斷[9-12];這類模型在間斷處通常須引入一系列假定（如假定間斷處局部弗氏數為某一經驗取值[2，7]），且模型中尚須引入卷吸系數、層間阻力系數以反映流體層間的質量交換和動量交換[9-10]，同時，模型無法描述流體界面附近的物理摻混過程。相比單層深度或平均異重流數學模型，三維模型可自動反映流體層間的質量和動量交換，能描述交界面附近異重流的摻混過程，因而在理論上比前者能更好地描述異重流的運動，但三維數學模型依然要處理界面間斷問題。近30 a來，近似黎曼求解器被廣泛應用于計算界面間斷問題，如潰壩洪水波、近岸海域波浪破碎、空氣動力學中激波問題的模擬研究[13-16]，并取得了較大成功。然而在異重流三維（立面）二維數值模擬方面鮮有研究報道。本研究嘗試基于文獻[14]提出的HLLC近似黎曼求解器建立模擬異重流運動的立面二維數學模型，以期檢驗近似黎曼求解器在異重流運動模擬方面的性能。模型中垂向采用σ坐標變換。

1 數學模型

采用Boussinesq近似、忽略非靜壓和垂向加速度影響，并引入垂向σ坐標變換（以適應不規則床面地形和捕捉水面變化）

t′=t，x′=x，σ=z′=z-zbh（1）

式中：t為時間;x，z為笛卡爾直角坐標;zb，h分別為河床高程及水深。可得σ坐標下雷諾時間平均的立面二維異重流運動控制方程組：

ht′+qxx′+ωσ=0（2）

qxt′+uqx+12gh2x′+uωσ=-ghzbx′

-ghx1ρ∫ηzρ-ρ0dz+S′ν，u（3）

qct′+uqcx′+cωσ=S′ν，c（4）

式中：ω=hdσdt=hσt+uσx+wσz;u，w分別為x，z 2個方向的流速分量;c為溶質（體積）濃度;qx=hu，qc=hc;g為重力加速度;ρ=1-cρ0+cρc，這里ρ0，ρc，ρ分別為水、溶質及混合密度;S′ν，φ為擴散項，并用下式計算：

S′ν，φ=x′ν+νtoφhφx′+σν+νtoφ1h2qφσ（5）

式中：φ=u，c;ν與νt分別為分子與湍流運動黏性系數，本文研究中νt采用經浮力修正的非線性K-ε兩方程湍流模型計算：

qKt′+uqKx′+Kωσ=S′ν，K+hGs+Gb-ε（6）

qεt′+uqεx′+εωσ=S′ν，ε+

hεKC1εGs+C3εGb-C2εε（7）

式（5）～（7）中：νt=CμK2ε;qK=hK，qε=hε;S′ν，K和S′ν，ε采用式（5）計算（式（5）中φ=K，ε）;Gb=1hgρνtocρσ為浮力產生項;Gs=-u′iu′juixj為剪切力產生項，在非線性K-ε模型中，

u′iu′j=23Kδij-CdK2εuixj+ujxi

-C1×K3ε2uixlulxj+ujxlulxi-23ulxkukxlδij

-C2×K3ε2uixkujxk-13ulxkulxk

-C3×K3ε2ukxiukxj-13ulxkulxkδij（8）

式中經驗參數取值為：Cμ=0.09，C1ε=1.44，C2ε=1.92，ou=oν=oc=oK=1.0，oε=1.3，Cd=2317.4+Smax，C1=1185.2+D2max，C2=158.5+D2max，C3=1370.4+D2max，Smax=Kmaxuixi， Dmax=Kmaxuixj（Smax與Dmax計算中不采用Einstein求和約定）;式（8）中i，j，k，l=1，2。參考文獻[17，18]，C3ε取值為0。

2 數值離散

本文基于同位網格的Godunov型有限體積法進行數值求解。為方便數值離散，將式（1）～（7）寫成如下守恒型矢量形式：

Ut′+Fx′=-Hσ+S0+Sν，H+Sν，σ+Sρ+Se（9）

其中，

U=hqxqcqKqε，F=qxuqx+g2h2uqcuqKuqε，H=ωuωcωKωεω（10）

S0=0-ghzbx′000

Sν，H≈0x′ν+νtouhux′x′ν+νtochcx′00

Sν，σ≈0σν+νtou1h2quσσν+νtoc1h2qcσ00（11）

Sρ=0-ghx1ρ∫ηzρ-ρ0dz000

Se≈000hGs+Gb-εhεKC1εGs+C3εGb-C2εε（12）

式（9）中Sν，H，Sν，σ為σ坐標下水平與垂向黏性擴散項矢量。

對式（9）采用兩步二階Runge-Kutta方法進行時間步進離散：

UIi，k=UNi，k+UI-1i，k2+αIΔt2Ri，k（13）

式中：i，k為網格單元編號;N為時間層;Δt為時間步長;I為Runge-Kutta步數。α（1）=2，α（2）=1。

Ri，k=FI-1i-1/2，k-FI-1i+1/2，kΔx′i+HI-1i，k-1/2-HI-1i，k+1/2Δσk+

Sl0+Slν，H+Slν，σ+Slρ+Slei，k（14）

Δx′=x′i+1/2-x′i-1/2，Δσk=Δσk+1/2-Δσk-1/2;當I取1和2時，UI-1i，k分別取UNi，k和U1i，k;式（14）中當l取（I-1）或I時，表示該項顯式或隱式離散，本文模型除Sν，σ隱式離散外Ri，k中其他部分均顯式離散。

式（14）中水平方向界面數值通量Fi±1/2，k采用HLLC近似黎曼求解器計算如下[13-14]：

Fi-1/2，k=FLi-1/2，k ?0≤ξLi-1/2，kF*Li-1/2，k ξLi-1/2，k≤0≤ξMi-1/2，kF*Ri-1/2，k ξMi-1/2，k≤0≤ξRi-1/2，kFRi-1/2，k ξRi-1/2，k≤0（15）

其中，

FL=FUL （16-1）

FR=FUR（16-2）

F*L=F*1F*2wLF*1KLF*1εLF*1T（16-3）

F*R=F*1F*2wRF*1KRF*1εRF*1T（16-4）

F*=F*1F*2F*3F*4F*5T=

ξRFL-ξLFR+ξLξRUR-ULξR-ξL（16-5）

式（15）中，波速ξL，ξR，ξM采用下式計算：

ξL=uR-2 ghRhL=0minuL- ghL，u*- gh*hL>0ξR=uL+2 ghLhR=0maxuR+ ghR，u*+ gh*hR>0ξM=ξLhRuR-ξR-ξRhLuL-ξLhRuR-ξR-hLuL-ξL（17）

u*=12uL+uR+ ghL- ghRh*=1g12 ghL+ ghR+14uL-uR2（18）

式（14）中垂向方向界面數值通量Hi，k±1/2采用加權二階中心-一階迎風格式計算，如

Hφi，k-1/2=wi，k-1/2θφC+1-θφUPi，k-1/2（19）

這里φ可取任意變量，如u，w，K，;θ為加權系數，取值越小模型越穩定、越大精度越高，本文取1.0以提高模擬精度;φCi，k-1/2=Δσk-1Δσk-1+Δσkφi，k+ΔσkΔσk-1+Δσkφi，k-1，φUPi，k-1/2=12φi，k-1+φi，k+12signwi，k-1/2φi，k-1-φi，k。

計算中床面阻力項采用全隱式離散，式（12）Sρ中水平方向偏導數在笛卡爾直角坐標系下進行計算，模型中通量項與源項的離散能嚴格保證靜水平衡。模型中的水流模塊計算原理及性能檢驗詳見文獻[13]，本文異重流模型基于該水流模型擴展得到。

3 模型檢驗

本文采用3個經典的室內異重流試驗以檢驗模型性能，其中包括平坡和反坡異重流試驗。考慮到所選取的試驗中異重流均沿水槽底部流動，為提高模型模擬精度，模擬中對床面附近網格進行加密處理，垂向網格結點的σ坐標值采用下式計算：

σk=β+1λk-β-1λkβ+1λk-1+β-1λk-1（20）

式中：λk=k-1∕nz，k=1，2，…，nz+1，為垂向網格結點編號，nz為垂向網格單元數;β為控制垂向網格疏密的因子（β>1）。當β越趨近于1時底部網格越密集，網格分布越不均勻;當β越大時網格分布越均勻。本文計算中β取值1.3。

3.1 Adduce等平坡異重流試驗

Adduce等[9]在羅馬大學的水力學實驗室進行了一系列開閘式平坡異重流試驗，試驗水槽如圖1所示，水槽兩端封閉，長度為L、寬度為b、高度為H0。在距離水槽左端x0處有一閘門，閘門左側為鹽水、密度為ρ1，右側為淡水、密度為ρ0，閘門兩側初始水深相同，均為h0。試驗時，迅速開啟的閘門，因流體密度差異，在重力作用下密度較大的流體會沿水槽底部運動而形成異重流。選取其中兩組具有代表性的試驗工況進行模擬，這兩組工況分別命名為“工況1”和“工況2”。計算中，L=3 m，b=0.2 m，其余參數如表1所列。表1中，g′=ρ1-ρ0ρ1g為因密度差異造成的初始有效重力加速度，這里g=9.81 m/s2。

基于線性-對數律計算床面切應力。為保證計算精度并提高計算效率進行了網格無關解分析。研究中選取了nx×nz=100×5，200×10，300×20及300×30四套網格進行模擬，其中nx、nz分別為水平和垂向方向的網格單元數。圖2（a）、2（b）分別為2組工況下計算的異重流頭部位置隨時間變化與實測值的比較。從圖2可以看出，nx×nz=300×20與nx×nz=300×30兩套網格計算結果接近，考慮到該算例計算量不大、本算例取nx×nz=300×30。圖2表明，本文模擬結果與試驗結果吻合較好。

圖3為計算的典型工況（以工況1為例）不同時刻的濃度分布及流場。在異重流開始坍塌階段（見圖3（b）），左側鹽水下潛入侵右側淡水形成一個大尺度的渦漩，同時在水面上形成一個微小振幅的波浪向右傳播（圖3（b）水槽中部流場為該波傳播到此處引起）。此后，異重流沿水槽底部繼續向前運動（見圖3（c）），并逐步進入自相似階段（見圖3（d））。由圖3可以看出，模型能較好地模擬異重流在平整床面上的運動過程。

3.2 Hatcher等平坡異重流試驗

為進一步檢驗模型精度，本文選取Hatcher等[2]異重流試驗數據對該模型進行驗證。在長L=9.14 m、寬b=0.127 m的水槽中共進行了兩個組次的試驗，每個組次的試驗均重復進行了3次以保證試驗的可重復性，試驗原理同2.1節，試驗具體參數如表2所列。試驗中采用MicroADV測量了x*=6.68所在斷面的水槽中部垂線上距床面3個不同位置高度的異重流內部流速。監測點處的無量綱高度分別為h*A=hA/h0=0.188，0.125和0.063，hA為監測點到水槽底部的垂直距離。

類似于3.1節算例，此節算例經網格無關解分析后選取nx×nz=1 800×30。圖4（a）、4（b）分別計算了工況1和工況2下無量綱化異重流運動距離隨時間變化與實測值的比較，其中“實測值1”“實測值2”“實測值3”為每組工況下重復的3次試驗數據。圖4中無量綱異重流頭部位置定義為x*f=x/x0，其中x為異重流頭部到閘門處的距離;無量綱時間定義為t*=t/x0 g′h0。由圖4可知，模擬的異重流頭部位置隨時間變化的結果與試驗數據吻合良好。

圖5給出了2種工況下不同位置高度處監測點水平流速隨時間變化的計算值與實測值的比較。圖5同時給出了Hatcher等[2]采用雙層深度平均模型的計算結果。由于該模型求解上層水體與下層異重流流速的垂線平均值，因此，圖5中同一工況下各個垂線位置處該模型的流速計算值相同。由圖5可看出，2種工況下，在異重流頭部到達監測點之前（t*≤15），各監測點的流速值相對較小，而當異重流頭部到達監測點時（t*≈15）流速值迅速增大，之后各監測點流速隨時間有逐漸減小的趨勢，本文模型和Hatcher等人的雙層深度平均數學模型均能較好地反映異重流內流速的這一變化過程。由圖5亦可知，相比雙層深度平均數學模型，本文模型能模擬異重流流速在垂線上的分布，且與實測值較為吻合。

3.3 Marleau等反坡異重流試驗

為檢驗模型在非平整床面上異重流的模擬性能，本節選取Marleau等[21]開閘式反坡異重流試驗作為驗證算例，試驗水槽如圖6所示，水槽長L1=1.975 m，寬b1=0.176 m，高H1=0.485 m。閘門與水槽左端的距離為Ll，水槽右端插入一塊剛性塑料板作為底坡、其長度為Ls。試驗時，閘門左側水槽底部填充高度為D、密度為ρ1的鹽水，周圍淡水的密度為ρ0。本次模擬的試驗參數為：D=H=0.3 m，ρ1=1 001 kg∕m3，ρ0=998.5 kg/m3，Ll=0.284 m，Ls=1.2 m，坡度s=0.25。經網格無關解分析后，計算中選取nx×nz=600×30的網格。

不同時刻的濃度分布及流場如圖7所示，圖中t從異重流到達斜坡底端時起算。當異重流到達斜坡底端時（見圖7（a）），異重流厚度約為初始水深的一半，異重流與周圍水體交界面處產生一大尺度的渦漩;之后，異重流在慣性作用下沿斜坡向上運動（見圖7（b）），當異重流爬坡到某一高度時，異重流與清水交界面附近出現兩個大尺度的渦漩（見圖7（c）），該階段在重力與床面阻力影響下異重流運動速度降低，頭部厚度也逐漸減小。

圖8為異重流運動距離隨時間變化與實測值的比較，其中x′為異重流沿斜坡向上運動的距離。由圖8可知，異重流沿斜坡向上做減速運動，總體上數值模擬結果與試驗數據吻合良好。

4 結 論

本文基于HLLC近似黎曼求解器建立了異重流高分辨率立面二維數學模型，該模型采用σ坐標變換以適應不規則床面地形和捕捉水面變化。文中選用3個涉及到物理量間斷的開閘式平坡及反坡異重流試驗對模型性能進行檢驗，檢驗結果表明模型能較好地捕捉間斷問題，并在異重流模擬方面具有較高精度。

需要說明的是，本文建立的模型為立面二維數學模型，但較易拓展成三維數學模型以模擬復雜邊界條件下異重流的三維運動過程。另外，本文模型中湍流模型采用的是非線性K-ε兩方程模型、水平方向數值通量采用HLLC近似黎曼求解器，今后可對比研究其它湍流模型及近似黎曼求解器在異重流模擬方面的模擬性能。

參考文獻：

[1] ROTTMAN J W，SIMPSON J E.Gravity currents produced by instantaneous releases of a heavy fluid in a rectangular channel [J].Journal of Fluid Mechanics，1983（135）：95-110.

[2] HATCHER T M，VASCONCELOS J G.Finite-volume and shock-capturing shallow water equation model to simulate Boussinesq-type lock-exchange flows [J].Journal of Hydraulic Engineering，ASCE，2013，139（12）：1223-1233.

[3] 方春明，韓其為，何明民.異重流潛入條件分析及立面二維數值模擬 [J].泥沙研究，1997（4）：70-77.

[4] 彭楊，李義天，槐文信.異重流潛入運動的剖面二維數值模擬 [J].泥沙研究，2000（6）：25-30.

[5] 張芝永，楊元平，程文龍.異重流三維非靜壓數值模擬研究 [J].中國水運，2018，18（11）：74-76.

[6] 陸俊卿，張小峰，崔占峰.各向異性密度流模型及其驗證 [J].水科學進展，2010，21（1）：95-100.

[7] KLEMP J B，ROTUNNO R，SKAMAROCK W C.On the dynamics of gravity currents in a channel [J].Journal of Fluid Mechanics，1994（269）：169-198.

[8] UNGARISH M，ZEMACH T.On the slumping of high Reynolds number gravity currents in two-dimensional and axisymmetric configurations [J].European Journal of Mechanics-B/Fluids，2005，24（1）：71-90.

[9] ADDUCE C，SCIORTINO G，PROIETTI S.Gravity currents produced by lock exchanges：experiments and simulations with a two-layer shallow-water model with entrainment [J].Journal of Hydraulic Engineering，ASCE，2012，138（2）：111-121.

[10] HU P，CAO Z X，PENDER G，et al.Numerical modelling of turbidity currents in the Xiaolangdi reservoir，Yellow River，China [J].Journal of Hydrology，2012（464-465）：41-53.

[11] LU X H，DONG B J，MAO B，et al.A robust and well-balanced numerical model for solving the two-layer shallow water equations over uneven topography [J].Comptes Rendus Mécanique，2015，343（7-8）：429-442.

[12] LIN J，MAO B，LU X H.A two-layer hydrostatic-reconstruction method for high-resolution solving of the two-layer shallow-water equations over uneven bed topography [J].Mathematical Problems in Engineering，2019（2019）：1-14.

[13] LU X H，MAO B，ZHANG X F，et al.Well-balanced and shock-capturing solving of 3D shallow-water equations involving rapid wetting and drying with a local 2D transition approach [J].Computer Methods in Applied Mechanics and Engineering，2020（364）：112897.

[14] TORO E F.Riemann solvers and numerical methods for fluid dynamics [M].Newjersey：Springer，2009.

[15] CAO Z X，PENDER G，WALLIS S，et al.Computational dam-break hydraulics over erodible sediment bed [J].Journal of Hydraulic Engineering，ASCE，2004，130（7）：689-703.

[16] LU X H，XIE S B.Depth-averaged non-hydrostatic numerical modeling of nearshore wave propagations based on the FORCE scheme[J].Coastal Engineering，2016（114）：208-219.

[17] MA G F，KIRBY J T，SHI F Y.Numerical simulation of tsunami waves generated by deformable submarine landslides [J].Ocean Modelling，2013（69）：146-165.

[18] 盧新華.基于近似黎曼求解器的三維淺水方程組求解方法 [J].人民長江，2018，49（20）：74-80.

[19] CHOWDHURY M R，TESTIK F Y.Laboratory testing of mathematical models for high-concentration fluid mud turbidity currents [J].Ocean Engineering，2011，38（1）：256-270.

[20] HUPPERT H E，SIMPSON J E.The slumping of gravity currents [J].Journal of Fluid Mechanics，1980，99（4）：785-799.

[21] MARLEAU L J，FLYNN M R，SUTHERLAND B R.Gravity currents propagating up a slope [J].Physics of Fluids，2014，26（4）：213-234.

（編輯：李 慧）

A vertical 2D high-resolution numerical model for gravity current based on HLLC scheme

LU Xinhua，QIN Chao，ZHANG Xiaofeng

（State Key Laboratory of Water Resources and Hydropower Engineering Science，Wuhan University，Wuhan 430072，China）

Abstract：

Gravity currents are common phenomenon in nature.Discontinuities may exist near the interface during the movement of gravity currents.In order to capture this discontinuity well，we develop a vertical 2D high-resolution numerical model for gravity currents.The developed model uses the Godunov-type finite-volume method based on the isometric grid to solve the Reynolds time-mean Navier Stokes equations in σ coordinates.The horizontal inter-cell numerical flux is evaluated by the HLLC approximate Riemann solver，and the MUSCL scheme is employed for horizontal interface value reconstructions.A nonlinear k-ε model is employed for turbulence closure.Three classic lock-exchange experiments of gravity currents propagating on flat and adverseslope beds are employed to verify the performance of the model.Results show that the developed model simulates the movement of gravity currents well on flat or uneven bed with high accuracy.

Key words：

gravity current;HLLC;gravity current on adverse slope;σ coordinates;nonlinear k-ε model