利用Flow-3D軟件數值模擬堰塞壩潰決過程

彭 成，谷劉偉凱，鐘啟明

（1.南華大學，湖南 衡陽 421001；2.南京水利科學研究院，南京 210000；3.水利部土石壩破壞機理與防控技術重點實驗室，南京 210000）

堰塞壩是在一定地質和地貌條件下，由于地震、暴雨等因素引起山體發生滑坡、崩塌、泥石流而形成的阻塞河道、山谷的堆積體［1］。堰塞壩作為一種特殊的土石壩在全球分布較廣且易導致重大災害。 調查表明［2］，90%以上的堰塞壩是由地震或降水引起山體滑坡、崩塌和泥石流而形成的，約8%的堰塞壩是由火山噴發形成,其他引發因素如砍伐、河流下切和人工開挖等形成的堰塞壩約占2%，還有的則是在多種因素共同作用下形成的堰塞壩。 堰塞壩由于其特殊的幾何特征和材料組成， 在上游水流的不斷沖刷下易發生潰決。調查表明［1,3］，堰塞壩在形成后1 d 內潰決約占27%，50%的堰塞壩在形成10 d 內潰決，90%的堰塞壩在形成1 a 內潰決。從全球堰塞壩潰壩統計的資料可知，堰塞壩的潰決模式主要分管涌、漫頂和壩坡失穩［4］。研究資料表明［3］，堰塞壩潰決主要以漫頂潰決為主， 原因是堰塞壩缺少溢洪道和其他泄洪設施，隨著上游來水量增加，庫水位不斷上漲，水流將漫過壩頂開始溢流，水流不斷沖刷壩頂和壩坡，使堰塞壩高程降低、體積減小，隨著水流流量增加，加速堰塞壩的破壞，最終引起壩體突然潰決。

進行堰塞壩潰決過程模擬， 可為潰壩洪水風險分析和應急預案編制提供理論和技術支撐， 最大限度降低下游災害損失。目前，大量國內外學者針對堰塞壩的潰決機理展開深入研究［5-8］，研究重點主要包括壩高、壩寬、壩可蝕性對堰塞壩連續潰壩的影響［9］，泥沙沖蝕和潰口橫向發展對流量的影響［10］，記錄唐家山堰塞壩的泄流過程， 確定唐家山堰塞壩壩體材料的級配， 記錄潰決過程的峰值流量、 流速和潰口寬度， 驗證現有閾值方法計算泥沙侵蝕起始速度的可靠性［11］，堰塞壩沿河流運動方向的形態變化對漫頂潰壩的影響等［12］。堰塞壩的潰決過程模擬中存在潰壩流量的高度瞬態性質、結構空間變異性、潰口形態變化難以真實地描述等難題，涉及多學科交叉、多物理場耦合、多尺度模擬，目前對此領域研究的數值模型大多是二維模型，還不能真實地模擬潰壩過程。

唐家山堰塞壩現有研究沒有采用實際地形模擬潰壩過程，其采用的數值模型限制潰口變化，無法對潰口發展的復雜性進行更加真實的描述。 針對堰塞壩潰決過程模擬中的關鍵問題， 通過基于衛星遙感測繪的地形數據和堰塞壩的簡化模型建立唐家山堰塞壩的三維數值模型， 利用Flow-3D 軟件模擬唐家山堰塞壩漫頂潰壩潰決過程， 深入分析潰口處的流速特征及潰口沖蝕情況， 對2 個重要參數進行參數敏感性分析， 為今后堰塞壩的搶險方案和應急預案的制定提供技術參考。

1 唐家山堰塞湖概況

唐家山堰塞壩是2008 年“5·12”汶川地震造成對下游構成重大威脅的堰塞壩。 堰塞壩平面形態為沿河長約803 m，寬約611 m 的長條形，堰塞壩左右兩側高約90、124 m，堰塞壩體積約2037 萬m3。堰塞壩頂寬約300 m，地勢起伏較大。堰塞壩上游坡長約200 m， 坡比約1∶4； 下游壩坡長約300 m， 坡比約1∶2.4［13］。唐家山堰塞壩體由原山上的碎石土和寒武系硅質砂巖組成［7］。為盡快降低和解決唐家山堰塞壩的潰壩危險, 對唐家山堰塞壩進行開挖泄流渠導流排險。6 月7 日， 湖水開始溢出， 堰塞湖水位約740.3 m。6 月10 日6:00 左右，水位達742.10 m；6 月10 日12:30 出現峰值流量約6500 m3/s；6 月11 日庫水位由743.23 m 下降至715 m 以下， 相應庫容水量由2.466 億m3降至0.861 億m3以內［13］。

2 潰決過程模擬

