剛性與柔性邊界下砂離散元三軸試驗力學特性對比分析

陳愛平，劉凱文，邱睿哲，熊志鵬，蔡志清

(1.成都工業職業技術學院軌道交通學院,成都 610218; 2.西南交通大學高速鐵路線路工程教育部重點實驗室,成都 610031;3.西南交通大學土木工程學院,成都 610031; 4.貴陽中天巖土工程有限公司,貴陽 550001)

引言

在三軸試驗中,對具有代表性的土體材料試樣通過控制沿試樣邊界的力/壓力或位移,可獲得試樣相關的力學特性參數。其中,材料的位移變化響應與應力大小有關;而材料膨脹或收縮的趨勢則由其密度和應力水平決定[1]。目前,已有學者在室內三軸試驗方面進行了土體的特性研究[2-5]。然而,受限于目前已有的觀測技術水平,難以在細觀尺度上觀測土體試樣內部的位移變形特性[6]。

離散元方法(discrete element method,DEM) 不僅可以節約物力、人力的成本,還能多次重復同一模型試驗并克服傳統的宏觀連續性假設。大量學者開展基于離散元的數值模擬研究[7-10]。顆粒流模型細觀參數的確定是離散元數值模擬研究的關鍵環節,因此基于顆粒流的三軸試驗研究尤為重要[11]。CHENG等[12]較早地提出將DEM應用于三軸試驗的模擬分析中,發現這種數值模型的塑性行為與真實的砂土非常相似。隨后,一批學者開始選擇采用顆粒流軟件建立三軸模型對砂土[13-14]、黏土[15-16]的細觀參數進行分析。結果表明,顆粒粒徑大小和孔隙率是影響砂土壓縮特性的重要參數,孔隙率、黏結強度和摩擦系數等是影響黏土宏觀力學性質的重要參數,并基于這些結論建立了顆粒細觀參數與宏觀力學指標之間關系[17]。盡管這些研究在幫助理解土體材料細觀參數與宏觀力學特性關聯機制上取得進展,但其結果仍與真實三軸試驗存在差異。因此,如何采用離散元方法更為合理地再現三軸試驗成為重要話題。當前相關研究主要采用剛性邊界墻作為伺服邊界,并結合邊界控制系統來控制應力[18],這種方法可以簡單有效地用于控制沿邊界的整體應力狀態。然而,CHEUNG等[19]認為,如果模擬一個局部發展的三軸壓縮試驗,這些剛性邊界將抑制這個局部剪切帶的自然發展,并且沿邊界施加的應力將存在顯著的不均勻性。因此,其提出了使用一種“應力控制膜”來模擬包裹樣品的柔性乳膠膜伺服模擬方法。

目前采用圍壓柔性加載的三軸試驗離散元模擬還較少,不同邊界加載圍壓下三軸砂土試樣的力學特性、變形特性以及變形破壞規律也尚未探明。為此,利用PFC6.0程序開展不同密實度的HST95標準砂三軸模型試驗模擬,并結合試驗數據討論不同圍壓加載方式下DEM三軸試驗模型合理性,從顆粒細觀尺度上對比分析了不同邊界加載方式下砂土試樣剪切帶的形成、發展過程。

1 模型基本參數

1.1 砂土顆粒生成

本文采用的DEM模擬模型中砂土是基于英國鄧迪大學一系列巖土試驗采用的HST95標準砂土[20],該砂土基本物理參數如表1所示。

表1 HST95砂土物理參數Table 1 Physical parameters of HST95 sand

為使DEM中生成的顆粒與砂土實際級配一致,依據土體級配曲線將顆粒劃分為多組,后采用隨機分布顆粒生成法在區域范圍內隨機生成目標孔隙率的顆粒,得到顆粒級配如圖1所示。由圖1可知,其級配分布情況與實際情況基本吻合。

圖1 離散元顆粒級配與砂土級配Fig.1 Gradation of discrete element particle and sand

1.2 DEM模型建立

采用如下方法建立三維HST95砂土DEM模型,模型砂土體采用ball球體單元進行模擬,該方法主要包括以下步驟。

