跳转至

LavaBlock

user_LavaBlock1.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
%step1: packing the elements
clear;
fs.randSeed(1);%random model seed, 1,2,3...
B=obj_Box;%declare a box object
B.name='LavaBlock';
%--------------initial model------------
B.GPUstatus='auto';%program will test the CPU and GPU speed, and choose the quicker one
B.ballR=5;
B.isShear=0;
B.isClump=0;%if isClump=1, particles are composed of several balls
B.distriRate=0.25;%define distribution of ball radius, 
B.sampleW=1000;%width, length, height, average radius
B.sampleL=0;%when L is zero, it is a 2-dimensional model
B.sampleH=500;
B.boundaryRrate=0.999999;
B.BexpandRate=2;%boundary is 4-ball wider than 
B.PexpandRate=1;
B.isSample=1;
B.type='TriaxialCompression';
B.setType();
B.buildInitialModel();%B.show();
d=B.d;%d.breakGroup('sample');d.breakGroup('lefPlaten');
%you may change the size distribution of elements here, e.g. d.mo.aR=d.aR*0.95;
d.showB=1;
%--------------end initial model------------

%----------remove overlap platen elements
d.mo.setGPU('off');
delId=[d.GROUP.topPlaten(end-1:end);d.GROUP.botPlaten(end-1:end)];
d.delElement(delId);
%----------end remove overlap platen elements

d.mo.zeroBalance();
d.mo.setShear('off');%because d.mo.isShear will be reset to 1 in d.delElement
%---------- gravity sedimentation
B.gravitySediment(1);%you may use B.gravitySediment(10); to increase sedimentation time (10)
B.compactSample(5);%input is compaction time
mfs.reduceGravity(d,1);

%------------return and save result--------------
d.status.dispEnergy();%display the energy of the model
d.Rrate=1;
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();
d.show('aR');
user_LavaBlock2.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
%set the material of the model
clear
load('TempModel/LavaBlock1.mat');
B.setUIoutput();%set output of message
d=B.d;
d.calculateData();
d.mo.setGPU('off'); 
d.mo.aX(d.GROUP.topPlaten)=d.aX(d.GROUP.topPlaten);
d.mo.aY(d.GROUP.topPlaten)=d.aY(d.GROUP.topPlaten);
d.mo.aZ(d.GROUP.topPlaten)=d.aZ(d.GROUP.topPlaten);
platenId=[d.GROUP.lefPlaten;d.GROUP.rigPlaten;d.GROUP.botPlaten;d.GROUP.topPlaten];
d.addFixId('XYZ',platenId);
d.getModel();%get xyz from d.mo
d.resetStatus();
%----------remove top part
sX=d.mo.aX(d.GROUP.sample);
sZ=d.mo.aZ(d.GROUP.sample);
topFilter=sZ>0.9*B.sampleH;
%----------make a hole
cx=B.sampleW/2;
cz=0;
holeR=50;
cFilter=sqrt((sX-cx).^2+(sZ-cz).^2)<holeR;
delId=find(topFilter|cFilter);
d.delElement(delId);

%---------define weak layer
sZ=d.mo.aZ(d.GROUP.sample);
weakFilter=sZ>0.2*B.sampleH&sZ<0.25*B.sampleH;
d.addGroup('WeakLayer',find(weakFilter));


%----------set material of model
matTxt=load('Mats\LavaRock.txt');
rate=1;
matTxt(4)=matTxt(4)*rate;
matTxt(3)=matTxt(3)*rate;
matTxt(1)=matTxt(1)*1;%lower 
Mats{1,1}=material('RockHydro',matTxt,B.ballR);
Mats{1,1}.Id=1;

rate2=0.1;
matTxt(4)=matTxt(4)*rate2;
matTxt(3)=matTxt(3)*rate2;
Mats{2,1}=material('RockWeak',matTxt,B.ballR);
Mats{2,1}.Id=2;
d.Mats=Mats;
%----------end set material of model
d.showB=2;

%---------assign material to layers and balance the model
d.mo.setShear('on');
d.removeFixId('Z',[d.GROUP.lefPlaten;d.GROUP.rigPlaten]);
d.addFixId('Z',[d.GROUP.lefPlaten(end-1),d.GROUP.lefPlaten(end),d.GROUP.rigPlaten(end-1),d.GROUP.rigPlaten(end)]);
d.setGroupMat('sample','RockHydro');
d.setGroupMat('WeakLayer','RockWeak');

d.groupMat2Model({'sample','WeakLayer'},1);

platenIds=[d.GROUP.lefPlaten;d.GROUP.rigPlaten];
d.mo.aKN(platenIds)=d.mo.aKN(platenIds)*0.1;
d.mo.aMUp(platenIds)=0;

%d.show('aMatId');return
%---------start computing
d.mo.setGPU('auto');
d.breakGroup();
d.connectGroup('sample');
d.removePrestress(0.1);
d.balance('Standard',2);
d.connectGroup('sample');
d.removePrestress(0.1);
d.balance('Standard',2);

d.show();
%---------end assign material to layers and balance the model1. 

