沖瀉區(qū)形態(tài)動力學耦合模型研究Ⅱ：模型建立與驗證

鄧 斌，蔣昌波，陳 杰，楊樹清

（1.長沙理工大學 水利工程學院，湖南 長沙 410114；2.水沙科學與水災害防治湖南省重點實驗室，湖南 長沙 410114；3.長沙理工大學 水科學與環(huán)境工程國際研究中心，湖南 長沙 410114；4.School of Civil，Mining and Environmental Engineering，University of Wollongong，Wollongong，Australia 2522）

1 研究背景

沖瀉區(qū)是由短波和長波引起的水位波動交替覆蓋和暴露的海灘區(qū)域，充分了解該區(qū)域水沙動力特性對研究岸灘侵蝕機制和海岸洪水運動規(guī)律至關重要[1]。由于該區(qū)域進行水動力[2]和岸灘演變[3]測量非常困難，對沖瀉區(qū)高瞬態(tài)、薄層、湍流和多相流等水動力特征及其與床面相互作用機制的認識非常有限。在近期舉辦的第二屆沖瀉區(qū)國際研討會上，眾多學者一致認為關于滲流對沖瀉區(qū)沉積物輸運的影響以及形態(tài)變化的模擬仍然是研究的重點。

近十年，隨著計算機技術進步，沖瀉區(qū)水動力學數(shù)值模型也得到了廣泛發(fā)展[4]。學者們基于雷諾平均N-S 方程[5]、大渦模擬[6-7]、淺水模型[1，8-9]和邊界層模型[10]等，開展了沖瀉區(qū)水動力特性的研究，對水深與流速的時空分布、岸線軌跡、紊動、波浪爬高、壓力梯度、邊界層流動等方面有了一定的認識。這些研究主要集中在沖瀉區(qū)水動力學數(shù)值模型上的改進[1]，較少關注泥沙輸運和岸灘形態(tài)的數(shù)值研究[11]。

目前沖瀉區(qū)形態(tài)動力學的數(shù)值研究，主要采用基于淺水非線性方程（NSWEs）的水深平均模型和簡化泥沙輸運公式相結合來模擬[1]，為掌握岸灘形態(tài)演變規(guī)律提供了有價值的成果。但這些模型中，部分為定床模型，僅考慮其水動力特性，忽略床面摩阻、泥沙運動以及岸灘形態(tài)的影響；部分為動床模型，一定程度上考慮床面摩阻以及床面變形的影響，但忽略了沖瀉區(qū)內復雜的水沙運動特性，如泥沙輸運采用簡單的經驗公式、忽略岸灘水流出滲/入滲效應等；同時，這些模型存在數(shù)值處理格式相對簡單、激波捕捉格式精度較低的問題，難以很好揭示沖瀉區(qū)內復雜水沙運動的普遍規(guī)律和動力特性[12]。

另一方面，沖瀉區(qū)的水流、泥沙運動及床面變形均有相對獨立的特征時間尺度，屬于典型的多重時間尺度問題，其求解方法主要有非耦合數(shù)值解法和耦合數(shù)值解法兩種[13]。非耦合數(shù)值解法是在每個時間步驟內順序求解水流和床面變形方程[14-15]，由于沖瀉區(qū)內水沙運動的高瞬態(tài)性，該方法尚不能充分描述水流與泥沙運動之間的相互作用機制。而耦合數(shù)值解法是同步求解水流和床面變形方程，在每個時間步內都計算水流和泥沙運動的特征參數(shù)值[16-18]，可有效解決上述問題。因此，充分考慮沖瀉區(qū)物理過程并結合耦合數(shù)值解法是當前開展沖瀉區(qū)形態(tài)動力學模擬的有效方法[9，13]。本文將基于分段輸沙率公式[19]，建立沖瀉區(qū)水動力及其岸灘形態(tài)動力學耦合分析模型，進一步研究沖瀉區(qū)水流與泥沙運動之間的相互作用機制。

2 模型基礎

