基于PFC3D的加權(quán)平均數(shù)建模法在三軸數(shù)值試驗中的應(yīng)用*

張志華 張謝東 邱紅勝 武 林 鄧雅思

(武漢理工大學(xué)交通學(xué)院 武漢 430063)

基于PFC3D的加權(quán)平均數(shù)建模法在三軸數(shù)值試驗中的應(yīng)用*

張志華 張謝東 邱紅勝 武 林 鄧雅思

(武漢理工大學(xué)交通學(xué)院 武漢 430063)

基于室內(nèi)三軸試驗研究成果，以PFC3D為工具建立三維顆粒流模型，開展三軸數(shù)值試驗研究.針對大型三軸數(shù)值試驗大范圍粒徑顆粒生成困難問題，提出了加權(quán)平均數(shù)數(shù)值建模方法，并以轉(zhuǎn)動顆粒場分析手段為依據(jù)進行局部裂隙細(xì)觀分析.研究表明，加權(quán)平均數(shù)數(shù)值建模法得到的顆粒級配與室內(nèi)試驗顆粒級配及三軸試驗結(jié)果吻合度都較好，此法有利于提高計算效率；轉(zhuǎn)動顆粒場有利于三維立體模型剪切面的細(xì)觀分析，將其與顆粒位移場相結(jié)合進行細(xì)觀分析能更好的揭示粗粒土宏觀力學(xué)及變形特性的內(nèi)在機制.

顆粒流模型；三軸試驗；加權(quán)平均數(shù)顆粒生成法；轉(zhuǎn)動顆粒場

0 引 言

由于粗粒土中存在大粒徑顆粒，一般的小型三軸試驗很難對粗粒土進行試驗研究，而大型三軸試驗試樣的體積較大，試驗研究耗時耗力，且對于大型三軸試驗的數(shù)值模型的建立又是一大難題.對于建立成功的模型，細(xì)觀分析手段的選用是揭示粗粒土宏觀特性成敗的關(guān)鍵.Cundall等[1-2]提出了離散元理論，而顆粒流(PFC)作為離散元方法的一種，解決顆粒體非連續(xù)介質(zhì)數(shù)值模擬問題，并將顆粒流的思想應(yīng)用到土體中.

近幾年來，國內(nèi)外學(xué)者對顆粒流理論運用越來越多.Feng等[3-5]給出了二維區(qū)域內(nèi)部圓盤顆粒填充的前進邊法，而對于三維區(qū)域內(nèi)球體顆粒的前進面法并未給出具體的實現(xiàn)過程；耿慶東[6]提出了三維區(qū)域內(nèi)球體顆粒生成的前進面法，確定了前進面的選取順序，還通過改變半徑解決了漏填的情況，得到了三維區(qū)域內(nèi)部球體的高效填充方法；童朝霞等[7]采用截取法制樣，先重力沉積大量的密集顆粒體，再從中截取與室內(nèi)試樣尺寸相同的矩形試樣進行雙軸試驗；周劍等[8]對真實物體進行數(shù)字圖像二值處理，得到反映土石混合體結(jié)構(gòu)的圖像，將數(shù)字圖像轉(zhuǎn)換為PFC模型；Liu等[9]探討了不同黏結(jié)強度下的屈服應(yīng)力導(dǎo)致顆粒間裂隙的產(chǎn)生和發(fā)展不同；以細(xì)小顆粒簇模擬不規(guī)則非球狀的顆粒體進行數(shù)值模擬試驗，數(shù)值結(jié)果吻合程度較高，但此方法耗時耗力[10-12].

基于此，文中提出了加權(quán)平均數(shù)顆粒生成法進行大型三軸數(shù)值試驗?zāi)M，并根據(jù)擬合室內(nèi)三軸試驗結(jié)果為依據(jù)，以轉(zhuǎn)動顆粒場為細(xì)觀分析手段，并與顆粒位移場等分析手段相結(jié)合，進行試樣局部破壞過程分析，對今后粗粒土三軸顆粒流仿真試驗研究提供一定的參考價值.

1 顆粒轉(zhuǎn)動模型

