来自 UDEC 6.0 Manual →VERIFICATION PROBLEMS→Cyclic Loading of a Specimen with a Slipping Crack
1.问题描述
这个问题主要研究具有内部倾斜的“非贯通”节理的弹性体试样在轴向循环荷载作用下的力学响应。如下图所示：
涉及的问题：(1)非贯通节理的处理；(2)table的使用；(3)平均应力的计算
当然，这里的“非贯通"节理是要打上引号的，实际上UDEC是不能直接模拟非贯通节理的，生成的节理必须将块体完全切割，非贯通节理在运算开始时都会被删除。文中的做法是先设置一条贯通的构造节理，然后通过设置不同区域节理的参数来模拟非贯通节理。
在试样顶部施加轴向位移，试样底部固定，在轴向荷载的作用下试样会沿着内部的倾斜节理变形。在轴向循环加卸载作用下，轴向应力和轴向位移的关系可以用下图表示：
OA段加载：应力应变为弹性关系，试样发生弹性变形，沿着内部节理滑移；
OB段初始卸载：节理不产生滑移；
BO段最终卸载：弹性变形恢复。
试样材料的力学参数如下：
节理的部分力学参数如下：
使用UDEC模拟过程中，节理模型采用了库伦滑移模型以及连续屈服模型。
这个问题有相关研究者提出的概念模型，并推导了概念模型解答的表达式，感兴趣的可以去手册看，这里不做介绍。
2.UDEC模型
UDEC模型假设这个问题是平面问题，且只有中间的节理允许滑动，节理末端不允许滑动。试样的属性设置为线弹性，各向同性，试样模型使用常应变三角形有限差分网格划分，如下图：
主要研究两种节理本构模型下的问题：库伦滑移节理模型及连续屈服模型。模型参数如下：
库仑滑移模型是线弹性-全塑性本构关系，并且对应于概念模型中三个刚度的表达式。 连续屈服节理模型是非线性的，因此不符合用于概念模型的开发。 通过将初始摩擦角（jif）设置为等于固有摩擦角（jfric），并将节理粗糙度参数（jr）设置为较小的值，可以使模型近似于库仑滑移模型。 还必须将构造节理部分的弹性刚度设置为与节理滑动部分的弹性刚度一致：
set join_jkn 2e10 join_jks 2e10请注意，默认情况下，根据周围可变形块区域的刚度确定构造节理的弹性刚度（请参见 User’s Guide 中的3.2.3节）
FISH函数slip_load用于计算模型的轴向应力和轴向位移，并将结果存储在表1中（table 1）。 在FISH函数ana_slip中，使用等式中给出的刚度计算轴向应力与位移的关系，并存储在表2中（table 2），以便与表1的结果进行比较。FISH函数包含在文件“ SLIP.FIS”中。
3.结果分析
对于具有内部倾斜滑动裂纹的弹性体问题，没有任何解析解，因为裂纹两端的应力条件非常复杂。 但是，这里描述的简单概念模型体现了循环载荷中观察到的问题的基本特征（即三个明显不同的整体刚度）。
UDEC中两个节理本构模型的结果都与概念模型的结果进行了比较，如图所示。
由于概念模型的简化假设，因此不能期望计算出的斜率与从概念模型确定的斜率紧密一致。 因为概念模型假设正应力在裂纹和弹性延伸上均匀分布，并且随着滑移裂纹长度的增加，应力集中（尤其是节理的法向应力）变得更加明显。
用库仑摩擦定律建模的带有滑动裂纹的试样在循环荷载作用下轴向应力与轴向位移的关系图：
用连续屈服节理模型建模的带有滑动裂纹的试样在循环荷载作用下轴向应力与轴向位移的关系图：
命令流文件点击下方👇链接查看：
这里主要解释一下slip.dat命令流以及slip.fisFISH文件的命令含义，另外一个文件只是节理模型不同，可以对照着看：
在slip_dat文件中，首先定义了一个函数const来给模型的力学参数赋值，这也是常用到的一种方法，这些赋值后的变量可以直接嵌入命令或者FISH函数中，后期参数有变化直接更改这个函数中的变量值即可，而不用把所有涉及到这些参数的地方统统改一遍。定义了函数以后，使用函数名const来调用这个函数完成变量赋值，以便给slip.fis中函数的变量使用。
call slip.fis即调用写好的slip.fisFISH文件，这个文件中定义了两个函数，分别为ana_slip和slip_load。ana_slip是进行概念模型解的计算，具体可以参照手册中的公式，该函数中，degrad是预定义好的变量，其值为∏/180，即将读书转化为弧度以便进行三角运算。float(ntab2)是将ntab2进行强制类型转换，转换为浮点型进行计算。xtable(2,ntab2)=ax_dis和ytable(2,ntab2)=ax_str是把变量ax_dis和ax_str的值存到ID号为2，序列号为ntab2的表中。slip_load函数是类似的。
之后，slip.dat中开始建立块体模型，生成单一的贯穿节理，划分网格：
round 0.001
edge 0.002
block 0,0 0,2 1,2 1,0
crack (0,0.5) (1,1.5)
gen edge 0.15下面赋予节理模型的几条命令：
;设置所有节理的模型为摩尔库伦，并设置参数(可以看到粘聚力被设置为一个非常大的值，这也意味着节理将不能滑动即为构造节理)
joint model area jks 2E10 jkn 2E10 jcoh 1e20 jfric 45
;将矩形范围为（0.3，0.7）到（0.74，1.28）中包含的节理分组，命名为“joint”
group joint 'joint' range 0.3,0.7 0.74,1.28
;使用range关键字将“joint”组的节理重新赋予参数，允许其产生滑动
joint model area jks 2E10 jkn 2E10 jfriction 16 range group 'joint'
; new contact default，即设置运算过程中新接触的默认属性
set jcondf joint model area jks=2E10 jkn=2E10 jfriction=16之后设置记录的变量和边界条件，保存模型：
;记录靠近坐标(0.5,2.0)的点的y方向位移以及靠近坐标(0.5,2.0)的单元的syy应力
history ydisplace 0.5,2.0
history syy 0.5,2.0
;fix the bottom boundary，即固定模型底部
boundary yvelocity 0 range -0.1 1.1 -0.1 0.1
;保存模型
save slip.sav后面就是在模型顶部轴向循环加卸载的过程：
;load step 1
;加载，运行400步
boundary yvelocity -4.43E-2 range -0.1 1.1 1.9 2.1
cycle 400
;加载停止，求解平衡
boundary yvelocity -0 range -0.1 1.1 1.9 2.1
solve force 10.0
;调用slip.fis中的函数slip_load进行计算应力应变
slip_load以上所述均为个人理解，如有错误欢迎指出！