沖瀉區(qū)內水深較淺，采用非線性淺水方程（NSWEs）能較好地描述沖瀉區(qū)內的水流運動特性[9]。因此，本研究基于1D NSWEs并結合床面變形Exner方程建立耦合分析模型。采用無量綱量來描述岸灘形態(tài)動力學理論模型：

式中：x為水平坐標；t為時間；h為水深；l為參考長度；u為水深平均速度；q為輸沙率（為h和u的函數(shù)）；h0為參考長度；q0為輸沙率參考尺度；g為重力加速度。為方便書寫，無量綱變量x*（水平坐標），t*（時間），h*（水深），u*（水深平均速度），B*（床面高程），q*（推移質輸沙率，為h和u的函數(shù)），分別采用x，t，h，u，B表示，則控制方程為：

3 模型建立

3.1 輸沙率公式首先考慮兩類廣泛應用于沖瀉區(qū)的輸沙率公式，即分別為假定推移質輸沙率為流速的函數(shù)q（u）與推移質輸沙率為流速和水深的函數(shù)q（h，u）。其次，采用已建立沖瀉區(qū)分段輸沙率公式封閉床面變形方程[19]。為便于計算分析，將6種推移質輸沙率公式均統(tǒng)一按q=q（h，u）的形式進行處理。表1給出不同輸沙率公式qb的量綱形式和無量綱形式，以及對應控制方程（2）—（4）中的q0、qh、qu和σ。6種輸沙率公式中的床面可動性參數(shù)A的取值詳見文獻[21]。

3.2 床面剪切應力現(xiàn)有考慮床面摩阻影響的剪切應力公式認為摩阻的影響相當于等效阻力，常采用半經驗公式（如：曼寧阻力公式、謝才公式）描述，其中采用謝才公式計算床面剪切應力在具有振蕩流特性的破碎區(qū)和沖瀉區(qū)獲得了較好結果[25]。本文亦采用謝才公式概化計算床面剪切應力：

表1 不同輸沙公式表達式及相關變量

式中：CD為阻力系數(shù)；h為水深；u為水深平均流速。

則原控制方程式（3）可改寫為：

式（6）中，在岸線位置，由于h→0，則對于沖瀉區(qū)的水流上爬階段，參考Antuono等[26]對岸線附近床面剪切應力項的處理方法，該項可通過hx平衡，即在波前處hx→-∞。在水流回落階段不能通過hx平衡；當u＜0時但不滿足hx→∞。因此，僅當u=0時，在水流回落階段，求解方程在岸線位置才成立。在水流上爬和回落階段，由于岸線位置不同，因此在數(shù)值處理時，可參考Zhu[9]在水流上爬階段采用外推法計算波前速度，在水流回落階段令岸線處u=0。

3.3 滲流源項現(xiàn)有岸灘形態(tài)動力學模型的研究多數(shù)忽略了滲流（含水流出入滲和沖流加減速產生的垂向流速）的作用。參照Delestre等[18]在坡面流數(shù)學模型中考慮降雨入滲的處理方式，在控制方程中考慮滲流源項：

式中：w為滲流流速，w＞0時表示水流出滲，w＜0時表示水流入滲。同時，對應動量方程中，也應體現(xiàn)滲流存在而引起的動量變化，即式（6）可進一步改寫為：

3.4 控制方程及數(shù)值耦合求解方法進一步考慮底坡、床面剪切應力、滲流等源項可得到描述沖瀉區(qū)水沙運動物理過程的模型控制方程：

為實現(xiàn)對上述方程的耦合求解，需同時考慮水流與床面變形的求解。由于上述控制方程可寫成不同的守恒形式，參考Hudson[28]對四種耦合求解方式下的計算結果，當在（hB）x=hBx+Bhx的形式下能得到較精確的結果。本文運用Hudson的做法，模型耦合求解的矩陣形式可寫為：

式中：U為守恒變量；F為x方向的通量；S為源項，分別為：

