跳转至

PorePermeability

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

d.mo.isShear=0;

%---------- 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();
%d.show('-Id');
user_L3PorePermeability2.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
%set the material of the model
clear
load('TempModel/PorePermeability1.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',1);
%---------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_L3PorePermeability3.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
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
%test the permeability of sample
%in this example, there are two material with different permeability

clear
fs.randSeed(2);
load('TempModel/PorePermeability2.mat');
sampleH=0.04;%height of sample is 4 cm
dH=1;%water head is 1m
%---------define permeability
kRate=0.5;%change the kRate to change permeability
k=2e-10*kRate;%permeability factor
lowKRate=0.1;%percentage of low permeability elements
totalBalance=100000;%increase the value when highK:lowK increase
%---------end define permeability

%-----------initializing the model
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.resetStatus();
d.getModel();
d.mo.isCrack=1;
%-----------end initializing the model

%------------remove top and bottom elements of sample
sampleHcenter=mean(d.mo.aZ(d.GROUP.sample));
topModelFilter=d.mo.aZ>sampleHcenter+sampleH/2;
botModelFilter=d.mo.aZ<sampleHcenter-sampleH/2;
delId=find(topModelFilter|botModelFilter);
maxS=max(d.GROUP.sample);
delId=delId(delId<=maxS);
d.delElement(delId);
%------------end remove top and bottom elements of sample

%------------initilizing pore system
p=pore(d);
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
p.aWaterdR=d.mo.aR*0.025;%defind the water diameter
p.setWaterdR();

sLength=length(d.GROUP.sample);
lowKId=randperm(sLength,round(sLength*lowKRate));
p.aWaterdR(lowKId)=p.aWaterdR(lowKId)*0.1;
p.setWaterdR();
d.mo.SET.aWaterdR=p.aWaterdR;
%d.show('SETaWaterdR');
%return
%-----------------end set cracks in the model

%--------------setting of the simulation
p.dT=p.d.mo.dT/kRate;
fName=['data/step/' B.name  num2str(B.ballR) '-lowKRate' num2str(lowKRate) 'loopNum'];
topBallId=ceil(mean(d.GROUP.topPlaten));
botBallId=ceil(mean(d.GROUP.botPlaten));

pressureHigh=p.pPressure(2)+1e3*9.8*dH;
pressureLow=p.pPressure(2);
topPoreId=p.getBallConnectedPore(topBallId);
botPoreId=p.getBallConnectedPore(botBallId);
p.pPressure(topPoreId)=pressureHigh;
p.pPressure(botPoreId)=pressureLow;
p.setPressure();%set pore pressure of seawater
p.isCouple=0;%no fluid flow coupling
%--------------end setting of the simulation

save([fName '0.mat']);%return;
%-----------apply high pressure

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

pIndex=p.getPIndex();
topPIndex=pIndex(topPoreId,:);
botPIndex=pIndex(botPoreId,:);
balanceRates=[];
%-------------balancing the pore pressure
for i=1:totalBalance
    p.pPressure(topPoreId)=pressureHigh;
    p.pPressure(botPoreId)=pressureLow;
    p.setPressure();%set pore pressure of seawater
    p.balance();%rate defines the balance time
    toppMass=p.poreFlowMass{topPIndex(1)};
    botpMass=p.poreFlowMass{botPIndex(1)};
    if mod(i,1000)==0
        fs.disp(['Calculating ' num2str(i/1000) '/' num2str(totalBalance/1000)]);
        balanceRate=-sum(botpMass,2)/sum(toppMass,2);
        balanceRates=[balanceRates;balanceRate];
        fs.disp(['Balance rate is ' num2str(balanceRate*100) ' percent']);
    end
end
%-------------end balancing the pore pressure
fs.disp('Balance of pore pressure is finished');

fs.disp('Start the permeability test');
stableT=p.totalT;
massI=0;
totalCircle=20;
stepNum=500;
topPoreMass=zeros(totalCircle*stepNum,1);
botPoreMass=zeros(totalCircle*stepNum,1);
for i=1:totalCircle
    for j=1:stepNum
        p.pPressure(topPoreId)=pressureHigh;
        p.pPressure(botPoreId)=pressureLow;
        p.setPressure();%set pore pressure of seawater
        p.balance();%rate defines the balance time
        massI=massI+1;
        toppMass=p.poreFlowMass{topPIndex(1)};
        botpMass=p.poreFlowMass{botPIndex(1)};
        topPoreMass(massI)=sum(toppMass,2);
        botPoreMass(massI)=sum(botpMass,2);
    end
    fs.disp(['Calculating ' num2str(i) '/' num2str(totalCircle)]);
    %save([fName num2str(i) '.mat']);
end

topMassAll=sum(topPoreMass);
botMassAll=sum(botPoreMass);
flowT=p.totalT-stableT;

%k = Q*L /( A*△h)
%---------calculate the hydraulic conductivity K
Q=botMassAll/1e3/flowT;
L=sampleH;
A=B.sampleW*p.pThickness;
K=Q*L/(A*dH);
fs.disp(['Permeability coefficient of the sample is ' num2str(K)]);
%---------end calculate the hydraulic conductivity K

p.show('pPressure');
p.showData('poreFlowMass');
save([fName '0.mat']);%return;