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

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

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

引言

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

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

目前采用圍壓柔性加載的三軸試驗離散元模擬還較少,不同邊界加載圍壓下三軸砂土試樣的力學(xué)特性、變形特性以及變形破壞規(guī)律也尚未探明。為此,利用PFC6.0程序開展不同密實度的HST95標(biāo)準(zhǔn)砂三軸模型試驗?zāi)M,并結(jié)合試驗數(shù)據(jù)討論不同圍壓加載方式下DEM三軸試驗?zāi)Ｐ秃侠硇?從顆粒細(xì)觀尺度上對比分析了不同邊界加載方式下砂土試樣剪切帶的形成、發(fā)展過程。

1 模型基本參數(shù)

1.1 砂土顆粒生成

本文采用的DEM模擬模型中砂土是基于英國鄧迪大學(xué)一系列巖土試驗采用的HST95標(biāo)準(zhǔn)砂土[20],該砂土基本物理參數(shù)如表1所示。

表1 HST95砂土物理參數(shù)Table 1 Physical parameters of HST95 sand

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

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

1.2 DEM模型建立

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

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

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

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

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

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

2 DEM加載模擬原理

2.1 柔性邊界模擬

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

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

2.2 圍壓邊界原理

在剛性邊界三軸試驗?zāi)M過程中,圍壓通過剛性墻進行加載,限制試樣發(fā)生同步側(cè)向變形。在柔性邊界三軸試驗?zāi)M過程中,圍壓則需通過等效集中力的方式施加在球體膜顆粒上,從而無需強迫發(fā)生同步側(cè)向變形,在每一次邊界循環(huán)過程中進行集中力的計算判斷,根據(jù)膜顆粒的相對位置不斷調(diào)整至所有試樣球體顆粒-膜顆粒平均應(yīng)力達到目標(biāo)值。對于柔性膜顆粒而言,任意球體顆粒皆與周圍顆粒接觸在一起。圖4為膜顆粒集中力的計算示意。

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

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

F=Aσ

(1)

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

A=l1×l2

(2)

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

2.3 試樣顆粒接觸模型

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

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

Fc=Fh+Fd

(3)

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

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

表2 DEM模型接觸參數(shù)Table 2 Contact parameters of DEM model

3 DEM試驗?zāi)M分析

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

3.1 力學(xué)特性分析

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

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

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

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

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

3.2 變形特性分析

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

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

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

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

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

圖8 體積應(yīng)變對比Fig.8 Comparison of volume strain

3.3 剪切破壞過程分析

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

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

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

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

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

4 結(jié) 論

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

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

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

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