在顆粒流理論中，顆粒間的相互作用看成一個動態(tài)的過程，在這個過程中顆粒內(nèi)部之間的力是平衡的，系統(tǒng)處于基本平衡狀態(tài)在運行.顆粒流理論運用的位移法與顆粒-顆粒間接觸點的相對位移和接觸力有關(guān)，球與球的接觸，在平行黏結(jié)模型中，隨著黏結(jié)劑發(fā)生變形，接觸點處會產(chǎn)生額外的接觸力和力矩；在接觸黏結(jié)模型中，接觸點處只會產(chǎn)生接觸力，由于接觸黏結(jié)劑不能傳遞力矩而不會產(chǎn)生力矩.

但是顆粒間的相互作用力大于黏結(jié)劑強度時黏結(jié)劑破裂，顆粒間形成裂隙.此時的顆粒間相互作用力就發(fā)生了改變，顆粒轉(zhuǎn)動模型見圖1.

圖1 顆粒轉(zhuǎn)動模型

(1)

(2)

(3)

2 三軸試驗顆粒流模擬

2.1 粗粒土三軸試驗宏觀力學(xué)參數(shù)

試驗使用天水洪山公司制造的1 500kN微機控制電液伺服靜態(tài)三軸試驗機，試驗用料取自水布埡大壩工程當(dāng)年的施工料場，是水布埡茅口組灰?guī)r堆石料.參照《土工試驗規(guī)程》(SL237—1999)對土的飽和吸水率、飽和時極限抗壓強度、比重、軟化系數(shù)等基本物性指標(biāo)進行測試，得到其物理性能.結(jié)果見表1.

表1 粗粒土的基本物性指標(biāo)

試驗三軸試樣的尺寸直徑×高度為300 mm×600 mm，干密度為2 180 kg/m3，得到一個試樣的總重量為92.457 kg.采用篩分法和混合法[13]測試粗粒土的顆粒級配，得到粗粒土不同粒徑范圍內(nèi)顆粒的含量見表2.

表2 粗粒土室內(nèi)試驗顆粒級配

2.2 加權(quán)平均數(shù)顆粒生成法

在進行顆粒生成時，如果按照試驗室內(nèi)的土料級配進行生成顆粒時，會發(fā)現(xiàn)一個很嚴(yán)重的問題，就是當(dāng)生成顆粒粒徑小于1 mm的顆粒時，雖然在該粒徑下的土料所占比例較小，但是由于單個顆粒的質(zhì)量很小，導(dǎo)致在顆粒流模型中會生成大量的粉末狀的小粒徑顆粒，嚴(yán)重影響了模型的運算速度.

加權(quán)平均數(shù)法是根據(jù)室內(nèi)試驗材料的顆粒級配作為數(shù)據(jù)基礎(chǔ)，利用若干個粒徑范圍內(nèi)的顆粒含量作為權(quán)數(shù)并與相應(yīng)的顆粒粒徑的乘積之和，計算得到的該加權(quán)算術(shù)平均數(shù).即

(4)

式中:Ri為顆粒平均半徑；Mi為某一粒徑下的顆粒含量(質(zhì)量分?jǐn)?shù))；Wi為某一粒徑的加權(quán)數(shù)，Wr為顆粒粒徑的加權(quán)平均數(shù);∑Mi為各顆粒粒徑占總質(zhì)量的含量之和，其值為1.

根據(jù)室內(nèi)顆粒級配(見表2)和式(3)，計算得到室內(nèi)試驗一個試樣內(nèi)的顆粒平均半徑為11.066 25mm.在顆粒流模型中，如果按照顆粒最小半徑及最大最小顆粒半徑比進行顆粒生成時，根據(jù)系統(tǒng)中的Random口令，顆粒半徑從最小到最大的生成概率是相同的.依據(jù)室內(nèi)試驗的加權(quán)平均顆粒半徑，如果固定顆粒的最小半徑，即可算出顆粒的最大半徑.計算式為

(5)

式中：Ri為表2中的顆粒半徑；ratioRmax/Rmin為顆粒最大半徑與最小半徑之比；Rmin為最小顆粒半徑.根據(jù)式(4)算得的Wr，即可算ratioRmax/Rmin的值.

