方衛(wèi)華，王潤英，許珉凡

（1.水利部南京水利水文自動化研究所，江蘇南京，210012；2.河海大學(xué)水利水電學(xué)院，江蘇南京，210098）

1 概述

汾河二庫位于太原市尖草坪區(qū)與陽曲縣交界的懸泉寺，下游距太原市區(qū)30 km，是汾河上游干流上一座具有綜合效益的大（2）型水利樞紐工程。樞紐工程由攔河大壩、供水發(fā)電洞和水電站組成。其中河床壩基為寒武系鳳山組巖層，上部為竹葉狀白云巖和薄層及中厚層白云巖，其表層為灰色白云質(zhì)頁巖、竹葉狀白云巖互層與條帶狀含泥白云巖，該層層厚為9.03 m，其中灰色白云質(zhì)頁巖、竹葉狀白云巖互層發(fā)育在上部，厚約1.5～2.0 m。壩址巖層主要為寒武系上統(tǒng)鳳山組和奧陶系下統(tǒng)，透水性微弱。

2007年7 月汾河二庫主體工程驗收時，遺留的工程有：部分固結(jié)灌漿工程、部分帷幕灌漿工程、全部接觸灌漿工程及并縫灌漿工程、壩體排水工程。

2014年開始進行以壩肩帷幕灌漿、壩基固結(jié)灌漿和壩后連續(xù)墻等施工為主的應(yīng)急除險加固。

汾河二庫攔河大壩為碾壓混凝土重力壩，考慮到大壩除險加固實施中河床溢流壩段并縫灌漿、壩基地下連續(xù)墻及固結(jié)灌漿結(jié)構(gòu)空間特性，對河床溢流壩段建立三維有限元模型，分析其加固后正常蓄水位工況條件下的結(jié)構(gòu)性態(tài)，評價除險加固工程的效果，為水庫正常運行提供決策依據(jù)。

2 數(shù)值計算方法及計算軟件選擇

2.1 應(yīng)力場、位移場求解的基本理論

考慮滲流場、應(yīng)力場、位移場三場耦合，鑒于規(guī)范[1]中混凝土重力壩的應(yīng)力控制標(biāo)準(zhǔn)基于考慮靜水壓力、動水壓力、揚壓力等水荷載作為外力的情況，對應(yīng)的方式為以材料和水的混合體作為研究對象，考慮總應(yīng)力狀態(tài)。這種情況下以應(yīng)力分量和位移分量表示的滲流場、應(yīng)力場、位移場耦合的微分方程分別如式（1）和式（2）所示：

式中：σij為總應(yīng)力張量，以拉應(yīng)力為正；γc為壩體材料或基巖的飽和容重；Fi為自重以外、包含水壓荷載的外荷載i=x,y,z；為剪切模量 ；為 Laplace算子 ；為體積應(yīng)變。

2.2 大型方程組求解方法及軟件選用

為滿足計算效率和計算精度的要求，方程組求解采用改進迭代格式的逐步超松弛共軛梯度法，該迭代計算方法的系數(shù)矩陣只需要存儲非零元素，節(jié)省計算機存儲容量，而且求解效率高，能保證多自由度有限元計算的有效進行[2]。

2.3 基于COMSOL Multiphysics的滲流計算

2.3.1 滲流方程選擇

基于Darcy定律及水流連續(xù)性方程的一般滲流偏微分方程為：

式中：ρ為流體密度；κ為滲透率；μ為流體黏度；p為壓力；D為位置水頭；Qm為源匯項；ε為孔隙率。

基于Rechards方程的滲流偏微分方程為：

其中：Cm為容水度；Se為有效飽和度；S為儲水系數(shù)；κS為飽和滲透率；κr為相對滲透率?？紤]穩(wěn)定滲流場，則式（4）的偏微分方程中不考慮這一項。

2.3.2 邊界條件

計算穩(wěn)定滲流場時，在COMSOL Multiphysics邊界條件中考慮定水頭邊界條件和用透水層邊界模擬的混合邊界條件。

上下游水位以下的定水頭邊界條件可表達為：

式中：H0為給定的水頭；D為高程。

混合邊界條件可表達為：

式中：Rb為傳導(dǎo)率；Hb為外部壓力水頭；n為外法線方向；u為滲透壓力。

由于滲流自由面及溢出點是未知的，溢出邊界屬于混合邊界，對于復(fù)雜滲控系統(tǒng)處理起來較困難。非飽和滲流計算中溢出邊界是一個特殊的邊界條件。

