鄭達, 周賽, 周鴻軻
(成都理工大學環境與土木工程學院, 成都 610059)
隨著中國道路工程建設逐步向西部高山峽谷地區深入,大型橋梁的橋位選址不可避免的會遇到不良地質體問題。例如,一些橋梁線位處存在的第四系堆積體,其承載力與穩定性問題就直接制約了橋墩與重力錨等工程設施的布置。正確認識這些堆積體的成因機制,準確獲取堆積體的強度參數是合理評價橋墩等工程建設穩定性的關鍵。
目前,對于山區大型堆積體工程特性的研究主要還是依賴于常規的室內外試驗手段,研究人員通過室外原位試驗可以掌握堆積體的變形特性與強度特性。但由于山區試驗場地的限制,此類試驗的條件常常無法滿足。許多堆積體的強度參數與力學機制還需要通過室內大型直剪試驗與三軸剪切試驗獲取。先有學者通過大型直剪試驗系統地研究了含石量對試樣強度參數的影響,并擬合多種物理、數學模型加以描述[1-2]。后因大型直剪試驗存在荷載偏心、破壞面固定等問題,更多的學者開始使用大型三軸試驗展開相關研究:如基于含石量指標,發現含石量35%的土樣在加載高圍壓時,試樣整體變形表現出剪縮-剪脹-二次剪縮-二次剪脹的特征,并建立表述土樣力學特征的強度準則[3-5];或是在試驗基礎上,總結應力與抗剪強度參數關系,并擬合出試驗關系式等[6]。以上研究均獲得了較為卓越的成果。然而近年來,考慮到室內試驗樣本數量有限,不能全面認識堆積體在外部加載等作用下的變形特性,部分學者開始依靠計算機數值仿真技術模擬堆積體材料的剪切變形過程,剖析此過程中堆積體內部的微細觀變化與能量分布。目前已有大部分學者使用顆粒流軟件PFC(particle flow code)的二維、三維模塊進行三軸試驗的仿真模擬[7-9],從細觀尺度對室內試驗進行補足或是拓展,發現在合適的參數標定下,室內試驗與數值模擬有相當高的吻合度。文獻[10-11]通過數值模擬進行邊坡穩定性與地質災害演化的分析,獲得不錯的成效。可見,將室內外試驗成果與數值分析的結論相結合是認識堆積體工程特性的有效途徑。
以往學者在對三軸試驗進行仿真模擬計算時,多用剛性邊界進行模擬,考慮到室內三軸試驗試樣外層為橡皮膜,在試驗過程中能較自由的發生變形。而剛性墻施加圍壓并不能很好地體現三軸試樣的實際變形情況,在反映三軸試驗的特性上也會有一定偏差[12]。鑒于此,本次數值模擬采用PFC-FLAC(fast Lagrangian analysis of continua)耦合方法,建立柔性邊界條件下的三軸試驗模型,可以更好地體現堆積體試樣在壓縮過程中的變形特征。
以擬建的西香高速公路瀘沽湖特大橋橋位處的碎石土堆積體為地質原型。針對該碎石土堆積體的強度特性,擬通過室內大型三軸試驗測得該處碎石土堆積體的抗剪強度參數,探究不同粒料含量堆積體的強度參數的變化規律,進而結合PFC-FLAC耦合數值仿真技術,建立柔性邊界伺服模型,獲得碎石土內部剪切帶隨圍壓變化的區域化特征,在此基礎上,研究把控研究區碎石土試樣的變形破壞特征與力學機制,為該處橋梁與類似工程的選址、施工方案提供一定依據。
堆積體位于滇西北橫斷山高山峽谷與云貴高原接壤地帶東側、四川盆地西南,地勢高差懸殊,屬構造剝蝕、侵蝕深切高山地貌。場區發育有多處斷裂帶,西以八二橋-金棉斷裂為界,與新村-長柏斷裂帶相鄰,東以金河-箐河斷裂為界,與中生代楚雄前陸盆地三級構造單元相鄰。該區域在三疊紀發生凹陷,沉積了一套厚約5 000 m以上的地層。

