郭 磊,張金鳳,郭利霞,申偉平,郭宇航
(1.華北水利水電大學水利學院,鄭州 450046;2.河南水谷研究院,鄭州 450046;3.河南省水環境模擬與治理重點實驗室,鄭州 450002)
極端暴雨使河道水量增多,河水流速變大,會在短時間內形成較大的洪峰流量,使庫水位急劇上升。目前,國內外關于水庫水位變化對混凝土壩的影響分析已有很多,BAYRAK T[1,2]采用一種新的動態變形分析方法,研究了亞穆拉大壩的位移和庫水位之間的關系,結果表明庫水位變化對壩體變形影響顯著。PYTHAROULI S I 等[3]利用譜分析技術對拉東大壩30 余年的大地檢測記錄進行分析,發現壩頂中部的垂直位移和水平位移最大,水位荷載與壩體變形為非線性關系。SUN K M 等[4]采用強度折減方法對混凝土壩進行了穩定分析,結果表明隨著水位荷載和降雨強度的增大,壩體內應力應變變化也越大,越容易發生失穩破壞。GUO L等[5]通過數值仿真手段研究了混凝土壩對極端暴雨洪水的響應,研究表明,水溫變化對壩體應力影響較小,而水位變化對壩體應力影響顯著。以上研究均是對靜動荷載下的混凝土大壩受力特性進行分析。本文以膠凝砂礫石壩為研究對象,對極端暴雨作用下水庫水位的影響進行分析,并在有限元分析中引入了膠凝砂礫石材料的彈塑性損傷模型,分析了暴雨歷時過程中壩體的應力、位移和損傷破壞規律,可為應對氣候變化的結構設計提供理論依據。
降雨徑流是一個復雜的過程,降雨通過產流計算扣除損失量之后剩下凈雨,再通過匯流作用經河網匯集到流域河道形成斷面流量(流域凈流量),水庫入庫流量通常是通過測量庫前河道斷面流速經計算得到的[6-8]。通過河道的全斷面流量Q為:

式中:A為河道斷面面積,m2;dA為A內的單元面積(其寬為db,高為dh),m2;v(h,b)為垂直于dA的流速,m/s;B為水面寬度,m;h為水深,m;hb為水邊到水面寬為b處的水深,m。
由上可知,水位-流量存在非線性關系,但受河道的斷面形態、彎曲狀態、底坡以及河床表面糙率等因素的影響,很難精確獲取其數學模型[9,10]。可利用有限元流體軟件模擬水流運動過程[11,12],其一般微分方程為:
連續方程:

動量方程:


式中:ρ為水的密度;u為沿水流方向的水流流速;w為垂直方向上的水流流速;μ為動力黏度系數;p為壓強;Su、Sw為廣義源項;為雷諾應力項,通用形式為:

式中:μt為湍動系數;k為湍動能;ui為時均流速;δij克羅內克符號(當i=j時,δij= 1;當i≠j,δij= 0)。
對于空間結構的有限元分析計算,是將建立的實體模型劃分為有限個具有塊狀結構離散的實體單元,相鄰單元之間通過節點以鉸接的方式連接,利用相互作用的連續單元運用數學方法去模擬物體真實的物理系統,實現利用有限個未知量求得無限個未知量的物理系統的近似解。在空間問題中,彈性力學物理方程可表示為:

式中:εx、εy、εz分別為X、Y、Z方向上的垂直應變;σx、σy、σz分別為X、Y、Z方向上的垂直應力;γxy、γyz、γzx分別為X-Y、Y-Z、Z-X平面上的剪切應變;τxy、τyz、τzx分別為X-Y、Y-Z、Z-X平面上的剪切應力;E為彈性模量,可將E換成考慮損傷后的彈性模量=(1 -d)E,其中d為損傷參數,可根據試驗得到。
選取土耳其的Oyuk 壩為研究對象,該壩位于在峽谷之間,兩岸壩肩陡峭,壩基以片麻巖為主。壩體高度為100 m,壩頂寬度為7.5 m,上、下游坡比為1∶0.7。在建立有限元模型時,為了提高仿真計算的可靠度,需適當擴大計算區域,壩踵處的基巖向上游延伸100 m,壩趾處的基巖向下游延伸100 m,壩底基巖豎直向下延伸100 m。采用SOLID185 單元劃分網格,模型共有56 896 個節點,計算單元共51 060 個,有限元模型如圖1 所示。應力計算時在地基底部施加全約束,在地基四周施加法向約束。

