跳转至

BoxTunnelNew

user_BoxTunnelNew1 - R.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
%build the geometrical model
clear;
fs.randSeed(1);%build random model
B=obj_Box;%build a box object
B.name='PoreTunnel';
B.GPUstatus='auto';
B.ballR=0.2;
B.isShear=0
B.isClump=0
B.distriRate=0.2
B.sampleW=40
B.sampleL=0
B.sampleH=40
B.boundaryRrate=0.999999;
B.BexpandRate=2;
B.PexpandRate=1
B.isSample=1;
B.setType('Fluid');
B.buildInitialModel();
d=B.d;

%----------remove overlap platen elements

%----------end remove overlap platen elements
%--------------end initial model------------
B.gravitySediment();

%------------return and save result--------------

d.status.dispEnergy();%display the energy of the model
d.clearData(1);%clear dependent data
d.recordCalHour('PoreTunnelNew1Finish');
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_BoxTunnelNew2 - R.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
%set the material of the model
clear
load('TempModel/PoreTunnel1.mat');
%------------initialize model-------------------
B.setUIoutput();
d=B.d;
d.calculateData();
d.mo.setGPU('off');



d.getModel();%get xyz from d.m

%------------end initialize model-------------------


%-------------set new material----------------
matTxt=load('Mats\SandYanjie.txt');
Mats{1,1}=material('SandYanjie',matTxt,B.ballR);
Mats{1,1}.Id=1;
matTxt2=load('Mats\Tunnel.txt');
Mats{2,1}=material('Tunnel',matTxt2,B.ballR);
Mats{2,1}.Id=2;
d.Mats=Mats;
%-------------end set new material----------------


%--------------------make tunnel-----------------------
sampleId=d.getGroupId('sample');
sX=d.mo.aX(sampleId);sY=d.mo.aY(sampleId);sZ=d.mo.aZ(sampleId);sR=d.mo.aR(sampleId);
dipD=0;dipA=90;radius=2.75;height=30;
columnFilter=f.run('fun/getColumnFilter.m',sX,sY,sZ,dipD,dipA,radius+B.ballR,height);
d.addGroup('Tunnel',find(columnFilter));
tunnelId=d.getGroupId('Tunnel');
d.delElement(tunnelId);


%--------------------end make tunnel-----------------------


B.name='PoreTunnel';
innerR=radius;
layerNum=2;
minBallR=0.1;
Rrate=0.8;


ringObj=f.run('fun/makeRing.m',innerR,layerNum,minBallR,Rrate);
ringObj=mfs.rotate(ringObj,'YZ',90);%rotate the group along XZ plane
ringId=d.addElement(1,ringObj);%add a slope boundary
d.addGroup('ring',ringId);%add a new group
d.setClump('ring');%set the pile clump
d.moveGroup('ring',(max(sX)+min(sX))/2,0,(max(sZ)+min(sZ))/2);
d.minusGroup('sample','ring',0.4);%remove overlap elements from sample


d.setGroupMat('sample','SandYanjie');
d.setGroupMat('ring','Tunnel');
d.groupMat2Model({'sample','ring'},2);


R1=max(d.mo.aX(d.GROUP.ring))-min(d.mo.aX(d.GROUP.ring))+2*max(ringObj.R)


tp=d.GROUP.topPlaten
blockWidth=10;
tpX=d.mo.aX(tp);
tpXCenter=(max(tpX)+min(tpX))/2;
blockFilter=tpX>(tpXCenter-blockWidth/2)&tpX<(tpXCenter+blockWidth/2);
blockId=tp(blockFilter);
d.addGroup('block',blockId);
gpuStatus=d.mo.setGPU('auto');


d.breakGroup();
d.mo.zeroBalance();
d.mo.setShear('off');
d.mo.setGPU('auto');
d.balance('Standard',3)


R2=max(d.mo.aX(d.GROUP.ring))-min(d.mo.aX(d.GROUP.ring))+2*max(ringObj.R)


d.mo.mGZ(blockId)=d.mo.mGZ(blockId)+(-6e+05);

d.breakGroup();
d.mo.zeroBalance();
d.mo.setShear('off');
d.mo.setGPU('auto');
d.balance('Standard',3)

%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

R3=max(d.mo.aX(d.GROUP.ring))-min(d.mo.aX(d.GROUP.ring))+2*max(ringObj.R)


%--------------------save data-----------------------
d.mo.setGPU('off');
d.clearData(1);
d.recordCalHour('BoxTunnel2Finish');
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-----------------------
user_BoxTunnelNew3 - R.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
%-------------------user_mxSlope3.m;
clear;
load('TempModel/PoreTunnel2.mat');
B.setUIoutput();
%------------initialize model-------------------
d=B.d;
d.calculateData();
d.getModel();
d.mo.setGPU('off');
d.breakGroup();
d.mo.setShear('on');
d.Rrate=1;
d.resetStatus();
d.mo.isCrack=1;
d.mo.mVis=d.mo.mVis./d.vRate*0.1;
top_pId=d.findNearestId(B.sampleW/2,0,B.sampleH);

%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
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.00001;%permeability factor
%---------end calcualte connection diameter and flow K
%k=k*5000;%change the permeability

pPressureHigh=3.25e+07;%use greater pressure
pPressure0=0.1e6;%use greater pressure
%---------setting of the simulation
totalCircle=10;%default value is 100, use 10 to increase speed
dPressure=pPressureHigh/totalCircle;
stepNum=20;
fName=['data/step/2' B.name num2str(B.ballR) '-' num2str(B.distriRate) 'loopNum'];
save([fName '0.mat']);%return;
d.tic(totalCircle*stepNum);%record the initial time of loop6547145258/*963
R4=zeros(totalCircle,1);
SUM=zeros(totalCircle,1);

tunnel_cx0=mean(d.mo.aX(d.GROUP.ring));
tunnel_cz0=mean(d.mo.aZ(d.GROUP.ring));




for i=1:totalCircle
   for j=1:stepNum
        d.toc();
        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)/2;%K of intacted bond is smaller
        %set pressure
        tunnel_cx=mean(d.mo.aX(d.GROUP.ring));
        tunnel_cz=mean(d.mo.aZ(d.GROUP.ring));
        x1=tunnel_cx-6.7;
        x2=tunnel_cx+6.7;
        z1=tunnel_cz-2.5;
        z2=tunnel_cz+2.5;
        border=0.6;
        pIdLeft=pfs.getRectAreaPoreId(p,x1-border,z1,x1+border,z2);
        p.setPressure(pIdLeft,pPressureHigh); 
       pIdRight=pfs.getRectAreaPoreId(p,x2-border,z1,x2+border,z2);

        p.setPressure(pIdRight,pPressureHigh);
        p.pDensity(pIdLeft)=1700;
        p.pDensity(pIdRight)=1700;
        %end set pressure
        p.balance();
       %d.mo.setShear('off');
        d.balance();
        d.recordStatus();
    end
    cla;%clear the previous figure
    p.show('pPressure');
    save([fName num2str(i) '.mat']);%save data
    R4(i)=max(d.mo.aX(d.GROUP.ring))-min(d.mo.aX(d.GROUP.ring))+2*max(d.mo.aR(d.GROUP.ring));


    save([fName num2str(i) '.mat']);%save data
    d.recordStatus();

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