高心墻堆石壩材料本構(gòu)模型計算的適用性研究

程瑞林，汪涇周，范欽煜，湛正剛，周偉,3，馬剛,3

(1.武漢大學 水資源工程與調(diào)度全國重點實驗室，湖北 武漢，430072；2.中國電建集團 貴陽勘測設計研究院有限公司，貴州 貴陽，550081；3.武漢大學 水工程科學研究院，湖北 武漢，430072)

近年來，隨著堆石壩設計理論的發(fā)展與高壩填筑技術(shù)的革新，我國高堆石壩的建設水平逐漸從200 m級上升到300 m級。相比于早期建設的中低土石壩，已建高堆石壩監(jiān)測結(jié)果表明，目前的堆石壩設計理論在高壩壩體的應力、變形控制方面仍存在不足，高堆石壩數(shù)值計算精度較低問題突出。為解決這一問題，學者們主要從材料本構(gòu)關(guān)系和模型參數(shù)確定2個方面展開了研究。在對材料本構(gòu)的研究方面，ALONSO等[1-4]提出了BBM(barcelona basic model)模型和RM(rock-fill model)模型，并將其用于預測Beliche壩的變形和心墻孔隙壓力。姚仰平等[5-6]在修正劍橋模型基礎上通過引入統(tǒng)一硬化參數(shù)，對修正劍橋模型進行改進，使其能夠統(tǒng)一反映土體的剪脹和剪縮現(xiàn)象，并建立了考慮顆粒破碎的粗粒土本構(gòu)模型。PASTOR等[7]提出了P-Z廣義塑性模型，該模型直接基于剪脹方程定義塑性流動方向和加載方向，并能靈活地定義不同加載條件下的塑性模量。董威信[8]在P-Z模型的基礎上建立了MPZG模型，將臨界狀態(tài)理論和狀態(tài)參數(shù)引入P-Z模型中。在對模型參數(shù)的研究方面，一些學者致力于改進堆石料材料參數(shù)的獲取方法?？讘椌┑萚9-10]研發(fā)了超大型三軸儀，通過降低縮尺程度提高了參數(shù)的準確性。周偉等[11-13]采用連續(xù)離散耦合分析方法對堆石料進行數(shù)值試驗，分析了筑壩堆石料的縮尺效應，提高了參數(shù)的合理性。另一些學者則致力于改進反演方法以便從監(jiān)測資料中獲得更精準的材料參數(shù)。朱晟等[14]提出了基于增量分析的堆石壩瞬變-流變聯(lián)合反演方法，并結(jié)合壩內(nèi)水管式沉降儀和測斜沉降管對南歐江7級面板壩進行了聯(lián)合反演分析。楊荷等[15-16]提出了基于改進響應面法的反演分析方法，并運用該方法進行了水布埡面板壩的材料參數(shù)的反演計算。盡管學者們進行了大量關(guān)于單一本構(gòu)模型及其參數(shù)取值方法的研究，但對于多種本構(gòu)模型在實際高堆石壩計算中的優(yōu)劣與適用性評價研究較少，新建高堆石壩工程究竟應采用何種模型進行計算評估仍有待考慮。本文以此為出發(fā)點，選擇非線性彈性的鄧肯E-B模型、傳統(tǒng)彈塑性的沈珠江雙屈服面模型和基于廣義塑性框架的MPZG模型這3類高堆石壩計算中較具代表性的模型，采用室內(nèi)試驗方法得到的參數(shù)，對筑壩堆石料在三軸排水條件下宏觀力學響應進行驗證，并對壩體應力、變形進行計算，探討這3種模型在高心墻堆石壩應力、變形計算中的適用性問題，以期為其他高心墻堆石壩工程的應力、變形計算提供借鑒與參考。

1 計算理論

1.1 鄧肯E-B模型

DUNCAN等[17]基于雙曲線模型進行了應力-應變關(guān)系的非線性分析，并提出了E-ν模型。由于E-ν模型在實際應用時計算得到的泊松比通常大于試驗測定結(jié)果，DUNCAN在E-ν模型基礎上采用切線體積模量代替切線泊松比，得到了E-B模型。在該模型中，砂土的彈性模量和體積模量采用切線值進行計算，切線模量Et和切線體積模量Bt計算式分別如式(1)和式(2)所示。

