章若茵，吳保生

（清華大學 水沙科學與水利水電工程國家重點實驗室，北京 100084）

1 研究背景

異重流是兩種比重相差不大的流體，因為比重差異而發(fā)生的相對運動。渾水入庫形成的泥沙異重流是水庫中常見的現(xiàn)象，其潛入過程以及流速與含沙量沿垂向的分布規(guī)律等是水庫異重流研究的重要內容，對優(yōu)化水庫調度、提高水庫運用效益具有重要意義。但由于異重流潛入?yún)^(qū)域的變動以及異重流沿庫底流動的特點，在天然水庫中觀測完整的異重流潛入以及流速和含沙量的垂線分布較為困難，前人多以水槽試驗或數(shù)值模擬的結果來探討其分布的規(guī)律。

范家驊［1］、蘆田和男［2］、曹如軒等［3］、方春明等［4］、焦恩澤［5］、李濤等［6］通過水槽試驗、水力分析和數(shù)值計算等手段對不同比降、不同粒徑、不同水沙條件下的異重流潛入進行研究，分析了潛入點Froude 數(shù)的大小以及潛入水深與流量、含沙量、底坡的關系。但多數(shù)只分析了潛入水深與潛入點所在位置的流速、含沙量關系，而董炳江等［7］利用進口密度Froude 數(shù)表征入口水流的慣性力與浮力作用之比，受到入口流量和含沙量的影響，并根據(jù)數(shù)值模擬結果得到潛入水深與進口水深之比會隨著入口Froude 數(shù)的增加而增加的認識。這一變化趨勢雖然在前人［8］根據(jù)實測資料、水槽試驗和數(shù)值模型得到的結論范圍內，但缺乏定量描述，且入口Froude 數(shù)變化在0.5 ～8 的范圍內，當入口Froude 數(shù)在更大范圍內變化時這一關系是否存在仍不確定。

異重流潛入后的流速和含沙量分布與明渠水流完全不同，異重流的流速分布中存在最大流速點，一般認為其在異重流高度的0.2 ～0.3 倍之間，這與異重流水沙交界面的確定方法有關，也與床面粗糙程度有關，當床面粗糙程度大時最大流速點的位置會提高，因為受到下邊界的阻力作用增強［9］。而含沙量的分布較為復雜，垂直分布與否及垂直分布的厚度等與泥沙粒徑、濃度以及沖淤狀態(tài)均有關系，并沒有較為一致的結論。焦恩澤［10］認為交界面以下的流速分布近乎對數(shù)分布，而含沙量在一般含沙和高含沙量時有不同的分布型式。Peakall 等［11］、Kneller 等［9］列舉了高含沙量、低含沙量、沖刷狀態(tài)以及非均勻沙情況下4 種不同形式的異重流含沙量分布。Lee 等［12］根據(jù)水槽試驗結果認為無論是急流還是緩流情況，異重流的流速和含沙量的分布都較一致，與水沙條件無關，但由于只在每組試驗中觀測了兩個斷面，斷面Froude 數(shù)只在1.05 ～1.40 范圍內，此結論存在局限性。而Sequeiros［13］認為異重流流速和含沙量的分布與水沙條件、床面形態(tài)等均有關，當斷面Froude 數(shù)減小時，最大流速的位置提高，含沙量在靠近河床處垂直分布的區(qū)間更厚，其Froude 數(shù)變化范圍在0.41 ～2.21 之間。由此可見，對于異重流流速和含沙量分布隨水沙條件的變化規(guī)律仍然存在一定的爭議，值得進一步研究。

