基于多種材料模型的混凝土/基巖膠結面現場抗剪試驗的數值模擬研究

王義鵬，賴國偉，程 亮，余 巍，朱 東

(1.中南電力設計院有限公司，武漢 430071；2.武漢大學 水資源與水電工程國家重點實驗室，武漢 430072)

0 引 言

重力壩作為人類最早使用的一種水壩壩型，其抗滑穩定一直是學者比較關心的問題，同時，重力壩與基巖的膠結面作為一種特殊的結構面，對重力壩的抗滑穩定起著控制作用。對于沿壩基面的抗滑穩定分析，《混凝土重力壩設計規范(SL319-2018)》[1]推薦了兩種計算公式：摩擦公式和抗剪強度公式。這兩種公式中的抗剪摩擦系數f和抗剪斷摩擦系數f′、抗剪斷凝聚力c′的取值對抗滑穩定有很大的影響，如果選取的數值偏大，則壩基的穩定性得不到保證；反之，則偏于保守，造成不必要的浪費。

《混凝土重力壩設計規范(SL319-2018)》[1]規定壩體混凝土/基巖膠結面的f、f′、c′的取值在可行性研究階段及其后的各設計階段應經現場抗剪試驗確定。李仕勝(1997)[2]認為巖性和試驗正應力大小是影響膠結面抗剪強度的主要因素，林偉平等[3]認為膠結面的粗糙度對抗剪強度有顯著影響，Haberfield和Seidel[4]通過總結室內試驗的結果提出了基于分形幾何概念的模擬剪切行為的理論模型。但物理模型試驗不能反映剪切過程中的應力狀態，不能讓人很直觀地了解整個試驗的剪切破壞過程，而且現場抗剪試驗是在假設剪切面的應力分布為均勻分布的條件下得出的，與實際的破壞過程不相符，其結果必然存在誤差。如今，隨著有限元計算方法的快速發展以及其方法的通用性和有效性，有限元數值計算在工程技術界得到了廣泛的應用，王宏碩等[5]通過平面彈性有限元計算，指出混凝土試塊與基巖膠結面上的應力分布是極不均勻的；鄧建文[6]在王宏碩等[5]的成果上采用三維彈性有限元進一步驗證了剪切面上的應力不均勻性及彈模比對膠結面應力分布的影響。前人的研究大多只使用單一的材料模型，缺乏不同材料模型之間的對比，有失客觀性。本文采用多種不同的材料模型模擬混凝土/基巖膠結面，通過數值計算分析現場抗剪試驗中膠結面的剪切破壞過程及應力、變形狀態，進而對現場抗剪試驗得出的參數進行評價。

1 材料模型簡介

1.1 Mohr-Coulomb 材料模型

Mohr-Coulomb強度準則在巖土工程中的應用非常廣泛，大量的巖土工程設計計算都采用了Mohr-Coulomb強度準則[7]。

(1)屈服函數。Mohr-Coulomb材料模型的屈服軌跡如圖1和圖2所示，表達式為：

圖1 Mohr-Coulomb材料模型

圖2 Mohr-Coulomb材料模型在子午面和π平面上的屈服面

(1)

(2)流動勢函數。G為流動勢函數，表達式如下所示：

(2)

式中：ψ是剪脹角；c|0是初始黏聚力，即沒有塑性應變時的黏聚力；∈為子午面上的偏心率，它控制了G在子午面上的形狀與函數漸近線之間的相似度。Rmw則控制了其在 π 平面上的形狀。

1.2 節理材料模型

ABAQUS/Standard提供了一種節理材料模型[7]，該模型可以模擬其裂隙成組出現且有規律的節理巖體。

節理材料模型主要應用于材料受壓的情形，但也規定了節理的張開能力，當垂直節理面即將變為拉伸時，此時，假定材料裂隙的法向沒有彈性剛度，只有應變，由于節理面的張開，導致了節理面的各向異性彈性。

節理材料模型的塑性行為被單獨的分為兩部分，一部分是節理體系的塑性行為，另一部分是巖塊材料的塑性行為，兩者的屈服準則各不相同。

1.2.1 節理體系的塑性行為

節理體系a的滑動破壞面定義如下：

fa=τa-patanβa-da=0

(3)

式中：βa為節理體系的內摩擦角；da為節理體系的黏聚力，如圖3所示；τa為節理面上的剪應力；pa是作用在節理面上的法向壓應力。當fa<0時，節理體系不發生滑動。當fa=0時，節理體系滑動。此時，塑性應變為：

圖3 節理材料模型

(4)

(5)

