【PFC6.0.30】三维Cluster模拟GBM矿物晶粒岩石单轴

lobby 2023年5月14日浏览：1992 评论：2 收藏：7

0 引言

首先祝大家劳动节快乐，劳动人民最光荣！

这应该是cluster目前最后一个应用场景了。

离散元的模拟思路是从微观力学行为去反映宏观特性，在这个过程中，能够实现的现实的因素越多，得到的力学行为也就越准确。所以我们很多人去做柔性三轴，并不是为了去研究力学特性，而是为了得到更加准确的宏观特性与现实相比对。

岩石也是一样，一个完整的岩石应当包括矿物晶粒和胶结物，并且除了致密的花岗岩这种岩浆岩，沉积岩变质岩在内部或多或少都会存在微小裂纹，甚至有一些碳酸岩体内部还存在微小孔洞。

本文主要是利用cluster的概念，使用一个个cluster来模拟矿物晶粒，从模拟的思路来看是能够更好的反映岩石行为的。但是需要注意的是，矿物晶粒的尺度对岩石而言是相当微小的，考虑矿物晶粒的破坏对岩石而言是否有必要还应当得到进一步的考量。当然本文是一个纯技术层次的探讨，不在模拟假定方面进行深入探讨。

1 生成晶粒颗粒

这里指定了晶粒的尺寸大小，当然各位可以根据晶粒名称去指定更加复杂的晶粒级配。需要注意的是，后面使用的rBlock构建方式是from-ball，这里采用的是ball的边界进行计算的，所以在第一步颗粒和墙体不能有过大的重叠，墙体和颗粒的刚度设置的大很多。

model new
def par    width=0.2    height=width*2        poro_jingti=0.3    jingti_radmin=0.015    jingti_radmax=0.02end@par
model domain extent [-width] [width] [-width] [width] [-height] [height]



wall generate box [-width*0.5] [width*0.5] [-width*0.5] [width*0.5] [-height*0.5] [height*0.5] expand 1.5

ball distribute porosity @poro_jingti radius @jingti_radmin @jingti_radmax box [-width*0.5] [width*0.5] ...            [-width*0.5] [width*0.5] [-height*0.5] [height*0.5]cmat default type ball-facet model linear method deform emod 100e9 kratio 1.5  cmat default type ball-ball model linear method deform emod 100e6 kratio 1.5   

ball attribute density 2e3 damp 0.7model cycle 2000 calm 50model solve model save "SampleFirst"

【PFC6.0.30】三维Cluster模拟GBM矿物晶粒岩石单轴的图1

2 构建rblock

使用rblock构建的命令自动进行GBM范围划分

model restore "SampleFirst"

rblock construct from-balls polydisperse true  ball delete
model save "rblock_Sample"

【PFC6.0.30】三维Cluster模拟GBM矿物晶粒岩石单轴的图2

3 生成新的晶粒内颗粒

在这里首先对已有的rblock的形状进行导出，然后删除rblock，后续使用这些范围对ball的归属进行划分

model restore "rblock_Sample"
def ball_par    ball_radmin=jingti_radmin*0.2    ball_radmax=jingti_radmin*0.3    ball_poro=0.3end@ball_par
def GetGeo    geoname_count=1    loop foreach rb rblock.list        groupName=string.build("geo_%1",geoname_count)        idRblock=rblock.id(rb)        command            rblock export to-geometry @groupName range id @idRblock        endcommand        geoname_count+=1    endloopend@GetGeorblock delete


cmat default type ball-facet model linear method deform emod 200e6 kratio 1.5  ball distribute porosity @ball_poro radius @ball_radmin @ball_radmax box [-width*0.5] [width*0.5] ...            [-width*0.5] [width*0.5] [-height*0.5] [height*0.5]ball attribute density 2e3 damp 0.7model cycle 2000 calm 50model solve model save "sample"

【PFC6.0.30】三维Cluster模拟GBM矿物晶粒岩石单轴的图3

4 确定颗粒归属

首先对每个形状进行随机定义组名，这里的组名中包含一个sin_count，这个指代形状的当前数目，所以实际上每个组的组名都是不同的。但是每个组名都有一个前缀，比如Rock1_1和Rock1_2这两个是不同的组，但是根据前缀可以判定为同一种矿物成分。这样可以避免后续加胶结的时候，同种矿物间被添加矿物内的胶结属性。

第二点是这里每个组的形状都导入进一个墙体，很多学者会用sj模型来模拟滑裂面的特性。实际上我对sj模型一致都是一种不信任的状态。将滑动面颗粒平整化，其实一样可以减小切向的刚度。

model restore "sample"

def GetTemplateName(sin_count)    rd=math.random.uniform    if rd<0.25 then        GetTemplateName="Rock1_"+string(sin_count)    else if rd<0.5 then         GetTemplateName="Rock2_"+string(sin_count)    else if rd<0.75 then         GetTemplateName="Rock3_"+string(sin_count)    else         GetTemplateName="Rock4_"+string(sin_count)    endifend
def GroupBallIn(geo_count_single)    groupName=GetTemplateName(geo_count_single)    geoName="geo_"+string(geo_count_single)    command        ball group @groupName range geometry-space @geoName count 1     endcommandend

