含结构面岩石直剪试验

  岩块和岩体的区别在于，岩块我们一般指的是完整的均匀的岩石，而岩体则包含岩石在现实中的各种结构面。目前对于完整岩石的研究基本上趋于饱和了，而复合岩石以及结构面岩石的研究还处于发展中状态。

    完整岩石的直剪试验可以参考一下单元试验，基本顺序为：成样-预压-加胶结-加剪切盒-加围压-加载。

    数值模拟时刻要以现实为准绳，这也是我一直强调的，所以结构面在什么时候生成也就显而易见了。岩石在形成后才会因为地层运动产生结构面，于是我们需要在加胶结形成完整岩石后，给岩石加结构面。

    岩石加胶结生成后的状态如图所示：

    可以看到图中显示的为岩石的胶结状态，其胶结已经完整的加进去了，下面我们需要进行的是结构面的生成。

    我们这里采用的是外部导入的结构面形状，我是用adobe illustrate绘制的结构面形状，没有注意尺寸，只绘制了大概的形状，注意这里一定需要使用多段线，如果是用CAD画的圆弧也一定要打断成多段线。

    导入后使用前文(使用fish操作导入的geometry)讲述的方法进行geometry的处理，并将其调整到适合的位置与大小。

    处理完后的geometry如图：

可以发现已经放在了中间位置。

    之后我们使用导入的geometry生成DFN，这里使用dfn gimport geometry 命令就可以完成。

    后面我们需要给DFN处的接触赋属性。

    我这里用的是SJ模型去模拟结构面，SJ模型我对其初步理解为面模型，也就是说一般模型是以接触法向为力的法向作用方向，但是SJ模型是以DFN位置的法向位置作为力的法向方向。由于这个原因，SJ模型经常出现颗粒有比较大的重叠量，不过这个是没有力的，如图所示：

    可以看到在geometry位置上的接触都被赋予了SJ模型，而且部分颗粒出现了比较大的重叠量，后面有时间会做一个双球试验去验证一下这个模型的力学效应。

    结构面有的是无胶结的，这种一般是由于剪切产生的构造面，有的是有胶结的，也就是一些胶结物，例如石英、铁锰物质充填形成的。SJ模型可以加胶结也可以不加胶结，这里只模拟了最简单的情况，也就是光滑的构造面，这里SJ模型只赋予了法向刚度，而切向几乎为0，而且也没有胶结。

    给试样加围压后，由于颗粒的离散性，会在局部形成一系列力的集中点，这个现象随着颗粒数的增加而慢慢减弱。本文的模型采用了4w6的颗粒数，已经足够多，但还是不太能使这个线性消失。

下面展示随着应变的加载，其变形破坏特性：

应变0.002

应变0.004

应变0.012

应变0.018

应变0.026

可以看一下破坏后的块体图：

可以发现上部岩体在右边，而下部岩体在左边形成了破坏区。

力链图：

总的应力应变曲线为：

    可以比较一下之前的单元试验结果(我的PFC岩土颗粒流离散元分析攻略（附赠学习资料）)，相比于完整岩石的9.8MPa的强度，这里不到8MPa的强度已经减小了很多了，说明结构面是起作用的。

    当然结构面的形成方式有很多，SJ模型算是比较好的一个方法了，因为可以考虑面移动，如果采用简单的线性模型其实也可以，但是对颗粒数的要求会比较高。如果是含胶结的话建议还是采用改参数的方式就可以了。

restore jiajiaojie

geometry import jiegoumian.dxfdef get_min_max(id)    global x_min=1e100    global x_max=-1e100        global y_min=1e100    global y_max=-1e100    local gs = geom.set.find(id)    loop foreach local gn geom.node.list(gs)        local pos = geom.node.pos(gn)        if x_min > comp.x(pos)            x_min = comp.x(pos)        endif        if y_min > comp.y(pos)            y_min = comp.y(pos)        endif           if x_max < comp.x(pos)            x_max = comp.x(pos)        endif        if y_max < comp.y(pos)            y_max = comp.y(pos)        endif    endloopenddef moveToOrigin(id)    get_min_max(id)    local gs = geom.set.find(id)    loop foreach local gn geom.node.list(gs)        geom.node.pos.x(gn)=geom.node.pos.x(gn)-(x_min x_max)*0.5        geom.node.pos.y(gn)=geom.node.pos.y(gn)-(y_min y_max)*0.5    endloop    get_min_max(id)enddef sacle(id,n)    local gs = geom.set.find(id)    loop foreach local gn geom.node.list(gs)        geom.node.pos.x(gn)=geom.node.pos.x(gn)*n        geom.node.pos.y(gn)=geom.node.pos.y(gn)*n    endloopend@moveToOrigin(1)@sacle(1,[wlx/(x_max-x_min)])

dfn gimport geometry jiegoumian

dfn model name smoothjoint  installdfn property sj_kn 5e9 sj_ks 1e-4 

cycle 1solvesave addjiegoumian