UDEC中的直接剪切以及循环剪切模拟案例

1.剪切模型

单元本构：ZONE model elastic（弹性）

节理本构：JOINT model area/residual （库伦滑移节理模型）或者 JOINT model cy （连续屈服节理模型）

法向应力：10MPa

2.直剪模拟

2.1库伦滑移节理模型

点击查看整个模型的源码
小提示：按住Ctrl再点击链接可在新标签页打开代码 😃

节理模型：

joint model area jks 4E4 jkn 4E4 jfriction 30 jdilation 6 zdilation 4E-4 &
 range group 'joint'
set jcondf joint model area jks=4E4 jkn=4E4 jfriction=30 jdilation=6 &
 zdilation=4E-4

这里节理使用的是(面接触的)库伦滑移模型(模型可选的参数可参考上篇:库伦节理模型以及连续屈服模型的参数介绍)，剪切刚度、法向刚度、内摩擦角、剪胀角等就不多说了，这里主要说一下 zdilation 这个参数，它不是代表 z方向的剪胀角，而是指 zero dilation ,即指定产生零剪胀时的剪切位移。

平均剪应力与剪切位移呈线性关系如图，最后极限剪应力大约在剪切位移为0.15mm时达到6MPa。而这时候，开始产生剪胀，法向的剪胀位移与剪切位移大致也呈线性关系，直到剪切位移达到约0.4mm时剪胀为零，即达到设置的 zdilation = 4E-4 m 这个条件，也就是法向位移不随剪切位移增加而变化了。最后最大的剪胀位移大约达到了 2.7E-5 m

库伦滑移模型(JOINT model area )只定义了节理的极限剪切强度值，在节理滑动之后的剪胀位移也是剪胀角的线性函数，具体的库伦滑移模型的原理以后有时间再整理叭。

平均剪应力随剪切位移的变化关系如下：

平均法向位移随剪切位移的变化关系如下：

2.2"强度软化"的库伦滑移节理模型

点击查看整个模型的源码

整个模型与2.1节中的几乎完全相同，只在节理模型定义以后增加了一行代码（set add_dil on）：

joint model area jks 4E4 jkn 4E4 jfriction 30 jdilation 6 zdilation 4E-4 &
 range group 'joint'
set jcondf joint model area jks=4E4 jkn=4E4 jfriction=30 jdilation=6 &
 zdilation=4E-4
set add_dil on

这行命令通过将剪胀角加到设置的内摩擦角中，形成有效内摩擦角。当达到零剪胀时的极限剪切位移时，极限剪应力会产生突变"软化"。与上面的模型相比，通过增加这行命令，仅仅只改变了模型的摩擦角而已。

平均剪应力随剪切位移的变化关系如下：

平均法向位移随剪切位移的变化关系如下：

2.3"残余强度"的库伦滑移节理模型

点击查看整个模型的源码

UDEC提供另一种带有残余强度的库伦滑移模型(JOINT model residual)，以此来模拟节理的峰值剪切强度以及残余剪切强度。

注意：这里为了在初始平衡模型的时候不减少节理的强度，平衡时仍要使用普通的(面接触)库伦滑移节理模型(JOINT model area)

joint model area jks 4E4 jkn 4E4 jfriction 30 jdilation 6 zdilation 4E-4 &
 range group 'joint'
; new contact default
set jcondf joint model area jks=4E4 jkn=4E4 jfriction=30 jdilation=6 &
 zdilation=4E-4

在最后再设置节理模型为 JOINT model residual，进行求解。

group joint 'residual joint'
joint model residual jks 4E4 jkn 4E4 jfriction 35 jcohesion 2 &
                      jdilation 6 zdilation 4E-4 jrescoh 0 jrfric 30 &
                      range group 'residual joint'
new contact default
set jcondf joint model residual jks=4E4 jkn=4E4 jfriction=35 &
                        jcohesion=2 jdilation=6 zdilation=4E-4 &
                        jrescoh=0 jrfric=30

在 JOINT model residual 这个模型中，设置了残余粘聚力与残余内摩擦角，都比初始的粘聚力与内摩擦角小，所以峰值剪应力出现断崖式下跌。最终的剪胀位移大约为 1.8E-5m。

平均剪应力随剪切位移的变化关系如下：

平均法向位移随剪切位移的变化关系如下：

2.4基于用户自定义调整的库伦滑移节理模型

点击查看整个模型的源码

这个说实话，这个还没搞太懂，等以后搞懂了再写吧。。。就放两张图。如果有看懂的也欢迎交流哦

