盧新華 秦超 張小峰

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

關(guān) 鍵 詞：

異重流; HLLC; 反坡異重流; σ坐標(biāo); 非線性K-ε模型

中圖法分類號： TV131.2

文獻(xiàn)標(biāo)志碼： A

DOI：10.16232/j.cnki.1001-4179.2021.06.021

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

1 數(shù)學(xué)模型

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

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

式中：t為時間;x，z為笛卡爾直角坐標(biāo);zb，h分別為河床高程及水深?？傻忙易鴺?biāo)下雷諾時間平均的立面二維異重流運(yùn)動控制方程組：

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為溶質(zhì)（體積）濃度;qx=hu，qc=hc;g為重力加速度;ρ=1-cρ0+cρc，這里ρ0，ρc，ρ分別為水、溶質(zhì)及混合密度;S′ν，φ為擴(kuò)散項，并用下式計算：

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

式中：φ=u，c;ν與νt分別為分子與湍流運(yùn)動黏性系數(shù)，本文研究中νt采用經(jīng)浮力修正的非線性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ρσ為浮力產(chǎn)生項;Gs=-u′iu′juixj為剪切力產(chǎn)生項，在非線性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）

式中經(jīng)驗參數(shù)取值為：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。參考文獻(xiàn)[17，18]，C3ε取值為0。

2 數(shù)值離散

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

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

其中，

U=hqxqcqKqε，F(xiàn)=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ν，σ為σ坐標(biāo)下水平與垂向黏性擴(kuò)散項矢量。

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

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

式中：i，k為網(wǎng)格單元編號;N為時間層;Δt為時間步長;I為Runge-Kutta步數(shù)。α（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;當(dāng)I取1和2時，UI-1i，k分別取UNi，k和U1i，k;式（14）中當(dāng)l取（I-1）或I時，表示該項顯式或隱式離散，本文模型除Sν，σ隱式離散外Ri，k中其他部分均顯式離散。

式（14）中水平方向界面數(shù)值通量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）中垂向方向界面數(shù)值通量Hi，k±1/2采用加權(quán)二階中心-一階迎風(fēng)格式計算，如

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

這里φ可取任意變量，如u，w，K，;θ為加權(quán)系數(shù)，取值越小模型越穩(wěn)定、越大精度越高，本文取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ρ中水平方向偏導(dǎo)數(shù)在笛卡爾直角坐標(biāo)系下進(jìn)行計算，模型中通量項與源項的離散能嚴(yán)格保證靜水平衡。模型中的水流模塊計算原理及性能檢驗詳見文獻(xiàn)[13]，本文異重流模型基于該水流模型擴(kuò)展得到。

3 模型檢驗

本文采用3個經(jīng)典的室內(nèi)異重流試驗以檢驗?zāi)Ｐ托阅?，其中包括平坡和反坡異重流試驗。考慮到所選取的試驗中異重流均沿水槽底部流動，為提高模型模擬精度，模擬中對床面附近網(wǎng)格進(jìn)行加密處理，垂向網(wǎng)格結(jié)點(diǎn)的σ坐標(biāo)值采用下式計算：

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

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

3.1 Adduce等平坡異重流試驗

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

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

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

3.2 Hatcher等平坡異重流試驗

為進(jìn)一步檢驗?zāi)Ｐ途龋疚倪x取Hatcher等[2]異重流試驗數(shù)據(jù)對該模型進(jìn)行驗證。在長L=9.14 m、寬b=0.127 m的水槽中共進(jìn)行了兩個組次的試驗，每個組次的試驗均重復(fù)進(jìn)行了3次以保證試驗的可重復(fù)性，試驗原理同2.1節(jié)，試驗具體參數(shù)如表2所列。試驗中采用MicroADV測量了x*=6.68所在斷面的水槽中部垂線上距床面3個不同位置高度的異重流內(nèi)部流速。監(jiān)測點(diǎn)處的無量綱高度分別為h*A=hA/h0=0.188，0.125和0.063，hA為監(jiān)測點(diǎn)到水槽底部的垂直距離。

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

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

3.3 Marleau等反坡異重流試驗

為檢驗?zāi)Ｐ驮诜瞧秸裁嫔袭愔亓鞯哪M性能，本節(jié)選取Marleau等[21]開閘式反坡異重流試驗作為驗證算例，試驗水槽如圖6所示，水槽長L1=1.975 m，寬b1=0.176 m，高H1=0.485 m。閘門與水槽左端的距離為Ll，水槽右端插入一塊剛性塑料板作為底坡、其長度為Ls。試驗時，閘門左側(cè)水槽底部填充高度為D、密度為ρ1的鹽水，周圍淡水的密度為ρ0。本次模擬的試驗參數(shù)為：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。經(jīng)網(wǎng)格無關(guān)解分析后，計算中選取nx×nz=600×30的網(wǎng)格。

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

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

4 結(jié) 論

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

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

參考文獻(xiàn)：

[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] 方春明，韓其為，何明民.異重流潛入條件分析及立面二維數(shù)值模擬 [J].泥沙研究，1997（4）：70-77.

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

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

[6] 陸俊卿，張小峰，崔占峰.各向異性密度流模型及其驗證 [J].水科學(xué)進(jìn)展，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，F(xiàn)LYNN 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