圖1 壩體有限元模型Fig.1 Finite element model of dam body
依據工程資料對壩體膠凝砂礫石材料進行了室內試驗,壩基為巖體,參數選取參考文獻[13],膠凝砂礫石壩仿真模型計算參數見表1所示。

表1 計算參數Tab.1 Calculation parameters
膠凝砂礫石是一種彈塑性材料[14],選取《混凝土結構設計規范》[15]的分段曲線模型作為膠凝砂礫石材料本構模型,通過單軸抗拉和單軸抗壓試驗得到其損傷演化公式為:
受壓損傷公式為:

式中:dc為膠凝砂礫石單軸受壓損傷演化參數;x為膠凝砂礫石壓應變與其單軸抗壓強度代表值相應的峰值壓應變的比值。
受拉損傷公式為:

式中:dt為膠凝砂礫石單軸受拉損傷演化參數;x為膠凝砂礫石拉應變與其單軸抗拉強度代表值相應的峰值拉應變的比值。
由于極端強降雨天氣的發生,導致水庫水位上升較快,水庫蓄滿后泄水建筑物無法宣泄多余洪水,使得水位持續升高出現漫頂甚至潰壩事故。本文選取當地某次典型極端暴雨洪水過程為研究依據,因缺乏該庫區實測資料,故假設庫區相對規則,入庫洪水過程線如圖2所示。

圖2 入庫凈流速過程線Fig.2 Process line of net inflow velocity
為了獲取極端暴雨最不利工況下的水位變化特征,利用Fluent 軟件采用VOF(Volume of Fluid,VOF)法和標準k-ε 模型對壩前水流進行模擬。模型使用非結構化四邊形網格劃分,共44 583 個網格。流體域從壩體分別向上游、下游延伸150 m,豎直向上延伸50 m。邊界條件設置:模型底部為默認的固體壁面邊界,上部設置為壓力型入口;模型上游分別設置速度型入口和壓力型入口,壩前庫水位為97 m,斷面平均流速設置為0.68 m/s,氣壓為一個工程大氣壓,下游設置為壓力型出口。計算模型如圖3 所示,對流場進行初始化,初始速度均為0 m/s,計算時長為100 s。

圖3 流體域網格剖分示意圖Fig.3 Schematic Diagram of Grid Subdivision in Fluid Domain
根據仿真計算結果繪制出時間-水位變化曲線,隨著匯流作用水位呈現先急速后緩慢上升趨勢,最高水位達到105 m,高出壩體5 m。
由圖4擬合出暴雨工況下時間-水位關系式為:

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

應力計算采用超載法,將漫頂水位荷載完全等效為作用在壩體上的水荷載。極端降雨工況下壩體應力計算結果如圖5~10所示。

圖5 膠凝砂礫石壩不同歷時應力云圖(t=0 s,單位:Pa)Fig.5 Stress nephogram of CSG dam with different durations(t=0 s)

圖6 膠凝砂礫石壩不同歷時應力云圖(t=20 s,單位:Pa)Fig.6 Stress nephogram of CSG dam with different durations(t=20 s)

圖7 膠凝砂礫石壩不同歷時應力云圖(t=40 s,單位:Pa)Fig.7 Stress nephogram of CSG dam with different durations(t=40 s)

圖8 膠凝砂礫石壩不同歷時應力云圖(t=60 s,單位:Pa)Fig.8 Stress nephogram of CSG dam with different durations(t=60 s)