基于有限體積法數(shù)值離散上述守恒形式的方程組（13），令空間和時間區(qū)域屬于（x，t）∈（I，[0，T]），空間步長設為Δx，單元中心節(jié)點為xi=i×h，i=0，…，N，所分割單元為Ii=[xi-1/2，xi+1/2]。時間步長為Δxt=tn+1-tn，時間離散為0=t0＜t1＜t2＜…＜tn+1＜…＜tN=T，則時間分割單元為[tn，tn+1]。根據上述網格劃分，則各單元變量在時間tn的平均值為：

對式（13）在空間時間區(qū)域上（x，t）∈（I，[0，T]）進行積分，得到有限體積的離散形式：

運用式（15），式（16）可寫成：

4 數(shù)值模擬驗證

4.1 床面形態(tài)變化驗證——沙壩算例為檢驗數(shù)值方法計算地形變化的有效性，選取經典的“沙壩算例”[30]進行模型驗證。初始條件設置如圖1所示，設置長2000 m的計算范圍，并在床面附近設置一個小沙壩，其無量綱初始地形表達如式（18）所示，計算時忽略床面剪切應力和滲流的影響。

圖1 沙壩算列

式中：Δ=B/h0為地形無量綱擾動振幅。其中，參考長度尺度取B0，入口流速尺度為u0，時間尺度取Ts=（B0/g）1/2。計算參數(shù)設置如下：Δ=0.1，xL=2000 m，u0=1 m/s，h0=10 m，其他參數(shù)設置與Hudson等[30]和Postacchini等[31]保持一致，計算網格取Δx=2 m。

圖2為t=186070 s的地形變化結果?？梢?，采用TVD-WAF格式計算，在沙壩迎水面和背水面與解析解吻合較好，特別是在沙壩波峰處（x=570 m），TVD-WAF格式計算的結果與解析解非常接近，計算值與解析解的均方根誤差RSME=0.019，可應用于后續(xù)的相關計算。

圖2 地形變化結果（t=186070 s）

圖3 不同輸沙公式下地形變化結果（t=186070 s）

為檢驗不同輸沙率公式對結果的影響，計算分別采用了表1中6種輸沙率公式，計算結果如圖3所示。從圖中可見，采用q（u）形式的輸沙率公式時，Grass、Van Rijn、Bagnold和Meyer-Peter Müller等輸沙率公式計算得到的結果基本一致，且均與近似值接近，對應的均方根誤差RSME值分別為：0.0505、0.0605、0.0565和0.0221；采用q（h，u）形式的輸沙率公式（如Pritchard&Hogg輸沙公式）計算得到的結果局部稍偏離采用q（u）形式的計算結果以及近似解，在x=400 m附近計算結果出現(xiàn)負值，在沙峰位置大于其他公式的計算值，且沙壩的整體移動要慢于其他公式，計算得到RSME值為0.0937。采用分段輸沙率公式[19]計算得到結果與q（u）型輸沙公式基本一致，但在沙壩背水面優(yōu)于q（u）型輸沙公式的計算結果，對應的RSME值為0.0296。

4.2 ZD13 岸灘剖面變化數(shù)值模擬Zhu 和 Dodd[32]（下文簡稱 ZD13）在 Peregrine和 Williams[8]（簡稱PW01）的基礎上，基于特征線理論建立了一維岸灘形態(tài)動力模型，數(shù)值模擬PW01問題，并與PW01解析解進行對比。由于PW01解析解是基于定床得到，ZD13通過假定床面可動系數(shù)σ=1×10-7進行數(shù)值驗證，并在此基礎上進行動床條件下的數(shù)值模擬。本研究以PW01問題為基礎，考慮分段輸沙率公式及有無床面剪切應力作用，數(shù)值模擬單次沖流事件下沖瀉區(qū)內岸灘剖面的變化，并與ZD13的研究進行對比，驗證模型的可靠性。

圖4 單次沖流事件下沖瀉區(qū)概化計算模型（ZD13）

