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

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

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

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

1 計算理論

1.1 鄧肯E-B模型

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

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

式中：Kur為量綱一模型參數。

1.2 沈珠江雙屈服面模型

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

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

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

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

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

1.3 MPZG模型

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

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

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

對于塑性部分，有

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

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

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

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

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

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

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

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

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

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

2 計算模型與參數

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層反濾層，上下游反濾層與堆石體之間設置過渡區。心墻與巖石基礎面之間設置接觸黏土層和混凝土蓋板，以適應岸坡剪切變形，防止接觸沖刷。心墻堆石壩典型斷面圖見圖1。

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

2.1 計算模型及工況

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

圖2 三維網格模型Fig.2 Three-dimensional grid model

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

2.2 計算參數

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

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

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

表3 MPZG模型參數Table 3 MPZG model parameters

3 不同本構模型適用性分析

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

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

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

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

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

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

3.2 壩體應力分析

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

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

3.3 壩體變形分析

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

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

3.4 模型適用性討論

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

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

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

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

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

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

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

表4 中國典型高心墻壩沉降監測值統計Table 4 Settlement monitoring statistics of typical high core wall dams in China

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

4 結論

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

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

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