(1)土體成樣。于DEM模型尺寸邊界生成剛性墻體,在設定的剛性墻體內采用隨機分布顆粒生成法生成目標孔隙率的級配顆粒,然后將接觸模型及各類參數賦值在顆粒與墻體間,再通過循環命令使得顆粒之間逐漸降低不平衡力,使顆粒穩定。

(2)邊界模擬。對平衡完畢的試樣顆粒進行邊界伺服,對于剛性邊界條件下:采用循環伺服機制對側向墻體的移動速度進行調整,以達到目標所需圍壓值。對于柔性邊界條件下:將側向剛性墻替換為由球體顆粒組成的“柔性墻體”,再通過邊界命令來控制調整邊界球體顆粒上的集中力,從而實現與實際物理三軸試驗相同的柔性圍壓。

(3)模型加載。清除上下邊界墻體的邊界狀態,保持側向圍壓的邊界狀態,對上下墻體施加一個恒定壓縮速度進行加載。隨加載過程記錄試樣的應力-應變曲線、孔隙率、體應變和位移變化云圖等宏、細觀參數,以供后續進行參數分析。

圖2為剛性和柔性邊界條件下顆粒試樣模型邊界狀態時的剖面,即邊墻與顆粒的分布情況。

圖2 不同邊界條件下模型剖面Fig.2 Model section under various boundary conditions

2 DEM加載模擬原理

2.1 柔性邊界模擬

在DEM三軸試驗模擬中,剛性邊界模擬通常只需建立剛性墻體即可快速實現。而為討論柔性邊界對三軸試驗仿真模擬的影響,文中柔性邊界采用大小相同的球體顆粒組成球體膜顆粒來模擬三軸物理試驗中的柔性膜。模型中柔性膜顆粒是由每層相互平行排列的顆粒圓組成,顆粒間彼此接觸黏結,相鄰顆粒間形心距離為L(兩倍的膜顆粒半徑r)。因三軸試驗中橡皮膜的特殊柔軟力學特性,故將膜顆粒間設置為不傳遞力矩的接觸-黏結模型進行模擬。具體模型和接觸示意如圖3所示。

圖3 柔性膜顆粒示意Fig.3 Schematic diagram of flexible membrane particles

2.2 圍壓邊界原理

在剛性邊界三軸試驗模擬過程中,圍壓通過剛性墻進行加載,限制試樣發生同步側向變形。在柔性邊界三軸試驗模擬過程中,圍壓則需通過等效集中力的方式施加在球體膜顆粒上,從而無需強迫發生同步側向變形,在每一次邊界循環過程中進行集中力的計算判斷,根據膜顆粒的相對位置不斷調整至所有試樣球體顆粒-膜顆粒平均應力達到目標值。對于柔性膜顆粒而言,任意球體顆粒皆與周圍顆粒接觸在一起。圖4為膜顆粒集中力的計算示意。

圖4 膜顆粒等效力計算示意Fig.4 Calculation of equivalent force for membrane particles

對于任意一個膜顆粒,施加在該顆粒的等效集中力計算公式如下

F=Aσ

(1)

式中,F為集中力;A為等效面積;σ為控制圍壓。其中,等效面積的計算公式如下

A=l1×l2

(2)

式中,l1,l2分別為與球體相切四邊形相鄰兩邊長的向量。

2.3 試樣顆粒接觸模型

根據上述過程建立三維DEM三軸試驗模型,對土體顆粒進行參數賦值,參數接觸采用赫茲非線性模型。該接觸模型適用于光滑的彈性散體材料,也是DEM軟件中適用于分析砂土顆粒特性的接觸模型。赫茲接觸模型基于Mindlin和Deresiewicz(1953)[21]的近似理論得到的非線性公式組成。赫茲接觸模型產生法向力和剪切力是基于對光滑彈性球體在摩擦接觸中的變形理論分析。

赫茲模型中,接觸力可分解為赫茲分量和阻尼分量,赫茲分量部分提供線彈性行為,阻尼分量部分提供黏性行為。接觸力的計算關系公式為

Fc=Fh+Fd

(3)

式中,Fh為非線性赫茲力;Fd為阻尼力。