本文選擇SCHISM 模型對Lee 和Yu［12］的水槽試驗進行三維數(shù)值模擬，不僅能表現(xiàn)異重流強烈的三維特性，明顯優(yōu)于以往采用一維模型［14］和立面二維模型［15］對此進行的模擬，而且能根據(jù)三維模擬結果分析不同時間和斷面的異重流運動情況，補充原試驗分析的不足。SCHISM 是基于雷諾時均方程（RANS）的三維數(shù)值模型［16-18］，相比于直接數(shù)值模擬（DNS）和大渦模擬（LES）等僅適用于實驗室尺度的模型，SCHISM 能夠應用于實驗室和實際工程等各個尺度的模擬。SCHISM 在水平方向上采用無結構化網(wǎng)格，能夠靈活適應復雜的邊界條件，優(yōu)于FLOW-3D 等模型中的結構化網(wǎng)格［19］；垂直上則包括SZ坐標［20］和LSC2坐標［21］兩種方式，兩種方式均有良好的垂向邊界跟隨性，其中SZ 坐標結合了z 坐標和張芝永等［22］、Pérez-Díaz 等［23］在三維模型中采用的σ坐標，一定程度上消除了σ坐標在底坡突變的情況時產(chǎn)生的不合理的跨密度面混合和壓力梯度誤差，同時LSC2坐標也可消除這種誤差［24］，因此均適合運用于發(fā)生在底坡陡降、水深突增情況下的水庫異重流。此外，在引用隱式TVD2格式計算物質運輸后，SCHISM 在模擬鹽水入侵、密度分層等現(xiàn)象時更加精確［25-26］。但目前來看，SCHISM 模型多用于河口處鹽度或溫度異重流模擬，而缺乏模擬水庫泥沙異重流的研究。

本文將首先通過對水庫異重流水槽試驗的模擬，驗證該模型模擬水庫異重流的可靠性和準確性；然后根據(jù)模擬結果，分析不同入口水沙條件下異重流潛入點的變化規(guī)律和潛入后含沙量、流速分布的規(guī)律，探討入口水沙條件與潛入水深之間的定量關系。研究結果不僅為下一步運用SCHISM 模型模擬大型水庫泥沙異重流奠定基礎，而且加深了對異重流水沙運動和分布規(guī)律的認識。

2 三維水沙數(shù)學模型

2.1 控制方程與湍流模型SCHSIM 模型基于RANS 方程，并且滿足靜壓假定和Boussinesq 渦黏性假定，其在笛卡爾坐標下的控制方程如下：

連續(xù)方程：

動量方程：

物質輸移方程：

狀態(tài)方程：

式中：x、y 為水平笛卡爾坐標，m；z 為垂向坐標，向上為正，m；u、v、w 分別為x、y、z 三個方向的流速，m/s；t 為時間，s；h 為水深，m；η為自由水面水位，m；f 為柯氏力系數(shù)，s-1；ρ0、ρ分別為參考密度和水體的密度，kg/m3；g 為重力加速度，m/s2；Kmh、Kmv分別為水平與垂直渦黏性系數(shù)，m2/s；pA為自由水面的大氣壓強，N/m-2；C 為物質濃度，此處為含沙量SSC，kg/m3；ωs為泥沙顆粒沉速，m/s；Ksh、Ksv分別為水平和垂直方向的泥沙擴散系數(shù)，m2/s；q 為輸運物質的源匯項。

式（1）—式（3）采用有限元和有限體積法進行求解，式（4）采用有限體積法求解。當式（4）用于計算其他變量（鹽度S/PSU、溫度T/℃等）時，則沒有泥沙沉速這一項。本文中暫時只考慮泥沙造成的異重流現(xiàn)象，因此模型中的鹽度設置為0，而入流水體和環(huán)境水體的溫度保持一致。

為了求解式（1）—式（4）需要采用湍流模型閉合，模型中采用Umlauf 和Burchard［27］的通用長度模型GLS（Generic Length-scale model），包括紊動能K 和通用長度變量ψ的控制方程，其中通用長度ψ由下式定義：

2.2 邊界條件針對式（2）—式（3），水面上不考慮風應力的影響時，水平流速的垂向梯度為0；床面上根據(jù)水體底層的雷諾應力與床面摩擦剪應力的平衡關系給出：

式中τbx、τby為床面的摩擦剪應力，m2/s2，由湍流邊界層和阻力系數(shù)確定［31］。對式（4）水面和床面上的邊界條件為：