節理材料模型對同一個巖體材料提供3組獨立的節理面，即在同一個點可以存在3組獨立的滑動節理面，而且每個節理體系的滑動不會改變其他體系的破壞準則和膨脹角。

圖4 巖塊材料模型

1.2.2 巖塊材料的塑性行為

巖塊材料的破壞服從線性的Drucker-Prager破壞準則，即：

q-ptanβb-db=0

(6)

當材料達到屈服時，巖塊材料的塑性流動定義為：

(7)

式中：gb=q-ptanψb是流動勢函數；ψb為散體材料的膨脹角，該巖塊破壞模式是擴展Drucker-Prager模型的簡化形式。

巖塊材料和節理體系的塑性流動都相互獨立，其塑性流動都不會相互影響。

1.3 基于ABAQUS的混凝土損傷塑性模型

ABAQUS軟件提供一種連續的、基于塑性的混凝土損傷塑性模型[7]，該模型假定混凝土材料主要因拉神開裂和壓縮破碎而破壞。

ABAQUS混凝土損傷塑性模型假定混凝土的單軸拉伸和壓縮性狀由損傷塑性描述，如圖5和圖6，并自動將用戶自定義的應力-“非彈性”應變數據轉換成σ-εpl關系曲線來描述混凝土的應力-應變關系。混凝土在應力-應變關系曲線的軟化段上卸載時，從圖中可以看出卸載的彈性剛度小于初始彈性剛度，這說明材料發生了損傷，用兩個損傷變量dt和dc表示單軸拉伸和壓縮下的應力-應變關系為：

圖6 混凝土單軸拉伸的損傷描述

(8)

(9)

(10)

(11)

(1)屈服函數。ABAQUS混凝土損傷塑性模型采用修正Lubliner-Fenves損傷模型的屈服面，如圖7，可以考慮拉伸和壓縮情況下不同的強度演化，并用有效應力來表達：

圖7 修正Lubliner-Fenves損傷模型在π平面上的屈服面

(12)

(2)流動勢函數。混凝土損傷塑性模型的流動勢函數采用非相關聯勢塑性流動。模型中所采用的流動勢G為Drucker-Prager拋物線函數，即：

(13)

2 混凝土/基巖膠結面現場抗剪試驗的數值模擬方法

2.1 有限元模型

根據《水利水電工程巖石試驗規程(SL264-2001)》[8]中對混凝土/巖體膠結面直剪試驗的描述，建立該試驗的二維有限元模型(圖8)：有限元模型中膠結面沿剪切方向的長度為50 cm，厚度設為1 cm，不考慮其起伏程度；混凝土試件推力面的中心線與剪切面成15°且過剪切面的中點，試件高40 cm；基巖上下游的長度和深度均為兩倍膠結面剪切方向長度，即100 cm。

圖8 二維有限元網格模型

二維有限元模型選用四邊形平面應變單元(CPE4R)，有單元2 208個、節點2 279個。

根據《巖基抗剪強度參數》[9]中對503組混凝土/基巖膠結面抗剪試驗的統計分析，以及《混凝土重力壩設計規范(SL319-2018)》[1]中對混凝土/基巖膠結面抗剪強度參數的建議值，取混凝土試件的標號為C20(老編號200號)，基巖的材料參數參考E2類巖體，具體取值見表1。其中《混凝土結構設計規范(GB50010-2010)》[10]對于C20等級的混凝土只提供了軸心抗拉強度和軸心抗壓強度，因此表1中的抗剪強度參數需要通過以下公式換算而得：

表1 各部分材料參數取值

(14)

(15)

2.2 荷載施加方法

《水利水電工程巖石試驗規程(SL264-2001)》[8]中對混凝土/巖體膠結面直剪試驗的試驗加載規定：每個試件的法向荷載分為3～5級施加；剪切荷載按照預估的最大值分8～10級施加，采用斜推法時，應同步降低因施加剪切載荷而產生的法向分量的增量，保持剪切面上的法向正應力不變。由于本文是為了研究剪切過程中膠結面的應力-應變狀態，而在剪切后期膠結面的破壞十分迅速，因此需要細分剪切荷載。本文的具體荷載步驟如下：對混凝土試件頂面施加法向荷載保證剪切面上的平均法向正應力分為4個等級，分別為1.0、2.0、3.0、4.0 MPa；對于每一級法向荷載，其對應的剪切荷載分17級施加，剪切前期推力增量較大，隨剪切的進行，增量減小，在增加推力的同時減小混凝土試件頂面的法向荷載，從而保持剪切面上的平均正應力為定值。

(16)

(17)

公式可簡化為

(18)

(19)

圖9 系統受力示意圖

2.3 本文名詞解釋及相關說明

為了方便后文對計算結果的分析，本文特作以下規定：

