混凝土大壩潰壩動(dòng)力響應(yīng)規(guī)律數(shù)值研究

李旨洪

（江西省水利水電建設(shè)集團(tuán)有限公司，江西 南昌 330025）

1 引言

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

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

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

2 工程概況

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

3 數(shù)值計(jì)算方程

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

式中：u?和uj由流體速度轉(zhuǎn)換而來(lái)（m/s）；xi和xj為笛卡爾坐標(biāo)方向；ρ為流體密度（kg/m3）；-p為壓力平均值（Pa）；t為時(shí)間（s）；u為流體分子黏度（Pa·s）；τij為雷諾應(yīng)力（Pa）。

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

式中：μt是子網(wǎng)格黏度（Pa·s）；Slj是更大尺度或定場(chǎng)的應(yīng)變率，無(wú)量綱；δij為張量符號(hào)；其余變量含義同上。

渦流黏度建模如下：

式中：Cs為模型參數(shù)，取0.1；Δ為過(guò)濾器長(zhǎng)度刻度（m）；ρ為流體密度（kg/m3）；為應(yīng)變率。

K-E模型是計(jì)算流體力學(xué)（CFD）中最廣泛使用的模型，用于模擬湍流條件下的平均流質(zhì)量。這是一個(gè)雙條件模型，通過(guò)2種傳輸條件對(duì)湍流進(jìn)行了一般描述。最初傳輸?shù)淖兞繘Q定了湍流中的能量，稱(chēng)為湍流動(dòng)能（k）。第二個(gè)傳輸變量是湍流耗散（ε），決定了湍流動(dòng)能的擴(kuò)散速度。對(duì)于湍流動(dòng)能k和湍流耗散ε分別采用下式描述：

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

4 數(shù)值模型建立

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

表1 大壩建模尺寸

圖1 數(shù)值模型

圖2 監(jiān)測(cè)點(diǎn)位置

5 數(shù)值結(jié)果分析

5.1 K-E模型和LES模型計(jì)算壓力分析

本次的數(shù)值結(jié)果包括流體速度、水面和底部壓力隨時(shí)間的變化。計(jì)算時(shí)一旦水從閘門(mén)排出，沿下游中軸線(xiàn)擴(kuò)散，然后在下游端離開(kāi)。之后，再利用LES模型和K-E模型進(jìn)行湍流描述。流經(jīng)閘門(mén)的水沿下游中軸線(xiàn)擴(kuò)散，在水庫(kù)側(cè)匯合，在下游側(cè)分流，大壩下游的河道面積隨著時(shí)間的推移而增大。

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

圖3 LES模型和K-E模型計(jì)算的不同監(jiān)測(cè)位置的壓力

圖4 下游位置的壓力變化

5.2 K-E模型和LES模型計(jì)算流速分析

監(jiān)測(cè)點(diǎn)平均速度隨時(shí)間的變化規(guī)律如圖5所示，使用LES和K-E模型在上游位置（如A和B）進(jìn)行速度計(jì)算，并相互比較。由圖5可知，上游側(cè)流速隨時(shí)間增加，在潰壩模擬后近2.5 s達(dá)到峰值，之后隨著時(shí)間的增加而進(jìn)一步降低，并保持接近恒定值。然而，在閘門(mén)位置，速度在短時(shí)間內(nèi)達(dá)到峰值，并緩慢下降。對(duì)于下游速度，閘門(mén)打開(kāi)和速度增加之間存在時(shí)間滯后。在位置F處，速度幾乎為零，直到時(shí)間t=2 s，然后緩慢增加；時(shí)間t=5 s后，流速保持不變，因?yàn)榇髩紊嫌魏拖掠蔚乃幌嗤?。此外，本文還計(jì)算了沿大壩長(zhǎng)度的速度變化。

圖5 監(jiān)測(cè)點(diǎn)平均速度隨時(shí)間的變化規(guī)律

計(jì)算結(jié)果表明，水波的速度在決口后立即沿長(zhǎng)度增加，在涌浪后隨距離減少。在大壩閘門(mén)周?chē)捎跊Q口，流速很高，渠道的特性允許水流擴(kuò)散到一個(gè)寬的橫截面，從而使流速沿水流方向降低。運(yùn)行時(shí)觀察到的監(jiān)測(cè)點(diǎn)流速在1.25和2.5 s左右大致接近，但隨著運(yùn)行時(shí)間增加到3.5 s，產(chǎn)生的速度會(huì)低于之前。這主要是由于潰壩后上游儲(chǔ)水量減少以及下游橫截面增加。

6 結(jié)論

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