郭 磊，張金鳳，郭利霞，申偉平，郭宇航

（1.華北水利水電大學(xué)水利學(xué)院，鄭州 450046；2.河南水谷研究院，鄭州 450046；3.河南省水環(huán)境模擬與治理重點實驗室，鄭州 450002）

0 引 言

極端暴雨使河道水量增多，河水流速變大，會在短時間內(nèi)形成較大的洪峰流量，使庫水位急劇上升。目前，國內(nèi)外關(guān)于水庫水位變化對混凝土壩的影響分析已有很多，BAYRAK T［1，2］采用一種新的動態(tài)變形分析方法，研究了亞穆拉大壩的位移和庫水位之間的關(guān)系，結(jié)果表明庫水位變化對壩體變形影響顯著。PYTHAROULI S I 等［3］利用譜分析技術(shù)對拉東大壩30 余年的大地檢測記錄進行分析，發(fā)現(xiàn)壩頂中部的垂直位移和水平位移最大，水位荷載與壩體變形為非線性關(guān)系。SUN K M 等［4］采用強度折減方法對混凝土壩進行了穩(wěn)定分析，結(jié)果表明隨著水位荷載和降雨強度的增大，壩體內(nèi)應(yīng)力應(yīng)變變化也越大，越容易發(fā)生失穩(wěn)破壞。GUO L等［5］通過數(shù)值仿真手段研究了混凝土壩對極端暴雨洪水的響應(yīng)，研究表明，水溫變化對壩體應(yīng)力影響較小，而水位變化對壩體應(yīng)力影響顯著。以上研究均是對靜動荷載下的混凝土大壩受力特性進行分析。本文以膠凝砂礫石壩為研究對象，對極端暴雨作用下水庫水位的影響進行分析，并在有限元分析中引入了膠凝砂礫石材料的彈塑性損傷模型，分析了暴雨歷時過程中壩體的應(yīng)力、位移和損傷破壞規(guī)律，可為應(yīng)對氣候變化的結(jié)構(gòu)設(shè)計提供理論依據(jù)。

1 計算理論

1.1 流量-水位關(guān)系

降雨徑流是一個復(fù)雜的過程，降雨通過產(chǎn)流計算扣除損失量之后剩下凈雨，再通過匯流作用經(jīng)河網(wǎng)匯集到流域河道形成斷面流量（流域凈流量），水庫入庫流量通常是通過測量庫前河道斷面流速經(jīng)計算得到的［6-8］。通過河道的全斷面流量Q為：

式中：A為河道斷面面積，m2；dA為A內(nèi)的單元面積（其寬為db，高為dh），m2；v(h，b)為垂直于dA的流速，m/s；B為水面寬度，m；h為水深，m；hb為水邊到水面寬為b處的水深，m。

由上可知，水位-流量存在非線性關(guān)系，但受河道的斷面形態(tài)、彎曲狀態(tài)、底坡以及河床表面糙率等因素的影響，很難精確獲取其數(shù)學(xué)模型［9，10］。可利用有限元流體軟件模擬水流運動過程［11，12］，其一般微分方程為：

連續(xù)方程：

動量方程：

式中：ρ為水的密度；u為沿水流方向的水流流速；w為垂直方向上的水流流速；μ為動力黏度系數(shù)；p為壓強；Su、Sw為廣義源項；為雷諾應(yīng)力項，通用形式為：

式中：μt為湍動系數(shù)；k為湍動能；ui為時均流速；δij克羅內(nèi)克符號（當(dāng)i=j時，δij= 1；當(dāng)i≠j，δij= 0）。

1.2 考慮損傷的有限元結(jié)構(gòu)分析

對于空間結(jié)構(gòu)的有限元分析計算，是將建立的實體模型劃分為有限個具有塊狀結(jié)構(gòu)離散的實體單元，相鄰單元之間通過節(jié)點以鉸接的方式連接，利用相互作用的連續(xù)單元運用數(shù)學(xué)方法去模擬物體真實的物理系統(tǒng)，實現(xiàn)利用有限個未知量求得無限個未知量的物理系統(tǒng)的近似解。在空間問題中，彈性力學(xué)物理方程可表示為：

式中：εx、εy、εz分別為X、Y、Z方向上的垂直應(yīng)變；σx、σy、σz分別為X、Y、Z方向上的垂直應(yīng)力；γxy、γyz、γzx分別為X-Y、Y-Z、Z-X平面上的剪切應(yīng)變；τxy、τyz、τzx分別為X-Y、Y-Z、Z-X平面上的剪切應(yīng)力；E為彈性模量，可將E換成考慮損傷后的彈性模量=(1 -d)E，其中d為損傷參數(shù)，可根據(jù)試驗得到。

2 計算模型及參數(shù)

2.1 有限元模型

