二维滚刀破岩模拟的是垂直于刀片转动方向的平面,所以不考虑刀片转动,只考虑刀片的垂向侵入。
这里的模拟步骤为:
1)成样
2)预压
3)加胶结
4)加围压
5)加刀片
6)加载切割
一步步的来讲,成样到加胶结和单元试验一模一样,所以这里可以认为是对煤岩进行切割。
一、成样
用户可以定义尺寸的大小,我这里没有过多的关注尺寸,大家可以自行修改。
newdef parwidth=2.0height=width*0.5rdmax=0.009rdmin=0.006poro=0.08end@pardomain extent [-height*5] [height*5]set random 10001wall generate box [-width*0.5] [width*0.5] [-height*0.5] [height*0.5] expand 1.5ball 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.5ball attribute density 2.7e3 damp 0.7cycle 2000 calm 50solvesave sample

二、预压
这部分原则上不需要修改
restore sample=-1e6]=-1e6]=0.5]=true]=true]=100]=global.step-1]def sevro_wallscompute_stressif timestepNow<global.step thenget_g(sevro_factor)=sevro_freqendifif do_xSevro=true thenXvel=gx*(wxss-txx)=-Xvel=Xvelendifif do_ySevro=true thenYvel=gy*(wyss-tyy)=-Yvel=Yvelendifenddef wp_iniwpDown=wall.find(1)wpRight=wall.find(2)wpUp=wall.find(3)wpLeft=wall.find(4)end@wp_inidef computer_chiCunwlx=wall.pos.x(wpRight)-wall.pos.x(wpLeft)wly=wall.pos.y(wpUp)-wall.pos.y(wpDown)enddef compute_stresscomputer_chiCunwxss=-(wall.force.contact.x(wpRight)-wall.force.contact.x(wpLeft))*0.5/wlywyss=-(wall.force.contact.y(wpUp)-wall.force.contact.y(wpDown))*0.5/wlxenddef get_g(fac)gx=0gy=0zongKNX=100e6*2*10zongKNY=100e6*2*10loop foreach ct wall.contactmap(wpLeft)=contact.prop(ct,"kn")endlooploop foreach ct wall.contactmap(wpRight)=contact.prop(ct,"kn")endlooploop foreach ct wall.contactmap(wpUp)=contact.prop(ct,"kn")endlooploop foreach ct wall.contactmap(wpDown)=contact.prop(ct,"kn")endloopgx=fac*wly/(zongKNX*global.timestep)gy=fac*wlx/(zongKNY*global.timestep)endset fish callback -1.0 @sevro_wallshistory id 1 @wxsshistory id 2 @wysshistory id 3 @gxhistory id 4 @gycycle 1solvesave yuya
三、加胶结
到这一步为止试样完全形成了
restore yuya[pb_coh=10e6][ten_coh=2.7]cmat default type ball-facet model linear method deformability emod 10e8 kratio 1.5cmat 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 0cmat applycleancycle 1solvecontact method bond gap [rdmin*0.2]cycle 1solvesave jiajiaojie

四、加围压
这里就和单元试验不一样了,我们要注意岩石经历过的应力状态。
首先岩石需要在围压下赋存,这时候岩石还没有被开挖。这里给的是6MPa的围压
restore jiajiaojie=-6e6]=-6e6]cycle 1solvesave addweiya
之后岩石开挖暴露,开挖面是没有力的,所以需要将上部的墙删除。按道理讲横向的围压还在,不过我觉得直接将两边墙固定也未尝不可,重新编写伺服函数也是比较麻烦的。这里可以参考***或仿真秀二维滚刀破岩的围压部分,那里是有围压的,但是并没有完全完善好。
wall delete walls range id 3set fish callback -1.0 remove @sevro_wallswall attribute vel 0cycle 1solvesave initial_state

这里是岩石面临开挖面的力链情况,可以发现力链总的还是横向为主。
五、加刀片
这里给出的是一个双刀模型,对于刀片形状如下图所示(随手画了个示意图,不好看见谅)。原则上需要四个参数:下底宽downwidth、下底高downheight、总高height、刀尖角度(这里给的10,可以在upwidth右边的math.tan中进行修改)。

后面两个函数就开始生成了,双刀模型多个参数是两个刀片的跨距(kuadu),这个可以用户自己修改。zhongxinX和zongxinY是刀片底部中心的坐标,这里的y坐标高一点,因为之前释放应力的时候试样膨胀过。这个1.01是放大系数,用户可以自己调整。
restore initial_statedef daopianxingzhuangdownwide=width*0.04downHeight=downwide*2.0height=downHeight*2.2upwide=downwide math.tan(10*math.pi/180.0)*downHeight*2kuadu=width*0.2end@daopianxingzhuangdef creatDaopianZuozhongxinX=-0.5*kuaduzhongxinY=0.5*wly*1.01commandwall 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 5endcommandenddef creatDaopianYouzhongxinX=0.5*kuaduzhongxinY=0.5*wly*1.01commandwall 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 6endcommandend@creatDaopianZuo@creatDaopianYousave jiadaopian
这里是代码运行后的样子。

五、加载切割
这部分主要是给刀片速度就可以了。主要定义刀片加载速度(jiazaisulv),切割深度(shendu),为了方便后处理,这里设定了每1s保存一个sav文件(baocunpinlv)。jiance函数中设定了位移和刀尖力。直接调用fracture文件就可以监测裂纹了。
restore jiadaopianset mech age 0ball attribute displacement multiply 0=1] ;;这个是多长时间保存一个save文件=height*0.02] ;单位是m/min=height*0.2]def jianceweiyi=jiazaisulv*mech.agewpzuodap=wall.find(5)wpyoudap=wall.find(6)zuodaoForce=wall.force.contact.y(wpzuodap)youdaoForce=wall.force.contact.y(wpyoudap)daoforce=zuodaoForce youdaoForceend=mech.age-1]=1]def savefileif mech.age-time_record > baocunpinlv thenfilename=string.build("jieguo%1",count)commandsave @filenameendcommandtime_record=mech.agecount =1endifendwall attribute yvelocity [-jiazaisulv] range id 5 6history deletecall fracture.p2fis@track_inithistory id 1 @daoforcehistory id 2 @weiyihistory nstep 50set fish callback -1.0 @jianceset fish callback -1.0 @savefilecalmsolve time [shendu/jiazaisulv]save zuizhongjieguo
下图为加载初期岩石中的力链图,可以看到刀尖力向周围的扩散规律

破坏后的力链图,这里因为速度比较大,岩石破碎形成,导致刀尖悬空,这个在刀尖力变化中也可以看出来

破坏的块体图,可以看出岩石的破碎区域

刀尖力的变化就很明显的看出来刀速度可能有点快了

下面给出裂纹图:
基本上是以剪裂纹为主,这个在裂纹数目中也可以看出来。


破碎体的数目基本上差不多,但是破碎体积占比有点问题,后续还需要再完善一下。