根據(jù)室內(nèi)試驗顆粒級配，選用模型內(nèi)最小顆粒半徑為0.5 mm，則根據(jù)式(4)得出，模型內(nèi)顆粒最大半徑與最小半徑之比為59.265，因此最大顆粒半徑為29.632 5 mm，接近室內(nèi)試驗試樣的最大顆粒尺寸.根據(jù)隨機生成概率相同的原理，得到顆粒流模型內(nèi)顆粒的級配，見表3，并進行室內(nèi)顆粒繼配及顆粒流模型中顆粒繼配曲線對比，見圖2.

表3 顆粒流模型內(nèi)顆粒級配

圖2 室內(nèi)試驗及數(shù)值試驗顆粒級配曲線

由圖2可知，加權(quán)平均數(shù)建模法生成的顆粒級配與室內(nèi)試驗顆粒級配基本吻合.

2.3 各細(xì)觀參數(shù)的確定

根據(jù)室內(nèi)試驗結(jié)果的數(shù)據(jù)處理，并結(jié)合顆粒流模型內(nèi)考慮到的各細(xì)觀參數(shù)如初始彈性模量Ec，顆?？辜魪姸萩，φ，顆粒剛度比kn/ks，顆粒最小半徑Rmin，顆粒最大、最小半徑比ratioRmax/Rmin；試樣的密度ρ，初始接觸點個數(shù)Nf，墻面的剛度Wstiff及墻面的摩擦系數(shù)Wfric等.文中根據(jù)室內(nèi)試驗結(jié)果中應(yīng)力應(yīng)變曲線的彈性階段(起始直線部分)計算初始彈性模量，根據(jù)摩爾應(yīng)力圓得到顆粒的抗剪強度指標(biāo)，根據(jù)試驗試樣確定材料的密度，基于加權(quán)平均數(shù)計算方法計算得到數(shù)值試驗顆粒級配曲線中各顆粒粒徑的含量.得到一系列細(xì)觀參數(shù)值，見表4.

表4 三軸試驗?zāi)Ｐ椭懈骷?xì)觀參數(shù)(圍壓0.8 MPa)

2.4 數(shù)值結(jié)果與試驗結(jié)果對比

根據(jù)表4中確定的模型內(nèi)的各細(xì)觀參數(shù)的值，并根據(jù)PFC3D手冊[14-15]中不同圍壓下模型內(nèi)顆粒的初始楊氏模量和顆粒剛度比不同.根據(jù)泊松比與顆粒剛度比之間的關(guān)系[16],在已知室內(nèi)試驗結(jié)束時測得試樣的泊松比，可計算各圍壓下的顆粒剛度比.圖3為不同圍壓下數(shù)值試驗應(yīng)力應(yīng)變曲線與室內(nèi)試驗結(jié)果的對比圖.

圖3 不同圍壓下室內(nèi)試驗與數(shù)值試驗應(yīng)力應(yīng)變曲線

3 以0.8 MPa圍壓為例的局部破壞過程分析

在室內(nèi)三軸試驗中，低圍壓下試樣容易發(fā)生剪切破壞，而在高圍壓下試樣則不容易發(fā)生剪切破壞.但由圖3中應(yīng)力應(yīng)變曲線可知，在不同圍壓下，試樣在軸應(yīng)變5%～10%之間發(fā)生了從彈性變形到塑性變形的過渡.在這個過程中試樣內(nèi)部是否發(fā)生了剪切破壞，這就需要從細(xì)觀角度進行分析.

為了分析數(shù)值試驗及室內(nèi)三軸試驗局部破壞過程，提出轉(zhuǎn)動顆粒場分析手段，并與顆粒位移場進行分析，對比模型剪切面的形成機理.

3.1 不同的細(xì)觀分析方法

3.1.1 轉(zhuǎn)動顆粒場

