二维滚刀破岩模拟的是垂直于刀片转动方向的平面,所以不考虑刀片转动,只考虑刀片的垂向侵入。
这里的模拟步骤为:
1)成样
2)预压
3)加胶结
4)加围压
5)加刀片
6)加载切割
一步步的来讲,成样到加胶结和单元试验一模一样,所以这里可以认为是对煤岩进行切割。
一、成样
用户可以定义尺寸的大小,我这里没有过多的关注尺寸,大家可以自行修改。
new
def par
width=2.0
height=width*0.5
rdmax=0.009
rdmin=0.006
poro=0.08
end
@par
domain extent [-height*5] [height*5]
set random 10001
wall generate box [-width*0.5] [width*0.5] [-height*0.5] [height*0.5] expand 1.5
ball distribute porosity @poro radius @rdmin @rdmax box [-width*0.5] [width*0.5] ...
[height*0.5]
cmat default model linear method deformability emod 10e9 kratio 1.5
ball attribute density 2.7e3 damp 0.7
cycle 2000 calm 50
solve
save sample
二、预压
这部分原则上不需要修改
restore sample
-1e6] =
-1e6] =
0.5] =
true] =
true] =
100] =
global.step-1] =
def sevro_walls
compute_stress
if timestepNow<global.step then
get_g(sevro_factor)
sevro_freq =
endif
if do_xSevro=true then
Xvel=gx*(wxss-txx)
-Xvel =
Xvel =
endif
if do_ySevro=true then
Yvel=gy*(wyss-tyy)
-Yvel =
Yvel =
endif
end
def wp_ini
wpDown=wall.find(1)
wpRight=wall.find(2)
wpUp=wall.find(3)
wpLeft=wall.find(4)
end
@wp_ini
def computer_chiCun
wlx=wall.pos.x(wpRight)-wall.pos.x(wpLeft)
wly=wall.pos.y(wpUp)-wall.pos.y(wpDown)
end
def compute_stress
computer_chiCun
wxss=-(wall.force.contact.x(wpRight)-wall.force.contact.x(wpLeft))*0.5/wly
wyss=-(wall.force.contact.y(wpUp)-wall.force.contact.y(wpDown))*0.5/wlx
end
def get_g(fac)
gx=0
gy=0
zongKNX=100e6*2*10
zongKNY=100e6*2*10
loop foreach ct wall.contactmap(wpLeft)
contact.prop(ct,"kn") =
endloop
loop foreach ct wall.contactmap(wpRight)
contact.prop(ct,"kn") =
endloop
loop foreach ct wall.contactmap(wpUp)
contact.prop(ct,"kn") =
endloop
loop foreach ct wall.contactmap(wpDown)
contact.prop(ct,"kn") =
endloop
gx=fac*wly/(zongKNX*global.timestep)
gy=fac*wlx/(zongKNY*global.timestep)
end
set fish callback -1.0 @sevro_walls
history id 1 @wxss
history id 2 @wyss
history id 3 @gx
history id 4 @gy
cycle 1
solve
save yuya
三、加胶结
到这一步为止试样完全形成了
restore yuya
[pb_coh=10e6]
[ten_coh=2.7]
cmat default type ball-facet model linear method deformability emod 10e8 kratio 1.5
cmat default type ball-ball model linearpbond method deformability ...
emod 12e8 kratio 1.5 pb_deformability emod 54e8 kratio 1.5 ...
property pb_coh [pb_coh] pb_ten [pb_coh*ten_coh] pb_fa 50 fric 0
cmat apply
clean
cycle 1
solve
contact method bond gap [rdmin*0.2]
cycle 1
solve
save jiajiaojie
四、加围压
这里就和单元试验不一样了,我们要注意岩石经历过的应力状态。
首先岩石需要在围压下赋存,这时候岩石还没有被开挖。这里给的是6MPa的围压
restore jiajiaojie
-6e6] =
-6e6] =
cycle 1
solve
save addweiya
之后岩石开挖暴露,开挖面是没有力的,所以需要将上部的墙删除。按道理讲横向的围压还在,不过我觉得直接将两边墙固定也未尝不可,重新编写伺服函数也是比较麻烦的。这里可以参考***或仿真秀二维滚刀破岩的围压部分,那里是有围压的,但是并没有完全完善好。
wall delete walls range id 3
set fish callback -1.0 remove @sevro_walls
wall attribute vel 0
cycle 1
solve
save initial_state
这里是岩石面临开挖面的力链情况,可以发现力链总的还是横向为主。
五、加刀片
这里给出的是一个双刀模型,对于刀片形状如下图所示(随手画了个示意图,不好看见谅)。原则上需要四个参数:下底宽downwidth、下底高downheight、总高height、刀尖角度(这里给的10,可以在upwidth右边的math.tan中进行修改)。
后面两个函数就开始生成了,双刀模型多个参数是两个刀片的跨距(kuadu),这个可以用户自己修改。zhongxinX和zongxinY是刀片底部中心的坐标,这里的y坐标高一点,因为之前释放应力的时候试样膨胀过。这个1.01是放大系数,用户可以自己调整。
restore initial_state
def daopianxingzhuang
downwide=width*0.04
downHeight=downwide*2.0
height=downHeight*2.2
upwide=downwide math.tan(10*math.pi/180.0)*downHeight*2
kuadu=width*0.2
end
@daopianxingzhuang
def creatDaopianZuo
zhongxinX=-0.5*kuadu
zhongxinY=0.5*wly*1.01
command
wall create vertices [zhongxinX 0.5*downwide] [zhongxinY] [zhongxinX 0.5*upwide] [zhongxinY downHeight] ...
[zhongxinY downHeight] [zhongxinX 0.5*upwide] [zhongxinY height] ...
[zhongxinY height] [zhongxinX-0.5*upwide] [zhongxinY height] ...
[zhongxinY height] [zhongxinX-0.5*upwide] [zhongxinY downHeight] ...
[zhongxinY downHeight] [zhongxinX-0.5*downwide] [zhongxinY] ...
[zhongxinY] [zhongxinX 0.5*downwide] [zhongxinY] id 5
endcommand
end
def creatDaopianYou
zhongxinX=0.5*kuadu
zhongxinY=0.5*wly*1.01
command
wall create vertices [zhongxinX 0.5*downwide] [zhongxinY] [zhongxinX 0.5*upwide] [zhongxinY downHeight] ...
[zhongxinY downHeight] [zhongxinX 0.5*upwide] [zhongxinY height] ...
[zhongxinY height] [zhongxinX-0.5*upwide] [zhongxinY height] ...
[zhongxinY height] [zhongxinX-0.5*upwide] [zhongxinY downHeight] ...
[zhongxinY downHeight] [zhongxinX-0.5*downwide] [zhongxinY] ...
[zhongxinY] [zhongxinX 0.5*downwide] [zhongxinY] id 6
endcommand
end
@creatDaopianZuo
@creatDaopianYou
save jiadaopian
这里是代码运行后的样子。
五、加载切割
这部分主要是给刀片速度就可以了。主要定义刀片加载速度(jiazaisulv),切割深度(shendu),为了方便后处理,这里设定了每1s保存一个sav文件(baocunpinlv)。jiance函数中设定了位移和刀尖力。直接调用fracture文件就可以监测裂纹了。
restore jiadaopian
set mech age 0
ball attribute displacement multiply 0
1] ;;这个是多长时间保存一个save文件 =
height*0.02] ;单位是m/min =
height*0.2] =
def jiance
weiyi=jiazaisulv*mech.age
wpzuodap=wall.find(5)
wpyoudap=wall.find(6)
zuodaoForce=wall.force.contact.y(wpzuodap)
youdaoForce=wall.force.contact.y(wpyoudap)
daoforce=zuodaoForce youdaoForce
end
mech.age-1] =
1] =
def savefile
if mech.age-time_record > baocunpinlv then
filename=string.build("jieguo%1",count)
command
save @filename
endcommand
time_record=mech.age
count =1
endif
end
wall attribute yvelocity [-jiazaisulv] range id 5 6
history delete
call fracture.p2fis
@track_init
history id 1 @daoforce
history id 2 @weiyi
history nstep 50
set fish callback -1.0 @jiance
set fish callback -1.0 @savefile
calm
solve time [shendu/jiazaisulv]
save zuizhongjieguo
下图为加载初期岩石中的力链图,可以看到刀尖力向周围的扩散规律
破坏后的力链图,这里因为速度比较大,岩石破碎形成,导致刀尖悬空,这个在刀尖力变化中也可以看出来
破坏的块体图,可以看出岩石的破碎区域
刀尖力的变化就很明显的看出来刀速度可能有点快了
下面给出裂纹图:
基本上是以剪裂纹为主,这个在裂纹数目中也可以看出来。
破碎体的数目基本上差不多,但是破碎体积占比有点问题,后续还需要再完善一下。