采用多物理場耦合軟件COMSOL Multiphysics進行滲流及多物理場耦合滲流分析。COMSOL Multiphysics是一款基于有限單元法的數(shù)值仿真軟件，以一般偏微分方程或偏微分方程組為基礎(chǔ)建立模型，具有較強的計算性能和多物理場直接耦合能力。COMSOL Multiphysics中提供的透水層邊界是一種混合邊界條件，并已通過若干算例驗證了技術(shù)路線的合理性和可行性，可以實現(xiàn)對這一類邊界條件的定義[3-4]。這類混合邊界實際上是水頭邊界與流量邊界的組合，求解時通過計算得到壓力p的分布進行判斷，基于此種情況，可以采用飽和-非飽和滲流模型及固定網(wǎng)格法進行迭代求解來確定混合邊界條件下的滲流場。

2.3.3 壩基排水管模擬

文獻[5]、[6]的計算結(jié)果表明，“以縫代井列”方法能模擬巖體排水孔幕處的流場性態(tài)，這種方法適用于排水孔間距小于6 m的排水井列。在水利水利水電工程中，排水孔間距一般都在5 m內(nèi)[5]，此方法滿足要求。在上述文獻研究的基礎(chǔ)上，在COM-SOL Multiphysics中將“以縫代井列”發(fā)展為“以面代井列”方法，建立模型時在排水管列的位置設(shè)置內(nèi)部邊界面，無需細分排水管單元，設(shè)置邊界條件時將壩基排水孔幕設(shè)置為內(nèi)部邊界，用透水層邊界條件來表達，其外部水頭值取為排水井（管）口高程。

2.3.4 應(yīng)力求解及標(biāo)準(zhǔn)

在COMSOL Multiphysics求解時，先用固體力學(xué)模塊計算基巖在自重作用下的應(yīng)力場。在壩體和壩基整體計算時，也采用固體力學(xué)模塊，考慮滲流的作用，并且將采用固體力學(xué)模塊計算得到的壩基應(yīng)力作為整體計算時的壩基的預(yù)應(yīng)力。在此基礎(chǔ)上建立三維有限元模型應(yīng)當(dāng)滿足SL 319-2005《混凝土重力壩設(shè)計規(guī)范》6.3.4條規(guī)定。

3 溢流壩段三維有限元分析

3.1 三維有限元模型建立

三維有限元模型中，壩基上、下游和鉛直方向各取1.5倍壩高范圍。針對溢流壩段除險加固后的結(jié)構(gòu)特點和幾何拓撲性質(zhì)，模型中考慮溢流壩段的結(jié)構(gòu)型式、壩體材料分區(qū)、基巖材料分區(qū)及帷幕灌漿、排水、固結(jié)灌漿的布置等因素。模型中考慮混凝土地下連續(xù)墻結(jié)構(gòu)，壩基灌漿帷幕的厚度取為1.5 m。

3.2 三維有限元邊界條件

滲流場計算時，上、下游水位以下取為等水頭邊界，其水位以上部分為混合邊界，所有廊道表面為可能溢出邊界，其余邊界為不透水邊界。

位移場、應(yīng)力場計算時，壩基底部視為固定約束，壩基兩側(cè)x方向取為連桿約束，壩基及并縫灌漿高程（855.0 m）以下壩體左右岸方向（y方向）兩側(cè)取為連桿約束。

3.3 三維有限元計算參數(shù)

應(yīng)力場、位移場計算時，壩體和壩基均采用塑性本構(gòu)關(guān)系，屈服準(zhǔn)則采用匹配Mohr-Coulomb準(zhǔn)則的Drucker-Prager屈服準(zhǔn)則，具體力學(xué)參數(shù)見表1。

下文以河床溢流典型壩段（0+099.0～0+147.2）在上游水位高程為905.7 m、下游水位高程為855.7 m、壩前淤沙高程為875.0 m的工況下為例，進行應(yīng)急除險加固后的大壩安全性態(tài)分析。

3.4 溢流壩段計算結(jié)果及分析

3.4.1 滲流場

表1 壩基材料的計算參數(shù)Table 1 Calculation parameters of dam foundation materials

表2 壩基強度計算參數(shù)Table 2 Parameters of dam foundation strength calculation

除險加固之后，壩基抽排設(shè)施正常發(fā)揮作用，壩體和壩基滲流場產(chǎn)生變化，壩體內(nèi)壓力水頭和壩基面滲壓水頭與實際值相比得到一定削減。三維有限元計算得到的滲流場等值面圖見圖1。

圖1 滲流場壓力水頭等值面圖Fig.1 Isopleth map of seepage field

3.4.2 位移場

圖2～4為壩體x方向位移、z方向位移和總位移云圖，其最大值和最小值及它們發(fā)生的位置列于表3中。結(jié)合圖、表分析可知：壩體在x方向產(chǎn)生向下游的位移，最大值為3.92 mm，發(fā)生在堰頂，最小值為0.095 5 mm，發(fā)生在基槽下游側(cè)壩趾處。壩體在z方向發(fā)生向下的沉降，絕對值最大為3.17 mm，發(fā)生在挑流鼻坎處；絕對值最小為0.840 mm，發(fā)生在壩體底部混凝土墊層上游側(cè)壩踵處。總位移最大值為4.92 mm，發(fā)生在堰頂略偏下游側(cè)；最小值為0.845 mm，發(fā)生在基槽下游側(cè)壩趾處。