選取土耳其的Oyuk 壩為研究對象，該壩位于在峽谷之間，兩岸壩肩陡峭，壩基以片麻巖為主。壩體高度為100 m，壩頂寬度為7.5 m，上、下游坡比為1∶0.7。在建立有限元模型時，為了提高仿真計算的可靠度，需適當(dāng)擴大計算區(qū)域，壩踵處的基巖向上游延伸100 m，壩趾處的基巖向下游延伸100 m，壩底基巖豎直向下延伸100 m。采用SOLID185 單元劃分網(wǎng)格，模型共有56 896 個節(jié)點，計算單元共51 060 個，有限元模型如圖1 所示。應(yīng)力計算時在地基底部施加全約束，在地基四周施加法向約束。

圖1 壩體有限元模型Fig.1 Finite element model of dam body

2.2 計算參數(shù)

依據(jù)工程資料對壩體膠凝砂礫石材料進行了室內(nèi)試驗，壩基為巖體，參數(shù)選取參考文獻［13］，膠凝砂礫石壩仿真模型計算參數(shù)見表1所示。

表1 計算參數(shù)Tab.1 Calculation parameters

膠凝砂礫石是一種彈塑性材料［14］，選取《混凝土結(jié)構(gòu)設(shè)計規(guī)范》［15］的分段曲線模型作為膠凝砂礫石材料本構(gòu)模型，通過單軸抗拉和單軸抗壓試驗得到其損傷演化公式為：

受壓損傷公式為：

式中：dc為膠凝砂礫石單軸受壓損傷演化參數(shù)；x為膠凝砂礫石壓應(yīng)變與其單軸抗壓強度代表值相應(yīng)的峰值壓應(yīng)變的比值。

受拉損傷公式為：

式中：dt為膠凝砂礫石單軸受拉損傷演化參數(shù)；x為膠凝砂礫石拉應(yīng)變與其單軸抗拉強度代表值相應(yīng)的峰值拉應(yīng)變的比值。

2.3 水位邊界條件的獲取

由于極端強降雨天氣的發(fā)生，導(dǎo)致水庫水位上升較快，水庫蓄滿后泄水建筑物無法宣泄多余洪水，使得水位持續(xù)升高出現(xiàn)漫頂甚至潰壩事故。本文選取當(dāng)?shù)啬炒蔚湫蜆O端暴雨洪水過程為研究依據(jù)，因缺乏該庫區(qū)實測資料，故假設(shè)庫區(qū)相對規(guī)則，入庫洪水過程線如圖2所示。

圖2 入庫凈流速過程線Fig.2 Process line of net inflow velocity

為了獲取極端暴雨最不利工況下的水位變化特征，利用Fluent 軟件采用VOF（Volume of Fluid，VOF）法和標(biāo)準(zhǔn)k-ε 模型對壩前水流進行模擬。模型使用非結(jié)構(gòu)化四邊形網(wǎng)格劃分，共44 583 個網(wǎng)格。流體域從壩體分別向上游、下游延伸150 m，豎直向上延伸50 m。邊界條件設(shè)置：模型底部為默認(rèn)的固體壁面邊界，上部設(shè)置為壓力型入口；模型上游分別設(shè)置速度型入口和壓力型入口，壩前庫水位為97 m，斷面平均流速設(shè)置為0.68 m/s，氣壓為一個工程大氣壓，下游設(shè)置為壓力型出口。計算模型如圖3 所示，對流場進行初始化，初始速度均為0 m/s，計算時長為100 s。

圖3 流體域網(wǎng)格剖分示意圖Fig.3 Schematic Diagram of Grid Subdivision in Fluid Domain

根據(jù)仿真計算結(jié)果繪制出時間-水位變化曲線，隨著匯流作用水位呈現(xiàn)先急速后緩慢上升趨勢，最高水位達到105 m，高出壩體5 m。

由圖4擬合出暴雨工況下時間-水位關(guān)系式為：

圖4 暴雨工況下庫水位變化曲線圖Fig.4 Changes of water level under heavy rain conditions

3 計算結(jié)果分析

應(yīng)力計算采用超載法，將漫頂水位荷載完全等效為作用在壩體上的水荷載。極端降雨工況下壩體應(yīng)力計算結(jié)果如圖5～10所示。

圖5 膠凝砂礫石壩不同歷時應(yīng)力云圖（t=0 s，單位：Pa）Fig.5 Stress nephogram of CSG dam with different durations（t=0 s）

圖6 膠凝砂礫石壩不同歷時應(yīng)力云圖（t=20 s，單位：Pa）Fig.6 Stress nephogram of CSG dam with different durations（t=20 s）

圖7 膠凝砂礫石壩不同歷時應(yīng)力云圖（t=40 s，單位：Pa）Fig.7 Stress nephogram of CSG dam with different durations（t=40 s）

