跳转至

3DJointStress

user_Box3DJointStress1.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
%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='3DJointStress';
B.GPUstatus='auto';
B.ballR=0.005;
B.isClump=0;
B.distriRate=0.2;
B.sampleW=0.1;
B.sampleL=0.1;
B.sampleH=0.15;
B.BexpandRate=B.sampleW*0.1/B.ballR;
B.PexpandRate=B.sampleW*0.1/B.ballR;
B.setType('TriaxialCompression');
B.buildInitialModel();%B.show();
B.setUIoutput();

d=B.d;
%--------------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.showFilter('Group',{'sample'});
d.show('mGZ');

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();%because data is clear, it will be re-calculated
user_Box3DJointStress2.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
%set the material of the model
clear
load('TempModel/3DJointStress1.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\StrongRock.txt');
Mats{1,1}=material('StrongRock',matTxt,B.ballR);
Mats{1,1}.Id=1;
d.Mats=Mats;
d.groupMat2Model({'sample'},1);
%-------------end set new material----------------

%-------------------apply stress, and balance model------------------
B.SET.stressXX=-6e6;
B.SET.stressYY=-6e6;
B.SET.stressZZ=-10e6;
B.setPlatenFixId();
d.resetStatus();
B.setPlatenStress(B.SET.stressXX,B.SET.stressYY,B.SET.stressZZ,B.ballR*5);

d.balanceBondedModel0(4);
d.mo.dT=d.mo.dT*4;%increase the dT to increase the speed of balance
aMUp=d.mo.aMUp;d.mo.aMUp(:)=0;%remove friction
d.balance('Standard',5);
d.mo.aMUp=aMUp;%restore friction
d.mo.dT=d.mo.dT/4;
%-------------------end apply stress, and balance model------------------

%--------------------save data-----------------------
d.mo.setGPU('off');
d.clearData(1);
d.recordCalHour('Step3Finish');
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();
d.show('ZDisplacement');
user_Box3DJointStress3.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
%-------------------user_mxSlope3.m;
clear;
load('TempModel/3DJointStress2.mat');
B.setUIoutput();
%------------initialize model-------------------
d=B.d;
d.calculateData();
d.mo.setGPU('off');
d.mo.bFilter(:)=true;
d.deleteConnection('boundary');
d.mo.zeroBalance();
d.getModel();%d.setModel();%reset the initial status of the model
d.resetStatus();
d.mo.isHeat=1;%calculate heat in the model
d.setStandarddT();
d.mo.isCrack=1;
d.moveBoundary('right',B.sampleW*0.2,0,0);
d.moveBoundary('back',0,B.sampleL*0.2,0);
d.mo.isCrack=1;%will record crack data

%------------end initialize model-------------------
%--------------------apply initial stress-----------------------
type='break';%change type to 'break' to create crack
if strcmp(type,'glue')
    d.mo.bFilter(:)=false;
end

%------------crack-------------
TriX=[0.01,0.08,0.02];TriY=[0.05,0.1,0.05];TriZ=[0,0.05,0.1];
bondFilter=mfs.setBondByTriangle(d,TriX,TriY,TriZ,type);%make triangle crack, run d.mo.zeroBalance() after the function
d.mo.zeroBalance();
d.showFilter('Group',{'sample'});
%------------end crack-------------

bondFilter=bondFilter&(d.mo.cFilter|d.mo.bFilter);
connectFilter=sum(bondFilter,2)>0;
connectId=find(connectFilter);
d.addGroup('JointLayer',connectId);%add a new group
figure;
d.showFilter('Group','JointLayer','aR');

%------------weak or strong joint-------------
TriX=[0.09,0.09,0.1,0.01];TriY=[0,0.1,0,0.1];TriZ=[0,0,0.1,0.1];
bondFilter=mfs.setBondByPolygon(d,TriX,TriY,TriZ,'glue');%make triangle crack
d.mo.nBondRate(bondFilter)=2;%make strong joint
%------------end weak or strong joint------------

d.show('Crack','Stereo');

%-------------using Tool_Cut Triangle-------------
TriX2=[TriX(1:3);TriX([1,3,4])]-0.05;
TriY2=[TriY(1:3);TriY([1,3,4])];
TriZ2=[TriZ(1:3);TriZ([1,3,4])];
C=Tool_Cut(d);
C.setTriangle(TriX2,TriY2,TriZ2);
bondFilter=C.setBondByTriangle(type);
d.mo.nBondRate(bondFilter)=0.5;%make weak joint, BF and FS0 will reduce
d.show('BondRate');
C.showTriangle();
%-------------end using Tool_Cut Triangle-------------

%added in MatDEM 2.01
%------------nBondRate for breaking force (BF), friction coefficient (MUp), initial shear resistance (FS0)
d.mo.isBondRate=1;
nBondRateMUp=ones(size(d.mo.nBall));
nBondRateMUp(bondFilter)=0;
d.mo.SET.nBondRate.MUp=nBondRateMUp;%MUp could be replaced by 'BF', 'FS0'
d.show('BondRateMUp');
%------------end nBondRate for breaking force (BF), friction coefficient (MUp), initial shear resistance (FS0)

%---------------make triangle by scattered points----------------
%we modify Surf SurfTri or TriangleX,Y,Z to build new surface
fs.randSeed(2);
d.mo.bFilter(:)=false;
d.mo.zeroBalance();
randN=10;
PX=rand(randN,1)*0.1;PY=rand(randN,1)*0.1;PZ=0.05+rand(randN,1)*0.04;
C.addSurf(PX,PY,PZ);
C.getSurfTri(1,1);
C.getTriangle(1);
bondFilter=C.setBondByTriangle(type);
d.mo.bFilter=bondFilter;
d.showFilter('Group',{'sample'});
d.show('--');
C.showTriangle();
%---------------end make triangle by scattered points----------------

%--------------------start calculation-----------------------
totalCircle=5;
d.tic(totalCircle);%record the initial time of loop
fName=['data/step/' B.name  num2str(B.ballR) '-' num2str(B.distriRate) 'loopNum'];
save([fName '0.mat']);%return;
gpuStatus=d.mo.setGPU('auto');
for i=1:totalCircle
    B.setPlatenStress(B.SET.stressXX,B.SET.stressYY,B.SET.stressZZ*(1+i/totalCircle),B.ballR*5);
    d.balance('Standard',1);%standard balance
    d.clearData(1);
    d.mo.setGPU('off');
    save([fName num2str(i) '.mat']);
    d.calculateData();
    d.mo.setGPU(gpuStatus);
    d.toc();%show the note of time;
end
d.recordCalHour('Step3Finish');
d.mo.setGPU('off');
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.showFilter('Group',{'sample'},'StressZZ');