圖1 堆積體工程地質剖面圖Fig.1 Engineering geological profile of accumulation body
根據現場地質調查成果,堆積體斜坡巖層產狀為N28°~60°E/NW∠20°~30°,堆積體長800 m、寬300 m,總方量約2.4×106m3。該堆積體后緣呈內凹型,中前緣外突,由于擬建主墩位于該堆積體邊緣,覆蓋層厚度大,對擬建大橋的潛在影響較大。因此對該碎石土堆積體的變形破壞特征與力學特性進行研究分析有一定的必要性。
本次試驗所使用的材料為堆積體原狀碎石土樣,其天然容重為22.68 kN/m3,密度為2.37 g/cm3。為減小原狀土樣顆粒大小對試驗的影響,試料粒徑須控制在60 mm以下,因此需對原狀土樣進行縮尺[13-14]。試樣超徑顆粒含量小于40%,直接使用等量替代法縮制試樣,便能得到與原料性狀差異較小的試樣。等量替代法計算公式為
(1)
式(1)中:Pi為等量替代后i粒組的百分含量;Poi為原級配某粒組的百分含量;P5為粒徑大于5 mm的百分含量;Pdmax為超粒徑的百分含量。
根據式(1)可確定重塑土顆粒級配參數,由此獲得其試驗級配曲線,如圖2所示。

圖2 級配曲線Fig.2 Gradation curve
試驗所用儀器為圖3所示的英國大型粗粒土動靜三軸試驗系統(geotechnical data system,GDS),主要包括控制試驗軟件和三軸試驗的硬件。該儀器為全自動試驗儀器,可在高圍壓條件下對粗粒土進行三軸壓縮、動三軸等試驗。

圖3 大型三軸試驗儀示意圖Fig.3 Diagram of large-scale triaxial test instrument
根據試驗條件,試驗采用等側壓常規三軸試驗。選定圍壓100、200、300、400 kPa,試樣制備成φ300×600 mm徑高比為1∶2的標準試件。本次試樣取最大粒徑為60 mm。通過控制軸向位移進行試驗,加載速率為3.0 mm/min。
對試驗數據進行整理,得到碎石土試樣的應力-應變曲線(圖4)。可以看出,試樣在初始變形期(軸向應變0~3%)存在一定的彈性增長階段,應力-應變圖中表現為直線,此線性階段的斜率為碎石土的初始彈性模量Ei。具體數值如表1所示。

表1 不同圍壓條件下的初始彈性模量Table 1 Initial elastic modulus under different confining pressures

σ1為軸向應力;σ3為圍壓
根據表1可知,Ei與圍壓呈正相關關系。初始彈性模量與應力應變關系密切相關,而圍壓影響應力應變關系最直接的因素在于試樣的壓密程度,即圍壓通過影響試樣的密度從而影響初始彈性模量Ei。分析其原因,在圍壓較大的情況下,試樣的壓密程度更高,此時試樣內部顆粒間的咬合作用較強,顆粒難以產生錯動的同時,試樣的抗變形能力也得到增強,因此表現出較大的彈性模量。
按照《土工試驗方法標準》(GB/T 50123—2019)中粗粒土三軸試驗強度取值規定獲取試驗中各圍壓下的破壞強度σcre。由此可得到該碎石土試樣的強度包絡線(圖5),對包絡線擬合切線可得到抗剪強度指標內聚力c為48.36 kPa,內摩擦角φ為31.42°。對于碎石土的抗剪強度參數,與顆粒級配有較大關系,文獻[10,15]在粗顆粒料占比50%~70%基礎上測得c為60~99.7 kPa,φ為20.8°~36.65°。本次試驗粗粒粒組的質量分數大于80%,相較以往學者內聚力偏小而內摩擦角偏大的原因在于,試驗縮尺后粗顆粒料較多,顆粒之間閉鎖、咬合等作用較強,而細顆粒料間的接觸相對減少,降低了試樣的黏性。因此得到結論:隨著粗顆粒料增加,碎石土的內聚力會變小而內摩擦角變大。

圖5 碎石土試樣莫爾包絡線Fig.5 Mohr envelope of gravel soil specimen
由破壞強度與圍壓關系曲線(圖6)可以看出,該處堆積體所取碎石土的破壞強度與圍壓的關系曲線可以用一次或二次函數曲線方程進行擬合。比較發現,一次函數相關性系數更高,更便于在實際工程中的使用。因此,所用碎石土試樣的破壞強度與圍壓關系可近似表示為

