Exploision

Step 1Step 2Step 3

%build the geometrical model
clear;
fs.randSeed(2);%build random model
B=obj_Box;%build a box object
B.name='BoxModel';
B.GPUstatus='auto';
B.ballR=1;
B.isClump=0;
B.distriRate=0.2;
B.sampleW=100;
B.sampleL=0;
B.sampleH=80;
%B.BexpandRate=4;
%B.PexpandRate=4;
B.type='topPlaten';
%B.type='TriaxialCompression';
B.setType();
B.buildInitialModel();%B.show();

B.setUIoutput();

d=B.d;%d.breakGroup('sample');d.breakGroup('lefPlaten');

%--------------end initial model------------
B.gravitySediment();
B.compactSample(1);%input is compaction time
%------------return and save result--------------
d.status.dispEnergy();%display the energy of the model

d.clearData(1);%clear dependent data
d.recordCalHour('BoxModel1Finish');
save(['TempModel/' B.name '1.mat'],'B','d');
save(['TempModel/' B.name '1R' num2str(B.ballR) '-distri' num2str(B.distriRate)  'aNum' num2str(d.aNum) '.mat']);
d.calculateData();
d.show('aR');

%set the material of the model
clear
load('TempModel/BoxModel1.mat');
B.setUIoutput();%set output of message
d=B.d;
d.calculateData();
d.mo.setGPU('off');
d.getModel();%get xyz from d.mo

%---------cut the model to make slope
C=Tool_Cut(d);%cut the model
lSurf=load('slope/layer surface.txt');%load the surface data
C.addSurf(lSurf);%add the surfaces to the cut
C.setLayer({'sample'},[1,2,3,4]);%set layers according geometrical data
gNames={'lefPlaten';'rigPlaten';'botPlaten';'layer1';'layer2';'layer3'};
d.makeModelByGroups(gNames);
%---------end cut the model to make slope

%----------set material of model
matTxt=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;
%----------end set material of model

%---------assign material to layers and balance the model
d.setGroupMat('layer2','Soil2');
d.groupMat2Model({'sample','layer2'});
d.balanceBondedModel();
%---------end assign material to layers and balance the model

%---------save the data
d.mo.setGPU('off');
d.clearData(1);
d.recordCalHour('BoxModel2Finish');
save(['TempModel/' B.name '2.mat'],'B','d');
save(['TempModel/' B.name '2R' num2str(B.ballR) '-distri' num2str(B.distriRate)  'aNum' num2str(d.aNum) '.mat']);
d.calculateData();

%note: run user_BoxModel1 and 2 before runing the code
clear;
load('TempModel/BoxModel2.mat');
B.setUIoutput();
d=B.d;
d.calculateData();
d.mo.setGPU('off');
B.name=[B.name 'Exploision'];
%---------enhance the strength of the slope
d.mo.aBF=d.mo.aBF*10;
d.mo.aFS0=d.mo.aFS0*10;

%-------make an exploisive point
centerX=15;centerZ=20;
bombR=2;
dX=d.mo.aX-centerX;
dZ=d.mo.aZ-centerZ;
bombId=find((dX.*dX+dZ.*dZ)<bombR.*bombR);%get the Id of bomb
d.addGroup('Bomb1',bombId);%add a new group
d.showFilter('Group',{'Bomb1'},'aR');
d.mo.zeroBalance();

d.recordStatus();
%-------end make an exploisive point
%-------increase the bomb element size
oldKe=d.status.elasticEs(end);%record the original energy
bombExpandRate=1.4;
d.mo.aR(bombId)=d.mo.aR(bombId)*bombExpandRate;%increase bomb element size
d.mo.zeroBalance();
d.recordStatus();
newKe=d.status.elasticEs(end);
dKe=newKe-oldKe;%calculate the energy increment

%calculate the TNT equivalent 
fs.disp(['Energy of the bomb is ' num2str(dKe) 'J', ' ~=' num2str(dKe/4.2e6) ' Kg TNT']);
%-------end increase the bomb element size

d.getModel();%reset the initial status of the model
d.resetStatus();%initialize model status, which records running information
%d.show();return;

%---------initialize the simulation
d.mo.isHeat=1;%calculate heat in the model
visRate=0.00001;
d.mo.mVis=d.mo.mVis*visRate;
gpuStatus=d.mo.setGPU('auto');
d.setStandarddT();

%----------start simulation
totalCircle=10;
d.tic(totalCircle);
fName=['data/step/' B.name  num2str(B.ballR) '-' num2str(B.distriRate) 'loopNum'];
save([fName '0.mat']);%return;
for i=1:totalCircle
    d.mo.setGPU(gpuStatus);
    d.balance('Standard',0.1);
    d.clearData(1);
    save([fName num2str(i) '.mat']);
    d.calculateData();
    d.toc();%show the note of time
end
%----------end simulation

%--------save data
d.mo.setGPU('off');
d.clearData(1);
d.recordCalHour('BoxModel3Finish');
save(['TempModel/' B.name '3.mat'],'B','d');
save(['TempModel/' B.name '3R' num2str(B.ballR) '-distri' num2str(B.distriRate)  'aNum' num2str(d.aNum) '.mat']);
d.calculateData();