圖9 膠凝砂礫石壩不同歷時應力云圖(t=80 s,單位:Pa)Fig.9 Stress nephogram of CSG dam with different durations(t=80 s)

圖10 膠凝砂礫石壩不同歷時應力云圖(t=100 s,單位:Pa)Fig.10 Stress nephogram of CSG dam with different durations(t=100 s)
(1)應力響應分析。由以上計算結果可知,在上游水位荷載的作用下,壩踵處所受拉應力最大,壩趾所受壓應力最大。由于壩體與壩基材料不同,因變形不協調致使二者接觸面較為薄弱,存在較大的不利應力。隨著水位的上升,壩趾處壓應力有所增大,而壩踵處拉應力明顯降低了很多,這是因為壩踵區域處的材料發生了損傷,并且隨著損傷的累計,現有破壞部位便會在應力集中的作用下,沿著順水流方向和壩體縱深方向進行進一步的擴展,最后形成較大的宏觀裂縫致使建筑物失效。故在極端暴雨工況下,膠凝砂礫石壩可能會在壩踵處發生損傷破壞,且為拉裂破壞;壩體與壩基接觸面處存在較大不利應力,但影響不大。
(2)位移響應分析。膠凝砂礫石壩順水流方向上最大變形面的位移量隨暴雨歷時變化如圖11所示。

圖11 暴雨工況下膠凝砂礫石壩順水流方向位移圖Fig.11 Displacement diagram of downstream flow direction on top of CSG dam under rainstorm condition
由圖11 可知,壩體位移量隨著上游水位的上升而逐漸增大,其增長速率由大變小,增長趨勢與壩體上游水位荷載增長趨勢相同,為正相關性。從大壩整體來看,壩體中部位移量較大,壩頂次之,壩底位移量最小。
(3)損傷單元分析。為了更好表述壩體的受力狀況和損傷演化規律,提取壩體若干損傷特征單元進行分析對比,特征單元的具體位置詳見圖1,其損傷變化如圖12所示。

圖12 暴雨工況下膠凝砂礫石壩特征點的彈性模量損傷圖Fig.12 Elastic modulus damage diagram of characteristic points of CSG dam under rainstorm conditions
由圖12 分析可知,隨著上游水位不斷的升高,作用在上游壩坡的荷載逐漸增大,內部的39 961 單元無損傷,其他特征單元均產生損傷,且內部39 751 單元的損傷值小于表層39 301 單元,可見壩體開始由無損狀態向損傷狀態過渡,損傷范圍逐漸擴大。由于損傷不斷累積,特征單元的彈性模量有減小的趨勢,其趨勢由大逐漸變小直至彈性模量為零,損傷量達到極值,即特征點處結構完全失效,壩體材料被破壞。單元損傷值隨時間的變化同時也驗證了壩體應力的變化。總之,隨著水位的升高壩踵部位的損傷量會優先達到極值,壩踵區域處的淺層地基會優先發生塑性屈服,塑性屈服區會分別沿著豎直方向向下擴展和建基面方向向下游擴展。
以典型工程為例,采用有限元法對暴雨洪水入庫時的水位變化過程進行模擬,得到漫頂時期的穩定水位,即暴雨時期的最不利工況水位,以此水位為邊界條件對壩體進行應力仿真分析,得出以下結論。
(1)在水位荷載的作用下,在壩踵處壩體所受拉應力最大,在壩趾處壩體所受壓應力最大;壩體與地基接觸面處承受較大的不利應力,此部位為壩體薄弱區域。
(2)隨著水位持續升高,壩體位移響應值不斷增大,且壩中部位的響應值最大,壩頂次之,壩底最小,說明壩中部位變形量最大。
(3)對于膠凝砂礫石壩體而言,受極端強降雨的作用,損傷部位主要集中在壩踵區域處,沿著順水流方向和壩體縱深方向進行進一步的擴展,壩體與地基接觸面處存在較大不利應力,但影響不大,說明外荷載越大,損傷范圍越大,壩體越容易失穩破壞。