GeoThermalBox

Step 1Step 2Step 3

%step1: packing the elements
clear;
fs.randSeed(1);%random model seed, 1,2,3...
B=obj_Box;%declare a box object
B.name='GeoThermalBox';
%--------------initial model------------
B.GPUstatus='auto';%program will test the CPU and GPU speed, and choose the quicker one
B.ballR=0.05;
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=3;%width, length, height, average radius
B.sampleL=0;%when L is zero, it is a 2-dimensional model
B.sampleH=1.5;
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)
B.compactSample(2);%input is compaction time
%------------return and save result--------------
d.status.dispEnergy();%display the energy of the model
d.show('-aR');
d.mo.bFilter(:)=1;
d.mo.zeroBalance();
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();

%set the material of the model
clear
load('TempModel/GeoThermalBox1.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.mo.bFilter(:)=false;
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();

%this code is used simulate the moisture-heat test in Suzhou
%set the material of the model
clear
fs.randSeed(2);
load('TempModel/GeoThermalBox2.mat');

%---------------regular setting
B.setUIoutput();
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

%---------------set different layers of the box model
interface=[0.0714;0.7964];%interface between layers
maxZ=max(d.mo.aZ);
botLayerFilter=d.mo.aZ<maxZ*interface(1);
midLayerFilter=d.mo.aZ>=maxZ*interface(1)&d.mo.aZ<=maxZ*interface(2);
topLayerFilter=d.mo.aZ>maxZ*interface(2);
d.addGroup('botLayer',find(botLayerFilter));%define the bottom layer of the model
d.addGroup('midLayer',find(midLayerFilter));
d.addGroup('topLayer',find(topLayerFilter));
d.mo.zeroBalance();
d.setGroupId();%distinguish different groups
d.show('groupId');%show groupId
%---------------end set different layers of the box model

%-------------initializing the pore network
p=pore(d);%make pore object
p.dT=p.d.mo.dT;%use the same step time, may be modified later
p.pathLimitRate=0.3;%path diameter<pathLimitRate*ballR will be connection
p.isCouple=0;%no fluid-solid coupling
p.setInitialPores();
p.setPlaten('fix');%fix the coordinates of platens
%-------------end initializing the pore network

%-----------set the fluid flow parameters (permeability)
p.aWaterdR=d.mo.aR*0.1;%%water radius deviation of model elements
p.aWaterdR(d.GROUP.midLayer)=p.aWaterdR(d.GROUP.midLayer)/5;
%elements of midLayer use lower value, i.e. lower permeability
p.setWaterdR();%calculate the p.addDiameter based on aWaterdR
%d.mo.SET.aWaterdR=p.aWaterdR;
%d.show('SETaWaterdR');%show the water radius of element
%return
%-----------end set the fluid flow parameters (permeability)

%---------calculate connection diameter and flow K
kFlow=0.00000001;%permeability factor
kT=1000;%heat conductivity factor
%---------end calcualte connection diameter and flow K

%---------find four corners of the box model
sX=d.mo.aX(d.GROUP.sample);sZ=d.mo.aZ(d.GROUP.sample);
[~,lefBotId]=min(sX+sZ);
[~,rigBotId]=min(-sX+sZ);
[~,lefTopId]=min(sX-sZ);
[~,rigTopId]=min(-sX-sZ);
pressureHigh=p.pPressure(1)*10000;%use great pressure to increase the speed
pressureLow=p.pPressure(1)*1;
%---------end find four corners of the box model

%---------------add temperature "solute"
TPara.Id=1;
TPara.name='T';
TPara.initialValue=12;%initial temperature
p.addSolutePara(TPara);
%---------------end add temperature "solute"

dNum=1000;%save the data every dNum steps
fName=['data/step/' B.name  num2str(B.ballR) '-' num2str(B.distriRate) 'loopNum'];
save([fName '0.mat']);%return;

for step=1:10000
    %----------heat conduct
    p.SET.cKT=p.cLength*kT;
    p.setBallPara('T',lefBotId,3);%3 degrees
    p.setBallPara('T',lefTopId,30);%30 degrees
    %----------end heat conduct

    %------------fluid flow
    cDiameterFlow=p.cDiameter+p.cDiameterAdd;%calculate the diameter of
    cDiameterFlow(cDiameterFlow<0)=0;
    p.cKFlow=cDiameterFlow*kFlow./p.cPathLength;%default K of throat is determined by diameter and path length
    p.setBallPressure(lefBotId,pressureHigh);%set the pore pressure around the ball
    p.setBallPressure(lefTopId,pressureHigh);
    p.setBallPressure(rigBotId,pressureLow);
    %------------end fluid flow

    p.balance();%calculation
    if mod(step,dNum)==0
        saveIndex=ceil(step/dNum);
        save([fName num2str(saveIndex) '.mat']);
        %figure;p.show('pPressure');
    end
end
p.show('SETpT');
%p.show('pPressure');
p.showData('poreFlowMass');

%---------save the data
d.mo.zeroBalance();
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();