式中：S為應力；Ei為初始切線模量；σ1和σ3分別為大、小主應力；Rf為破壞比，表示砂土實際破壞強度與極限破壞強度的比值；Pa為大氣壓強；K、n、Kb和m均為量綱一模型參數(shù)；c和?分別為Mohr-Coulomb強度破壞準則中的黏聚力和內(nèi)摩擦角，其中，黏聚力一般認為是常數(shù)，內(nèi)摩擦角則可利用材料參數(shù)?0和Δ?及圍壓σ3采用非線性模型?=?0-Δ?lg(σ3/Pa)計算得到。此外，為模擬砂土在加卸載時宏觀響應的差異性，砂土卸載時的切線模量可采用卸荷模量Eur進行計算。

式中：Kur為量綱一模型參數(shù)。

1.2 沈珠江雙屈服面模型

沈珠江雙屈服面模型屬于彈塑性模型[18]，分別采用橢圓曲線和冪曲線作為其屈服面，屈服函數(shù)定義為

式中：p和q分別為平均主應力和廣義剪應力；r和s為根據(jù)土性特點調(diào)整的屈服面參數(shù)，對于粗粒土，二者均可取值為2，對細粒土可分別取2和3。根據(jù)彈塑性力學的基本原理，砂土的變形Δε即為彈性部分變形Δεe和塑性部分變形Δεp之和，二者可分別按照胡克定律和正交流動法則計算。

式中：C為彈性剛度張量，可由切線模量Et和切線泊松比μt計算得到；Δσ為應力張量增量；Δf1和Δf2分別為屈服函數(shù)f1和f2的增量；A1和A2分別為相應塑性系數(shù)，反映2個屈服面各自產(chǎn)生的塑性應變的大小，可按式(7)計算。

式中：η=q/p，為應力比；Kur和Gur分別為按照卸荷模量Eur和卸荷泊松比μur計算得到的體積模量和剪切模量；切線模量Et和卸荷模量Eur可參照E-B模型分別按照式(1)和式(3)進行計算，卸荷泊松比μur取經(jīng)驗值0.3，切線泊松比μt可以依照式(8)計算。

式中：cd、nd、Rd和Rf均為模型參數(shù)。

1.3 MPZG模型

PASTOR等[7]提出的P-Z模型是一種廣義塑性模型，MPZG模型則為在其基礎上引入臨界狀態(tài)理論和狀態(tài)參數(shù)后的改進模型。在MPZG模型中，同樣將砂土的應變增量分解成彈性和塑性應變增量2個部分。

對于彈性部分，可按下式計算體積模量Kev和剪切模量Ges：

式中：K0和G0分別為彈性階段的體積模量和剪切模量；e為孔隙比；Pa為大氣壓強。

對于塑性部分，有

式中：Δλ為塑性乘子；De為彈性張量；ngL/U為塑性流動方向；nf為加載方向；HLU為塑性模量。

MPZG模型的剪脹方程采用劍橋模型剪脹方程的修正形式，即

式中：η為應力比；αg和mg為模型參數(shù)；Mg為臨界狀態(tài)的應力比；Md為Mg的修正值；ψ為狀態(tài)參數(shù)，能夠反映土體狀態(tài)對應力、應變的影響。

式中：ec為當前有效應力對應的臨界孔隙比；eΓ、λc和ζ均為模型參數(shù)，對于砂土材料ζ，可依據(jù)經(jīng)驗取0.7。

根據(jù)剪脹方程可得p-q應力空間中的塑性流動方向ngL為

加載方向nf采用與塑性流動方向相似的格式，

式中：αf為模型參數(shù)；Mf為峰值應力比。塑性模量HLU分別考慮加載和卸載2種情況。加載時，塑性模量HL可按式(16)計算。

式中：H0為模型參數(shù)；Hf、Hv、HDM分別為描述偏應變映射法則、相變轉(zhuǎn)換線相關(guān)特性及邊界面塑性映射規(guī)則的系數(shù)，三者可分別按照式(18)、式(19)、式(20)進行計算。

