离散元方法与实践慕课matdem-第九章2段9.2矿山斜坡

离散元方法与实践刘春南京大学9.1矿山斜坡爆破0概述本节将介绍矿山爆破示例。该模拟基于BoxModel示例代码，这个代码最初用于二维和三维滑坡数值模拟。本示例需先运行BoxModel1和BoxModel2代码来生成斜坡，再利用Model3Exploisive代码实现矿山爆破过程的数值模拟。爆破模拟的基本原理为：搜索出一定区域内的单元，定义为爆炸点，增加这些单元的半径，使单元间相互重叠，以产生巨大压力和弹性应变能。在数值模拟过程中，爆破点的单元在内压力作用下，迅速向外运动，释放应变能，并模拟爆破作用。1建立斜坡模型BoxModel1与其他示例的第一步相同，在此不作介绍。本节直接使用第一步得到的二维堆积模型。在加载第一步生成的堆积模型数据后，BoxModel2的主要建模代码分为导入数据并切割模型、以及设置材料和平衡模型两个部分1.1导入数据并切割模型首先创建Tool_Cut对象C（第1行）；并采用Tool_Cut对象C和addSurf命令，读入记事本文档中的高程数据（第2、3行）；C=Tool_Cut(d);%cutthemodellSurf=load('slope/layersurface.txt');%loadthesurfacedataC.addSurf(lSurf);%addthesurfacestothecutC.setLayer({'sample'},[1,2,3,4]);%setlayersaccordinggeometricaldatagNames={'lefPlaten';'rigPlaten';'botPlaten';'layer1';'layer2';'layer3'};d.makeModelByGroups(gNames);为方便查看层面，这些数据记录于“滑坡高程.xls”的BoxModel数据表中，具体数据和格式如表中所示，相应的excel折线表如图所示。在Excel中设定完高程数据后，再将其复制到txt文档中，以方便MatDEM读取。折线号surface0surface1surface2surface3坐标xyzxyzxyzxyz高程数据1001029103310434004026403040427007022702670421000100161002010038………………………………高程数据.xls文件中的折线表Excel表中用于生成折线图的矩阵数据在图中，我们只定义了水平方向0-80米处的高程，MatDEM会自动根据已有的高程数据向外延伸，切分出80-100米的地层（图b中）。a高程数据.xls文件中的折线表b由折线表数据生成的斜坡模型然后，通过直接设定材料性质的方法来生成两种材料Soil1和Soil2，得到材料数组Mats（第1-8行）。%----------setmaterialofmodelmatTxt=load('Mats\Soil1.txt');Mats{1,1}=material('Soil1',matTxt,B.ballR);Mats{1,1}.Id=1;matTxt2=load('Mats\Soil2.txt');Mats{2,1}=material('Soil2',matTxt2,B.ballR);Mats{2,1}.Id=2;d.Mats=Mats;%---------assignmaterialtolayersandbalancethemodeld.setGroupMat('layer2','Soil2');d.groupMat2Model({'sample','layer2'});d.balanceBondedModel();1.2设置材料和平衡模型如表所示，Soil2的强度只有Soil1的十分之一，且内摩擦系数也较小，将作为软弱层。由于单元默认的材料号为1，只需将软弱层（layer2）的材料设为Soil2，并使用groupMat2Model将材料性质赋给单元。最后平衡模型，得到如图所示的离散元模型，其中layer2为软弱层。材料名称杨氏模量/MPa泊松比单轴拉伸强度/kPa单轴压缩强度/kPa内摩擦系数密度/kg/m3Soil1200.14202000.81900Soil2100.182200.62000Soil1和Soil2两种材料的输入力学性质

2设置爆破点和爆破能量模型建好后，进入第三步代码Model3Exploision。由于第二步默认的材料力学性质较弱，为保证边坡稳定，此处先将单元的断裂力（aBF）和初始抗剪力（aFS0）增加为10倍（第2、3行）。…B.name=[B.name'Exploision'];d.mo.aBF=d.mo.aBF*10;d.mo.aFS0=d.mo.aFS0*10;然后通过代码定义爆破点在斜坡中的位置和半径，搜索出相应的单元编号bombId，并将其定义为Bomb1组。centerX=15;centerZ=20;bombR=2;dX=d.mo.aX-centerX;dZ=d.mo.aZ-centerZ;bombId=find((dX.*dX+dZ.*dZ)<bombR.*bombR);%gettheIdofbombd.addGroup('Bomb1',bombId);%addanewgroupd.mo.zeroBalance();d.recordStatus();我们将爆破点的单元组号定义为1，利用d.show函数显示出爆破点的位置（图）；然后计算平均应力Stress，并保存于d.data.Stress中，通过d.show(‘Stress’)显示结果，实现了自定义后处理显示。d.setData();d.data.groupId(d.GROUP.Bomb1)=1;d.show('groupId');d.data.Stress=(abs(d.data.StressXX)+abs(d.data.StressZZ));d.show('Stress');爆破点单元示意图

（a）组号图；

（b）自定义平均应力图

在获得爆破单元的基础上，通过增加其半径来产生爆破能量。为此，先将单元半径增加40%，并记录新的系统状态。oldKe=d.status.elasticEs(end);%recordtheoriginalenergybombExpandRate=1.4;d.mo.aR(bombId)=d.mo.aR(bombId)*bombExpandRate;%increasebombelementsized.mo.zeroBalance();d.recordStatus();newKe=d.status.elasticEs(end);dKe=newKe-oldKe;%calculatetheenergyincrementfs.disp(['Energyofthebombis'num2str(dKe)'J','~='num2str(dKe/4.2e6)'KgTNT']);可以看到，爆破点处单元的应力非常大。为了计算爆破能量，利用oldKe记录原始的弹性应变能（第1行），用newKe记录单元半径增加后的弹性应变能（第6行），计算得到由于爆破点单元半径增加而增加的弹性应变能。最后一条命令，将能量转化为TNT当量，以量化爆破能量。运行代码后，会得到以下信息Energyofthebombis17973262.3911J~=4.2793KgTNT即模拟的爆破能量相当于4.3公斤的TNT炸药。需要注意的是，在本示例中，为了较清晰地显示爆破点的位置，采用了较大的爆破点半径。在实际应用中，可根据需要设置较小的区域半径，并通过增加爆破点单元的刚度来增加爆破能量。同时，也可以通过给单元施加较大初速度的方法来增加能量，等。3迭代计算和模拟结果进一步，可以开始迭代计算。由于本例为动力作用，需要将活动单元的阻尼系数（d.mo.mVis）设小一些，相当于模拟空气阻力（1、2行）。同时，运行d.mo.setGPU(‘auto’)命令，测试GPU和CPU的速度，选择更快的一个，用于后续计算（第3行）

。visRate=0.00001;d.mo.mVis=d.mo.mVis*visRate;gpuStatus=d.mo.setGPU('auto');在设置好模型后，通过循环完成数值模拟：totalCircle=20;d.tic(totalCircle);fName=['data/step/'B.namenum2str(B.ballR)'-'num2str(B.distriRate)'loopNum'];save([fName'0.mat']);fori=1:totalCircled.mo.setGPU(gpuStatus);d.balance('Standard',0.1);d.clearData(1);save([fNamenum2str(i)'.mat']);d.calculateData();d.toc();%showthenoteoftimeendd.show('mV');这个模拟分为20次循环（totalCircle），每次循环进行0.1次标准平衡（第7行）共模拟真实世界时间5s。图a为0.1s时的单元速度场图，在图中可以清晰地看到，应力波以爆破点为中心向外传播，并在坡面形成较大的单元速度场（40m/s）。单元速度场分布