跳转至

BoxShear

user_BoxShear0.m
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
%make a box model, which will be put on a slope
clear;
fs.randSeed(1);%build random model
B=obj_Box;%build a box object
%B.name='BoxShear';
B.name='BoxShear';
B.SET.shearType='torsional';%may change it to shear
B.GPUstatus='auto';
B.ballR=0.005;
B.isClump=0;
B.distriRate=0.2;
B.SET.sampleR=0.0309;
B.SET.sampleH=0.02;
if strcmp(B.SET.shearType,'torsional')
    B.SET.sampleR=0.0309;
    B.SET.sampleH=0.06;
end

B.sampleW=B.SET.sampleR*2.2;
B.sampleL=B.SET.sampleR*2.2;
B.sampleH=B.SET.sampleH*1.5;
%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(2);%input is compaction time
mfs.reduceGravity(d,10);%reduce the gravity of element

%------------return and save result--------------
d.status.dispEnergy();%display the energy of the model
d.show('aR');

d.mo.setGPU('off');
d.clearData(1);%clear dependent data
d.recordCalHour('Step0Finish');
save(['TempModel/' B.name '0.mat'],'B','d');
d.calculateData();%because data is clear, it will be re-calculated
user_BoxShear1.m
  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
clear
load('TempModel/BoxShear0.mat');
sampleR=B.SET.sampleR;
sampleH=B.SET.sampleH;
ballR=B.ballR;
Rrate=0.7;
tubeR=sampleR+ballR;
ringWidth=0.005;
innerWidth=sampleR*0.9;
innerHeight=sampleH/6;
botTubeH=sampleH/2+ballR*2;
boxType='wall';%boxType could be 'model' or 'wall', see element type in help

boxSampleId=d.getGroupId('sample');
sX=d.mo.aX(boxSampleId);sY=d.mo.aY(boxSampleId);sZ=d.mo.aZ(boxSampleId);
dipD=0;dipA=0;radius=sampleR-ballR;height=sampleH-ballR;
columnFilter=f.run('fun/getColumnFilter.m',sX,sY,sZ,dipD,dipA,radius+B.ballR,height);
d.addGroup('column',find(columnFilter));
sampleObj=d.group2Obj('column');

sampleObj=mfs.moveObj2Origin(sampleObj);
%------------------------make object of shearBox-----------------

botTubeObj=mfs.denseModel(Rrate,@mfs.makeTube,tubeR+(1-Rrate)*ballR*2,botTubeH,ballR);
discObj=mfs.denseModel(Rrate,@mfs.makeDisc,sampleR+(1-Rrate)*ballR*1,ballR);
botTubeObj=mfs.moveObj2Origin(botTubeObj);
discObj=mfs.moveObj2Origin(discObj);
[botTubeObj,botDiscObj]=mfs.alignObj('bottom',botTubeObj,discObj);
botBoxObj=mfs.combineObj(botTubeObj,botDiscObj);

botBoxObj=mfs.align2Value('bottom',botBoxObj,0);
botBoxTopZ=mfs.getObjEdge('top',botBoxObj);
botRingObj=mfs.makeRing2(tubeR+ballR-ballR*(1-Rrate)*2,tubeR+ballR+ringWidth,ballR,Rrate);
botRingObj=mfs.align2Value('top',botRingObj,botBoxTopZ);
botBoxObj=mfs.combineObj(botBoxObj,botRingObj);

if strcmp(B.SET.shearType,'torsional')
    planeObj=mfs.makeBox(innerWidth,innerHeight,ballR,ballR);
    planeObj=mfs.rotate(planeObj,'YZ',90);
    planeObj=mfs.move(planeObj,sampleR-innerWidth,0,ballR*2);
    planeObj=mfs.rotateCopy(planeObj,60,6);
    botBoxObj=mfs.combineObj(botBoxObj,planeObj);
end

topTubeObj=mfs.align2Value('bottom',botTubeObj,botBoxTopZ);
[topTubeObj,topPlatenObj]=mfs.alignObj('top',topTubeObj,botDiscObj);
topTubeBotZ=mfs.getObjEdge('bottom',topTubeObj);
topRingObj=mfs.align2Value('bottom',botRingObj,topTubeBotZ);
topTubeRingObj=mfs.combineObj(topTubeObj,topRingObj);

%------------------------end make object of shearBox-----------------
%make a big box for the simulation
fs.randSeed(1);%random model seed, 1,2,3...
B=obj_Box;%declare a box object
B.name='BoxShear';
%--------------initial model------------
B.GPUstatus='auto';%program will test the CPU and GPU speed, and choose the quicker one
B.ballR=ballR;
B.isShear=0;
B.isClump=0;%if isClump=1, particles are composed of several balls
B.distriRate=0.2;%define distribution of ball radius, 
B.sampleW=(sampleR+ringWidth)*2.5;
B.sampleL=(sampleR+ringWidth)*2.5;
B.sampleH=sampleH*1.2;
B.BexpandRate=2;%boundary is 4-ball wider than 
B.PexpandRate=0;
B.type='botPlaten';
B.isSample=0;
%B.type='TriaxialCompression';
B.setType();
B.SET.boxType=boxType;
B.buildInitialModel();
d=B.d;
d.mo.setGPU('off');
B.SET.sampleR=sampleR;
B.SET.sampleH=sampleH;

botBoxId=d.addElement(1,botBoxObj,boxType);
d.addGroup('botBox',botBoxId);%add a new group
d.setClump('botBox');
topTubeRingId=d.addElement(1,topTubeRingObj,boxType);
d.addGroup('topTubeRing',topTubeRingId);%add a new group
d.setClump('topTubeRing');
topPlatenId=d.addElement(1,topPlatenObj);
d.addGroup('topPlaten',topPlatenId);%add a new group
d.setClump('topPlaten');
d.addGroup('shearBox',[d.GROUP.botBox;d.GROUP.topTubeRing;d.GROUP.topPlaten]);
d.delElement('botPlaten');

