PoreTunnel

Step 1Step 2Step 3

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

%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.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
%------------end initialize model-------------------

%-------------set new material----------------
matTxt=load('Mats\StrongSoil.txt');
Mats{1,1}=material('StrongSoil',matTxt,B.ballR);
Mats{1,1}.Id=1;
matTxt=load('Mats\Tunnel.txt');
Mats{1,2}=material('TunnelMat',matTxt,B.ballR);
Mats{1,2}.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=4;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=3;
minBallR=min(sR);
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','StrongSoil');
d.setGroupMat('ring','TunnelMat');

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

d.breakGroup();
d.mo.zeroBalance();
mfs.randomV(d,1:d.mNum);

d.mo.setShear('off');
d.mo.setGPU('auto');
d.balance('Standard',2);

d.balanceBondedModel();
for i=1:50
d.connectGroup('sample');
d.removePrestress(0.1);%remove internal stress of sample
d.balance('Standard',0.5);
end
d.mo.setShear('on');
d.balance('Standard',5);

%--------------------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();
d.show('ZDisplacement');

clear
fs.randSeed(2);
load('TempModel/PoreTunnel2.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,0.01~0.1

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
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*5000;%change the permeability

pPressureHigh=2e6;%use greater pressure
pPressure0=0.1e6;%use greater pressure
%---------setting of the simulation
totalCircle=100;
dPressure=pPressureHigh/totalCircle;
stepNum=200;
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 loop6547145258/*963

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)/2;%K of intacted bond is smaller

        p.setBallPressure(top_pId,pPressure0);%set the pressure around the sample

        %set pressure
        tunnel_cx=mean(d.mo.aX(d.GROUP.ring));
        tunnel_cz=mean(d.mo.aZ(d.GROUP.ring));
        x1=tunnel_cx-10;
        x2=tunnel_cx+10;
        z1=tunnel_cz-2.5;
        z2=tunnel_cz+2.5;
        border=0.3;
        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);
        %end set pressure

        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