平均剪应力随剪切位移的变化关系如下：

平均法向位移随剪切位移的变化关系如下：

2.5连续屈服节理模型

点击查看整个模型的源码

joint model cy jfriction 30 jif 60 jrough 1E-4 jks 4E4 jkn 4E4 range group &
 'cyjoint'
; new contact default
set jcondf joint model cy jfriction=30 jif=60 jrough=1E-4 jks=4E4 jkn=4E4

针对"强度软化"的节理模型，UDEC 中的连续屈服节理模型( JOINT model cy )就有这种效果。(模型可选的参数可参考上篇:库伦节理模型以及连续屈服模型的参数介绍)。C-Y节理模型模拟了剪切过程中的渐进破坏，具体的C-Y模型中的原理以后再与库伦滑移模型一起整理。与前面模型中的不同之处在于，C-Y 模型的模拟结果从峰值剪应力到残余剪应力之间的过渡比较光滑，剪胀位移也是非线性的光滑的曲线。

平均剪应力随剪切位移的变化关系如下：

平均法向位移随剪切位移的变化关系如下：

其实不管是库伦滑移节理模型还是连续屈服节理模型，还是另外的Barton-Bandis节理模型，在选择节理模型时，要选择最合适的能符合预期实验结果的模型而不是说选最复杂的模型。一般来说，在模拟更复杂的剪切之前，先以简单的模型模拟为基础，来评估各个参数对模型的影响。

3.循环剪切模拟

点击查看整个模型的源码

循环剪切模拟时，仍使用简单的库伦滑移模型，通过与2.1节中一样设置 set add_dil on ，来达到下面的效果：在朝同一个方向剪切时，将剪胀角加到内摩擦角中去，而剪切反向时则减去剪胀角。

joint model area jks 1E6 jkn 1E6 jfriction 30 jdilation 6 range group &
 'cycle joint'
; new contact default
set jcondf joint model area jks=1E6 jkn=1E6 jfriction=30 jdilation=6
set add_dil on

平均剪应力随剪切位移的变化关系如下：

平均法向位移随剪切位移的变化关系如下：

该模型模拟的节理行为可能足以评估循环载荷的影响。但是如果要模拟在载荷作用下节理产生渐进性损伤，则库仑滑移模型可能会不合适，因为这个模型不能模拟这种行为。 连续屈服模型可以模拟节理的渐进性损伤，但不能模拟节理在反向剪切时剪胀的恢复， 可以考虑使用Barton-Bandis模型试试。

4.主要FISH函数说明

FISH函数中主要涉及到两个函数， ini_jdisp 和 av_str ，前者用来初始化剪切位移，后者用来计算平均的剪应力、法向应力、剪位移、法向位移。

注意：在UDEC 项目文件的存储路径或者命令流、注释中尽量不要出现中文，有时候会产生一些错误或者乱码，下面注释的中文只是为了说明。

def ini_jdisp
   njdisp0 = 0.0
   sjdisp0 = 0.0
   ic = contact_head ;获取接触索引表头
   loop while ic # 0  ;判断接触索引非空时循环累加，遍历接触
      njdisp0 = njdisp0 + c_ndis(ic)
      sjdisp0 = sjdisp0 + c_sdis(ic)
      ic = c_next(ic)
   endloop
end
ini_jdisp
;
;Name:av_str
def av_str
   whilestepping    ;每次cycle都运行这个函数
   sstav = 0.0      ;初始化各个值
   nstav = 0.0
   njdisp = 0.0
   sjdisp = 0.0
   ncon = 0
   jl   = 0.2              ;设置节理的长度
   ic = contact_head        ;每次循环中都获取接触索引的表头
   loop while ic # 0        ;遍历接触
      ncon = ncon+1         ;计数
      sstav = sstav + c_sforce(ic)  ;累加接触的剪切力、法向力、剪切位移、法向位移
      nstav = nstav + c_nforce(ic)
      njdisp = njdisp + c_ndis(ic)
      sjdisp = sjdisp + c_sdis(ic)
      ic = c_next(ic)
   endloop
   if ncon # 0
      sstav = sstav / jl        ;计算应力平均值（总的力/节理长）
      nstav = nstav / jl
      njdisp = (njdisp-njdisp0) / ncon  ;计算位移平均值（总位移/总数）
      sjdisp = (sjdisp-sjdisp0) / ncon
   endif
end

长风破浪会有时，直挂云帆济沧海。