def WallIn    loop gen_cur(1,geoname_count-1)        GroupBallIn(gen_cur)        geoName="geo_"+string(gen_cur)        id_count=10000+gen_cur        command              wall import  from-geometry @geoName id @id_count                  endcommand          endloopend@WallInmodel cycle 2000 calm 50model solve model save "sample_group"

【PFC6.0.30】三维Cluster模拟GBM矿物晶粒岩石单轴的图4

5 加胶结

这一步略微有一点麻烦，主要的麻烦的点在于要实现三点：

1）同种矿物给到同种属性

2）组间颗粒应当给到一个比较弱的属性

3）同种矿物不同组间也是一个比较弱的属性

针对这种比较复杂的接触定义，range fish无疑是能够更方便一点，所以各位如果模型种牵扯到复杂接触定义，建议都需要认真学一下range fish的概念。【Range fish的使用】

model restore "sample_group"
wall delete walls range id 3 10000000 

def GetSubGroupName(ct)    end1=contact.end1(ct)    end2=contact.end2(ct)    if type.pointer.name(end1)#"Ball" then        GetSubGroupName=""        exit    endif    if type.pointer.name(end2)#"Ball" then         GetSubGroupName=""         exit    endif    if ball.group(end1)# ball.group(end2) then         GetSubGroupName=""         exit    endif    groupName=ball.group(end1)    str=string.sub(groupName,9,5)    GetSubGroupName = strend
def isRock1In(pos,ct)    sr= GetSubGroupName(ct)        if sr=="Rock1" then        isRock1In=true        exit    endif    isRock1In=falseenddef isRock2In(pos,ct)    sr= GetSubGroupName(ct)        if sr=="Rock2" then        isRock2In=true        exit    endif    isRock2In=falseenddef isRock3In(pos,ct)    sr= GetSubGroupName(ct)        if sr=="Rock3" then        isRock3In=true        exit    endif    isRock3In=falseenddef isRock4In(pos,ct)    sr= GetSubGroupName(ct)        if sr=="Rock4" then        isRock4In=true        exit    endif    isRock4In=falseend
cmat default model linearpbond method deform emod 1e9 kratio 1.5 pb_deform emod 1e9 kratio 1.5 ...        property pb_coh 10e6 pb_ten 10e6 pb_fa 50 fric 0.1 lin_mode 1cmat add 1 model linear method deform emod 10e9 kratio 1.5 ...        property  lin_mode 1  range contact type "ball-facet" cmat add 2 model linearpbond method deform emod 5e9 kratio 1.5 pb_deform emod 5e9 kratio 1.5 ...        property pb_coh 15e6 pb_ten 15e6 pb_fa 50 fric 0.1 lin_mode 1  range fish @isRock1In cmat add 3 model linearpbond method deform emod 6e9 kratio 1.5 pb_deform emod 6e9 kratio 1.5 ...        property pb_coh 18e6 pb_ten 18e6 pb_fa 50 fric 0.1 lin_mode 1  range fish @isRock2In cmat add 4 model linearpbond method deform emod 7e9 kratio 1.5 pb_deform emod 7e9 kratio 1.5 ...        property pb_coh 20e6 pb_ten 20e6 pb_fa 50 fric 0.1 lin_mode 1  range fish @isRock3In cmat add 5 model linearpbond method deform emod 8e9 kratio 1.5 pb_deform emod 8e9 kratio 1.5 ...        property pb_coh 22e6 pb_ten 22e6 pb_fa 50 fric 0.1 lin_mode 1  range fish @isRock4In 
      cmat apply
model clean
contact method bond gap [ball_radmin*1.8]
contact property lin_force 0 0 0
ball attribute force-contact 0 0 0 moment-contact 0 0 0model calm
model cycle 1
model solve
model save "jiaojie"

【PFC6.0.30】三维Cluster模拟GBM矿物晶粒岩石单轴的图5

6 加载

加载这步和之前都一样。

model restore "jiaojie"model mechanical time-total 0wall delete walls range id 3 6 ball attribute displacement multiply 0def wall_init    wpUp=wall.find(2)    wpDown=wall.find(1)end@wall_init

[IZ0=wall.pos.z(wpUp)-wall.pos.z(wpDown)]
[strainrate=0.1]wall attribute vel-z [-IZ0*strainrate*0.5] range id 2wall attribute vel-z [IZ0*strainrate*0.5] range id 1
def jiance    whilestepping    wzss=(wall.force.contact.z(wpUp)-wall.force.contact.z(wpDown))*0.5/(width*width)    wlz=wall.pos.z(wpUp)-wall.pos.z(wpDown)end
[final_time=1e-2/strainrate][baocunpinlv=final_time/40.0][time_record=-100][count=0]def savefile    time=mech.time.total    wezz=(wlz-IZ0)/IZ0    if time-time_record >= baocunpinlv then        filename=string.build("jieguo_%1",count)        command           model save @filename        endcommand        time_record=time        count +=1    endifend
fish callback add @savefile -1.0
program call "fracture.p3fis"
history deletehistory id 1 @wzsshistory id 2 @wezz@track_initmodel solve time @final_time
model save "result"