%---------save the 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']);
d.calculateData();
user_LavaBlock3.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
%apply initial stress on the model
clear
load('TempModel/LavaBlock2.mat');
B.setUIoutput();%set output of message
d=B.d;
d.calculateData();
d.mo.setGPU('off'); 
d.getModel();
d.resetStatus();
d.mo.isCrack=1;


sZ=d.mo.aZ(d.GROUP.sample);
topFilter=sZ>max(d.mo.aZ(d.GROUP.sample))-B.sampleH*0.05;
d.addGroup('topBlock',find(topFilter));
area=B.sampleW*B.ballR*2;
stresszz=-5e6;%additional vertical pressure on the model
gZ=area*stresszz/sum(d.mo.mM(d.GROUP.topBlock));
d.mo.mGZ(d.GROUP.topBlock)=d.mo.mM(d.GROUP.topBlock)*gZ;

d.mo.aMUp(d.GROUP.topBlock)=2;
d.mo.aKS(d.GROUP.topBlock)=d.mo.aKS(d.GROUP.topBlock)*5;
platens=[d.GROUP.lefPlaten;d.GROUP.lefPlaten;d.GROUP.lefPlaten;d.GROUP.lefPlaten;];
d.addGroup('platens',platens);
d.setClump('platens');

d.mo.setGPU('auto'); 
d.balance('Standard',4);
%---------end assign material to layers and balance the model1. 

%---------save the data
d.mo.setGPU('off');
d.show();
d.clearData(1);
d.recordCalHour('Step2Finish');
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_LavaBlock4.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
clear
fs.randSeed(2);
load('TempModel/LavaBlock3.mat');
B.setUIoutput();
%--------initializing the model
d=B.d;
d.calculateData();
d.mo.setGPU('off');

d.showB=0;
d.Rrate=1;
d.getModel();
d.resetStatus();
d.mo.isCrack=1;
d.mo.mVis=d.mo.mVis./d.vRate*0.1;%use uniform viscosity

hole_pId=d.findNearestId(B.sampleW/2,0,0);
top_pId=d.findNearestId(B.sampleW/2,0,B.sampleH);
d.Rrate=0.2;

%start building pore system
p=pore(d);
p.fluid_k=0.4615e-4;p.fluid_c=870;%use oil property, default value is water

d.mo.dT=d.mo.dT;%
p.dT=p.d.mo.dT;
p.pathLimitRate=0.3;%path diameter<pathLimitRate*ballR will be connection
p.isCouple=1;%fluid-solid coupling
p.setInitialPores();
p.setPlaten('fix');%fix the coordinates of platens
%---------end initializing the model

%----------set the drawing of result during iterations
%setappdata(0,'simpleFigure',1);%use simpleFigure to increase drawing speed
%setappdata(0,'ballRate',0.01);%use small ballRate to increase drawing speed
showType='*pPressure';
d.figureNumber=1;
%----------end set the drawing of result during iterations

%d.show('Crack');hold all;d.show('--');

%---------calculate connection diameter and flow K
k=0.00000001;%permeability factor
%---------end calcualte connection diameter and flow K
k=k*10000;%change the permeability

pPressureHigh=100e6;%use greater pressure
pPressure0=0.1e6;%use greater pressure
%---------setting of the simulation
totalCircle=10;%default value is 40, use 10 to increase speed
dPressure=pPressureHigh/totalCircle;
stepNum=50;
fName=['data/step/2' B.name num2str(B.ballR) '-' num2str(B.distriRate) 'loopNum'];
save([fName '0.mat']);%return;
d.tic(totalCircle);%record the initial time of loop
%d.mo.setGPU('on');
%p.setGPU('on');

for i=1:totalCircle
    pPressure=dPressure*totalCircle;
    for j=1:stepNum
        cDiameterFlow=p.cDiameter+p.cDiameterAdd;%calculate the diameter of
        cDiameterFlow(cDiameterFlow<0)=0;
        p.cKFlow=cDiameterFlow*k./p.cPathLength;%default K of throat is determined by diameter and path length
        cbFilter=d.mo.bFilter(p.cnIndex);%bonded filter of cList
        p.cKFlow(cbFilter)=p.cKFlow(cbFilter)/100;%K of intacted bond is very small
        p.setBallPressure(hole_pId,pPressure);%set the pressure in the crack
        p.setBallPressure(top_pId,pPressure0);%set the pressure around the sample

        p.balance();
        d.balance(10);
        d.recordStatus();
        if j==0
            d.show('mV');
            return
        end
    end
    cla;%clear the previous figure
    p.show('pPressure');

    %p.show(showType);
    %pause(0.1);%pause and show the figure
    save([fName num2str(i) '.mat']);%save data
    d.recordStatus();
    d.toc();
end
%---------save the data
d.mo.setGPU('off');
d.clearData(1);
d.recordCalHour('Step3Finish');
save(['TempModel/' B.name '4.mat'],'B','d');
save(['TempModel/' B.name '4R' num2str(B.ballR) '-distri' num2str(B.distriRate)  'aNum' num2str(d.aNum) '.mat']);
d.calculateData();
%user_makeGIF