式中：Cb為床面含沙量，kg/m3，可根據(jù)文獻［32］的方法計算；α*為非平衡輸沙響應系數(shù)，表示床面含沙量與垂向平均含沙量之比，可隨水沙條件動態(tài)調整；S*為垂向平均挾沙力，kg/m3。與方程（9）相對應的河床變形方程為：

式中：ρ′為河床組成物質的干密度，kg/m3；Z 為床面厚度，m。

3 水槽異重流試驗模擬

3.1 水槽試驗條件為了研究異重流潛入后垂向上流速和含沙量分布的特點，Lee 和Yu［12］針對水庫異重流進行了一系列水槽試驗。試驗所用的水槽長20 m，寬0.2 m，坡度為0.02，通過設置出口水位保持明渠段長4 m，壅水段長16 m，以及斜坡尾部的蓄水池長1 m，其側面示意圖和高程如圖1 所示。選用的泥沙是粒徑為0.0068 mm 的均勻沙，粒徑較細，泥沙比重2.65。

選擇水槽試驗中的TC01—TC18 組進行模擬，每組試驗的入流水沙條件、溫度如表1。不同的入口水沙條件可用進口密度Froude 數(shù)來表示：

式中：Fr0為入口密度Froude 數(shù)；U0、h0和q0分別為入口的流速、水深和單寬流量；g′為有效重力；Δρ為入口渾水（ρin）與清水（ρa）的密度差，各組對應的Fr0數(shù)也列于表1。

圖1 水槽設置示意圖及網(wǎng)格高程

表1 TC01—TC18 組水槽試驗入庫水沙條件

3.2 模型設置與選擇考慮到矩形水槽較為規(guī)則，本文在水平向采用四邊形網(wǎng)格，尺度為2.5 cm×2 cm，則網(wǎng)格數(shù)8400，節(jié)點數(shù)8251。模擬時間步長為0.1 s，每1 s 輸出結果，模擬時長810 s。選用清華大學“探索100”集群進行計算，由于模擬尺度較小，每次計算使用12 個核［33］。由于垂直分層的方式和層數(shù)對于數(shù)值模擬的計算精度影響較大，所以根據(jù)SCHISM 兩種垂向網(wǎng)格方式，設置3 種不同的垂直分層方式如表2 所示，分別模擬TC14 組的異重流過程，從而確定最終的垂直分層方法。

表2 不同垂直分層方法的層數(shù)要求

本文以垂線上流速為0 的點為清渾水交界面的位置，則TC14 組交界面以下實測及模擬的流速和含沙量分布如圖2 所示，沿程含沙量分布如圖3 所示。由圖2（a）可知，3 組模擬的流速分布均出現(xiàn)了最大流速點，且最大流速點上、下區(qū)域的流速分布規(guī)律不同。A 組模擬的流速分布與實測情況最為吻合，在距水槽底部4.8 cm 處達到最大流速點10.32 cm/s，略大于實測的10.27 cm/s；B 組最大流速點的位置在距底3.8 cm 處達到10.05 cm/s，但最大流速點以下的部分與實測點相差較大；C 組的流速分布明顯與實測分布差別較大。由圖2（a）可以看到，實測最大流速點以下只有一個實測數(shù)據(jù)點，為了彌補實測數(shù)據(jù)點較少的問題，利用Geza 和Bogich［34］提出的異重流最大流速點以下的流速分布公式進行了補充計算，其具體計算公式如下：

式中：ue為高度z 處對應的流速；u′em為最大流速；h1e為最大流速點到河底距離； u*為摩阻流速；κ為卡門常數(shù)，與最大流速點以下的平均流速有關，經(jīng)計算此處為0.27。

根據(jù)式（12）得到異重流最大流速點以下的流速分布曲線，見圖2（a）?？梢钥吹剑碚摴角€具有與A 組相近的分布，說明采用A 組模擬的流速分布是較為可靠的。需要說明的是，根據(jù)邊界無滑動條件假定，垂向流速在河底處均為0，而文中所有的流速分布均只繪制到河底上一層流速，因此并不為0。

圖2 不同垂直分層方法模擬的流速和含沙量分布（TC14）

圖3 不同垂直分層方法模擬的TC14 組沿程含沙量分布（t=150 s）