依照位移法提出了轉(zhuǎn)動顆粒場細(xì)觀分析手段，根據(jù)局部破壞過程中的顆粒轉(zhuǎn)角大小分布情況，找出了二維顆粒流雙軸試驗剪切帶的分布情況，表明在模型加載過程中顆粒材料發(fā)生剪切破壞時局部顆粒伴隨有較大的轉(zhuǎn)動.采用X射線及光學(xué)測量手段分析了砂土的微觀結(jié)構(gòu)變形，發(fā)現(xiàn)在局部化砂土破壞形成過程中，局部區(qū)域的顆粒顯示出較大的顆粒旋轉(zhuǎn).

圖4 不同軸應(yīng)變時的轉(zhuǎn)動顆粒場

圖5給出了不同軸應(yīng)變時的顆粒旋轉(zhuǎn)速度場，當(dāng)顆粒的瞬時轉(zhuǎn)動速度達(dá)到110 rad/s時(以0.8 MPa為例，粘結(jié)顆粒發(fā)生破裂時，接觸力均達(dá)到兆帕，顆粒間接觸力很大，由式(3)可知，質(zhì)量很小的顆粒在較大力矩作用下產(chǎn)生的瞬時角加速度很大)，顆粒被標(biāo)定為黑色，因此可以明顯的看出隨著軸應(yīng)變的增大顆粒在局部區(qū)域內(nèi)的旋轉(zhuǎn)速度情況，尤其是在剪切面附近的顆粒.因此，隨著剪切面的形成，剪切面附近的顆粒伴隨有較大的轉(zhuǎn)動.

圖5 不同軸應(yīng)變時的顆粒位移場

由圖4可知，軸應(yīng)變在5%～10%內(nèi)模型中轉(zhuǎn)動速度較大的顆粒數(shù)目較少，在軸應(yīng)變達(dá)到15%時，從外觀上可以很明顯的看到黑色顆粒的斜面分布規(guī)律，該斜面為剪切面.且隨著軸應(yīng)變不斷增大，分布在剪切面附近的顆粒數(shù)目越來越多.通過轉(zhuǎn)動速度較大的顆粒可以很好的找出模型中的剪切面的分布情況.比較圖3～4可知，試樣在壓縮剪切過程中達(dá)到應(yīng)力峰值之后的應(yīng)變軟化階段，剪切面才會完全形成.

3.1.2 顆粒位移場

觀察不同軸應(yīng)變時的顆粒位移場發(fā)現(xiàn)，見圖5a)，從壓縮剪切開始，由于上層加載板向下運動，導(dǎo)致模型上半部分顆粒集體向下運動，下半部分顆粒呈現(xiàn)兩側(cè)運動.在軸應(yīng)變5%左右時模型上半部分顆粒的位移分布比較均勻，上層顆粒均勻向下移動，基本沒有分塊移動；當(dāng)軸應(yīng)變達(dá)到10%時，上半部分顆粒開始出現(xiàn)一個三角區(qū)域，見圖5b)，該三角區(qū)域是試樣剪切破壞之后的塊體，剪切面分布在三角區(qū)域的下邊緣.也就是說當(dāng)軸應(yīng)變在10%左右，試樣出現(xiàn)剪切面初始結(jié)構(gòu)，試樣處于彈性到塑性階段的過渡過程，應(yīng)力應(yīng)變曲線出現(xiàn)應(yīng)力峰值，與圖3中室內(nèi)試驗結(jié)果吻合.在軸應(yīng)變?yōu)?5%～25%階段時，見圖5c)～圖5e)，可以很明顯的看到模型右上角三角區(qū)域顆粒位移整體呈現(xiàn)右下方向，剪切面徹底形成，試樣處于塑性階段，應(yīng)力應(yīng)變曲線出現(xiàn)應(yīng)變軟化現(xiàn)象.因此剪切面的完整形成是在試樣應(yīng)力應(yīng)變曲線達(dá)到峰值強度之后的軟化階段.

3.2 轉(zhuǎn)動顆粒與空間剪切面

由于三軸試驗三維立體的內(nèi)部不可視性，導(dǎo)致形成的空間剪切面也無法看見.文中提出的轉(zhuǎn)動顆粒場以轉(zhuǎn)動顆粒跟蹤試樣空間剪切破壞面，比較圖5～6可以確定轉(zhuǎn)動顆粒場的可用性.基于此，文中在提出轉(zhuǎn)動顆粒場的基礎(chǔ)上，根據(jù)數(shù)值模型顆粒的空間處理，得到去除轉(zhuǎn)動顆粒的試樣空間剪切帶和轉(zhuǎn)動顆粒的空間分布圖，見圖7.

