离散元方法与实践慕课matdem第九章2段91陨石撞击地面

离散元方法与实践刘春南京大学9.1陨石撞击地面0概述本节主要讲解撞击作用的离散元模拟，以陨石撞击地面为例（BoxCrash），建立了二维数值模型。陨石撞击作用模拟的基本思路是先堆积地层，然后在地层之上一定距离建立圆形的陨石模型，赋予其一定的初速度，并撞击地面。通过少量的修改，这个示例可以模拟很多问题。如减小模型尺寸，提高陨石单元的刚度和胶结力（或使用clump），其可以模拟子弹击穿水泥板的过程；通过切割模型和定义新部件，其可以改为霍普金森杆试验模拟器，等。1堆积地层模型首先，在BoxCrash1代码中，通过常规的堆积建模操作，建立如图所示的地层模型。…B.name='BoxCrash';B.ballR=0.5;B.distriRate=0.2;B.sampleW=500;B.sampleL=0;B.sampleH=300;B.BexpandRate=4;B.type='topPlaten';B.setType();B.buildInitialModel();…B.gravitySediment();d.mo.aMUp(:)=0;B.compactSample(2);%inputiscompactiontime…堆积地层模型

2建立陨石模型在BoxCrash2代码中加载第一步建立的堆积模型，并利用过滤器将堆积模型上半部分单元（代码1-3行）。最后得到图示地层模型。mZ=d.mo.aZ(1:d.mNum);topLayerFilter=mZ>max(mZ)*0.5;d.delElement(find(topLayerFilter));删除模型上半部分单元后的地层模型

由于在模型中需要对陨石和地层设置不同的材料，共导入两种材料（1-7行）：

材料1为WeakRock，作为地层模型的材料；材料2为StrongRock，作为陨石模型的材料。…matTxt=load('Mats\WeakRock.txt');Mats{1,1}=material('WeakRock',matTxt,B.ballR);Mats{1,1}.Id=1;matTxt2=load('Mats\StrongRock.txt');Mats{2,1}=material('StrongRock',matTxt2,B.ballR);Mats{2,1}.Id=2;d.Mats=Mats;d.groupMat2Model({'sample'},1);…所加载的两种材料的性质如下表所示。注意，这里所采用的是直接输入的材料性质，基于这些性质设定单元参数，其堆积体力学性质可能与设定值有一定误差。最后，通过d.groupMat2Model({‘sample’},1)命令设置地层模型单元的材料（接上页8行）。材料名称杨氏模量/GPa泊松比单轴拉伸强度/MPa单轴压缩强度/MPa内摩擦系数密度/kg/m3WeakRock50.22200.62600StrongRock100.15101000.62850WeakRock和StrongRock两种材料的输入力学性质

完成以上两个步骤后，我们来建立陨石结构体，并将其导入模型中。其基本思想为：先在地层模型中间筛选出一个圆形区域的单元（第1~5行）；将其添加为一个组（第6行），并导出为陨石结构体（第7行）。sampleId=d.GROUP.sample;sX=d.aX(sampleId);sZ=d.aZ(sampleId);sR=d.aR(sampleId);discCX=mean(sX);discCZ=mean(sZ);discR=20;discFilter=(d.aX-discCX).^2+(d.aZ-discCZ).^2<discR^2;d.addGroup('Disc0',find(discFilter));discObj=d.group2Obj('Disc0');具体来看，首先得到样品单元（地层）的Id，并将其储存在sampleId矩阵中（第1行）；然后，将中地层单元X、Z坐标与单元半斤aR分别保存在变量sX、sZ与aR中（第2行）；并求出地层单元的X、Z坐标平均值保存在discCX和discCZ中，则坐标（discCX，discCZ）可以表示长方形地层的中心点（第3行）；sampleId=d.GROUP.sample;sX=d.aX(sampleId);sZ=d.aZ(sampleId);sR=d.aR(sampleId);discCX=mean(sX);discCZ=mean(sZ);discR=20;discFilter=(d.aX-discCX).^2+(d.aZ-discCZ).^2<discR^2;d.addGroup('Disc0',find(discFilter));discObj=d.group2Obj('Disc0');discR为选取圆形区域的半径，即陨石的半径为20m

（第4行）；以（discCX，discCZ）为圆心选取了半径20m的圆形区域，建立了过滤矩阵discFilter（第5行）；最后，通过addGroup命令，将discFilter选中的单元添加为Disc0组（第6行）；并通过group2Obj命令将Disc0组转化为结构体discObj（第7行），以供进一步使用。创建的陨石结构体如图所示。切割出的陨石结构体

进一步，通过addElement命令将陨石结构体discObj导入模型，材料为StrongRock，并通过addGroup命令将其添加为Disc组（第2行）；最后将Disc刚好移动到地面之上。discId=d.addElement('StrongRock',discObj);d.addGroup('Disc',discId);disZ=max(sZ+sR)-min(discObj.Z-discObj.R);d.moveGroup('Disc',0,0,disZ);d.balanceBondedModel0();得到如图所示的最终模型。地层与陨石模型

3陨石撞击过程模拟在第三步代码BoxCrash3中，实现了陨石撞击地面过程的模拟。首先加载第二步的计算数据，并初始化模型，然后设置相关的模拟参数d.mo.isHeat=1;visRate=0.0001;d.mo.mVis=d.mo.mVis*visRate;discId=d.GROUP.Disc;d.mo.mVZ(discId)=-1000;d.setStandarddT();d.mo.dT=d.mo.dT*0.05;由于该模拟为动力作用数值模拟，阻尼系数需按实际情况取较小值，为此，该处将d.mo.mVis设为最优阻尼的1/10000（2、3行）；将陨石单元的速度设为-1000m/s（4、5行），即以1km/s的速度向下撞击地面；最后，适当减小计算时间步，设为标准时间步的1/20（6、7行）。当模拟参数设置完成后，可通过for循环进行撞击作用的数值计算：这里总共设置了40次循环，并保存40个中间过程文件（totalCircle）；每次循环计算进行0.4次标准平衡（第4行）；在关闭GPU计算并压缩数据后，通过save命令保存计算数据。totalCircle=40;fori=1:totalCircled.mo.setGPU(gpuStatus);d.balance('Standard',0.4);d.mo.setGPU('off');d.clearData(1);save([fNamenum2str(i)'.mat']);d.calculateData();end计算完成后，可通过后处理得到撞击过程的各类场图及动画，该模型包含41万单元，可精细地展现应力波的产生和传播。当前图示显示了在0.42秒时间内应力分布StressZZ的变化。从图中可以看到，当陨石接触地面时，二者间产生巨大的应力作用，应力波向四周传播，并在遇到下边界时发生反射。同时，在地面形成巨大陨石坑。图撞击作用过程图（Stre