本文DEM赫茲模型具體細觀參數以Micheal(2019)模型所標定出的接觸參數為標準[22],如表2所示,其主要包含的參數有:剪切模量S,摩擦系數fc,泊松比υ。

表2 DEM模型接觸參數Table 2 Contact parameters of DEM model

3 DEM試驗模擬分析

利用前文所述圍壓加載原理,分別對不同邊界方式和密實度的HST95砂土建立60 kPa圍壓下三軸DEM模型,與文獻[22]中試驗結果數據對比分析不同邊界加載方式的合理性,進而對HST95砂土試樣力學特性、變形特性和剪切過程進行分析。

3.1 力學特性分析

圖5給出了不同邊界條件下模擬HST95砂土剪切破壞過程的應力應變曲線,分析圖5可得到如下結論。

圖5 不同密實度下應力應變曲線Fig.5 Stress-strain curves under different relative densities

(1)相較于柔性邊界,剛性邊界情況下應力曲線增長速度更快,應力峰值也更大。

(2)在4%左右應變對應的峰值應力上,30%密實度剛性邊界、柔性邊界的峰值應力為180.53 kPa和160.06 kPa,分別為試驗數據的115.4%和102.6%;70%密實度剛性邊界、柔性邊界的峰值應力為255.26 kPa和236.92 kPa,分別為試驗數據的107.8%和100.5%。

(3)在20%應變對應的殘余應力數值上,30%密實度剛性邊界、柔性邊界的殘余應力為148.15 kPa和134.15 kPa,分別為試驗數據的110.9%和101.6%;70%密實度剛性邊界、柔性邊界的殘余應力為220.55 kPa和154.52 kPa,分別為試驗數據的153%和107.6%,剛性邊界條件下的殘余應力較大。由此表明,該DEM仿真模型能有效地模擬還原HST95砂土的室內三軸試驗,柔性邊界模型應力應變的發展過程與砂土實際的壓縮過程更為接近,其能夠更為有效地模擬出室內試驗的效果。

3.2 變形特性分析

為準確分析三軸試樣的變形特性,分別在柔性和剛性邊界三軸模型中設置高度值為0.2、0.3、0.5、0.6和0.7倍試樣高度,直徑為0.9倍試樣寬的5個測量球,如圖6所示。通過測量球監測試樣顆粒孔隙率的變化,分析不同邊界條件下試樣的變形特性。

圖6 測量球分布示意(柔性左,剛性右)Fig.6 Distribution diagram of measurement ball (Left: flexible; Right: rigid)

圖7給出了不同邊界條件下三維模型的孔隙率變化曲線對比情況,可以得出如下結論。(1)在30%密實度條件下,孔隙率變化呈現逐漸減小至穩定的趨勢。在數值上,20%應變時剛性邊界和柔性邊界的孔隙率分別為0.381和0.383。(2)在70%密實度條件下,不同邊界條件對孔隙率變化的影響差別較大。剛性邊界孔隙率變化呈現“先小幅度減小,再逐漸增大,后再次減小”的變化規律;柔性邊界孔隙率變化與密砂的體應變變化規律一致,為“先剪縮,后剪脹”的特點。在數值上,20%應變時剛性邊界和柔性邊界的孔隙率分別為0.363和0.369。

圖7 孔隙率變化對比Fig.7 Comparison of porosity changes

結合孔隙率對比情況,進一步分析試樣的體積變形特性,圖8給出了三維模型的體積應變曲線與模型試驗對比情況,可以得出如下結論。(1)變形趨勢上,柔性邊界模型較好地體現出了“松砂剪縮,密砂剪脹”的砂土剪切變形規律;而剛性邊界在70%密實度情況下與試驗數據存在較大差異。(2)在數值上,20%剪切應變時30%密實度剛性邊界和柔性邊界的體應變為-6.4%和-6.2%,分別為試驗數據的106.7%和100.6%;70%密實度剛性邊界和柔性邊界的體應變為1.9%和2.9%,分別為試驗數據的67.9%和100.3%。(3)無論在何種密實度工況下,柔性邊界模型的體積應變都與模型試驗結果更為吻合。具體而言,30%低密實度條件下,體積應變均與試驗一樣,呈現“不斷減小,逐步放緩”的特點,且柔性邊界數據與試驗更為接近。70%高密實度條件下,柔性邊界體應變與試驗數據吻合較好,表現為“先剪縮,后剪脹”的變化特點;而剛性邊界條件下,體積應變則表現為“先剪縮,后剪脹,再剪縮”的趨勢,其最后階段的“剪縮”特征與實際試驗存在較大誤差。

