考慮異重流的紫坪鋪水庫一維水沙數值模擬研究

由麗華，歐陽莉，趙 陽，王 強，楊 飛

（1.四川省紫坪鋪開發有限責任公司，四川 成都 610091；2.黃河水利委員會 黃河水利科學研究院，河南 鄭州 450003；3.河南省湖庫功能恢復與維持工程技術研究中心，河南 鄭州 450003）

河道內挾帶大量細沙的水流進入水庫后，含沙水流因與庫內清水存在密度差而潛入水庫底部，形成水庫異重流。水庫異重流對于水庫的淤積總量以及淤積形態影響深遠［1］，因此科學把握水庫異重流運動過程可有效判別水庫的淤積形態，同時可指導多沙水庫泥沙管理。目前，對于水庫異重流的研究主要分為理論研究、現場觀測以及數學模型試驗等方面。基礎理論研究主要集中在異重流形成條件、潛入點位置及持續運行條件3個方面。不少學者通過水槽試驗對異重流的潛入過程進行研究，代表學者包括范家驊［2］、曹如軒等［3］、焦恩澤［4］。Xia等［5］為了研究異重流潛入點判別公式，重點考慮了清水弗勞德數以及含沙量影響，并通過前人的水槽試驗結果擬合了新的潛入點判別公式。范家驊等［6－7］全面回顧了異重流潛入現象的各種成果，并對不同進口水沙條件下異重流潛入點的水沙特性、潛入點與水沙因子的關系進行了探討。現階段對異重流現場觀測研究多圍繞小浪底水庫開展，自2002年首次開展調水調沙試驗后［8］，黃河水利委員會多次開展黃河調水調沙運行。經過多年的調水調沙，黃河下游河道主槽不斷萎縮的態勢得到了初步遏制，在防洪減淤方面取得了重要進展［9］。此后，也有不少學者針對黃河調水調沙期間小浪底水庫異重流運動特性進行研究，包括潛入點及流速分布［10］、壩前渾液面沉降［11］、壩前異重流厚度［12］以及排沙效率［13］等。由于異重流測量需要把握時機，且現階段其他水庫較少有現場觀測的資料，因此構建數學模型研究水庫異重流運動過程也是目前較為主流的研究手段之一。熊正偉［14］構建了同時適用于異重流及溯源沖刷的水沙輸移與河床沖淤耦合的一維數學模型，并對1975年三門峽水庫溯源沖刷過程進行了模擬，分析了溯源沖刷發展特點及規律。常云華等［15］提出采用一維非飽和不平衡輸沙異重流數學模型進行水庫泥沙淤積計算，在黃河潼關至三門峽庫區段進行計算的結果良好。王增輝［1］綜合考慮干支流倒灌因素，對水庫溯源沖刷以及異重流過程進行模擬，同時改進淺水二維水沙耦合模型對不同邊界條件下的水庫支流口溯源沖刷進行模擬研究。筆者針對近年來紫坪鋪水庫來沙偏多、水庫淤積日漸嚴重的現狀，開展異重流與明流耦合水沙數學模型構建工作，并利用2019年實測水沙數據進行驗證，以期為紫坪鋪水庫開展水庫清淤與泥沙調度提供指導。

1 研究區域與數據

紫坪鋪水利樞紐工程位于四川省岷江上游映秀鎮至都江堰河段，是一座以灌溉和供水為主，兼有防洪、發電、環境保護和旅游等綜合效益的一等大（I）型水利水電工程，也是都江堰灌區和成都市的重要水源調節工程。壩址以上流域面積2.26萬km2，多年平均徑流量148億m3，占岷江上游徑流總量的97%。水庫正常蓄水位877.0 m，相應庫容9.98億m3；校核洪水位883.1 m，相應庫容11.12億m3；死水位817.0 m，相應庫容2.24億m3；水庫調節庫容7.74億m3。水庫設計多年平均懸移質沙量792萬t，含沙量0.572 kg／m3。多年汛期（6—10月）平均沙量778.9萬t，占全年沙量的98.3%；5—9月沙量為760.7萬t，占全年沙量的96%；7—9月沙量為533.43萬t，占全年沙量的67.4%。多年平均入庫推移質70萬t，最大粒徑為500 mm，中值粒徑為113 mm。紫坪鋪水庫入庫控制站為映秀水文站，出庫控制站為紫坪鋪水文站。

紫坪鋪水庫2019年來水偏多，形成幾次較大的洪峰，最大洪峰為1 792 m3／s。受降雨影響，上游發生多起泥石流以及滑坡現象，造成罕見的入庫沙峰，最大沙峰為63.4 kg／m3。圖1為紫坪鋪水庫2019年汛期進出庫流量與水位變化過程，由圖1可知，2019年汛期（6—10月）水量總體偏大，形成幾次不小的洪峰，水位按照調度規程，圍繞850 m調度。然而，如圖2所示，入庫泥沙只出現了一次較大的峰值（8月20日），出庫含沙量出現了歷史罕見的63.4 kg／m3。由于2019年有異重流到達壩前，因此選用2019年的數據進行模型驗證。

