PoreHydraulic

Step 1Step 2Step 3Step 4

%step1: packing the elements
clear;
fs.randSeed(1);%random model seed, 1,2,3...
B=obj_Box;%declare a box object
B.name='PoreHydraulic';
%--------------initial model------------
B.GPUstatus='auto';%program will test the CPU and GPU speed, and choose the quicker one
B.ballR=0.002;
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=0.2;%width, length, height, average radius
B.sampleL=0;%when L is zero, it is a 2-dimensional model
B.sampleH=0.2;
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------------

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

%---------- gravity sedimentation
B.gravitySediment(1);%you may use B.gravitySediment(10); to increase sedimentation time (10)
%------------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('-Id');

%set the material of the model
clear
load('TempModel/PoreHydraulic1.mat');
B.setUIoutput();%set output of message
d=B.d;
d.calculateData();
d.mo.setGPU('off');
d.getModel();%get xyz from d.mo

%----------set material of model
matTxt=load('Mats\RockHydro.txt');
Mats{1,1}=material('RockHydro',matTxt,B.ballR);
Mats{1,1}.Id=1;
d.Mats=Mats;
%----------end set material of model

%---------assign material to layers and balance the model
B.setPlatenFixId();
d.setGroupMat('sample','RockHydro');
d.groupMat2Model({'sample'});
d.balanceBondedModel0();
d.balance('Standard',2);
%---------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();

clear
fs.randSeed(2);
load('TempModel/PoreHydraulic2.mat');
B.setUIoutput();

%--------initializing the model
d=B.d;
d.calculateData();
d.mo.setGPU('off');
d.getModel();%get xyz from d.mo
d.showB=2;
d.deleteConnection('boundary');
d.Rrate=1;
d.getModel();
d.mo.isCrack=1;
d.mo.bFilter(:)=true;%break all connections
d.mo.zeroBalance();
d.mo.mVis=d.mo.mVis./d.vRate*0.001;%use uniform viscosity

p=pore(d);
d.mo.dT=d.mo.dT/5;%use small step time
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';
figureNumber=d.show('mV');%define figureNumber, figure will shown in one form during iterations
d.figureNumber=figureNumber;
%----------end set the drawing of result during iterations

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

pPressureHigh=p.pPressure(1)*5000;%use greater pressure
%---------setting of the simulation
totalCircle=10;
stepNum=100;
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

for i=1: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
        p.setBallPressure(210,pPressureHigh);%you may fix the element
        p.balance();
        d.balance();
    end
    %cla;%clear the previous figure
    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 '3.mat'],'B','d');
save(['TempModel/' B.name '3R' num2str(B.ballR) '-distri' num2str(B.distriRate)  'aNum' num2str(d.aNum) '.mat']);
d.calculateData();

%set the material of the model
clear
fs.randSeed(2);
load('TempModel/PoreHydraulic2.mat');
B.setUIoutput();
%---------------regular setting
d=B.d;
d.calculateData();
d.mo.setGPU('off');
d.getModel();%get xyz from d.mo
d.showB=2;
d.deleteConnection('boundary');
d.Rrate=1;
d.getModel();
d.mo.isCrack=1;
%---------------end regular setting

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

%-----------------set cracks in the model
C=Tool_Cut(d);%cut the model
lSurf=[0.3,0,0;10,0,5;20,0,20]/100;%load the surface data
lSurf2=[0.5,0,0;5,0,15]/100;%load the surface data
C.addSurf(lSurf);%add the surfaces to the cut
C.addSurf(lSurf2);%add the surfaces to the cut
s1Filter=d.setSurfBond(C.Surf(1),'break');
s2Filter=d.setSurfBond(C.Surf(2),'break');
%-----------------end set cracks in the model
d.show('--');
C.showSurf();

%---------calculate connection diameter and flow K
k=0.00000001;%permeability factor
%---------end calcualte connection diameter and flow K
k=k/10;
pPressureHigh=p.pPressure(1)*10000;%use greater pressure

for step=1:2000
    cDiameterFlow=p.cDiameter+p.cDiameterAdd;%calculate the diameter of
    cDiameterFlow(cDiameterFlow<0)=0;
    cbFilter=d.mo.bFilter(p.cnIndex);%bonded filter of cList
    cs1Filter=s1Filter(p.cnIndex);%s2Filter, select connection of crack 2
    cs2Filter=s2Filter(p.cnIndex);%s2Filter, select connection of crack 2

    p.cKFlow=cDiameterFlow*k./p.cPathLength;%default K of throat is determined by diameter and path length
    p.cKFlow(cbFilter)=p.cKFlow(cbFilter)*0.01;%K of intacted bond is very small
    p.cKFlow(cs1Filter)=p.cKFlow(cs1Filter)*2;%bottom crack K is greater
    p.cKFlow(cs2Filter)=p.cKFlow(cs2Filter)*1;%top crack K is greater
    p.setBallPressure(1,pPressureHigh);%you may fix the element
    p.balance();
end
p.show('pPressure');

%---------save the data
d.mo.setGPU('off');
d.clearData(1);
d.recordCalHour('Step3Finish');
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();