圖8 體積應變對比Fig.8 Comparison of volume strain

3.3 剪切破壞過程分析

柔性和剛性邊界下試樣顆粒的位移矢量如圖9所示。當采用柔性邊界時,70%密實度試樣兩側出現明顯的腰型鼓脹變形,變形特性與實際的三軸試驗更為接近;而采用剛性邊界時,試樣的側向鼓脹變形受到側向剛性墻體的約束,其發展受限,側向未發生明顯的鼓脹變形。因此,柔性邊界能夠更為清楚地反映出三軸剪切過程中的側向變形特性。

圖9 試樣中顆粒的位移矢量圖Fig.9 Diagram of displacement vector for particles in samples

進一步選取不同密實度情況下,對各模型的加載過程進行切片,分析各邊界條件下試樣的剪切破壞特性與變形規律。相關研究中,顆粒材料破壞往往伴隨著剪切帶的形成。通過對加載過程中的HST95砂土模型試樣進行切片后,根據試樣外部的變形特性和內部顆粒的轉動特征,分析HST95砂土試樣破壞后形成的剪切帶。

圖10給出了各邊界條件在5%、10%和20%相對密實度下,試樣顆粒軸向應變的位移剖面云圖。從圖10中可以看出,各組模型均體現出三軸加載過程中剪切帶的形成,說明模型能有效地還原土體試樣的剪切破壞特性。隨著軸向應變增加,試樣的剪切帶開始由局部形成,之后隨著軸向應變繼續發展。對比柔性和剛性邊界的試樣位移云圖,柔性邊界試樣在剪切應變作用下開始出現側向鼓脹變形,使得剪試樣內位移云圖進一步變化。由此可知,試樣在軸向位移加載下,隨著軸向應變增大,試樣的剪切帶會不斷發展,且柔性邊界會使剪切帶進一步變化。對比30%和70%密實度的工況,可以發現在高密實度狀態下,剪切破壞帶發展更快。

圖10 試樣中顆粒軸向應變的位移云圖Fig.10 Displacement nephogram of particles axial strain in samples

4 結 論

采用DEM軟件PFC,基于赫茲非線性接觸模型、剛性墻圍壓邊界和柔性黏結顆粒膜的圍壓邊界方法,開展了不同密實度下HST95標準砂試樣的DEM三軸試驗模擬,分析了不同邊界條件對DEM模擬三軸試樣的力學特性、變形特性和變形破壞規律,得出本試驗特定條件下的結論。

(1)30%和70%密實度試樣在剛性邊界加載下,其結果與試驗結果對比,誤差在7%～53%區間,而柔性邊界加載結果與試驗結果對比誤差均在10%以內, 對于相對密實度越大的試樣,柔性邊界加載對于三軸試驗力學特性的還原程度越優于剛性加載。

(2)柔性邊界加載“松砂剪縮硬化、密砂剪脹軟化”的變化趨勢與數值及室內試驗結果基本一致,而剛性邊界加載下密砂剪切應變超過12.5%后體應變呈現增大的變化趨勢與試驗不符。柔性邊界加載下對于30%和70%密實度試樣其20%剪切應變下的體應變與試驗結果誤差均在1%以內,而剛性邊界加載下誤差分別為6.7%和32.1%。對于相對密實度越大的試樣,柔性邊界加載對于三軸試驗變形特性的還原程度越優于剛性加載。

(3)通過對剪切帶形成過程分析發現,剛性和柔性邊界條件下,模型均體現出了剪切破壞的特性,且在柔性邊界狀態下,隨著應變繼續增加,試樣側向不均勻鼓脹變形加劇,試樣內部顆粒位移進一步變化,能更好反映剪切帶的形成過程。