陳 平，朱家奇，王慶祥

(1.中國電建集團中南勘測設計研究院有限公司，湖南 長沙 410014；2.河海大學水利水電學院，江蘇 南京 210024)

近年來，我國水利水電建設發(fā)展迅速，在土石壩填筑設計理論和筑壩技術等方面均取得了長足的進步，面板堆石壩正成為當今壩工建設中的主要壩型之一。由于水電工程的復雜性和建設運行經(jīng)驗的欠缺，DL/T 5016—2011《混凝土面板堆石壩設計規(guī)范》規(guī)定對高混凝土面板堆石壩應進行變形和應力有限元計算復核；同時，考慮到堆石料本構模型尚未成熟，對于特高混凝土面板堆石壩要求采用不同本構模型進行對比驗證，以提高計算成果的精度和適用性。因此，選擇合適的本構模型對特高混凝土面板堆石壩的有限元分析十分重要。

學者們一直致力于堆石料本構關系的研究并提出了相應模型[1]。在眾多堆石體計算的本構模型當中，以鄧肯-張模型[2]為代表的彈性非線性模型，由于其參數(shù)簡單易懂、物理意義明確、適用性強、并且能夠能充分反映出土體的非線性特性而得到廣泛應用，但該模型不能反映堆石料高緊密度條件下的剪脹性[3]。沈珠江[4]采用鄧肯-張模型應力應變雙曲線將硬化規(guī)律通過顯式表示，提高了雙屈服面模型的實用性，應用較廣；但該模型高估了堆石料的剪脹性，使得堆石體的計算變形一般偏小。以往傳統(tǒng)的塑性理論常常需要提前假設出土體的塑性勢面和關聯(lián)流動法則，并不考慮其應力主軸的旋轉(zhuǎn)，而廣義塑性模型[5- 6]不需要采用塑性勢面的顯式表達式及屈服面，可以通過定義微分形式的加載方向和塑性流動方向，并簡單方便地定義出不同加載情況下的塑性模量，以完成各種加載條件下土體應力變形的預測任務。陳生水等[7]提出了基于廣義塑性理論的本構模型，模型通過改變摩擦角反映堆石料顆粒破碎導致峰值應力隨平均應力減小。朱晟等[8]采用更具普適性的剪脹方程，提出了適應于復雜應力路徑的靜動力耦合統(tǒng)一廣義塑性模型。

拉哇水電站混凝土面板堆石壩目前為世界上在建的最高的混凝土面板堆石壩，壩頂總長398.0 m，壩頂寬15.0 m，壩頂高程2 709.0 m，最大壩高239.0 m。河床深槽部位趾板設置基座混凝土結(jié)構，擬最大厚度為5 m，基座頂部高程2 475.0 m，上游面板底部壩坡坡比1∶1.45，下游壩坡在2 590.0 m高程以上布置“之”字形上壩道路，綜合坡比為1∶1.628。本文以拉哇水電站的特高混凝土面板堆石壩為例，通過采用鄧肯E-B模型和統(tǒng)一廣義塑性模型兩種計算方法，計算其壩體變形及應力分布情況，以此進一步分析堆石體和面板的變形協(xié)調(diào)特性，為壩體分區(qū)結(jié)構設計提供相應依據(jù)。

1 堆石體計算本構模型

1.1 鄧肯E-B模型[2]

堆石體在自重荷載及蓄水期水壓力的作用下往往會表現(xiàn)出明顯的非線性特征，因此可以采用充分反映土體非線性的鄧肯E-B模型，切線彈性模量Et和體積變形模量Bt的計算公式為

(1)

(2)

(3)

(4)

式中，pa為大氣壓力；K、n、Rf、Kb、m、φ0、Δφ為模型參數(shù)。

1.2 河海統(tǒng)一廣義塑性模型[8]

統(tǒng)一廣義塑性模型的彈塑性矩陣為

(5)

式中，De為彈性矩陣；Dp為塑性矩陣；ngL/U為加卸載過程中對應的塑性流動方向；HL/U為加卸載過程中對應的塑性模量；n為加載時的方向。

根據(jù)等向壓縮特性，堆石體的彈性模量為

(6)

式中，υ為泊松比，一般取0.3。

堆石體的塑性流動方向，在p-q平面內(nèi)也就是它對應的剪脹方程為

(7)

式中，β，α，Mg為粗粒土參數(shù)。

加載時的塑性流動方向為

(8)

為了模擬堆石體的“卸載體縮”現(xiàn)象，將卸載時的塑性流動方向定義為

(9)

