真三轴试验和节理建模

离散元方法与实践刘春南京大学7.2真三轴试验和节理建模本节课程我们将学习MatDEM中真三轴试验的实现和节理建模。概述

高应力作用下的真三轴试验非常困难，离散元数值模拟提供一个简单高效的测试方法。我们提供的示例代码3DJointStress演示了如何建立三轴试验模拟器，并进行高压条件下裂隙岩石变形和破坏的数值模拟。这个示例的第一和第二步通过常规方法堆积了一个岩石块体；第三步演示了如何向这个块体中增加复杂的节理和裂隙，并施加真三轴压力。运行前两步代码，得到图示模型定义裂隙面模型成功建立后，可通过以下多种方式来添加节理和裂隙面（3DJointStress3

），主要包括。1）利用三点坐标，2）利用多个点坐标，（不一定要在同一个面上）3）利用Tool_Cut工具，这三种方式需要注意的是，由于设置材料后单元力学性质突变模型需要重新平衡，因此设置节理和裂隙步骤应在第二步模型平衡之后，或在第三步施加新荷载前。通常在第三步开始时执行节理定义操作会较方便。利用三角形函数定义裂隙面在第三步代码中（3DJointStress3

），首先来看如何利用三角面顶点坐标来生成裂隙面。这种方法的主要思想为：通过指定的三个点的坐标，在空间中生成一个三角形的平面；寻找与该平面相接触的单元连接，并进行断开或胶结操作。TriX=[0.01,0.08,0.02];TriY=[0.05,0.1,0.05];TriZ=[0,0.05,0.1];bondFilter=mfs.setBondByTriangle(d,TriX,TriY,TriZ,type);d.mo.zeroBalance();d.showFilter('Group',{'sample'});d.show('Crack');首先，在空间中定义三角面三个顶点的X，Y和Z坐标(第1行)。利用mfs.setBondByTriangle函数，根据三个点的坐标，在模型对象d中产生一个三角面，并返回与三角面接触的单元连接的布尔矩阵bondFilter(第2行)

。TriX=[0.01,0.08,0.02];TriY=[0.05,0.1,0.05];TriZ=[0,0.05,0.1];bondFilter=mfs.setBondByTriangle(d,TriX,TriY,TriZ,type);d.mo.zeroBalance();d.showFilter('Group',{'sample'});d.show('Crack');在bondFilter中，与三角面相接触的单元连接值被记为true（1），其余值为false（0）（给出bondFilter的截图@黄）。参数type为计算操作类型，可声明为进行胶结（‘glue’）或者生成裂隙（‘break’）。通过运行这些代码，将生成图中所示的裂隙面。TriX=[0.01,0.08,0.02];TriY=[0.05,0.1,0.05];TriZ=[0,0.05,0.1];bondFilter=mfs.setBondByTriangle(d,TriX,TriY,TriZ,type);d.mo.zeroBalance();d.showFilter('Group',{'sample'});d.show('Crack');三角形切割出的微裂隙面；基于bondFilter可以获得这些连接对应的单元。首先，将连接压缩状态矩阵cFilter和连接胶结状态矩阵bFilter取并集，得到连接接触矩阵再与bondFilter取交集，则得到与三角面接触的连接connectFilter（第2行）；bondFilter=bondFilter&(d.mo.cFilter|d.mo.bFilter);connectFilter=sum(bondFilter,2)>0;connectId=find(connectFilter);d.addGroup('JointLayer',connectId);%addanewgroupfigure;d.showFilter('Group','JointLayer','aR');

然后利用find命令得到与三角面接触的单元编号，定义为'JointLayer'组，并显示（第3-5行）

三角面连接的单元。bondFilter=bondFilter&(d.mo.cFilter|d.mo.bFilter);connectFilter=sum(bondFilter,2)>0;connectId=find(connectFilter);d.addGroup('JointLayer',connectId);%addanewgroupfigure;d.showFilter('Group','JointLayer','aR');

裂隙面连接的单元

裂隙面的多边形定义函数的原理与三角面定义函数类似。函数setBondByPolygon将多个点构成的空间多边形分解成若干个三角形这些三角面均以第一个顶点为起点，顺序生成相邻三角面。利用多边形函数定义裂隙面TriX=[0.09,0.09,0.1,0.01];TriY=[0,0.1,0,0.1];TriZ=[0,0,0.1,0.1];bondFilter=mfs.setBondByPolygon(d,TriX,TriY,TriZ,'glue');d.mo.nBondRate(bondFilter)=2;%makestrongjoint运行这些代码，即可生成所需节理面（图右侧节理面）。连接强度系数BondRate

利用Tool_Cut定义复杂的节理面通过上述两种方法可以在模型中生成简单的节理面。但对于复杂三维空间裂隙和节理面，需要用一系列的三角面来表示显然无法使用函数setBondByTriangle或setBondByPolygon来逐一设置三角面的连接。对此，MatDEM提供了更加高效模型的切割工具Tool_Cut。Tool_Cut的对象通常命名为C，主要用于利用三维面来定义组，以及节理面相关操作等。该工具能利用输入的一系列三角面来定义三维节理面；也可以根据离散的点自动生成一系列三角面，以拟合出复杂的三维节理面。以下命令实现了利用Tool_Cut来定义三角形节理面首先，利用已赋值的参数定义了两个三角面的顶点（1-3行）。TriX2=[TriX(1:3);TriX([1,3,4])]-0.05;TriY2=[TriY(1:3);TriY([1,3,4])];TriZ2=[TriZ(1:3);TriZ([1,3,4])];C=Tool_Cut(d);C.setTriangle(TriX2,TriY2,TriZ2);bondFilter=C.setBondByTriangle(type);d.mo.nBondRate(bondFilter)=0.5;d.show('BondRate');C.showTriangle();基于模型对象d创建切割对象C

