顆粒流PFC3D數值模擬中標定細觀參數的方法——擬合直剪試驗法

湯保新，徐亮，李琪

（1.揚州大學，江蘇 揚州 225127；2.信息產業電子第十一設計研究院科技工程股份有限公司華東分院，江蘇 無錫 214000）

顆粒流PFC3D數值模擬中標定細觀參數的方法——擬合直剪試驗法

湯保新1，徐亮2，李琪1

（1.揚州大學，江蘇 揚州 225127；2.信息產業電子第十一設計研究院科技工程股份有限公司華東分院，江蘇 無錫 214000）

建立在傳統的連續介質力學基礎上的有限元法等數值計算方法在關于研究非連續介質的問題上難以直接用于計算和模擬材料具體的破壞形式和破壞的整個過程，而離散元法在這一方面顯示出巨大的優勢。但是，對于顆粒模型的細觀參數——法向剛度與切向剛度等參數的取值問題成為創建精確的顆粒模型的一道障礙。文章利用顆粒流離散元軟件PFC3D模擬試驗室直剪試驗，通過不斷調整法向剛度與切向剛度的比值，即剛度比，擬合直剪試驗的結果，最終達到標定顆粒流模型細觀參數的目的。

離散元法；PFC3D；細觀參數；剛度比

1 緒 論

在運用顆粒流離散元軟件 PFC3D進行數值模擬時，需要定義顆粒流模型的細觀參數。對于非動力學問題的分析，需要定義的參數一般包括：顆粒流的摩擦系數、顆粒密度、法向剛度、切向剛度等其它參數。其中摩擦系數和顆粒密度可通過一般的資料查詢即獲取，或者通過簡單的試驗測定得到。而法向剛度與切向剛度是顆粒流的特有參數，一般不能直接簡單地與材料一系列的宏觀物理力學參數聯系，這和連續介質的模型是有本質差別的。沒有確定的剛度比，將使數值模擬的結果顯得毫無依據。本文通過擬合直剪試驗中不同法向應力σ作用下土體的τ-L曲線，并當兩者的τ-L曲線所確定的土體不同正應力σ作用下的抗剪強度無限接近時，此顆粒模型的細觀參數即可用以標定該種土體在其他顆粒模型中的一系列參數。

2 試驗室直剪試驗概述

直剪試驗是測定土體抗剪強度的一種常用方法。原理主要根據庫侖定律，土的內摩擦力與剪切面上的法向壓力成正比。將試驗用土制備成幾個土樣（本次試驗為4個），分別在不同的法向壓力下，沿固定的剪切面直接施加水平力進行剪切，在剪切過程中實時記錄一定階段的手輪圈數或剪切位移量L和其所對應的用以剪切土體的水平力的大小，換算后得到τ-L曲線，根據曲線得其剪壞時的剪應力，即抗剪強度fτ。

本次試驗取4組土樣，分別在100kPa、200kPa、300kPa和400kPa的法向應力作用下測定該土體的抗剪強度fτ 。另，本次試驗所采用的土體為砂土，含水量w=15.67%，密度為1.57g/cm3。試驗室直剪儀如圖1所示。

圖1 試驗室直剪儀

3 PFC3D模擬

顆粒流模型中一般只包含顆粒（ball）與剛性墻體（wall）。墻體不能被賦予密度這一物理屬性，即不服從牛頓第二定律，致使顆粒流模型中的墻體成為一個沒有厚度、沒有“從動”、且不產生變形的完全剛性的墻體。這種墻體可以用來限制顆粒的流動，也可以被定義在某一方向上產生給定的位移。

對于直剪試驗的顆粒流模型，可僅建立直剪儀剪切盒與盒內的土樣。通過定義，使剪切盒發生一定量的位移并使土樣產生剪切破壞。建模的過程可描述為：使用wall命令（關鍵詞為cylinder）首先創建2個底部接觸并且完全吻合的圓柱面，2個圓柱貫通且形成一個大的圓柱，大圓柱的高和底半徑同試驗室的剪切盒，然后為大圓柱創建上下兩個底面，再定義函數“assemble”并使用generate命令在剪切盒中隨機生成足夠量的顆粒；輸入 solve，讓顆粒在無其他外力作用的情況下自由調整各自的位置，最終達到“自平衡”狀態。然后就可以定義剪切位移并實時記錄剪應力。