唐家山堰塞壩根據衛星遙感測繪圖， 三維實體模型按實際地形1∶1 比例建立， 堰塞壩則基于堰塞壩的簡化模型和現場實測結果， 考慮壩體形態特征及筑壩材料的物理力學特性， 建立一個可模擬堰塞壩漫頂潰壩過程的三維模型并導入Flow-3D 軟件中，如圖1。為提升計算效率且避免影響研究區域大小， 唐家山堰塞壩三維計算模型計算區域為長約1100 m，寬約700 m，高約150 m，包含堰塞壩的計算范圍。對計算區域進行網格設置，為更好地觀察潰口處情況在潰口部分區域采用更精細的網格單元，細分處網格大小為2.5 m×2.5 m×2.5 m，其他網格設置為5 m×5 m×5 m，總有效網格為390 萬個。

圖1 唐家山衛星遙感圖像及堰塞壩三維計算模型

通過嚴格控制上游水位來消除堰塞壩的庫容、上游來水量對潰決流量過程的影響。 模型邊界設定如下，將上游端指定為壓力邊界，下游出口處指定為外流，側面指定為對稱，底面采用墻面邊界，頂部邊界設定為一個大氣壓的壓力邊界。

通過研究唐家山堰塞壩上選取的5 個鉆孔所得的級配曲線［14］（如圖2）， 三維模型的輸入參數［11］由此得出。堰塞壩區域設置為泥砂模型，泥砂顆粒干容重為2200 kg/m3，攜帶系數為0.018， 臨界希爾茲數為0.05，其余模型輸入參數設置如表1。

表1 泥砂模型輸入參數

圖2 唐家山堰塞壩級配

推移質起動時相對臨界牽引力稱為臨界希爾茲數，可由Soulsby-Whitehouse 方程［15］求解。

式中θcr為臨界希爾茲數，無量綱；d*為無量綱參數；sn=ρn/ρ，vf為流體黏度；τcr為泥砂的臨界剪應力；ρn為推移質的密度；ρ 為水流密度。

沉積的泥砂顆粒在水流的沖蝕和自重共同作用下在沉積狀態和懸浮狀態之間轉換的過程稱為泥砂挾帶。 軟件通過推移質的挾帶過程通過Mastbergen和Van den Berg 經驗公式［16］計算，挾帶速度公式為：

式中d*,i為無量綱顆粒當量粒徑；di為泥砂顆粒當量粒徑；θi為基于河床局部剪應力τ 的希爾茲數；αi為攜帶系數，ns為河床的外法向量。

懸移質控制方程為：

式中Cs,i為含沙流體的質量濃度；D 為擴散系數；us,i為含沙流體速度。

運用邁耶-皮特泥砂輸運模型［17］推算泥沙輸移率，推移質泥沙厚度通過范瑞恩公式［18］計算，河床輸沙率方程為：

式中βi為推移質系數。

3 計算結果分析

3.1 模擬結果與實測結果對比分析

在堰塞壩潰口處設立流量監測斷面， 具體情況如圖1（b）。計算結果顯示，6 月10 日11：50 達到峰值流量6937 m3/s，現場實測數據顯示潰口在6 月10日12：30 達到峰值流量，Flow-3D 模擬結果與實測峰值流量差異較小，相對誤差小于15%，到達峰值流量時間較實際時間提前40 min， 峰值流量模擬值較實測值偏大7.05% ， 潰口最終寬度模擬值較實測值偏大4.16%， 表明通過現場實測數據和模擬的峰值流量和潰口最終寬度的比較驗證模擬的有效性和準確性。實測數據與模擬數據對比如表2 和圖3，圖4。

表2 潰決峰值流量比較

圖3 潰口流量模擬值與實測值比較

圖4 潰口寬度模擬值與實測值比較

3.2 地面高程變化分析

唐家山堰塞壩漫頂潰決過程復雜，如圖5。按地面高程變化率和過流量大小分3 個時間段， 分別是潰口快速發展階段、洪峰階段和潰口發展減緩階段。地面高程變化率是指泄流渠底部高程隨時間的變化［19］，高程增加為正，高程下降為負，反應泄流渠受水流沖蝕作用的強弱程度。從潰決初期到洪峰階段，這段時間水流沖蝕強度高，河床泥砂沖蝕強度大，河床高程快速變化，潰口橫向發展迅速。洪峰結束后，隨著庫水位和潰口流量的降低，泥沙沖蝕程度降低，潰口形態變化速度減緩，地形變化隨之減緩。

圖5 唐家山堰塞壩潰決過程中地形變化模擬結果

3.3 流速變化規律

漫壩水流不同階段的流速模擬結果如圖6，開始溢流后流速一直處于上升趨勢，流速最大是泄流渠中陡坎的附近區域， 流速峰值約18 m/s。 由于水流流速較大，沖蝕程度較強， 導致陡坎高度不斷增加，陡坎坡度不斷變陡，陡坎底部被劇烈沖蝕， 水流變換劇烈，出現旋流，泥砂易被水流卷起帶向下游。洪峰結束后，河床下切和坡降都明顯變小，同時，隨著潰決過程的不斷發展，上游水位下降迅速， 流速趨于平緩，在該階段流速最高達16 m/s 左右，河道中陡坎處和末段出現流速最大值。

圖6 唐家山堰塞壩潰決過程水流流速變化計算結果