（第4行），并通过C.setTriangle函数将三角面顶点坐标赋给C的三角面属性TriangleX、Y和Z（第5行）。查看TriX2的数据，可以看到此矩阵包括两行数据，其中每一行数据表示一个三角面的三个顶点的X坐标。因此将生成两个三角面。TriX2=[TriX(1:3);TriX([1,3,4])]-0.05;TriY2=[TriY(1:3);TriY([1,3,4])];TriZ2=[TriZ(1:3);TriZ([1,3,4])];C=Tool_Cut(d);C.setTriangle(TriX2,TriY2,TriZ2);bondFilter=C.setBondByTriangle(type);d.mo.nBondRate(bondFilter)=0.5;%makeweakjointd.show('BondRate');C.showTriangle();C.setTriangle函数的输入三角面坐标数据

在执行C.setBondByTriangle函数时，程序会逐行读取C.TriangleX、Y、Z矩阵中记录的三角面坐标，根据输入参数type，来确定对连接进行胶结还是断开操作。TriX2=[TriX(1:3);TriX([1,3,4])]-0.05;TriY2=[TriY(1:3);TriY([1,3,4])];TriZ2=[TriZ(1:3);TriZ([1,3,4])];C=Tool_Cut(d);C.setTriangle(TriX2,TriY2,TriZ2);bondFilter=C.setBondByTriangle(type);d.mo.nBondRate(bondFilter)=0.5;%makeweakjointd.show('BondRate');C.showTriangle();进一步，将函数返回参数bondFilter用来设置连接强度矩阵nBondRate，将其三角面对应的连接的强度设为原始值的0.5倍（第7行）。TriX2=[TriX(1:3);TriX([1,3,4])]-0.05;TriY2=[TriY(1:3);TriY([1,3,4])];TriZ2=[TriZ(1:3);TriZ([1,3,4])];C=Tool_Cut(d);C.setTriangle(TriX2,TriY2,TriZ2);bondFilter=C.setBondByTriangle(type);d.mo.nBondRate(bondFilter)=0.5;%makeweakjointd.show('BondRate');C.showTriangle();随后，我们可以通过命令C.showTriangle可以查看C中存储的三角面也可以通过后处理命令d.show('BondRate')显示连接强度系数图。TriX2=[TriX(1:3);TriX([1,3,4])]-0.05;TriY2=[TriY(1:3);TriY([1,3,4])];TriZ2=[TriZ(1:3);TriZ([1,3,4])];C=Tool_Cut(d);C.setTriangle(TriX2,TriY2,TriZ2);bondFilter=C.setBondByTriangle(type);d.mo.nBondRate(bondFilter)=0.5;%makeweakjointd.show('BondRate');C.showTriangle();Tool_Cut工具能利用离散的点来生成复杂的面，包括数字高程数据或GPS记录的地面坐标。利用离散元生成裂隙面以下代码随机生成10个离散坐标，并建立三角面构建的曲面。其中，fs.randSeed(seed)首先确定了随机数种子，如果输入参数不变，程序每次都会生成相同的三角面，修改seed值则可以获得不同的三角面（第1行）；fs.randSeed(2);d.mo.bFilter(:)=false;d.mo.zeroBalance();randN=10;PX=rand(randN,1)*0.1;PY=rand(randN,1)*0.1;PZ=0.05+rand(randN,1)*0.04;C.addSurf(PX,PY,PZ);C.getSurfTri(1,1);C.getTriangle(1);然后，利用rand函数生成了一系列随机点坐标，包括X、Y、Z（第5行）。并基于这些离散点来生成三角面，具体包括为：1.通过C.addSurf命令导入离散点数据（第6行）；2.通过C.getSurfTri命令获得离散元三角面信息（第7行）；3.通过C.getTriangle命令获得三角面顶点坐标（第8行）三个步骤。fs.randSeed(2);d.mo.bFilter(:)=false;d.mo.zeroBalance();randN=10;PX=rand(randN,1)*0.1;PY=rand(randN,1)*0.1;PZ=0.05+rand(randN,1)*0.04;C.addSurf(PX,PY,PZ);C.getSurfTri(1,1);C.getTriangle(1);其中，C.addSurf命令基于Matlab内置函数scatteredInterpolant函数，对散点数据进行插值。运行C.addSurf(PX,PY,PZ)命令后，离散点的坐标信息将会被储存在C.Surf{1,1}对象中（图）。可以看到，C.Surf是一个cell数组，C对象可以多次运行C.addSurf命令，并导入多个面。离散点的X、Y坐标被储存在Surf{1,1}.Points里，Z坐标被储存在Surf{1,1}.Values里（图）。但是，C.Surf中存储的只是离散点的坐标，需要通过命令C.getSurfTri(m,rate)函数将离散点坐标转化成三角面信息。这个函数有两个输入参数，第一个输入参数指明使用C.Surf数组中的第m组数据来生成三角面信息，并将其存储在C.SurfTri的相应位置。SurfTri记录了导入点的X、Y、Z坐标信息（图），而其中的Tri矩阵则逐行记录了生成的三角面的顶点的编号。如图所示，Tri矩阵大小为13×3，表示共生成了13个三角面。。C.SurfTri{1,1}中的三角面信息

最后，利用已有的三角面信息SurfTri，通过C.getTriangle(1)命令获得第一组三角面的顶点坐标数据，即C.TriangleX、Y、Z。最后，使用后处理命令显示试样中的胶结连接，以及使用C.showTriangle()命令来显示三角面，得到图示的节理面和相关胶结图。由离