由圖2（b）可知，A 組和C 組的含沙量與實測分布較為接近，含沙量在交界面以下迅速變大，且在靠近底部的一定厚度呈現(xiàn)垂直分布的情況，與三門峽水庫實測高含沙量異重流的含沙量垂線分布十分相似，主要是由于顆粒較細、含沙量較高產(chǎn)生的黏性作用，導致含沙量呈現(xiàn)幾乎垂直分布的特點［10］，一般稱這個垂直分布的區(qū)間為致密層。而B 組的含沙量分布明顯較為均勻，圖3（b）的沿程含沙量分布沒有明顯的異重流潛入點，這可能是因為垂向分層過少導致模型計算時垂向擴散劇烈，異重流潛入過程不明顯。綜合流速和含沙量分布來看，A 組的分層方法具有較明顯的優(yōu)勢，因此最終選擇A 組方式進行其余組次的模擬。

3.3 模型驗證采用上述配置模擬TC01—TC18 組水槽試驗后，以TC14 組為例，渾水進入壅水段后形成異重流并穩(wěn)定運動的完整過程如圖4 所示。渾水從明渠段進入壅水區(qū)后，會在滿足異重流潛入條件的位置處下潛形成異重流。圖中顯示約40 s 時，接近河床位置的含沙量變大，而表層含沙量較小。異重流有開始向下潛入的趨勢，隨著時間變化，異重流潛入的位置逐漸向下游移動，直到一個固定的位置，形成穩(wěn)定的異重流潛入點［35-36］。圖5 是160 s 時潛入點附近的流場分布圖，可以看到，異重流潛入點及清渾水交界面處流速接近為0，而上層清水流速較小，且有反向流速，同時異重流內部的流速都較大。

為了分析異重流潛入以后垂線上流速和含沙量的分布特點，可提取潛入點下游斷面垂線上的流速和含沙量，根據(jù)式（13）得到該垂線上的異重流平均厚度、平均流速和平均含沙量［37］：

圖4 異重流潛入及運動時的沿程含沙量分布（TC14）

圖5 潛入點附近流速（t=160s）

式中：U 為不同高度的流速；Ugc為異重流平均流速；hgc為異重流平均厚度；C 為不同高度的含沙量； Cgc為異重流平均含沙量；δ表示垂線上流速為0 的位置與水槽底部之間的距離。

利用計算的異重流厚度hgc、平均流速Ugc和平均含沙量Cgc代入式（11）可以得到此時該斷面上異重流的Froude 數(shù)Frd。

原試驗在潛入點下游測量并計算了兩個斷面的相關數(shù)據(jù)，本文從模擬結果選取相同斷面的流速和含沙量得到hgc、Ugc和Cgc，實測與模擬結果的對比如圖6 所示。由圖6 可知，hgc、Ugc和Cgc的R2分別為0.95、0.92 和0.94，三者的符合程度都較高，說明SCHISM 模擬的異重流潛入后的運動過程與實測過程較為一致，模擬結果基本可靠。

圖6 實測與模擬異重流平均流速、厚度和平均含沙量對比（TC01—TC18）

4 異重流影響因素分析

4.1 入口含沙量的影響由表1 中TC14—TC17 的入流水沙條件可知，這4 組流量接近，但入流的含沙量逐漸增加，可以對比分析入流含沙量對異重流流速和含沙量分布的影響。因160 s 時異重流運動受到水槽尾部的回流影響較小，所以以該時刻距入口14 m 的斷面為例，據(jù)式（11）和式（13）計算不同組次的hgc、Ugc、Cgc和對應的Frd數(shù)列于表3，該時刻沿程的含沙量分布如圖7，分別以表3 中的hgc、Ugc和Cgc為參考值，可以得到在交界面以下流速和含沙量隨深度的無量綱分布，見圖8。

表3 不同入流含沙量對異重流的影響

圖7 入口含沙量不同時異重流的含沙量沿程分布（t=160 s）