(1)平均正應力、平均剪應力：是指在假設試驗系統為剛體且膠結面上的應力均勻分布的情況下，由混凝土試件頂面的法向荷載和推力面的剪切荷載通過試件幾何尺寸關系換算出的膠結面上的法向正應力和切向剪應力，其大小與數值計算結果無關，在概念上等同于外荷載。

(2)法向正應力、切向剪應力：是指通過數值計算得出的膠結面上每個點的真實應力，其大小與數值計算結果相關。

(5)計算結果分析中的符號規定：法向正應力以壓為負、拉為正；切向剪應力以剪切方向為正；法向位移以張開為正、壓縮為負；切向位移以剪切方向為正。

(6)膠結面各節點及測點的位置說明(圖10)。

圖10 膠結面各節點及測點位置說明

根據《水利水電工程巖石試驗規程(SL264-2001)》[8]中對現場抗剪試驗裝置的說明，選取圖10中x=-20、x=0和x=20處的3對節點分別作為上游測點、中部測點和下游測點(每處有上下兩個節點，共同組成一個測點)。

(7)后文中的切向位移和法向位移均為測點位置處上下兩節點的相對位移。

3 基于多種材料模型的膠結面應力變形及破壞分析

3.1 膠結面破壞判別依據

對于數值模擬破壞的判別依據有很多，其中使用非常廣泛的穩定判別依據主要有以下幾種[11]：收斂性判據；屈服區貫通法；位移突變法；能量法。

現場抗剪試驗的加載過程類似于拱壩的超載試驗，都是在逐步增加外荷載的過程中尋找模型的破壞時刻，因此本文膠結面破壞的判別依據參考拱壩破壞的判別依據：以屈服區貫通法作為膠結面破壞的主要判別依據，膠結面測點的位移突變作為參考依據，其中將等效塑性應變大于10-5作為屈服標準。

3.2 計算結果分析

混凝土試件和基巖采用線彈性模型，膠結面分別采用節理材料模型(下文簡稱JOINT)、Mohr-Coulomb材料模型(下文簡稱M-C)、理想彈塑性的損傷塑性模型(下文簡稱CDP-S)和考慮軟化的損傷塑性模型(下文簡稱CDP)進行破壞分析。

表2 計算抗剪強度參數

圖11 平均正應力與平均剪應力散點圖

3.2.2 膠結面應力分布

圖12至圖19所示的膠結面應力分布曲線均在平均正應力為1 MPa時平均剪應力分別為4、13、15 MPa三個荷載步下獲得的。

圖12 JOINT不同平均剪應力下膠結面的法向正應力分布

圖13 JOINT不同平均剪應力下膠結面的切向剪應力分布

圖14 M-C不同平均剪應力下膠結面的法向正應力分布

圖15 M-C不同平均剪應力下膠結面的切向剪應力分布

圖16 CDP-S不同平均剪應力下膠結面的法向正應力分布

圖17 CDP-S不同平均剪應力下膠結面的切向剪應力分布

圖18 CDP不同平均剪應力下膠結面的法向正應力分布

圖19 CDP不同平均剪應力下膠結面的切向剪應力分布

當平均剪應力為4MPa時，各材料模型的膠結面仍處于彈性階段，不同材料模型膠結面的法向正應力分布和切向剪應力分布規律幾乎相同，膠結面大約80%的區域應力均勻分布，僅在膠結面首尾兩端稍微有些應力集中。

隨著平均剪應力逐漸增大，各材料模型的膠結面上的應力分布均越來越不均勻，兩端應力集中現象愈發明顯。隨著推力荷載的增加，JOINT膠結面法向正應力不變的區域其切向剪應力也不變，說明了JOINT破壞方向的唯一性；M-C和CDP-S法向正應力變小而切向剪應力變大，說明了該點處的破壞方向與剪切方向不一致，也體現了M-C和CDP-S各向同性的性質；CDP應力集中處的應力數值反而減小的規律體現了CDP材料應變軟化的性質。

3.2.3 剪力位移曲線

圖20和圖21分別是不同材料模型膠結面上、中、下3個測點的平均剪應力-切向位移關系曲線和平均剪應力-法向位移關系曲線。

從圖20中可得，不同材料模型在相同測點處的平均剪應力-切向位移關系曲線大致相同，所有測點的切向位移隨外荷載的增加先呈線性緩慢增加，說明在該階段該測點還未屈服，隨后切向位移的增長速率增大，切向位移快速發展。不同材料模型之間在平均剪應力-切向位移關系曲線上的差別僅體現在切向位移發展的速率上，容易發現，M-C和CDP-S測點的切向位移隨外荷載呈曲線光滑增長，而JOINT和CDP測點的切向位移曲線有十分明顯的拐點，切向位移近似成折線增長。