式中：HL0、βH、mb、γ均為模型參數(shù)。卸載時，塑性模量可按式(21)計算。

式中：HU0為模型參數(shù)；ηU為卸載時的應力比；Hden為表征粗粒土壓密和硬化效應的硬化項，可根據(jù)體積應變εv和模型參數(shù)rden按式(22)計算。

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

RM心墻堆石壩壩頂高程為2 907.00 m，壩頂寬為18 m，河床段心墻建基面高程為2 592.00 m，最大壩高315.00 m，防浪墻頂高程為2 908.20 m。壩頂長為650.20 m，壩體最大底寬約1 250 m。大壩上游坡比為1:2.1，在高程2 860.00 m和2 810.00 m處分別設置5 m寬馬道，下游綜合坡比為1:2.0。防滲心墻采用直心墻，心墻頂寬5 m，上下游坡比為1:0.23，心墻上下游分別設置2層反濾層，上下游反濾層與堆石體之間設置過渡區(qū)。心墻與巖石基礎面之間設置接觸黏土層和混凝土蓋板，以適應岸坡剪切變形，防止接觸沖刷。心墻堆石壩典型斷面圖見圖1。

圖1 心墻堆石壩典型斷面圖Fig.1 Typical sectional view of core wall rockfill dam

2.1 計算模型及工況

三維網(wǎng)格模型沿橫河方向單元邊長均為20 m，共包括16 703個結(jié)點和18 623個單元。圖2和圖3所示分別為整體三維計算網(wǎng)格及其順河向最大斷面剖面圖。計算工況為壩體分級填筑至壩頂高程2 907 m，并蓄水至正常蓄水位2 895 m，模擬的填筑施工和蓄水過程如圖4所示。

圖2 三維網(wǎng)格模型Fig.2 Three-dimensional grid model

圖3 三維網(wǎng)格模型順河向最大斷面剖面圖Fig.3 The maximum section view of the three-dimensional mesh model along the river

2.2 計算參數(shù)

常規(guī)三軸試驗所得壩體堆石料Ⅰ、堆石料Ⅱ和過渡料關(guān)于鄧肯E-B模型、沈珠江雙屈服面模型和MPZG模型這3種模型的計算參數(shù)分別如表1～3所示。

表1 鄧肯-張E-B模型參數(shù)Table 1 Duncan-chang E-B model parameters

表2 沈珠江雙屈服面模型參數(shù)Table 2 Shen Zhujiang double yield surface model parameters

表3 MPZG模型參數(shù)Table 3 MPZG model parameters

3 不同本構(gòu)模型適用性分析

3.1 三軸應力路徑的宏觀力學響應分析

分別采用鄧肯E-B模型、沈珠江雙屈服面模型和MPZG模型及前文所述的模型參數(shù)對堆石料Ⅰ、堆石料II、過渡料進行三軸排水試驗的數(shù)值模擬，試驗圍壓σ3分別設置為400、1 200和2 000 kPa，得到宏觀響應曲線，分別如圖5～7所示，其中沈珠江雙屈服面模型簡稱為SZJ模型。

圖5 堆石料I常規(guī)三軸排水試驗結(jié)果對比Fig.5 Comparison of conventional triaxial drainage test results of rockfill I

圖6 堆石料II常規(guī)三軸排水試驗結(jié)果對比Fig.6 Comparison of conventional triaxial drainage test results of rockfill II

圖7 過渡料常規(guī)三軸排水試驗結(jié)果對比Fig.7 Comparison of conventional triaxial drainage test results for transition materials

從試驗結(jié)果來看，不同試驗堆石料的宏觀力學響應規(guī)律一致。在應力響應上，3種模型均能較好地反映堆石料的應力響應，在較高圍壓下，MPZG模型的擬合程度稍差于鄧肯E-B模型和沈珠江雙屈服面模型的擬合程度。在體積應變響應上，鄧肯E-B模型不能反映材料的剪脹特性，在低圍壓下擬合效果很差；沈珠江雙屈服面模型和MPZG模型雖然能反映材料的剪脹性，但在低圍壓下，二者剪脹程度也存在差異。在高圍壓下，堆石料的體積響應以體縮為主，3種模型均能較好地吻合。綜合而言，較高圍壓下的鄧肯E-B模型和沈珠江模型宏觀力學響應比MPZG模型的略優(yōu)，在較低圍壓下，沈珠江雙屈服面模型和MPZG模型的宏觀力學表現(xiàn)明顯優(yōu)于鄧肯E-B模型的宏觀力學表現(xiàn)。