d.moveGroup('shearBox',B.sampleW/2,B.sampleL/2,0);

boxSampleId=d.addElement(1,sampleObj);
d.addGroup('sample',boxSampleId);%add a new group
boxZ=d.mo.aZ(d.GROUP.shearBox);
d.moveGroup('sample',B.sampleW/2,B.sampleL/2,mean(boxZ));
d.minusGroup('sample','shearBox',0.5);

%d.showFilter('Group',{'shearBox'},'aR');return;

d.removeGroupForce(d.GROUP.topTubeRing,d.GROUP.topPlaten);
d.mo.zeroBalance();
fixId=[d.GROUP.botBox;d.GROUP.topTubeRing];
d.addFixId('X',[fixId;d.GROUP.topPlaten]);
d.addFixId('Y',[fixId;d.GROUP.topPlaten]);
d.addFixId('Z',fixId);

%------------return and save result--------------
B.gravitySediment(0.5);
B.compactSample(2);%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('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();
d.showFilter('SlideY',0.3,1,'StressXX');
user_BoxShear2.m
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
%set the material of the model
clear
load('TempModel/BoxShear1.mat');
%------------initialize model-------------------
B.setUIoutput();
d=B.d;
d.calculateData();
d.mo.setGPU('off');
d.getModel();%get xyz from d.mo
%------------end initialize model-------------------

%-------------set new material----------------
matTxt=load('Mats\Soil1.txt');
Mats{1,1}=material('Soil1',matTxt,B.ballR);
Mats{1,1}.Id=1;
d.Mats=Mats;
d.groupMat2Model({'sample'},1);
%-------------end set new material----------------

%-------------------apply stress, and balance model------------------
d.resetStatus();
B.SET.platenArea=pi*B.SET.sampleR*B.SET.sampleR;
verticalStress=50e3;
verticalForce=verticalStress*B.SET.platenArea/length(d.GROUP.topPlaten);
d.mo.mGZ(d.GROUP.topPlaten)=-verticalForce;
d.breakGroup();
d.mo.zeroBalance();
B.compactSample(2);%input is compaction time
d.balance('Standard',2);
d.showFilter('SlideY',0.3,1,'mV');


%--------------------save data-----------------------
d.mo.setGPU('off');
d.clearData(1);
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']);
%--------------------end save data-----------------------
d.calculateData();
user_BoxShear3.m
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
clear;
load('TempModel/BoxShear2.mat');
B.setUIoutput();
d=B.d;
d.calculateData();
d.mo.setGPU('off');
%initializing
d.getModel();%d.setModel();%reset the initial status of the model
d.status=modelStatus(d);%initialize model status, which records running information
if strcmp(B.SET.boxType,'model')
d.status.SET.FX=[];
d.status.recordCommand='gId=d.GROUP.topTubeRing;FX=sum(sum(d.mo.nFnX(gId,:)));obj.SET.FX=[obj.SET.FX;FX];';
end

d.mo.bFilter(:)=0;%break all bonds, no tensile force
d.mo.isHeat=1;%calculate heat in the model
d.mo.setGPU('auto');
d.setStandarddT();

totalCircle=20;
stepNum=10;
dis=0.005;%total distance
dDis=dis/totalCircle/stepNum;%distance of each step
d.tic(totalCircle*stepNum);
%.mat files will be saved in the folder data/step
fName=['data/step/' B.name num2str(B.ballR) '-' num2str(B.distriRate) 'loopNum'];
save([fName '0.mat']);%return;

for i=1:totalCircle
    for j=1:stepNum
        d.toc();%show the note of time
        d.moveGroup('botBox',dDis,0,0,'mo');
        d.balance('Standard',0.1);%
    end
    d.clearData(1);%clear data in d.mo
    save([fName num2str(i) '.mat']);
    d.calculateData();
end

d.mo.setGPU('off');
d.clearData(1);
d.recordCalHour('BoxPile3Finish');
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();
user_BoxShearTorsional3.m
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
clear;
load('TempModel/BoxShear2.mat');
B.setUIoutput();
d=B.d;
d.calculateData();
d.mo.setGPU('off');
B.name='BoxShearTorsional';
%initializing
d.getModel();%d.setModel();%reset the initial status of the model
d.status=modelStatus(d);%initialize model status, which records running information

d.mo.bFilter(:)=0;%break all bonds, no tensile force
d.mo.isHeat=1;%calculate heat in the model
d.mo.setGPU('auto');
d.setStandarddT();

d.addGroup('topBox',[d.GROUP.topTubeRing;d.GROUP.topPlaten]);

totalCircle=20;
stepNum=10;
angle=60;%total rotation
dAngle=angle/totalCircle/stepNum;%distance of each step
d.tic(totalCircle*stepNum);
%.mat files will be saved in the folder data/step
fName=['data/step/' B.name num2str(B.ballR) '-' num2str(B.distriRate) 'loopNum'];
save([fName '0.mat']);%return;

for i=1:totalCircle
    for j=1:stepNum
        d.toc();%show the note of time
        d.rotateGroup('botBox','XY',dAngle,B.sampleW/2,B.sampleL/2,0,'mo');%only d.mo data will be rotate
        d.balance('Standard',0.1);%
    end
    d.clearData(1);%clear data in d.mo
    save([fName num2str(i) '.mat']);
    d.calculateData();
end

d.mo.setGPU('off');
d.clearData(1);
d.recordCalHour('BoxPile3Finish');
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();