圖6 去除轉(zhuǎn)動顆粒的空間剪切帶空隙

圖7 不同軸應(yīng)變時的轉(zhuǎn)動顆粒空間分布

圖6是去掉轉(zhuǎn)動顆粒之后的顆粒模型，可以很直觀的看出圍壓0.8 MPa下試樣的剪切破壞“面”，該破壞面并不是一個簡簡單單的單顆粒厚度的凸凹面，而是具有一定厚度，且該剪切面的厚度在模型邊緣較薄，在模型內(nèi)部較厚.

圖7表示轉(zhuǎn)動顆粒在模型剪切壓縮過程中的生成與分布，由圖7a)～圖7b)可知，在初始壓縮階段，試樣內(nèi)部出現(xiàn)了大量的轉(zhuǎn)動顆粒，這是試樣內(nèi)產(chǎn)生局部裂隙導(dǎo)致的.圖7c)可知，轉(zhuǎn)動顆粒的分布呈斜面分布，與圖5c)比較得知，該斜面是剪切面.因此，在軸應(yīng)變10%～15%(含15%)，某一時刻試樣被剪破，這與模型宏觀力學(xué)反應(yīng)-應(yīng)力應(yīng)變曲線是吻合的.

將轉(zhuǎn)動顆粒與其位移場結(jié)合見圖8，可以看出轉(zhuǎn)動顆粒的位移場與顆粒的總位移場(圖6)是一致的.轉(zhuǎn)動顆粒在轉(zhuǎn)動過程中沿著剪切面的走勢滾動，滾動的顆粒更加促使了剪切面上相對摩擦力的減小，使得宏觀力學(xué)特性表現(xiàn)出抗剪強度下降的現(xiàn)象，從而形成了峰值過后應(yīng)變軟化的現(xiàn)象.

圖9表示在三軸數(shù)值試驗過程中轉(zhuǎn)動顆粒的數(shù)量，得到轉(zhuǎn)動顆粒數(shù)量與軸應(yīng)變的關(guān)系式為一元二次方程，由于顆粒的轉(zhuǎn)動是顆粒間黏結(jié)劑破裂導(dǎo)致的，表明在剪切破壞初期，模型內(nèi)裂隙數(shù)量呈線性增加，當(dāng)試樣剪切破壞以后，模型內(nèi)部顆粒間裂隙的數(shù)量增加減緩，這是模型內(nèi)形成了剪切面造成的.剪切面的形成使得分布在剪切面上下的顆粒體會沿著最小抗剪強度的位置移動，也導(dǎo)致了模型整體抗剪強度的下降.

圖8 軸應(yīng)變?yōu)?0%時的轉(zhuǎn)動顆粒位移場

圖9 轉(zhuǎn)動顆粒的數(shù)量與軸應(yīng)變的關(guān)系

4 結(jié) 論

1) 以加權(quán)平均數(shù)顆粒生成法建立粗粒土三軸顆粒流模型，顆粒級配與室內(nèi)試樣級配大致吻合，且得到的應(yīng)力應(yīng)變曲線與室內(nèi)試驗結(jié)果基本吻合，驗證了該方法在三軸顆粒流數(shù)值試驗中的可用性.

2) 轉(zhuǎn)動顆粒場分析手段能很好的揭示低圍壓下粗粒土應(yīng)變軟化現(xiàn)象機理，并將其與顆粒位移場相結(jié)合進行細(xì)觀分析，能更好的揭示粗粒土三軸試驗宏觀力學(xué)特性機理.尤其是以轉(zhuǎn)動顆粒為分析對象，充分解釋了低圍壓下剪切破壞初期模型內(nèi)部裂隙數(shù)量線性增加，裂隙貫穿形成剪切面的過程，而且根據(jù)轉(zhuǎn)動顆粒空間位置分布，得出剪切破壞面是一個厚度不均勻，形狀不規(guī)則的“斜面板”.