圖4給出了ZD13的概化計算模型，數(shù)學模型設置均與ZD13問題完全一致，令岸灘坡度λ=0.1，則床面高程B=γx，在x=0處把岸灘左側定為碎浪區(qū)，右側為沖瀉區(qū)。左側初始水深為h（x＜0，t=0）=1，速度u（x＜0，t=0）=0，右側初始水深與速度均設為0，并假定海側邊界x=-250 m，h（-250，t）=1，B（-250，t）=-250γ，u（-250，0）=0，u（-250，t）=-γt。計算首先考慮無床面剪切應力作用，即令CD=0（ZD13假定曼寧系數(shù)n=0），網格取Δx=0.002 m，時間步長取Δt=0.0002 s。

圖5為計算得到的無量綱水深、流速和床面變化值的時空分布與ZD13的結果對比。從圖中可見，當采用與ZD13一樣的輸沙公式時（即q=u3），本研究計算得到的h、u和ΔB與ZD13的計算結果基本一致。結果表明泥沙在水流上爬階段是向岸運動，在高沖瀉區(qū)產生淤積，而在回落階段更多的泥沙離岸運動，導致中沖瀉區(qū)和低沖瀉區(qū)被沖刷，與本研究實驗觀測到的現(xiàn)象以及與Jiang等[7]通過對沖瀉區(qū)水動力特性進行數(shù)值計算預測得到的沖淤趨勢相吻合。

圖5 無量綱水動力參數(shù)計算結果的時空分布與ZD13的結果對比

為研究床面剪切應力對岸灘形態(tài)變化的影響，采用3種不同的阻力系數(shù)（CD分別取0、1×10-3和1×10-2），其他初始條件設置保持不變。圖6為3種不同阻力系數(shù)下，計算得到的地形變化值與ZD13的計算值對比。其中ZD13的結果為2個工況（輸沙率公式分別為q=u3和q=hu3）計算值，ZD13的計算結果顯示，采用輸沙公式q=u3時，岸線附近的地形變化值較采用q=hu3得到的計算值小，高沖瀉區(qū)無明顯的淤積。而采用q=hu3計算得到高沖瀉區(qū)會產生少量淤積，這是由于公式q=hu3包含h和u的影響；在水流回落階段，高沖瀉區(qū)內水深和流速的逐漸減小會導致該區(qū)域較少的泥沙離岸輸運。從圖6可以看出，當阻力系數(shù)CD增大、對應床面阻力增大時，各計算工況得到的地形變化值隨之降低，即沖刷量均隨床面阻力的增大而變小，同時最大爬坡高度也隨之減小。值得注意的是，在不考慮床面剪切應力時（CD=0），采用q=hu3計算得到的最大爬坡高度達到26.34，即在岸線后出現(xiàn)一條狹長的薄層水流，而考慮床面剪切應力作用（CD≠0）會限制薄層水流的發(fā)展，從而導致水深增加，爬坡高度減小。ZD13采用兩種不同類型的輸沙率公式計算得到的結果具有明顯差別，表明在沖瀉區(qū)不同位置考慮不同的輸沙率公式，是精確模擬沖瀉區(qū)岸灘形態(tài)變化的有效途徑。

圖6同樣給出了采用本文分段輸沙率式在不同阻力系數(shù)下的岸灘地形變化值。從圖中可見，在岸線附近，計算值與ZD13采用q=u3計算得到結果接近，在高沖瀉區(qū)的計算值與ZD13采用q=hu3計算得到的結果基本一致，表明了本模型計算結果能充分反映低、中沖瀉區(qū)輸沙強度由流速主導、高沖瀉區(qū)輸沙強度由水深主導的內在機制。

圖6 不同阻力系數(shù)下本文計算得到的地形變化值與ZD13的計算值對比