圖20 不同材料模型膠結面上、中、下3個測點的平均剪應力-切向位移關系曲線

從圖21中可得，不同材料模型下相同測點的平均剪應力-法向位移關系曲線大致相同。在曲線拐點之前(測點屈服之前)，測點的法向位移與膠結面平均剪應力近似呈線性關系，法向位移增長緩慢；超過拐點之后(測點屈服之后)，位移發展迅速，CDP的拐點最為突出。

圖21 不同材料模型膠結面上、中、下3個測點的平均剪應力-法向位移關系曲線

不同測點的平均剪應力-切向位移曲線和平均剪應力-法向位移曲線也不完全相同，上游測點的位移(切向位移和法向位移)發展相對更加平緩且在相同的外荷載下位移更大。表3是不同測點平均剪應力-位移關系曲線拐點所對應的平均剪應力，從表3也可以看出，除CDP外，其他材料模型的上游測點拐點所對應的平均剪應力最小，下游測點拐點所對應的平均剪應力最大，體現了不同測點的屈服并不同步，也說明了膠結面的屈服區從上游向下游發展的規律。CDP材料模型各測點的拐點對應的平均剪應力相同，超過該拐點，位移迅速增加，材料應變軟化的性質導致了膠結面屈服區和位移的突變。

表3 不同測點平均剪應力-位移關系曲線拐點所對應的平均剪應力 MPa

圖22給出了重慶魚嘴長江大橋錨碇混凝土/基巖膠結面抗剪強度試驗τN1組應力與位移關系曲線[12]。將數值計算結果(圖20，圖21)與實際工程試驗結果(圖22)相比，兩者的切向位移發展規律和法向位移發展規律大致相同，說明數值計算結果比較合理。

圖22 重慶魚嘴長江大橋錨碇混凝土/基巖膠結面抗剪強度試驗τN1組應力與位移關系圖

3.2.4 膠結面破壞分析

圖23至圖27為平均正應力為1 MPa時的計算結果。圖23為不同材料模型膠結面屈服區比例隨外荷載(平均剪應力)增長的發展曲線。圖24-圖27為不同材料模型膠結面的屈服區發展云圖(紅色為屈服區，藍色為未屈服區)。

從圖23可得，M-C、CDP-S和CDP膠結面屈服區比例的發展曲線大致相同，當平均剪應力為0.64 MPa時，膠結面開始出現破壞，隨外荷載的增加，屈服區發展相對平緩，最終在平均剪應力為2.5 MPa屈服區比例達到100%，其中CDP膠結面屈服區比例在破壞后期發生突變，從72%直接增長到100%；JOINT膠結面屈服區比例的發展曲線較其他材料模型膠結面屈服區比例的發展曲線差別較大，在平均剪應力為1.28 MPa時，膠結面才開始發生破壞，且曲線存在拐點，平均剪應力超過該點后，屈服區快速發展，最終也在平均剪應力為2.5 MPa左右達到100%。

圖23 膠結面屈服區比例

從圖24至圖27可得，不同材料模型下的膠結面屈服區的發展規律大致相同，均從上游逐步向下游發展，直至貫通整個膠結面。

圖24 JOINT膠結面的屈服區發展云圖

圖25 M-C膠結面的屈服區發展云圖

圖26 CDP-S膠結面的屈服區發展云圖

圖27 CDP膠結面的屈服區發展云圖

4 結 論

本文通過有限元法模擬混凝土/基巖膠結面的現場抗剪試驗，分析了不同材料模型下膠結面的應力、位移及屈服區的分布及發展規律。研究結果表明：

(2)通過不同材料模型計算得出的膠結面應力分布規律大致相同，膠結面應力均呈現不均勻分布狀態，上下游端部應力集中明顯。

(3)通過不同材料模型計算得出的膠結面屈服區發展規律大致相同，隨外荷載增加，膠結面屈服區均從上游逐漸向下游發展，直至貫通整個膠結面。

(4)《水利水電工程巖石試驗規程(SL264-2001)》[8]中規定通過繪制測點在各法向應力下的剪應力與剪切位移及法向位移關系曲線來確定混凝土/基巖膠結面在各法向應力下的抗剪斷強度峰值，這種方法假設了測點處的破壞代表了整個膠結面的破壞，即假設了膠結面的應力分布為均勻分布。由于膠結面的應力分布并不均勻，所以實際上測點處的破壞并不能代表整體破壞，有可能膠結面并沒有完全破壞，因此筆者認為通過平均應力整理出來的抗剪強度是偏保守的。

□