采用非關聯(lián)流動法則，加載方向為

(10)

加載時的塑性模量為

(11)

由于顆粒破碎的影響，堆石料的峰值應力比Mf并不是一個常數(shù)，而是與平均壓力或者是圍壓相關的，將堆石料的峰值強度統(tǒng)一描述為

(12)

圖1 河床典型剖面的剖分網(wǎng)格示意

表1 壩體填筑分級加載過程

式中，Mf0與nf為材料參數(shù)，可通過常規(guī)三軸試驗峰值強度擬合得到。

1.3 流變本構模型[9]

目前混凝土面板堆石壩常采用碾壓的施工方式，在分級填筑碾壓的過程中，堆石體的流變往往會呈現(xiàn)出“繼效特性”。如果考慮到流變過程中的遺傳效應，可以采用朱晟等[8]提出的堆石料增量流變模型來進行分析。

假設混凝土面板堆石壩在ζi時刻加載第i級的應力增量，則在該級應力增量下大壩的累積流變量可以表示為(t>ζl)

(13)

(14)

式中，εv(t)、εs(t)分別為i級應力增量作用下的累積體積流變量、累計剪切流變量；Δεvfi、Δεsfi分別為第i級應力增量作用下的最終體積流變量和最終剪切流變量；c、α均為試驗常數(shù)。

切線流變體積模量Kt和剪切模量Gt分別為

(15)

(16)

式中，kv、nv、ks、ns、Rsf為模型參數(shù)，可根據(jù)室內(nèi)三軸流變試驗確定。

2 特高混凝土面板堆石壩的三維有限元計算

2.1 壩體網(wǎng)格剖分及施工加載模擬

針對拉哇特高混凝土面板堆石壩，按照前期地質(zhì)勘察所得到的地質(zhì)水文條件，可以將其地基作為剛性來考慮。因此依據(jù)河床斷面的典型特征及壩體分區(qū)，采用自動剖分程序進行有限元前處理的過程，同時結(jié)合面板分縫特點，沿壩軸線方向共計設置55個剖面，其中河床典型剖面0+188如圖1a所示，大壩三維整體剖分網(wǎng)格如圖1b所示。坐標系按照常規(guī)方向進行設置，以上游指向下游為x軸正方向，以鉛直向上的高程增加方向為y軸正方向，以沿壩軸線左岸指向右岸的方向為z軸正方向。拉哇混凝土面板堆石壩總計剖分得到18 518個結(jié)點，17 102個單元，其中包括14 154個堆石體單元，964個接觸面、面板及豎縫單元，56個周邊縫單元。

由于采用考慮流變遺傳效應的增量流變模型，故按照大壩實際填筑的分級施工順序，計算中將荷載施加的過程與設計的填筑分期過程相對應，按逐級加載的方式共分為45個加載級來進行模擬。壩體填筑加載過程及填筑程序分別見表1和圖2。

圖2 壩體填筑分期(單位：m)

2.2 本構模型參數(shù)

結(jié)合室內(nèi)三軸試驗結(jié)果及擬合過程，確定壩料廣義塑性模型計算參數(shù)如表2所示。利用上述壩料的室內(nèi)三軸試驗資料，整理鄧肯E-B模型參數(shù)如表3所示。同時根據(jù)室內(nèi)三軸流變試驗結(jié)果，整理壩料的增量流變模型參數(shù)于表4。

表2 壩料河海統(tǒng)一廣義塑性模型參數(shù)

表3 壩料鄧肯E-B模型參數(shù)

表4 壩料增量流變模型計算參數(shù)

表5 不同本構模型的堆石體有限元計算極值

2.3 不同本構模型大壩計算結(jié)果對比分析

2.3.1 堆石體的應力和變形

選取河床典型0+188縱向剖面進行不同模型計算結(jié)果的對比分析，為便于比較，將相同試驗資料的應力變形有限元計算極值列于表5。圖3為滿蓄期該典型剖面的主應力等值線，由圖3可知，兩種模型計算得到的壩體應力結(jié)果數(shù)值相近，第一主應力的等值線與壩坡平行，最大值分別為3.82 MPa和3.78 MPa，產(chǎn)生在壩底的軸線附近。統(tǒng)一廣義塑性模型計算的第三主應力最大值較大為1.70 MPa，是鄧肯E-B模型的1.4倍。