通過調整顆粒半徑、孔隙率及剛度比等相關參數，得到不同法向應力作用下的τ-L曲線對比圖，如圖2所示。

圖2 不同法向應力作用下的τ-L擬合曲線

圖中第一組S表示法向應力為100kPa時的試驗數據，M則代表模擬數據。第二組為 200kPa，以此類推。從圖中可以看出：第一組與第二組曲線吻合的非常好，幾乎沒有相差很大的值，第三組與第四組的擬合曲線相差稍大，但整體來看，數據的偏差可以接受，尤其根據模擬值與試驗值所得到的對應的抗剪強度相對較接近（見表1與圖3），所以可以采用擬合結果。若在模擬中有較高的精度要求，可繼續調整模型的細觀參數值，直至得到滿意的擬合結果。

最終確定的該砂土顆粒流離散元模型的細觀參數見表2。

4 結 論

應用PFC3D擬合試驗室直剪試驗為顆粒流模型標定細觀參數提供了一種實用的方法，這種方法可以快速地在顆粒流模型細觀參數與材料一系列宏觀物理力學參數之間搭建一座橋梁，但這種方法也存在一定缺陷，主要是在不斷嘗試擬合的過程中幾乎不可能完全與試驗結果一致，這時應根據模擬對象的精度要求判斷擬合結果的可行性。

試驗室直剪試驗的抗剪強度與模擬值 表1

圖3 直剪試驗擬合結果

模型的細觀參數 表2

附：部分函數的定義

def assemble；定義assemble函數

s_stiff=7e3；切向剛度預賦值

n_stiff=1e5；法向剛度預賦值

tot_vol=0.02*pi*0.0309^2.0；圓柱的體積

rbar=0.5*（rlo + rhi）；顆粒的平均半徑

num=int（1.0-poros）*tot_vol/（4.0/3.0*pi*rbar^3）；顆粒的個數

mult=1.6；放大系數

rlo_0=rlo/mult；初始小半徑

rhi_0=rhi/mult；初始大半徑

command

gen id=1，num rad=rlo_0，rhi_0 x=-0.1 0.1 y=-0.1 0.1 z=-0.01，0.01 &

filter ff_cylinder tries 1000000；按要求生成顆粒

prop dens=1570 ks=s_stiff kn=n_stiff；賦予顆粒屬性

end_command

ii=out（string（num）+' particles were created'）

sum=0.0

bp=ball_head

loop while bp # null

sum=sum + 4.0/3.0*pi*b_rad（bp）^3

bp=b_next（bp）

end_loop

pmeas=1.0-sum/tot_vol

mult=（（1.0-poros）/（1.0-pmeas））^（1.0/3.0）

command

ini rad mul mult

plot create assembly

plot set back white

plot add ball orange

plot add wall lblue

pl show assembly

prop fric 2

set grav 0 0-9.81

set dt dscale

end_command

end

；

def ff_cylinder；定義函數ff_cylinder

ff_cylinder=0；命令關閉

_brad=fc_arg（0）

_bx=fc_arg（1）

_by=fc_arg（2）

_bz=fc_arg（3）

_rad=sqrt（_bx^2 + _by^2）；圓柱底面半徑

if _rad + _brad ＞ 0.0309 then

ff_cylinder=1；命令打開

end_if

end

［1］PFC3DVersion 4.0 USER'S GUIDE Section 3 ：Problem Solving with PFC3D.

［2］路培國.壓入式自旋樁頭在砂土中承載特性的研究［D］.揚州：揚州大學，2014.

［3］徐亮.壓入式自旋樁頭在砂土中試驗的顆粒流數值模擬分析［D］.揚州：揚州大學，2015.

［4］郭書魁.豎向加筋砂土直剪試驗的顆粒流數值分析［D］.杭州：浙江大學，2013.

［5］鄧益兵，周健，劉文白，等.螺旋擠土樁下旋成孔過程的顆粒流數值模擬［J］.巖土工程學報，2011，33（9）：1391-1398.

TU502

A

1007-7359（2016）02-0250-03

10.16330/j.cnki.1007-7359.2016.02.089

江蘇省建設系統科技指導項目（2013ZD16）。

湯保新（1967-），男，湖北襄陽人，碩士，副教授，國家注冊一級結構工程師，主要從事混凝土結構教學及抗震研究工作。