圖8 膠凝砂礫石壩不同歷時應(yīng)力云圖（t=60 s，單位：Pa）Fig.8 Stress nephogram of CSG dam with different durations（t=60 s）

圖9 膠凝砂礫石壩不同歷時應(yīng)力云圖（t=80 s，單位：Pa）Fig.9 Stress nephogram of CSG dam with different durations（t=80 s）

圖10 膠凝砂礫石壩不同歷時應(yīng)力云圖（t=100 s，單位：Pa）Fig.10 Stress nephogram of CSG dam with different durations（t=100 s）

（1）應(yīng)力響應(yīng)分析。由以上計算結(jié)果可知，在上游水位荷載的作用下，壩踵處所受拉應(yīng)力最大，壩趾所受壓應(yīng)力最大。由于壩體與壩基材料不同，因變形不協(xié)調(diào)致使二者接觸面較為薄弱，存在較大的不利應(yīng)力。隨著水位的上升，壩趾處壓應(yīng)力有所增大，而壩踵處拉應(yīng)力明顯降低了很多，這是因為壩踵區(qū)域處的材料發(fā)生了損傷，并且隨著損傷的累計，現(xiàn)有破壞部位便會在應(yīng)力集中的作用下，沿著順?biāo)鞣较蚝蛪误w縱深方向進行進一步的擴展，最后形成較大的宏觀裂縫致使建筑物失效。故在極端暴雨工況下，膠凝砂礫石壩可能會在壩踵處發(fā)生損傷破壞，且為拉裂破壞；壩體與壩基接觸面處存在較大不利應(yīng)力，但影響不大。

（2）位移響應(yīng)分析。膠凝砂礫石壩順?biāo)鞣较蛏献畲笞冃蚊娴奈灰屏侩S暴雨歷時變化如圖11所示。

圖11 暴雨工況下膠凝砂礫石壩順?biāo)鞣较蛭灰茍DFig.11 Displacement diagram of downstream flow direction on top of CSG dam under rainstorm condition

由圖11 可知，壩體位移量隨著上游水位的上升而逐漸增大，其增長速率由大變小，增長趨勢與壩體上游水位荷載增長趨勢相同，為正相關(guān)性。從大壩整體來看，壩體中部位移量較大，壩頂次之，壩底位移量最小。

（3）損傷單元分析。為了更好表述壩體的受力狀況和損傷演化規(guī)律，提取壩體若干損傷特征單元進行分析對比，特征單元的具體位置詳見圖1，其損傷變化如圖12所示。

圖12 暴雨工況下膠凝砂礫石壩特征點的彈性模量損傷圖Fig.12 Elastic modulus damage diagram of characteristic points of CSG dam under rainstorm conditions

由圖12 分析可知，隨著上游水位不斷的升高，作用在上游壩坡的荷載逐漸增大，內(nèi)部的39 961 單元無損傷，其他特征單元均產(chǎn)生損傷，且內(nèi)部39 751 單元的損傷值小于表層39 301 單元，可見壩體開始由無損狀態(tài)向損傷狀態(tài)過渡，損傷范圍逐漸擴大。由于損傷不斷累積，特征單元的彈性模量有減小的趨勢，其趨勢由大逐漸變小直至彈性模量為零，損傷量達到極值，即特征點處結(jié)構(gòu)完全失效，壩體材料被破壞。單元損傷值隨時間的變化同時也驗證了壩體應(yīng)力的變化?？傊S著水位的升高壩踵部位的損傷量會優(yōu)先達到極值，壩踵區(qū)域處的淺層地基會優(yōu)先發(fā)生塑性屈服，塑性屈服區(qū)會分別沿著豎直方向向下擴展和建基面方向向下游擴展。

4 結(jié) 論

以典型工程為例，采用有限元法對暴雨洪水入庫時的水位變化過程進行模擬，得到漫頂時期的穩(wěn)定水位，即暴雨時期的最不利工況水位，以此水位為邊界條件對壩體進行應(yīng)力仿真分析，得出以下結(jié)論。

（1）在水位荷載的作用下，在壩踵處壩體所受拉應(yīng)力最大，在壩趾處壩體所受壓應(yīng)力最大；壩體與地基接觸面處承受較大的不利應(yīng)力，此部位為壩體薄弱區(qū)域。

（2）隨著水位持續(xù)升高，壩體位移響應(yīng)值不斷增大，且壩中部位的響應(yīng)值最大，壩頂次之，壩底最小，說明壩中部位變形量最大。

（3）對于膠凝砂礫石壩體而言，受極端強降雨的作用，損傷部位主要集中在壩踵區(qū)域處，沿著順?biāo)鞣较蚝蛪误w縱深方向進行進一步的擴展，壩體與地基接觸面處存在較大不利應(yīng)力，但影響不大，說明外荷載越大，損傷范圍越大，壩體越容易失穩(wěn)破壞。