圖4為滿蓄期典型縱向剖面堆石體的豎向位移等值線。由圖4可知，兩種模型計算的最大沉降分別為209.9 cm和202.8 cm，約為壩高的0.88%和0.85%，位置均出現(xiàn)在2 600 m高程附近，且位移分布規(guī)律保持一致。由表5可得，無論是竣工期還是滿蓄期，鄧肯E-B模型計算的壩體水平向下游的變形明顯較大，滿蓄期水平向下游的位移分別為32.2 cm和21.7 cm，約為統(tǒng)一廣義塑性模型計算結(jié)果的1.5倍。

圖4 滿蓄期典型縱向剖面豎向位移等值線(單位：cm)

表6為不同工程滿蓄期的堆石體沉降極值情況。由表6可知，拉哇混凝土面板堆石壩采用兩種不同堆石體本構模型的計算結(jié)果與類似工程的實際監(jiān)測數(shù)據(jù)數(shù)值接近，壩高占比符合一般規(guī)律。從堆石體的沉降變形來看，兩種模型均可以較好地計算得到大壩變形的分布規(guī)律。

表6 不同工程滿蓄期的堆石體沉降極值比較

2.3.2 混凝土面板的應力和變形

不同本構模型面板有限元計算極值如表7所示。

圖5為滿蓄期面板軸向應力的等值線分布。由圖5可知，兩種模型計算得到的軸向應力結(jié)果分布規(guī)律類似，面板在水荷載作用下基本處于受壓狀態(tài)，兩岸有局部區(qū)域存在較小的拉應力。從數(shù)值上來看，鄧肯E-B模型計算的應力最大值為5.66 MPa，而統(tǒng)一廣義塑性模型計算最大值為2.82 MPa，兩者均處于2 550 m高程附近?？傮w而言，統(tǒng)一廣義塑性模型計算的面板軸向應力均比鄧肯E-B模型小。

圖6為滿蓄期面板撓度的等值線分布。滿蓄期時，鄧肯E-B模型和廣義塑性模型計算出的面板整體分布規(guī)律相近，與類似工程的變形經(jīng)驗相一致[10]。鄧肯E-B模型計算的撓度最大值產(chǎn)生在河床典型剖面的2 620 m高程附近，為85.7 cm；而統(tǒng)一廣義塑性模型計算的面板中下部壩高部位數(shù)值較大，最大值為50.0 cm，可以看出兩者計算得到的面板撓度極值差異較大。

圖7為滿蓄期面板軸向位移的等值線分布。在蓄水后，兩種模型計算的面板軸向位移分布規(guī)律較為一致，面板左側(cè)向右岸移動，右側(cè)向左岸移動，整體呈現(xiàn)由外部向中心擠壓的趨勢。兩者的最大變形均出現(xiàn)在左右側(cè)中部位置，鄧肯E-B模型的軸向位移最大值分別為21.1、23.5 cm，而統(tǒng)一廣義塑性模型的軸向位移最大值分別為10.7、11.0 cm，量值上鄧肯E-B模型計算值約為統(tǒng)一廣義塑性模型的兩倍。

3 結(jié) 論

利用基于分形理論的堆石料室內(nèi)縮尺試驗成果，分別采用鄧肯E-B模型及統(tǒng)一廣義塑性模型，對拉哇特高混凝土面板堆石壩進行三維有限元分析計算，研究不同堆石體本構模型的壩體內(nèi)部應力變形分布規(guī)律，得到如下結(jié)論：

(1)采用鄧肯E-B模型和統(tǒng)一廣義塑性模型計算得到的滿蓄期堆石體最大沉降分別為209.9 cm和202.8 cm，占壩高的百分比分別為0.88%和0.85%，兩者與水布埡、三板溪及洪家渡等類似工程的實際監(jiān)測數(shù)據(jù)比較一致，兩種模型均可較好地描述面板堆石壩堆石體的變形規(guī)律。

表7 不同本構模型的面板有限元計算極值

圖5 滿蓄期面板軸向應力等值線分布(單位：MPa)

圖6 滿蓄期面板撓度等值線分布(單位：cm)

圖7 滿蓄期面板軸向位移等值線分布(單位：cm)

(2)兩種本構模型計算的面板變形與應力差異較大，如滿蓄時面板撓度最大值分別為85.7 cm和50.0 cm，統(tǒng)一廣義塑性模型計算面板的拉應力區(qū)域更小，分析原因主要是鄧肯模型不能考慮堆石的剪脹性，計算的水平變形更大引起。

(3)兩種本構模型計算的壩體應力變形規(guī)律相似，計算應力與變形分布規(guī)律良好，說明了兩種模型應用于特高混凝土面板堆石壩有限元計算分析的可靠性和合理性，計算結(jié)果可為后續(xù)特高壩工程的設計施工提供依據(jù)。