3.2 壩體應力分析

分別采用鄧肯E-B模型、沈珠江雙屈服面模型和MPZG模型及前文所述的模型參數(shù)對如美心墻壩三維有限元模型進行滿蓄期的有限元計算，得到壩體大主應力、小主應力等值線分別如圖8和圖9所示。

圖8 滿蓄期壩體大主應力分布Fig.8 Distribution diagram of major principal stress of dam body in full storage period

圖9 滿蓄期壩體小主應力分布Fig.9 Distribution diagram of minor principal stress of dam body in full storage period

從計算結(jié)果來看，3種模型的應力計算結(jié)果相當，應力分布規(guī)律基本一致。其中，3種模型的大主應力最大值均出現(xiàn)在心墻下游側(cè)和過渡料部位壩底與建基面交界處，分別為4.24、4.30和4.35 MPa；小主應力最大值均出現(xiàn)在心墻內(nèi)下游下部，分別為1.87、2.27和1.79 MPa。在應力分布上，應力等值線基本平行于壩坡，由外向內(nèi)逐漸增大，上游應力低于下游應力且表層出現(xiàn)小范圍的拉應力區(qū)，E-B模型計算所得拉應力區(qū)范圍略大于其他2種模型的拉應力區(qū)范圍。

3.3 壩體變形分析

分別采用鄧肯E-B模型、沈珠江雙屈服面模型和MPZG模型及前文所述的模型參數(shù)對如美心墻壩三維有限元模型進行滿蓄期的有限元計算，得到壩體順河向位移和沉降等值線分別如圖10和圖11所示。

圖10 滿蓄期壩體順河向位移等值線圖Fig.10 Contour map of displacement along river of dam body in full storage period

圖11 滿蓄期壩體沉降等值線圖Fig.11 Contour map of settlement of dam body in full storage period

從計算結(jié)果來看，3種模型的壩體變形計算總體規(guī)律一致，但具體數(shù)值尤其是順河向位移差異明顯。3種模型的順河向位移最大值均出現(xiàn)在壩體中上部，分別為91.4、47.3和69.3 cm；壩體沉降最大值同樣位于壩體中部，分別為207.4、199.7和183.9 cm。在變形分布上，順河向位移等值線呈“V”形與上下游壩坡斜交，在近壩頂處出現(xiàn)閉環(huán)，且上游壩坡位移梯度大于下游位移梯度，沉降變形呈現(xiàn)同心環(huán)狀分布，上下游近似對稱分布，自內(nèi)向外沉降逐漸減小。

3.4 模型適用性討論

首先，從宏觀力學方面來看，在三軸試驗條件下，3種模型的應力響應比較接近，均能反映堆石料的宏觀力學特征，因而，在壩體應力計算中，3種模型得到的最大和最小主應力分布均較接近。但是在變形響應上，鄧肯E-B模型不能考慮堆石料的剪脹特性，中低圍壓下的體積應變與室內(nèi)試驗結(jié)果相差較大，沈珠江雙屈服面模型和MPZG模型均能反映堆石料的剪脹特性，二者所得結(jié)果與室內(nèi)試驗結(jié)果更接近，因而，在壩體變形計算結(jié)果尤其是順河向位移上，沈珠江模型和MPZG模型計算結(jié)果更接近實際結(jié)果?？梢?，沈珠江雙屈服面模型和MPZG模型變形計算結(jié)果更具有說服力，更適用于堆石壩有限元計算分析。

其次，從應力、變形分布上看，盡管不同工程的壩體應力、變形有所差別，但具有大致相同的分布特點[19]。綜合參考文獻[19-23]中高心墻堆石壩壩體應力、變形結(jié)果，高心墻堆石壩應力、變形的一般規(guī)律為：