3.4.3 應(yīng)力

應(yīng)力以拉為正，以壓為負。壩體第三主應(yīng)力及垂直應(yīng)力云圖見圖5和圖6。壩體第三主應(yīng)力及垂直應(yīng)力最大值和最小值及它們發(fā)生的位置列于表4中。圖7為垂直應(yīng)力的拉應(yīng)力分布示意圖，圖8為垂直應(yīng)力局部位置拉應(yīng)力分布示意圖。

圖2 x位移云圖Fig.2 Cloud map of x-direction displacement

圖3 z位移云圖Fig.3 Cloud map of z-direction displacement

圖4 總位移云圖Fig.4 Cloud map of total displacement

表3 工況5壩體位移最大值和最小值及其位置Table 3 Maximum and minimum displacement of dam body and their locations

圖5 第三主應(yīng)力云圖Fig.5 Cloud map of minor principal stress

圖7 垂向拉應(yīng)力分布示意圖Fig.7 Distribution of vertical tensile stress

圖8 局部位置垂向拉應(yīng)力分布示意圖Fig.8 Distribution of vertical tensile stress in some components

表4 壩體應(yīng)力最大值和最小值及其位置Table 4 Maximum and minimum stress of dam body and their locations

計算顯示：壩體第三主應(yīng)力的最小值為-5.17MPa，即最大主壓應(yīng)力為5.17 MPa，發(fā)生在下游廊道底部的上游側(cè)。從圖8可以看出，在廊道表面附近、上游壩踵、挑流鼻坎、反弧段、堰頂都有拉應(yīng)力產(chǎn)生，其中，堰頂最大拉應(yīng)力為0.003 85 MPa，上游壩踵處的最大拉應(yīng)力為0.010 1 MPa，反弧段最大拉應(yīng)力為0.098 5 MPa。

根據(jù)調(diào)查，壩基寒武系上統(tǒng)白云巖的飽和抗壓強度為95.59～179.60 MPa，奧陶系下統(tǒng)白云巖的飽和抗壓強度為67～239 MPa[7]；C20混凝土的軸心抗壓強度標(biāo)準(zhǔn)為13.4 MPa，軸心抗拉強度標(biāo)準(zhǔn)值為1.54 MPa。壩踵部位拉應(yīng)力較小，且拉應(yīng)力區(qū)寬度小于壩底寬度的0.07倍且小于壩踵至帷幕中心線的距離，均滿足規(guī)范要求。

4 結(jié)語

汾河二庫除險加固后大壩安全性態(tài)的計算結(jié)果顯示，各物理場分布基本符合相關(guān)規(guī)律，表明COMSOL Multiphysics軟件可用于多場耦合分析。計算結(jié)果顯示，除險加固后大部分變形、滲流和應(yīng)力在合理區(qū)間內(nèi)，除險加固工程基本達到預(yù)期效果。但應(yīng)力場計算結(jié)果顯示，在廊道表面附近、上游壩踵、挑流鼻坎、反弧段、堰頂都有拉應(yīng)力產(chǎn)生，應(yīng)加強相關(guān)區(qū)域相關(guān)物理量的安全監(jiān)測（如廊道裂縫、滲流、鋼筋銹蝕等的監(jiān)測），進一步檢驗施工效果和計算的準(zhǔn)確性。

[1]SL 319-2005,混凝土重力壩設(shè)計規(guī)范[S].

[2]林紹忠.對稱逐步超松弛預(yù)處理共軛梯度法的改進迭代格式[J].數(shù)值計算與計算機應(yīng)用,1997,18(4):266-270.

[3]Chui T M,Freyberg D L.The use of COMSOL for integrated hydrological modeling[C].Proc.,COMSOL Conf.2007 Boston,COMSOL AB.Newton.Mass.,2007：217-223.

[4]徐軼,徐青.基于COMSOL Multiphysics的滲流有限元分析[J].武漢大學(xué)學(xué)報(工學(xué)版),2014,47(2):165-170.

[5]王恩志,王洪濤,王慧明.“以縫代井列”——排水孔幕模擬方法探討[J].巖石力學(xué)與工程學(xué)報,2002,21(1):98-101.

[6]杜京濃,宋漢周,霍吉祥,等.混凝土重力壩基礎(chǔ)排水孔模擬方法對比分析[J].勘察科學(xué)技術(shù),2015(2):9-13.

[7]水利部水利水電規(guī)劃設(shè)計總院.山西省汾河二庫水利樞紐工程蓄水安全鑒定報告[R].1993.