圖1 紫坪鋪水庫2019年汛期進出庫流量與水位變化過程

圖2 紫坪鋪水庫2019年汛期進出庫含沙量變化過程

2 一維水沙數學模型構建

基于水流連續性方程、水流動量方程、泥沙連續方程、河床變形方程以及異重流運動方程構建紫坪鋪水庫一維水沙數學模型。天然水庫的邊界條件十分復雜，為了便于求解，本研究對部分條件做如下簡化：①假定可以將最小計算時段內的非恒定流作為恒定流處理，具體是將進口斷面的實際流量過程線改為若干個不同流量級組成的梯級過程線進行計算，對于每一個梯級來說，流量為常數，水流為恒定流；②假定河床沖淤過程中每一個短時段內河床變形對水流條件影響不大；③在最小計算時段內水體中含沙量恒定。天然庫區的泥沙是非均勻沙，因此在實際的模擬計算中，需要將泥沙顆粒級配按照粒徑大小進行分級處理，最終按照分級級配計算水流的挾沙能力。其中，水流泥沙基本運動方程見文獻［15］，異重流模塊主要分為計算潛入點的水深、異重流均勻流水深、判斷異重流發生條件以及異重流的運動方程，具體計算過程如下。

（1）計算潛入點的水深h0。公式為

式中：h0為潛入點水深；Q為流量；ηg為重力修正系數；g為重力加速度；B為河寬；S為含沙量；γs為泥沙容重，取2 650 kg／m3；γ0為清水容重，取1 000 kg／m3。

（2）計算異重流均勻流水深hn。公式為

式中：hn為均勻流水深；λ為異重流阻力系數，取λ＝0.025～0.03；J為底坡（實際為能坡）。

（3）判斷水深h。當異重流潛入且能以均勻流形式運動時，其潛入點水深h＞max（h0，hn），若滿足，則可以認定異重流潛入，并能以均勻流方式運行，不上升到水面（式中：A為過水面積；i－1和i分別為前后兩個計算時間點）。

根據上述方程以及水流泥沙基本運動方程，將紫坪鋪水庫一維水沙數學模型的計算步驟總結如下：①輸入邊界條件、初始條件及有關參數；②計算各水力要素；③計算不同斷面的平均含沙量；④計算子斷面含沙量；⑤計算子斷面的沖淤面積；⑥修正子斷面的河床組成；⑦按照相鄰子斷面的沖淤厚度以及過水斷面的寬度所占的權重，修正節點上的河底高程；⑧輸出計算結果，進行下一斷面的計算。

3 模型結果及分析

3.1 出庫流量模擬與分析

由于紫坪鋪水庫于2018年才開始建設映秀水文站，因此在此之前并無實測的入庫含沙量資料。本文選取2019年的水沙、地形資料對一維異重流水沙數學模型進行驗證。具體的參數設置如下：①地形數據，選擇大斷面數據（不包括加測斷面數據），斷面地形經過內插后，精度達到1 m；②大斷面的距壩里程，按照河道深泓線進行量測；③不可沖刷斷面，由于紫坪鋪水庫屬于逐漸淤積型水庫，因此不可沖刷地形設置為2008年的初始地形，認為無論如何沖刷，地形均比2008年初始地形要高；④進口水沙數據，以2019年的進出庫水沙為例，進口控制站使用映秀水文站的數據，出庫站使用紫坪鋪水文站的數據；⑤壩前水位，按照紫坪鋪水庫2019年壩前水位進行計算；⑥支流庫容，采用紫坪鋪水庫2019年實測支流庫容數據；⑦床沙級配，采用2020年實測河床級配；⑧懸沙粒徑，以取樣分析結果為主。

圖3為紫坪鋪水庫實測出庫流量與計算出庫流量對比，圖中黑色線條表示計算出庫流量過程，紅色線條表示實測出庫流量過程。由圖3可知，兩條曲線的吻合度較高。流量計算結果與實測數據非常吻合，無論是峰值還是過程線都達到了較高的精度，含沙量計算結果與紫坪鋪水文站實測結果也較為接近，這也反映出模型的參數選擇基本上可以作為紫坪鋪水庫水沙運動情景計算的參數。

圖3 紫坪鋪水庫實測與計算出庫流量變化過程對比

3.2 出庫含沙量模擬與分析

