混凝土大壩潰壩動力響應規律數值研究

李旨洪

（江西省水利水電建設集團有限公司，江西 南昌 330025）

1 引言

混凝土壩以混凝土作為壩身材料，其結構簡單、施工方便，是一種整體性較好的剛性壩，中、高壩常用這種壩型[1，2]。此外，混凝土壩相對安全可靠，耐久性好，抵抗滲漏、洪水漫溢、地震和戰爭破壞能力都比較強，廣泛用于水利工程中。然而，由于地震等地質作用與材料強度之間的不匹配，混凝土大壩存在潛在的潰壩風險，會導致下游河谷洪水波快速傳動，造成嚴重的財產損失和人員傷亡[3，4]。

當前國內外通過數值模擬和室內試驗等手段對混凝土大壩漫頂潰壩進行了相關研究。胡文兵等[5]對該潰決方式進行水動力分析，得到混凝土壩傾倒過程的潰口流量計算公式，并結合概化典型流量過程線法得到了新的潰口流量計算方法；姚霄雯等[6]總結了混凝土壩的潰壩特點，統計了主要潰壩原因，進而提出了混凝土壩的主要潰壩模式及其潛在的潰壩路徑；王冰玲等[7]以某一混凝土壩為實例，利用開發的隨機網格生成算法建立了大壩模型，對爆炸載荷下混凝土壩的潰壩過程進行了模擬仿真；詹明強等[8]以某一面板堆石壩為工程實例模擬計算了大壩在不同潰決模式下的洪水傳播過程，得到了壩下河道沿程典型斷面的洪峰流量、最高洪水位及相應的出現時間。此外，還有部分學者通過統計分析、風險評價等方式對混凝土潰壩進行了相關研究[9-15]。

以位于江西省某簡化的大壩模型為例，采用計算流體力學CFD程序中的ANSYS FLUENT對其潰壩后的水動力演變過程進行模擬，并對比分析了LES和K-E湍流模型對計算結果的影響，研究成果可為相關工程提供參考。

2 工程概況

大壩位于江西省，是一座以發電為主，兼顧改善城區水環境、航運等綜合利用的水利水電工程。壩址以上控制流域面積16546 km2，總庫容1.52億m3，水庫正常蓄水位為103.8 m，電站裝機容量為32 MW，多年平均發電量13887×106kW·h，85%保證出力5300 kW，裝機年利用小時數3528 h。根據《防洪標準》（GB5021-2014）、《水利水電工程等級劃分及洪水標準》（SL252-2017）確定該工程屬大（2）型，工程等別Ⅱ等。主要建筑物級別為3級，次要建筑物級別為4級，臨時建筑物級別為5級，船閘級別為V級。擋水和泄水建筑物設計洪水標準為50 a一遇，校核洪水標準為300 a一遇；消能防沖洪水標準為50 a一遇。廠房非擋水部分設計洪水標準為20 a一遇，校核洪水標準為100 a一遇。大壩電站工程于2011年6月26日動工，主體工程于2011年10月12日正式開工，首臺機組（3#機組）于2013年11月12日完成72 h試運行，2014年2月21日全部3臺機組完成試運行。工程于2018年4月通過下閘蓄水驗收，2020年12月通過竣工驗收。

3 數值計算方程

本文通過計算流體力學CFD程序中的ANSYS FLUENT進行模擬。計算采用三維Navier-Stokes方程，并采用LES模型來模擬湍流。LES條件在空間和時間上都是離散的，在本研究中，PISO算法用于解決壓力場和速度場之間的轉換，使計算以更快的方式收斂。不可壓縮液體的質量和力的分離或雷諾平均守恒條件可以單獨由下式表達（總控制方程）：

式中：u?和uj由流體速度轉換而來（m/s）；xi和xj為笛卡爾坐標方向；ρ為流體密度（kg/m3）；-p為壓力平均值（Pa）；t為時間（s）；u為流體分子黏度（Pa·s）；τij為雷諾應力（Pa）。

對總控制方程進行建模后，模擬流體最重要的是對湍流方程進行建模分析。湍流是一種非對稱運動，當流體、液體或蒸汽流過堅固表面時，甚至當同一液體的相鄰流流過或超過彼此時，湍流都會出現在流體、液體或蒸汽中。目前進行湍流模擬的方法包括大渦模擬法（LES）、分離渦模擬法（DES）、比例-自適應模擬方法（SAS）、雷諾應力（7個方程式）、過渡SST方法、K-Omega湍流模型以及K-Epsilon湍流模型。其中，LES模型和K-Epsilon（以下簡稱K-E）模型是本文研究的重點。在大渦模擬方法中，先確定了較大的渦，并對較小或亞網格尺度的渦進行建模，從而比一般方法具有更好的穩定性。本研究中，采用了Smagorinsky提出的亞網格尺度模型：

式中：μt是子網格黏度（Pa·s）；Slj是更大尺度或定場的應變率，無量綱；δij為張量符號；其余變量含義同上。

渦流黏度建模如下：

式中：Cs為模型參數，取0.1；Δ為過濾器長度刻度（m）；ρ為流體密度（kg/m3）；為應變率。