圖6 破壞強度與圍壓關系曲線Fig.6 Relation curve of failure strength and confining pressure
σcre=2.95σ3+176.105,R2=0.995
(2)
考慮到僅憑室內試驗,缺乏對堆積體在加載過程中變形破壞特性,以及力鏈、能量等微細觀方面的認識,因此使用數值模擬軟件對室內試驗進行仿真模擬計算,對以上缺陷進行補足。并且,為了反映堆積體土樣在壓縮過程中的整體變形特征,本次數值模擬通過PFC-FLAC耦合,使用柔性邊界代替剛性側墻,在此基礎上建立模型模擬不同圍壓條件下的大型三軸試驗。
目前PFC-FLAC耦合模擬大型三軸試驗是通過邊界控制墻體實現[16]。這種情況下,在PFC模型的墻(wall)單元外側接入FLAC實體單元殼(shell),基于此,墻單元可以一同傳遞作用于球(ball)單元與殼單元間的力。在球單元運動過程中,作用于墻上的接觸力等效分配到墻單元的頂點,而后參與連續域的FLAC分析。同樣,連續域中殼單元的變形也以相同的方法傳遞給墻單元,進而將位置、速度等信息同步到離散域中的球單元[17]。這些值的不斷傳遞與更新,加之墻單元的頂點附著于殼單元的網格點上,在墻頂點與殼單元網格點同步運動的基礎上,給予柔性邊界不同于剛性墻邊界條件下的幾何參數變化。
在數值仿真計算過程中,模型的建立是第一步,尤為關鍵。本次基于PFC-FLAC耦合的模型建立主要考慮以下幾個方面。
2.2.1 顆粒級配
顆粒級配是模型建立成功與否的關鍵因素。若直接使用室內篩分試驗或縮尺后所獲顆粒級配進行顆粒的生成,粒徑小于1 mm的顆粒雖然在顆粒總量中占比較少,但因其顆粒質量較小,在計算過程中將產生極大的數量,會對模型的建立產生不利的影響,因此需重新分配粒徑以優化計算結果。為解決粒徑問題,同時保證所建模型的準確性,參考文獻[18],使用加權平均顆粒生成法,通過該方法設置相關顆粒級配,計算公式為
(3)
(4)
式中:Ri為1.2節中顆粒半徑;Mi為各粒徑范圍的質量占比;ra為最大粒徑Rmax與最小粒徑Rmin之比;Wr為加權平均半徑。
由1.2節中顆粒級配曲線,結合上式計算得到,Wr=6.012 5,最小顆粒半徑取1 mm,最大粒徑與最小粒徑之比為7.15。
在模型建立過程中,顆粒的生成有表面光滑的ball顆粒與表面凹凸不平的clump顆粒兩種形式。對于碎石土來說,在PFC中可以設定較大的摩擦系數使得模擬結果與室內試驗結果相吻合[19],故選用ball生成顆粒方法。
2.2.2 接觸模型
在離散元中,接觸模型是顆粒間、顆粒與墻之間相互作用的物理準則,決定了整個顆粒介質的宏觀、細觀的力學特性。在本次數值模擬過程中主要有顆粒-顆粒(ball-ball)、顆粒-三角形墻面(ball-fa-cet)兩種接觸類型,使用線性接觸模型模擬其間的接觸[20-21]。
2.2.3 柔性膜
在三軸試驗中,試樣使用橡膠膜包裹住以承受圍壓,上下底面用剛性板承受并傳遞荷載。而在數值模擬中,墻可以用來代替剛性加載板,但無法像橡膠膜一樣發生彈性變形,反而會給試樣側向約束,使得試樣無法在側向發生不均勻變形,從而影響顆粒與墻體的接觸,與實際不符。基于上述問題,引入柔性膜方法,側墻使用FLAC中可以發生彈性變形的殼單元,利用PFC與FLAC耦合方法創建側向的柔性邊界。
綜合以上考慮,最終確定模型的搭建是通過PFC-FLAC耦合方法建立柔性伺服邊界條件,內部顆粒粒徑在1~7.15 mm,并以線性接觸模型模擬顆粒、顆粒與墻之間的接觸。模型如圖7所示。