4.3 實驗條件下岸灘剖面變化數(shù)值模擬基于鄧斌[21]的實驗數(shù)據，開展相應實驗條件的泥沙運動與岸灘演變研究。主要分析以下兩種情況：一是考慮摩阻系數(shù)是否為常數(shù)對計算結果的影響；二是考慮有無滲流對計算結果的影響，進一步評估模型在源項改變時的適應性。

數(shù)值模擬工況與實驗布置完全一致。計算網格取Δx=0.002 m，時間步長取Δt=0.0002 s，摩阻系數(shù)取CD=0.01，無滲流作用（w=0 m/s）。此外，為簡化研究，輸沙率公式中泥沙代表粒徑均以中值粒徑表示。圖7為數(shù)值模擬得到的3個工況下不同斷面平均水深計算值與實驗值的對比。圖8為數(shù)值模擬得到的3個工況下岸線軌跡計算值與實驗值的對比?？梢姡瑪?shù)值結果與實驗值吻合較好，表明所建立的數(shù)學模型能適用于沖瀉區(qū)水動力特性的模擬。

圖7 3個工況下不同斷面平均水深計算值與實驗值對比

圖8 3個工況下岸線軌跡計算值與實驗值的對比

4.3.1 不同摩阻系數(shù)對岸灘剖面變化的影響 前人計算結果表明將阻力系數(shù)CD設為常數(shù)可獲得較好結果，但部分學者認為在沖瀉區(qū)中采用常數(shù)形式的床面摩阻系數(shù)CD計算床面剪切應力是不合適的[1]。Barnes等[33]通過現(xiàn)場觀測發(fā)現(xiàn)摩阻系數(shù)在水流上爬階段約為回落階段的2倍，即CDu≈2CDb，其中CDu和CDb分別表示水流上爬和回落階段的摩阻系數(shù)，因此本研究分別考慮水流上爬和回落階段摩阻系數(shù)采用不同的取值?？紤]到實驗基于定床得到水流上爬和回落階段的摩阻系數(shù)相比Barnes等[33]等通過現(xiàn)場研究得到的結果偏于保守，因此取CDb=0.5CDu。同時，為保證床面泥沙運動的一致性，床面可動性參數(shù)需滿足σb=（CDb/Cd）3/2σ，即q∝τ3/2[23]。

數(shù)值模擬工況與上述驗證工況完全一致，源項僅考慮底坡和床面剪切應力的影響。圖9為動床3個工況第一次沖流后，分別考慮CDu=CDb=0.01和CDu=2CDb=0.01兩種不同摩阻系數(shù)影響下的岸灘剖面改變的計算值與實驗值對比，Case2下計算得到RSME值分別為0.0003和0.0024，Case3下計算得到的RSME值分別為0.0017和0.0017，Case4下計算得到RSME值分別為0.0024和0.0023，可見數(shù)值模擬結果與實驗值吻合較好，模型考慮CDu=2CDb=0.01時，計算值與實驗值更為接近。在單次沖流事件后，岸灘剖面變化均較小，Case2下岸灘剖面呈沙壩形狀，在靜水面以下（x=0.45～0.75 m）出現(xiàn)淤積，在低沖瀉區(qū)和中沖瀉區(qū)均出現(xiàn)沖刷，高沖瀉區(qū)略有淤積，Case3和Case4下沖刷和淤積均較小，地形變化不夠明顯。

圖9 不同摩阻系數(shù)下岸灘剖面改變量的計算值與實驗值比較

計算中分別采用水流上爬階段和回落階段不同床面摩阻系數(shù)與同一摩阻系數(shù)，對比兩種計算結果表明，前者計算得到的淤積沙壩高度相對較小，在低、中沖瀉區(qū)的沖刷也相對較少，并在高沖瀉區(qū)呈現(xiàn)薄層的淤積，與實驗觀測到的現(xiàn)象一致。當水流回落階段床面摩阻小于水流上爬階段時，由于阻力的減少，回落水流向岸運動更快，將會導致更多的泥沙離岸運動，沖刷將會增大；然而，水流回落階段σ的減少將會使泥沙運動減少，從而導致較少泥沙離岸輸運，在整個周期內淤積增大。當采用CDb=0.5CDu時，沙壩淤積量有所減少，低、中沖刷量略有減少，在高沖瀉區(qū)淤積量有所增大，這表明CD和σ的聯(lián)合改變導致岸灘沖刷量減少，在高沖瀉區(qū)會產生凈淤積。因此在沖瀉區(qū)岸灘形態(tài)變化模擬時，考慮不同的床面剪切應力是非常有必要的。

