-
Notifications
You must be signed in to change notification settings - Fork 9
/
johansen3D.m
150 lines (118 loc) · 5.34 KB
/
johansen3D.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
142
143
144
145
146
147
148
149
150
mrstModule add ad-core;
mrstModule add ad-props;
%% Grid and rock
% Load Johansen model
[G, rock, bcIx] = makeJohansenVEgrid();
% Plot grid and rock properties
figure; plotCellData(G, rock.poro); view(-35,15); colorbar; % porosity
set(gcf, 'position', [531 337 923 356]); axis tight;
figure; plotCellData(G, rock.perm(:,1)); view(-55, 60); colorbar; % lateral permeability
set(gcf, 'position', [152 419 689 404]);
figure; plotCellData(G, rock.perm(:,3)); view(-55, 60); colorbar; % vertical permeability
set(gcf, 'position', [152 419 689 404]);
%% Initial state
gravity on; % tell MRST to turn on gravity
g = gravity; % get the gravity vector
rhow = 1000; % density of brine corresponding to 94 degrees C and 300 bar
initState.pressure = rhow * g(3) * G.cells.centroids(:,3);
initState.s = repmat([1, 0], G.cells.num, 1);
initState.sGmax = initState.s(:,2);
%% Fluid model
co2 = CO2props(); % load sampled tables of co2 fluid properties
p_ref = 30 * mega * Pascal; % choose reference pressure
t_ref = 94 + 273.15; % choose reference temperature, in Kelvin
rhoc = co2.rho(p_ref, t_ref); % co2 density at ref. press/temp
cf_co2 = co2.rhoDP(p_ref, t_ref) / rhoc; % co2 compressibility
cf_wat = 0; % brine compressibility (zero)
cf_rock = 4.35e-5 / barsa; % rock compressibility
muw = 8e-4 * Pascal * second; % brine viscosity
muco2 = co2.mu(p_ref, t_ref) * Pascal * second; % co2 viscosity
mrstModule add ad-props; % The module where initSimpleADIFluid is found
% Use function 'initSimpleADIFluid' to make a simple fluid object
fluid = initSimpleADIFluid('phases', 'WG' , ...
'mu' , [muw, muco2] , ...
'rho' , [rhow, rhoc] , ...
'pRef', p_ref , ...
'c' , [cf_wat, cf_co2] , ...
'cR' , cf_rock , ...
'n' , [2 2]);
% plot relperm curves
sw = linspace(0, 1, 200);
figure; hold on;
plot(sw, fluid.krW(sw), 'b', 'linewidth', 1.5);
plot(sw, fluid.krG(1-sw), 'r', 'linewidth', 1.5);
xlabel('brine saturation'); ylabel('relative permeability')
set(gca, 'fontsize', 14);
% Change relperm curves
srw = 0.27;
src = 0.20;
fluid.krW = @(s) fluid.krW(max((s-srw)./(1-srw), 0));
fluid.krG = @(s) fluid.krG(max((s-src)./(1-src), 0));
% Plot relperm curves
figure; hold on;
sw = linspace(srw, 1, 200);
plot(sw, fluid.krW(sw), 'b', 'linewidth', 1.5);
plot(sw, fluid.krG(1-sw), 'r', 'linewidth', 1.5);
line([srw, srw], [0 1], 'color', 'black', 'linestyle', ':', 'linewidth', 1);
xlabel('brine saturation'); ylabel('relative permeability')
set(gca, 'fontsize', 14, 'xlim', [0 1]);
% Add capillary pressure curve
pe = 5 * kilo * Pascal;
pcWG = @(sw) pe * sw.^(-1/2);
%fluid.pcWG = @(sg) pcWG(max((1-sg-srw)./(1-srw), 0));
fluid.pcWG = @(sg) pcWG(max((1-sg-srw)./(1-srw), 1e-5)); %@@
%% Wells
% Well cell indices in 'global' grid: 48, 48, 6:10
wc_global = false(G.cartDims);
wc_global(48, 48, 6:10) = true;
wc = find(wc_global(G.cells.indexMap));
% plot well cells
plotGrid(G, 'facecolor', 'none', 'edgealpha', 0.1);
plotGrid(extractSubgrid(G, wc), 'facecolor', 'red');
% Calculate the injection rate
inj_rate = 3.5 * mega * 1e3 / year / fluid.rhoGS;
% Start with empty set of wells
W = [];
% Add a well to the set
W = addWell([], G, rock, wc, ...
'type', 'rate', ... % inject at constant rate
'val', inj_rate, ... % volumetric injection rate
'comp_i', [0 1]); % inject CO2, not water
%% Boundary conditions
% Start with an empty set of boundary faces
bc = [];
% Computing hydrostatic pressure for boundary cells
p_bc = G.faces.centroids(bcIx, 3) * rhow * g(3);
% Add hydrostatic pressure conditions to open boundary faces
bc = addBC(bc, bcIx, 'pressure', p_bc);
%% Schedule
% Setting up two copies of the well and boundary specifications.
% Modifying the well in the second copy to have a zero flow rate.
schedule.control = struct('W', W, 'bc', bc);
schedule.control(2) = struct('W', W, 'bc', bc);
schedule.control(2).W.val = 0;
% Specifying length of simulation timesteps
schedule.step.val = [repmat(year, 100, 1); ...
repmat(10*year, 100, 1)];
% Specifying which control to use for each timestep.
% The first 100 timesteps will use control 1, the last 100
% timesteps will use control 2.
schedule.step.control = [ones(100, 1); ...
ones(100, 1) * 2];
%% Model
model = twoPhaseGasWaterModel(G, rock, fluid, 0, 0);
%% Simulate
[wellSol, states] = simulateScheduleAD(initState, model, schedule);
%% Plot results
figure; plotCellData(G, states{100}.s(:,2)); view(-63, 68); colorbar;
set(gcf, 'position', [531 337 923 356]); axis tight;
figure; plotCellData(G, states{200}.s(:,2)); view(-63, 68); colorbar;
set(gcf, 'position', [531 337 923 356]); axis tight;
% vertical plot
[i j k] = ind2sub(G.cartDims, G.cells.indexMap);
figure; plotCellData(extractSubgrid(G, j==48 & i>18 & i < 60), states{100}.s(j==48 & i>18 & i < 60,2));
view(0,0); axis tight;
%set(gcf, 'position', [500 500 760 300]);
figure; plotCellData(extractSubgrid(G, j==48 & i>18 & i < 60), states{200}.s(j==48 & i>18 & i < 60,2));
view(0,0); axis tight;
%set(gcf, 'position', [500 500 760 300]);