K-E模型是計算流體力學（CFD）中最廣泛使用的模型，用于模擬湍流條件下的平均流質量。這是一個雙條件模型，通過2種傳輸條件對湍流進行了一般描述。最初傳輸的變量決定了湍流中的能量，稱為湍流動能（k）。第二個傳輸變量是湍流耗散（ε），決定了湍流動能的擴散速度。對于湍流動能k和湍流耗散ε分別采用下式描述：

式中：C1?和C2?為常量；σk和σ?分別為k方程和ε方程的湍流Prandtl數；ε為耗散率；k為湍動能；μt是子網格黏度（Pa·s）；xi和xj為笛卡爾坐標方向；t為時間（s）；ρ為流體密度（kg/m3）；Eij表示變形率分量，無量綱。

4 數值模型建立

本文流體力學計算的前期階段是創建液流區的計算幾何條件，Z軸與潰壩流體流動方向有關，X軸與大壩橫向有關，Y軸與平行于大壩高度的方向有關。本次建立的模型上游水庫長寬2 m、高1 m，下游水庫長寬高分別為8、2、0.3 m，其中混凝土密度設置為2400 kg/m3，彈性模量為14 GPa，泊松比為0.2，下游無蓄水，大壩尺寸詳見表1，潰壩數值模型如圖1所示。數值研究的第二步也是最關鍵的一步是建立與幾何發展相關的網格。Navier-Stokes方程是非線性偏微分方程，將整個流體區域視為一個連續統一體，利用笛卡爾坐標框架，根據區域離散化方法，求解了液流控制方程，即動量條件、連續條件。在有限元法中，通過在適當的空間中加入加權單元來獲得數值排列。該技術適用于非結構化網格，采用非結構化網格劃分方法分別將流體和壩體劃分為3200和2178個網格。完成大壩模型的幾何和網格劃分后，必須給出邊界條件。為了創建實時潰壩模擬，將頂面和下游邊界定義為壓力出口。此外，水閘不指定任何邊界條件，出水口的所有側面部分全約束。對于頂部自由曲面，大部分采用對稱邊界條件，并且與無滑移邊界條件一樣，頂面上既沒有對流通量也沒有擴散通量。本文計算的監測位置，如圖2所示。

表1 大壩建模尺寸

圖1 數值模型

圖2 監測點位置

5 數值結果分析

5.1 K-E模型和LES模型計算壓力分析

本次的數值結果包括流體速度、水面和底部壓力隨時間的變化。計算時一旦水從閘門排出，沿下游中軸線擴散，然后在下游端離開。之后，再利用LES模型和K-E模型進行湍流描述。流經閘門的水沿下游中軸線擴散，在水庫側匯合，在下游側分流，大壩下游的河道面積隨著時間的推移而增大。

LES模型和K-E模型計算的不同監測位置的壓力如圖3所示，如A、B、C、D、E、F、G，并給出其隨時間變化的規律。上游水庫位置（如位置A和位置B）的底部壓力最初為1 m水深的壓力，隨著時間增加，壓力會緩慢降低。在時間t=0時，閘門位置C和位置G附近沒有壓力，模擬開始后，底部壓力突然增加，并瞬時達到峰值。隨著時間的進一步增加，水在下游擴散并到達下游區域的終點，底部壓力降低。沿著下游的位置E和F的壓力變化，如圖4所示。由圖4可知，下游位置的底部壓力變化與C點和G點相似，在一段時間后達到峰值，此外閘門開啟和壓力增加之間具有一定時間差，閘門開啟位置附近的峰值壓力明顯高于下游位置。

圖3 LES模型和K-E模型計算的不同監測位置的壓力

圖4 下游位置的壓力變化

5.2 K-E模型和LES模型計算流速分析

監測點平均速度隨時間的變化規律如圖5所示，使用LES和K-E模型在上游位置（如A和B）進行速度計算，并相互比較。由圖5可知，上游側流速隨時間增加，在潰壩模擬后近2.5 s達到峰值，之后隨著時間的增加而進一步降低，并保持接近恒定值。然而，在閘門位置，速度在短時間內達到峰值，并緩慢下降。對于下游速度，閘門打開和速度增加之間存在時間滯后。在位置F處，速度幾乎為零，直到時間t=2 s，然后緩慢增加；時間t=5 s后，流速保持不變，因為大壩上游和下游的水位相同。此外，本文還計算了沿大壩長度的速度變化。

圖5 監測點平均速度隨時間的變化規律

計算結果表明，水波的速度在決口后立即沿長度增加，在涌浪后隨距離減少。在大壩閘門周圍，由于決口，流速很高，渠道的特性允許水流擴散到一個寬的橫截面，從而使流速沿水流方向降低。運行時觀察到的監測點流速在1.25和2.5 s左右大致接近，但隨著運行時間增加到3.5 s，產生的速度會低于之前。這主要是由于潰壩后上游儲水量減少以及下游橫截面增加。

6 結論

本文以位于江西省某簡化的大壩模型為例，采用計算流體力學CFD程序中的ANSYS FLUENT對其潰壩后的水動力演變過程進行模擬，并對比分析了LES和K-E湍流模型對計算結果的影響。結果表明，用LES和K-E湍流模型預測的速度剖面顯示出顯著差異，尤其是在大壩下游開口附近。然而，使用LES和K-E湍流模型預測的底部壓力在上游區域監測點表現出相似性，但在下游監測點出現了顯著的變化。最后，水波的速度在決口后立即沿長度增加，在涌浪后隨距離減少，研究成果可為相關工程提供參考。