圖7 離散元模型Fig.7 Discrete element model
在三軸試驗中,本次使用的GDS試驗系統可直接測得偏應力σ、應變ε等數據。而在數值模擬中,則需通過公式推導,并使用history命令進行監測獲得數據,計算公式為
(5)
(6)
式中:σ為偏應力;σ1為軸向應力;σ3為圍壓;P1和P3分別為上下墻面承受的荷載;r為試樣的截面半徑;ε為應變;h0和h分別為試樣加載前后的高度。
基于監測結果,圖8給出了試驗與模擬所得應力應變曲線對比,兩者吻合度較好,應力應變曲線變化趨勢大體一致。其次,模擬試樣內部顆粒的位移矢量體現出了與試驗相同的剪脹特性(圖9)。可見,數值模擬在本次三軸試驗的應變范圍與圍壓范圍內很好地模擬了試樣本身的力學特征。

圖8 室內試驗與數值仿真應力-應變曲線對比Fig.8 Comparison of stress-strain curves between indoor test and numerical simulation

圖9 試樣剪脹與模擬位移矢量圖Fig.9 Dilation of specimen and displacement vector of numerical simulation
為展示試樣的整體變形特征,圖10列出了不同圍壓條件模型加載結束后柔性邊界上殼單元xoz與yoz平面的位移云圖。對比發現,基于PFC-FLAC耦合的柔性邊界模擬,模型很好地體現了在三軸壓縮過程中試樣的不均勻變形特征。其中主要變形區出現在試樣中部,且不同圍壓條件下試樣的整體形態也略有不同。

圖10 柔性邊界位移云圖Fig.10 Displacement cloud atlas of flexible boundary
顆粒類材料發生破壞時往往會形成剪切帶,且剪切帶內部顆粒的歐拉角會發生明顯的變化,外部則鮮有變化,即歐拉角的變化呈一定的區域化特征[22-23]。于是可以通過觀察模擬試件剖面上的歐拉角變化情況分析剪切帶的形成。在室內三軸試驗中,難以獲知試樣內部顆粒的變化情況,而通過離散元數值模擬計算,很好的補足了室內試驗無法可視化分析的缺陷。
由各圍壓條件下試樣內部歐拉角的變化情況(圖11)可知,隨著圍壓的變化,試樣內部歐拉角的變化也有不同:當試樣加載100 kPa低圍壓時,可發現試樣剖面形成明顯的“X”形剪切帶,兩端的三角形區域為彈性變形區,而中部為膨脹變形區,區域化特征較為明顯;對于加載200 kPa的試樣,僅有一條斜向剪切帶較為明顯;圍壓為300、400 kPa的試樣并未顯示有規則的剪切帶,區域化特征不明顯。

圖11 試樣內部顆粒歐拉角變化Fig.11 Eulerian angle variation diagram of the particles inside the specimen
結合400 kPa圍壓下試樣內部顆粒歐拉角隨應變變化情況(圖12),可以看出,在10%應變時試樣的歐拉角變化有一定區域化特征,內部出現有剪切帶,而在應變達到16%時剪切帶外部的歐拉角也出現變化,歐拉角的變化呈紊亂-規則-擴散趨勢。究其原因在于:當變形達到一定程度時,僅靠剪切帶內的力鏈已不足以承受高圍壓條件,剪切帶外的力鏈需要承受更多荷載,顆粒歐拉角也隨之變化,歐拉角的變化整體呈擴散趨勢。

圖12 顆粒歐拉角-應變變化情況Fig.12 Changes of Euler angle-strain
顆粒材料是由眾多離散單元組成的非連續介質,相鄰顆粒間以某種線性的鏈狀結構傳遞應力,該鏈狀結構稱為力鏈。在重力與外部荷載的作用下,力鏈間會交叉、融合,產生力鏈網絡,其復雜的空間分布與應力狀態決定了非連續性試樣的力學性能。
力鏈可分為強接觸力鏈和弱接觸力鏈兩種。其中,強接觸力鏈承受并傳遞大部分荷載,但其承受荷載存在上限,當傳遞的荷載大于力鏈的限值時,接觸力鏈就會斷裂;而弱接觸力鏈只承擔小部分荷載,并在強接觸力鏈斷裂后參與其重組[24]。因此可通過接觸力鏈的粗細來判斷傳遞的力的大小。
為觀察在歐拉角變化時力鏈的分布情況,圖13列出400 kPa圍壓條件下力鏈分布情況,取應變為10%與16%時的力鏈分布情況進行對比分析。對試樣的力鏈局部放大發現:試樣在達到應變10%時,內部力鏈指向在試件兩端呈豎直向而中部呈傾斜向,強接觸力鏈的分布整體呈“X”形,力鏈代表著力傳遞的大小與方向,應變10%時力鏈方向與剪切帶方向相似,正是此時的力鏈引導著剪切帶的形成;16%應變率時,強接觸力鏈明顯增多,分布區域并無規則,呈擴散趨勢,說明此時有更多的力鏈來承受荷載,與前文結論基本相似。

