SandPenetration

Step 1Step 2Step 3

%step1: packing the elements
clear;
fs.randSeed(1);%random model seed, 1,2,3...
B=obj_Box;%declare a box object
B.name='Sand_Penetration';

%--------------initial model------------
B.GPUstatus='auto';%program will test the CPU and GPU speed, and choose the quicker one
B.ballR=0.001;
B.isShear=1;
B.isClump=0;%if isClump=1, particles are composed of several balls
B.distriRate=0.2;%define distribution of ball radius, 
B.sampleW=0.25;%width, length, height, average radius
B.sampleL=0;%when L is zero, it is a 2-dimensional model
B.sampleH=0.06;
B.BexpandRate=4;%boundary is 4-ball wider than sample
B.PexpandRate=0;
% B.type='topPlaten';%add a top platen to compact model
B.isSample=1;
B.setType();
B.buildInitialModel();;
d=B.d;

%---------- gravity sedimentation
B.gravitySediment();
d.mo.aMUp(:)=1;   

%------------return and save result--------------
d.status.dispEnergy();
d.mo.setGPU('off');
d.clearData(1);%clear dependent data
d.recordCalHour('Step1Finish');
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();

%set the material of the model
clear
load('TempModel/Sand_Penetration1.mat');
B.setUIoutput();%set the output
d=B.d;
d.calculateData();%calculate data
d.mo.setGPU('off');%close the GPU calculation
d.getModel();%get xyz from d.mo

%---------------delele the box boundary and topelement
d.delElement(d.GROUP.lefB);
d.delElement(d.GROUP.rigB);
d.delElement(d.GROUP.botB);
d.delElement(d.GROUP.topB);
mZ=d.mo.aZ(1:d.mNum);%get the Z of elements
topLayerFilter=mZ>max(mZ)*0.8;
d.delElement(find(topLayerFilter));%delete elements according to id
mX=d.mo.aX(1:d.mNum);%get the Z of elements
rigLayerFilter=mX>max(mX)*0.8;
d.delElement(find(rigLayerFilter));%delete elements according to id
B=mfs.move(B,0,0,-0.024);

%--------------assign new material
matTxt=load('Mats\Soil1.txt');
Mats{1,1}=material('Soil1',matTxt,B.ballR);
Mats{1,1}.Id=1;
matTxt2=load('Mats\Brass.txt');
Mats{2,1}=material('Brass',matTxt2,B.ballR);
Mats{2,1}.Id=2;
d.Mats=Mats;
d.setGroupMat('sample','Soil1');
d.groupMat2Model({'sample'},2);

%----------make ProjectileObj------------
ProjectileObj=mfs.denseModel(0.8,@mfs.makeDisc,0.006,0.001);
ProjectileObj=mfs.rotate(ProjectileObj, 'YZ', 90);
ProjectileObj.Y(:)=0;
ProjectileObj=mfs.move(ProjectileObj,0.2,0,0.018);

%----------import bulletObj to model------------
ProjectileId=d.addElement(2,ProjectileObj);
d.addGroup('Projectile',ProjectileId);
d.setClump('Projectile');
d.mo.zeroBalance();
d.balance('Standard');
d.balanceBondedModel();%赋材料后，自动胶结和平衡模型，有单元间摩擦力

d.show('aR');
d.clearData(1);%clear dependent data
d.recordCalHour('Step2Finish');
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();

clear;
load('TempModel/Sand_Penetration2.mat');
d.calculateData();
d.mo.setGPU('off');
B.setUIoutput();
d=B.d;
d.getModel();
d.status=modelStatus(d);

d.mo.isHeat=1;
visRate=0.0001;
d.mo.mVis=d.mo.mVis*visRate;
ProjectileId=d.GROUP.Projectile;
d.mo.mVX(ProjectileId)=-200;
d.setStandarddT();

gpuStatus=d.mo.setGPU('auto');

totalCircle=20;%default value is 500, use 20 to increase speed
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(50,2);
    d.figureNumber=1;
    d.showB=0;
    d.show('mV');
    axis([0 0.25 -1 1 -0.05 0.1]);
    axis manual

    d.mo.setGPU('off');
    d.clearData(1);
    save([fName num2str(i) '.mat']);
    d.calculateData();
    d.toc()
end
d.mo.setGPU('off');
d.clearData(1);
d.recordCalHour('Step3Finish');
save(['TempModel/' B.name '3.mat'],'d');
save(['TempModel/' B.name '3R' num2str(B.ballR) '-' num2str(B.distriRate)  'aNum' num2str(d.aNum) '.mat']);
d.calculateData();
d.show('aR');