由圖6 可知， 流速變化規律在不同階段有一定共同點： 漫壩水流流入泄流渠前， 水流速度趨近于0，庫水壓處于穩定狀態；湖水進入泄流渠后，受到重力作用， 入口處的水流流速明顯增加。 從圖6 可看出， 隨著潰口起始位置泥砂顆粒不斷被水流帶向下游，在泄流渠的中部逐漸發展成兩個陡坎。隨著陡坎高程增加，水頭落差變大，陡坎的坡度變陡，水流對底部的沖蝕越來越強烈， 帶動陡坎面上的泥砂劇烈沖蝕，使得陡坎面不斷向上游推進。在泄流渠尾段，下游壩坡上的水流在重力作用下流速迅速增加，水流沖蝕程度增強，水流卷起的泥砂顆粒較多，形成高含沙水流對泄流渠不斷沖蝕下切， 直到潰口發展減緩時期，潰口橫向發展顯著減緩。如圖7，泄流渠在水流不斷沖蝕下，河床中出現明顯溯源現象，流速較快的區域位于泄流渠的中段和末段。

圖7 陡坎溯源沖蝕示意圖

3.4 參數敏感性分析

潰壩峰值流量取決于起潰水位、 潰口形狀和發展過程等因素。 為進一步研究模型參數對唐家山堰塞壩潰壩過程的影響， 特定選取起潰水位進行敏感性分析。此外，為驗證泄流渠對排險的作用，對沒有開挖泄流渠情況下的工況進行模擬， 研究泄流渠對潰壩過程的影響。

3.4.1 起潰水位敏感性分析

為探究起潰水位的敏感性， 分別模擬起潰水位為745 m 和748 m 兩種情況。 起潰水位敏感性分析計算結果如表3， 其中相對誤差為不同起潰水位的計算結果和初始起潰水位 （水位為742 m 時） 的比較；如圖8，不同起潰水位情況下的潰口流量過程。

表3 不同起潰水位參數敏感性分析

圖8 不同起潰水位工況下對潰口流量影響分析

從模擬結果發現， 堰塞壩潰決模擬計算得出的各項參數均對起潰水位較為敏感， 其中潰口底寬最為敏感。 如圖8 和表3， 當起潰水位由742 m 增至748 m 時， 相應潰壩峰值流量由6937 m3/s 增加至8542 m3/s，增幅23.1%，到達潰口峰值流量的時間也由5.83 h 提前到3.55 h，提前39.1%，其下泄水量也增加明顯。 說明起潰水位的精準測量對堰塞壩潰決過程模擬至關重要。

3.4.2 泄流渠開挖敏感性分析

對于大多數堰塞湖，受制于地理環境因素，常無法采取有效的人工干預措施。研究發現，在現實情況許可下，通過開挖泄流渠可使庫水位降低，從而降低潰口峰值流量和減小下泄水量， 以降低堰塞壩潰決造成的損失。為探究是否開挖泄流渠對潰決的影響，通過無泄流渠的三維模型模擬其潰決過程， 模型的其余參數保持不變。 是否開挖泄流渠的敏感性分析結果如表4。如圖9，依據湖水流入情況，在沒有泄流渠的工況下，唐家山堰塞壩在約為60 h 后發生漫頂潰壩， 在66.3 h 后潰口流量達到峰值且下泄水量遠大于開挖泄流渠時的下泄水量。 表4 和圖9 的泄流渠開挖敏感性分析結果表明， 對于庫容較大的堰塞湖，開挖泄流渠可有效降低潰口峰值流量，減小潰口寬度，從而減少對下游造成的損失。在現實情況許可下，可作為搶險計劃的首要措施。

表4 開挖泄流渠參數敏感性分析

圖9 開挖泄流渠對潰口流量影響分析

4 結語

（1）基于堰塞壩漫頂潰決機理，通過衛星遙感測繪獲取的高程數據、 現場實測數據和堰塞壩模型建立唐家山堰塞壩的三維模型， 并采用Flow－3D 軟件對唐家山堰塞壩漫頂潰決過程進行模擬分析。 通過模擬結果與實測數據對比，驗證三維模型的合理性。為探究模型參數對潰壩過程的影響， 開展重要參數的敏感性分析。

（2）由模擬結果可得出，唐家山堰塞壩漫頂沖刷是一個非常復雜的過程，可分為潰口快速發展階段、洪峰階段、 潰口發展減緩階段3 個階段。 潰決過程中， 泄流渠中段出現陡坎的區域和下游壩坡的水流流速大，泄流渠的兩側沖蝕深度大。唐家山堰塞壩潰決過程中有明顯的溯源沖刷現象。 從參數敏感性分析可得：起潰水位對于堰塞壩潰決過程影響顯著，對該參數的測量予以充分重視。此外，對具有較大庫容的堰塞壩，在現實情況許可下，開挖泄流渠是一種能夠降低堰塞湖潰決洪水帶來災害的有效措施， 可作為堰塞湖應急處置措施的首選方案。