1) 在應力方面，心墻區(qū)存在明顯的拱效應，由于心墻料比壩殼料軟，心墻部分的相對變形大，壩殼堆石料部分的相對變形小，心墻與壩殼料交界處由于變形不協(xié)調(diào)發(fā)生應力重分配，使得心墻部分的應力轉(zhuǎn)移向兩側(cè)壩殼料，顯著降低了心墻區(qū)的壓應力；同時，心墻兩側(cè)的過渡料附近則出現(xiàn)密集的應力等值線。

2) 在變形方面，由于心墻料比壩殼堆石料軟，更容易發(fā)生明顯的變形。在順河向位移上，受兩邊不對等的水壓作用，水平位移基本出現(xiàn)在下游，最大值出現(xiàn)在心墻中部或下游側(cè)壩殼料中上部；在豎直沉降上，其最大值出現(xiàn)在心墻內(nèi)，大約位于壩高中部。

結(jié)合本文3種模型的應力、變形計算結(jié)果可知，3種模型得到的應力、變形分布規(guī)律均滿足堆石壩應力、變形的一般性分布規(guī)律，表明這3種模型都能為高心墻堆石壩建設提供參考。

此外，由于E-B模型較為簡單，較早地被引入到堆石壩工程設計中，更容易通過工程類比等方法從已建工程中獲得材料參數(shù)。由此可見，在高心墻堆石壩建設初期，如果缺少試驗數(shù)據(jù)來確定材料參數(shù)，可以通過工程類比方法從已建工程中獲取，進行初步的應力、變形合理性分析，為設計方案比選提供計算依據(jù)。

表4所示為我國幾座典型心墻堆石壩的沉降監(jiān)測值統(tǒng)計結(jié)果，沉降最大值相對壩高的百分數(shù)均在1%以上，而本文3種模型的最大沉降計算結(jié)果為最大壩高的0.66%～0.58%，均小于已建高壩沉降監(jiān)測值，其原因如下：一是本文直接采用了由室內(nèi)縮尺級配試驗得到的材料參數(shù)，忽略了堆石料縮尺效應對材料參數(shù)的影響[11]；二是因為本文沒有考慮流變、濕化等因素引起的堆石料力學性能的劣化，高估了堆石料的力學性能[24-25]。在實際工程計算中，應考慮這兩者的影響才能得到更合理的計算結(jié)果。

表4 中國典型高心墻壩沉降監(jiān)測值統(tǒng)計Table 4 Settlement monitoring statistics of typical high core wall dams in China

綜上分析可知，本文所采用的3種模型計算得到的壩體應力、變形分布均符合一般性認識規(guī)律，均能為工程設計施工提供一定參考，其中E-B模型不能反映堆石料的剪脹特性，在壩體變形計算上有所欠缺，而沈珠江雙屈服面模型和MPZG模型等彈塑性模型能反映堆石料的剪脹特性。在工程建設初期，可通過工程類比獲取E-B模型參數(shù)，進行工程建設的早期評價，初步論證方案的可行性。

4 結(jié)論

1) 3種模型均能較好地模擬應力響應，但在體積響應的模擬上均存在不同程度的缺陷。整體上看，較高圍壓下的鄧肯E-B模型和沈珠江模型的宏觀力學表現(xiàn)稍優(yōu)于MPZG模型的宏觀力學表現(xiàn)，在較低圍壓下，沈珠江雙屈服面模型和MPZG模型的宏觀力學表現(xiàn)明顯優(yōu)于鄧肯E-B模型的宏觀力學表現(xiàn)。

2) 3種模型的應力、變形計算結(jié)果整體上均符合一般規(guī)律。3種模型中大、小主應力分布相似且數(shù)值相近。鄧肯E-B模型的變形計算值最大，沈珠江雙屈服面模型的順河向位移計算值較小，MPZG模型的沉降計算值較小。

3) 鄧肯E-B模型和沈珠江雙屈服面模型形式更簡單，在工程上應用較多；MPZG模型和沈珠江雙屈服面模型在理論上能反映堆石料的剪脹特性。3種模型均能反映高土石壩的應力、變形分布特點，建議在建設初期采用鄧肯E-B模型進行比選，在建設后期采用沈珠江雙屈服面模型或MPZG模型分析成果作為結(jié)構(gòu)設計的重要依據(jù)。