圖13 力鏈分布Fig.13 Distribution graph of force chain
試件的三軸壓縮過程也可解釋為能量變化的過程,能量變化是探究材料力學機制的重要手段[25]。在實驗過程中能量的變化主要包括被試樣系統吸收的總能量,記為WU,分為軸向和徑向,軸向輸入的能量,記為Wz;徑向可能是剪縮或剪脹,因此徑向能量存在正負,記為Wr;而摩擦能量消耗與溫度帶來的能量損耗均較小,可忽略不計。則外力所做功W均可看作被試樣吸收,即
W=Wz+Wr≈WU
(7)
假設在應力狀態σi,j下所產生的應力增量為dεi,j,外力做的功δW可表示為
δW=σi,jdεi,jdxdydz
(8)
對單位體積試件所做的功δw為
(9)

(10)
式(10)中:σz、σr分別為軸向、徑向應力;εz、εr分別為軸向、徑向應變。
根據上述推導能量計算公式,計算出試樣在各圍壓條件下單位體積各項能量變化,如圖14所示。

徑向本是釋放能量,為負值,為便于展示取其相反數
對比4組曲線(圖14)可以看出,不同圍壓下試樣能量變化大體趨勢相近。隨著軸向應變的增加,試樣總吸收能量逐漸增加;相同軸應變下,圍壓越大試樣各個方向吸收、釋放的能量越大。在100、200 kPa圍壓初始階段出現短暫的總能量為負的情況,應是低圍壓條件初始做功較少,忽略摩擦等能量導致的誤差。
由于能量守恒,試樣吸收的總能量不會憑空耗散,而是內部消耗,影響顆粒間力鏈、歐拉角的變化、剪切帶的形成等。結合應力-應變曲線,從以下3個階段分析能量的變化。
(1)彈性階段。試樣被壓縮,發生顆粒間彈性變形,能量主要轉化為彈性應變能。
(2)屈服階段。可以發現隨著圍壓增加,屈服階段隨之變長,在此階段吸收的能量也越多。此時能量主要用于試件內部微裂隙的發育、演化、滑移和強度劣化,直接導致試件強度降低。顆粒內部歐拉角發生變化,并有從紊亂變得規則的趨勢,剪切帶逐漸形成。
(3)變形破壞階段。在低圍壓條件下(100、200 kPa),試件內部形成完整的剪切帶,能量主要作用于剪切帶內部顆粒歐拉角的變化以及部分力鏈的斷裂,上述變化發生在剪切帶內。而在較高圍壓條件下(300、400 kPa),較低圍壓有更多能量需要釋放,多余能量向剪切帶外擴散,持續發育新的裂隙,對應前文歐拉角的擴散變化。
基于大型三軸試驗與PFC-FLAC耦合的柔性邊界數值模擬,綜合分析西香高速公路瀘沽湖特大橋橋位處碎石土堆積體的宏觀變形特征、微細觀力學機制,得到如下結論。
(1)圍壓通過影響試樣的壓密程度,與初始彈性模量呈現正相關的關系。
(2)碎石土的抗剪強度參數,體現出隨著粗顆粒料增加,內聚力減少而內摩擦角增大的規律。
(3)擬合出線性方程以描述該處碎石土試樣破壞強度與圍壓的關系:σcre=2.95σ3+176.105。
(4)從變形角度看,基于PFC-FLAC耦合加載,低圍壓條件下碎石土試樣內部出現呈“X”形的兩條剪切帶;由于試樣內部紊亂的強力鏈作用,剪切帶在圍壓達到200 kPa時變為一條,并隨著圍壓的增加而逐漸消失。
(5)從能量角度看,低圍壓下,試樣內部上下錐形區域持續壓密吸收能量,而在中部釋放能量導致了膨脹變形;而隨著圍壓的增加,圍壓做了更多的功造成能量冗余,多出的能量向剪切帶外擴散,產生新的裂隙,導致剪切帶逐漸減少直至消失。