由表3 和圖7 可知，不同的入口含沙量導致異重流的運動速度不同，含沙量越大，渾水與環(huán)境水體的密度差增大，異重流頭部的能量與動量越大，則異重流運動速度越快，14 m 斷面上的平均流速從0.090 m/s 增長到0.131 m/s。但同樣環(huán)境水體形成的反向流速也較大，對異重流的壓力增大，異重流環(huán)境水體的內外壓力差增大，導致異重流厚度越小，TC14 時該斷面的異重流平均厚度達0.145 m，而TC17 中平均厚度只有0.102 m，與圖7 中的分布一致，也與謝鑒衡［38］、董炳江等［7］的結論一致。此外，從圖7 還能看出同一時刻，入口含沙量越大，則異重流潛入點的位置越向上游，潛入水深越小。異重流潛入以后，基本保持一個較為穩(wěn)定的厚度持續(xù)向下游運動，但是頭部的厚度略大，主體部分的厚度則略小。

圖8 不同入口含沙量條件的無量綱化流速和含沙量分布

結合表3 中的Frd數(shù)和圖8 中的無量綱化的流速和含沙量分布來看，這4 組的Frd數(shù)從1.15 增長到1.28，流速和含沙量的分布有區(qū)別但不明顯，可能是因為Frd數(shù)的增長幅度較小。最大流速點的相對位置在0.38 ～0.42 附近波動，而最大流速從Ugc的1.20 倍增長到1.23 倍；致密層的厚度從hgc的0.47 倍降低到0.3 倍，對應的含沙量則從Cgc的1.15 倍增加到1.22 倍。雖然存在波動，但從變化趨勢來看，最大流速和致密層的含沙量均會隨Frd數(shù)的增大而增大，而最大流速點的位置和致密層厚度則會隨Frd數(shù)的增大而減小，且含沙量的變化比流速的變化更敏感。較大的Frd數(shù)說明水流慣性更大，則水流持續(xù)向下運動的能力越強，最大流速點越低。

4.2 入口流量的影響表1 中TC02、TC04、TC13 和TC18 組的流量逐漸增加，而含沙量接近，因此對比這4 組的結果可以分析入庫流量對異重流的影響。由于入流流量不同，這4 組異重流前進速度相差較大，因此保證異重流前峰運動到大致相同的位置時，每組的時間不同，分別選取該4 組試驗在268 s、213 s、154 s 和136 s 時斷面14 m 的結果，不同組次的異重流厚度、平均流速、平均含沙量和對應的Frd數(shù)列于表4，沿程的含沙量分布如圖9，分別以表4 中的hgc、Ugc和Cgc為參考值得到的流速和含沙量無量綱分布如圖10 所示。

表4 不同入庫流量對異重流的影響

圖9 入口流量不同時異重流的含沙量沿程分布

圖10 不同入口流量條件的無量綱化流速和含沙量分布

由表4 和圖9 可知，入庫流量越大，則異重流潛入位置越向下游移動，潛入水深越大，這也與前人發(fā)現(xiàn)的規(guī)律一致［4］。并且入流流量越大，動能越大，異重流厚度越大，平均厚度從0.079 m 增加到0.101 m，異重流前進速度也越大，平均流速從0.081 m/s 增加到0.103 m/s，但是TC18 的平均流速較小，這可能是由于該組異重流潛入點距入口較遠，異重流尚未充分發(fā)展至較穩(wěn)定的狀態(tài)即到達水槽尾部，因此14 m 處的水沙狀態(tài)可能受到異重流頭部強烈紊動作用的影響，而前三組的14 m 處基本已經(jīng)處于異重流主體較穩(wěn)定的狀態(tài)，因此TC18 組的結果不一定符合規(guī)律。

結合表4 中的Frd數(shù)和圖10 中的無量綱化流速和含沙量分布可知，如TC02 的Frd數(shù)較大為1.22，則最大流速點位置低，只達到hgc的0.28 倍，而最大流速為Ugc的1.33 倍，致密層的相對厚度越小，甚至在接近床底時才達到最大含沙量，但最大含沙量能達到Cgc的1.66 倍；而TC18 的Frd數(shù)為0.85，最大流速點能夠達到hgc的0.46 倍，而最大流速減小為Ugc的1.19 倍，致密層的相對厚度可達0.66，而致密層的含沙量則只有Cgc的1.05 倍。