2019年，紫坪鋪水庫遭遇了“8·20”特大山洪泥石流災害，紫坪鋪入庫流量大幅上漲，由20日3時的884 m3／s猛漲到4時的3 457 m3／s，洪水過程形成迅速。水庫水位也持續上漲，從3時的845.19 m上漲到8時的848.13 m。為應對洪水，水庫管理部門開始通過調度方式對洪水進行攔蓄，將最大洪峰3 457 m3／s削減至620 m3／s出庫，削峰率為82%。至20日16時，紫坪鋪水庫共攔蓄洪水6 266萬m3，有效緩解了下游河道防洪壓力，保護了下游城市安全。在本次調度過程中，水庫管理部門于20日8時50分在水位848.75 m時開啟沖沙放空洞泄洪，14時水位迫近汛限水位850 m時開啟1＃泄洪排沙洞進行泄洪，并控制最大出庫流量不超過1 800 m3／s。上游暴雨山洪泥石流引發高濁度洪水進入水庫，導致紫坪鋪泄水濁度達到30 000 NTU，16時沖沙放空洞全關。

為驗證2019年紫坪鋪水庫汛期含沙量變化過程，本文利用一維異重流水沙模型計算出庫含沙量，并與出庫站紫坪鋪水文站含沙量變化過程進行對比。圖4為計算與實測出庫含沙量變化過程對比，整體上看，兩條過程線吻合程度較高。實測最大出庫含沙量為7.73 kg／m3，計算最大出庫含沙量為7.673 kg／m3，吻合度較高，其他時段的計算出庫含沙量略大于實測出庫含沙量，究其原因，主要是低含沙量的泥沙即使以異重流的形式到達壩前，若排沙設施不打開，出庫含沙量依舊很低，其出庫過程與水庫調度過程相關。

圖4 紫坪鋪水庫實測與計算出庫含沙量變化過程對比

3.3 庫區深泓線模擬與分析

深泓線即沿河流方向最大水深處的連線，沿此線的剖面為河流的縱剖面。深泓線的高程基本可以反映河床縱向變化情況。圖5為模型模擬2019年汛后與2008年汛后實測、2019年汛前實測以及2020年汛后實測深泓線對比。由圖5可知，2019年汛期結束后的深泓線基本在2019年汛前與2020年汛后之間，因此本文構建的一維異重流水沙數學模型的精度相對較高，既可以對流量過程進行精確模擬，也可以對輸沙過程進行精確模擬。

圖5 紫坪鋪水庫實測與計算深泓線變化過程對比

3.4 庫區淤積分布模擬與分析

利用一維異重流水沙數學模型可以模擬水庫不同位置的淤積量，圖6為紫坪鋪水庫2019年汛期不同位置單位距離實測與計算淤積量結果對比，黑色實線為模型計算結果，紅色實線為實測結果，兩條曲線的趨勢大致相同，淤積的位置相對準確。然而，紅色實線較黑色實線總體偏高，也就是說，實測的淤積量比模擬計算的淤積量大，究其原因，主要是模型計算結果僅為2019年汛期的淤積量，實測淤積量為2019—2020年兩年汛期的淤積量，由于2020年水沙數據未獲取，因此并沒有復核2020年紫坪鋪水庫的淤積形態。但是通過實測與計算淤積量的對比可知，該模型可以計算庫區不同位置的淤積量，且具有較好的精度。由圖6可知，2019年模型計算結果為汛期淤積2 274萬m3，而2019—2020年實測淤積量為2 477萬m3。2020年汛期水量比較大，但是來沙量相對偏少，可以推測，2020年全年的淤積量較2019年少很多，因此該模型計算的結果較為合理。

圖6 紫坪鋪水庫2019年汛期不同位置單位距離實測與計算淤積量結果對比

4 結 論

（1）基于水流連續性方程、水流動量方程、泥沙連續方程、河床變形方程以及異重流運動方程構建紫坪鋪一維水沙數學模型。依據實測資料對該模型進行了率定與參數的確定，最終得到了可以計算異重流演進過程的水沙一維數學模型。

（2）汛期的流量數據驗證結果表明：流量計算結果與實測數據較為吻合，含沙量計算結果與紫坪鋪水文站實測結果也較為接近，這反映出本模型的參數選擇基本上可以作為紫坪鋪水庫水沙運動情景計算的參數。

（3）汛期的含沙量數據驗證結果表明：實測最大出庫含沙量誤差很小，吻合度較高，其他時段的計算出庫含沙量略大于實測出庫含沙量。究其原因，主要是因為低含沙量的泥沙即使以異重流的形式到達壩前，若排沙設施不打開，出庫含沙量依舊很低，其出庫過程與水庫調度過程相關。

（4）汛期的地形數據驗證結果表明：2019年汛期結束后的深泓線基本在2019年汛前與2020年汛后之間，水沙數學模型的精度相對較高。2019年模型計算結果為汛期淤積2 274萬m3，而2019—2020年實測淤積量為2 477萬m3。