[1]CUNDALL P A. A computer model for simulating progressive large scale movements in blocky rock systems[J]. Proceedings of the Symposium of the International Society of Rock Mechanics,1971(1):2-8.

[2]CUNDALL P A, STRACK O D L. A discrete numerical model for granular assemblies[J]. Geotechnique,1979,29(1):47-65.

[3]FENG Y T, HAN K, OWEN D R J. Filling domains with disks: an advancing front approach[J]. International Journal for Numerical Methods in Engineering,2003,56(3):699-713.

[4]FENG Y T, HAN K, OWEN D R J. Filling domains with disks[D]. Glasgow: University of Glasgow,2001.

[5]FENG Y T, HAN K, OWEN D R J. An advancing front packing of polygons, ellipses and spheres[C]. Proceedings of the 3rdinternational Conference on Diserete Element Methods: Numercial Modeling of Discontinua,Santa Fe, New Mexico, USA,2002.

[6]耿慶東.公共面搜索算法改進及三維顆粒體生成算法[D].大連：大連理工大學(xué),2009.

[7]童朝霞,周敏,張連衛(wèi),等.各項異性顆粒材料雙軸壓縮試驗的離散元數(shù)值模擬[J].巖石力學(xué)與工程學(xué)報,2014,33(增刊2)：4227-4232.

[8]周劍,張路青,戴福初,等.基于黏結(jié)顆粒模型某滑坡土石混合體直剪試驗數(shù)值模擬[J].巖石力學(xué)與工程學(xué)報,2013,32(增刊1)：2650-2659.

[9]LIU S H, SUN D A, WANG Y S. Numerical study of soil collapse behavior by discrete element modeling[J]. Computers and Geotechnics,2003,30:399-408.

[10]PRICE M, MURARIU V, MORRISON G. Sphere clump generation and trajectory comparison for real particles[J]. Proceedings of Discrete Element Modelling,2007,15:245-256.

[11]MATSUSHIMA T. Discrete element simulation of an assembly of irregularly shaped grains: quantitative comparison with experiments[C]. 16th ASCE Engineering Mechanics Conference, University of Washington, Seattle,2003.

[12]LEE Y. A packing algorithm for three-dimensional convex particles[J]. Granular Matter,2009,11(5):307-315.

[13]李廣信.高等土力學(xué)[M].北京：清華大學(xué)出版社,2004.

[14]Itasca Consulting Group Incorporation. PFC3Dtheory and background[R]. Minneapolis：Minnesota,2004.

[15]Itasca Consulting Group Incorporation. PFC3Duser’s guide[R]. Minneapolis：Minnesota,2004.

[16]張志華.基于PFC3D的粗粒土三軸試驗細(xì)觀數(shù)值模擬[D].宜昌:三峽大學(xué),2015.

Weighted Average Assembling Method Using in Triaxial Numerical Experiment Based on PFC3D

ZHANG Zhihua ZHANG Xiedong QIU Hongsheng WU Lin DENG Yasi

(SchoolofTransportation,WuhanUniversityofTechnology,Wuhan430063,China)

Based on the results of indoor triaxial experimental research, a three-dimensional particle flow model is established based on PFC3Dto conduct triaxial numerical experiments. Based on the problem of a wide range of particle size in large-scale triaxial numerical test, numerical modeling method of weighted average assembling is proposed. Based on the analysis method of rotating particles field, local fracture microstructure analysis has been conducted. The results show that the particle gradations from laboratory agree well with the simulation results generated by weighted average generate method. In addition, the results of macro responses have great agreement. This method is beneficial to improving the calculation efficiency. Besides, the field of rotating particles in the shear mesoscopic analysis of three-dimensional model is advantageous. By combining the displacement field, the particle mesoscopic analysis can better reveal the internal mechanism of mechanics and deformation characteristics of coarse grained soil.

particle flow model; triaxial test; weighted average generate method; field of rotating particles

2016-02-06

*國家自然科學(xué)基金項目資助(51308429)

TU43

10.3963/j.issn.2095-3844.2017.02.017

張志華(1991—)：男，博士生，主要研究領(lǐng)域為巖土體及隧道圍巖穩(wěn)定性離散元數(shù)值仿真