4.3.2 有無滲流對岸灘剖面變化的影響 為討論滲流對岸灘剖面變化的影響程度，利用實驗數(shù)據并結合數(shù)值模擬進行定性分析。數(shù)值模擬以Case2為例，源項完全考慮底坡、床面剪切應力和滲流的影響，其中床面剪切應力項中摩阻系數(shù)按CDb=0.5CDu計算。在原實驗條件下，在整個沖流周期內岸灘出現(xiàn)了水流入滲現(xiàn)象，但未出現(xiàn)明顯的出滲現(xiàn)象[19]。因此，在水流上爬和回落階段均只考慮水流的入滲作用，且其取值假定為常數(shù)，w的取值采用平均值，即w=-0.001 m/s。

圖10 有無滲流影響下岸灘剖面形態(tài)變化計算值與實驗值比較（Case2）

圖10分別為Case2第一次沖流后，分別考慮有無滲流影響下（w=0 m/s和w=-0.001 m/s）的岸灘剖面形態(tài)與地形沖淤量的計算與實驗值對比。從圖10可以看水流入滲下，岸灘剖面分布在高沖瀉區(qū)呈現(xiàn)較明顯的淤積，沙壩區(qū)淤積量則略有減少，在低、中沖瀉區(qū)的沖刷量稍有增加，表明水流入滲情況下，泥沙向岸輸運量相應增加，離岸輸運則有所減少，與Butt等[34]、Masselink和Li[35]和Li等[36]等的結果趨勢一致，進一步說明滲流對岸灘形態(tài)的變化存在影響。但圖中計算值與實驗值存在一定的差異，這是由于計算中僅概化滲流為入滲，未全面考慮實際滲流的沿程變化。因此，有必要進一步獲取整個沖瀉區(qū)岸灘滲流流速特征，為公式的完善提供更全面的基礎數(shù)據。

4 結論

本文基于一維非線性淺水方程和床面變形方程，進行了沖瀉區(qū)水動力特性和岸灘形態(tài)變化耦合求解模型的研究，主要結論如下：

（1）基于一維非線性淺水方程和床面變形方程，考慮兩類（共6種）不同輸沙率公式，增加了床面剪切應力源項和滲流源項，建立了描述沖瀉區(qū)復雜水沙動力特性的岸灘形態(tài)動力學耦合計算模型，并對模型可靠性進行了驗證。

（2）采用分段輸沙率公式，模擬分析了單次沖流事件下沖瀉區(qū)內岸灘剖面形態(tài)變化。模型計算得到的水深、流速和床面形態(tài)特征等與ZD13的結果基本一致。計算結果能充分反映低、中沖瀉區(qū)輸沙強度由流速主導、高沖瀉區(qū)輸沙強度由水深主導的內在機制，計算結果較ZD13的結果更為準確。

（3）數(shù)值模擬了潰壩波作用下沖瀉區(qū)水動力與岸灘剖面形態(tài)變化。計算得到的水深、岸線軌跡和岸灘剖面變化結果與實驗結果吻合良好。床面剪應力源項的改變對沖瀉區(qū)的水沙動力變化影響較大，沖瀉區(qū)內岸灘形態(tài)變化應充分考慮水流上爬和回落階段摩阻系數(shù)的不同以及滲流的影響。

研究建立的岸灘形態(tài)動力學模型為一維模型，下階段有必要考慮更高精度的數(shù)值求解格式和紊流模型，基于氣、液、固三相流的完全耦合數(shù)學模型是未來的研究方向。