4.3 Froude 數(shù)對異重流潛入的影響前文分別分析了入口的流量和含沙量對異重流潛入位置的影響，兩者的影響各不相同。下面以入口的Fr0數(shù)作為特征量綜合考慮入口水沙條件對潛入點造成的影響。

分析各組的異重流潛入位置、式（11）和表1 可知，入口含沙量越大時對應的Fr0數(shù)越小，而入口流量越大時對應的Fr0數(shù)越小，后者看似與式（11）矛盾，實際上若保持入口水深h0不變，則入口流量增加會導致Fr0數(shù)增加，但由于本文的模擬中沒有保持入口水深h0的不變，因此入口流量增加的同時會造成h0變大，而h0的增加又會造成Fr0數(shù)的減小，因此綜合考慮入口流量和水深的影響后，發(fā)現(xiàn)入口流量的增加最終會減小Fr0數(shù)。入口含沙量和流量的增加都會導致Fr0的減小，但前者導致潛入水深Hps減小，后者卻導致Hps增加，這種矛盾可能是由于h0對Fr0數(shù)的影響大于流量的影響產(chǎn)生的?？紤]到Fr0數(shù)為無量綱數(shù)，則以潛入水深和入口水深的比值Hps/h0作為無量綱參數(shù)，根據(jù)本文18 組模擬結果得到圖11 所示Hps/h0與Fr0的關系，可用式（14）來表示，且R2能夠達到0.90：

圖11 Hps/h0與入口Froude 數(shù)關系

圖11 還給出了前人所得Hps/h0與Fr0的關系［8］，雖然它們分別表示水槽試驗、數(shù)值模擬或現(xiàn)場實測等不同環(huán)境中的異重流，且底坡不同、入口水沙條件不同或產(chǎn)生異重流的原因不同，如溫度、鹽度異重流等，但是Hps/h0與Fr0始終遵循冪函數(shù)相關的關系，且本文結果在入口Froude 數(shù)大于8 時仍然滿足這一關系，說明這一關系在異重流運動中是普遍存在的，即可以率定式（14）在不同異重流的參數(shù)后根據(jù)入口流量、水深和含沙量確定異重流潛入的位置。但是是否存在一個普遍適用的統(tǒng)一公式可以包括不同情況的異重流仍需要更多的數(shù)據(jù)進行驗證和進一步的研究。

5 結論

（1）本文采用SCHISM 三維水沙數(shù)學模型模擬了水槽異重流試驗，選擇LSC2的垂直分層方法以及最多51 層的層數(shù)保證異重流模擬的精度及效率。模擬結果與實測異重流平均厚度、平均流速和平均含沙量吻合較好，對應R2均能達到0.9 以上，同時潛入點下游斷面的流速和含沙量的垂線分布也符合實測的分布情況。

（2）異重流潛入后的運動和水沙分布規(guī)律與入口水沙條件相關，入口含沙量越大，異重流運動速度越快，厚度越?。蝗肟诹髁吭酱?，異重流運動速度越快，厚度越大。流速和含沙量的無量綱分布與斷面Froude 數(shù)密切相關，斷面Froude 數(shù)從0.85 增長到1.22 時，最大流速從1.19 增長到1.33 倍平均流速，致密層含沙量從1.05 增長到1.66 倍平均含沙量，而最大流速點與河床距離從0.46 減少為0.28倍平均厚度，同時致密層厚度從0.66 倍平均厚度減小為0。

（3）異重流潛入點的變化與入口水沙條件相關，入口含沙量越大，入口流量越小，則異重流潛入位置越向上游，潛入水深越小。潛入水深與進口水深的比值Hps/h0隨入口Froude 數(shù)的增加而增加，建立了兩者之間的冪函數(shù)關系，可以用于確定異重流潛入點的水深。但是否存在一個普遍適用的統(tǒng)一公式來描述不同異重流的情況還需要進一步研究，為推廣到實際水庫中預測潛入點的水深和位置提供依據(jù)。