-
Notifications
You must be signed in to change notification settings - Fork 9
/
runJohansenSimulation.m
162 lines (142 loc) · 6.22 KB
/
runJohansenSimulation.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
151
152
153
154
155
156
157
158
159
160
161
162
%% Vertical-Averaged Simulation of the Johansen Formation
% The Johansen formation is a candidate site for large-scale CO2 storage
% offshore the south-west coast of Norway. In the following, we will use a
% simple vertically averaged model to simulate the early-stage migration of
% a CO2 plume injected from a single well positioned near the main fault in
% the formation. The formation is described by a geological model that has
% been developed based on available seismic and well data. A more thorough
% presentation of the geological model can be found in the script
% <matlab:edit('showJohansen.m') showJohansen.m>
%
% The data files
% necessary to run the example can be downloaded from the
% <http://www.sintef.no/Projectweb/MatMorA/Downloads/Johansen/ MatMoRA
% website>.
mrstModule add co2lab mimetic
%% Display header
clc;
disp('================================================================');
disp(' Vertical averaging applied to the Johansen formation');
disp('================================================================');
disp('');
%% Input data and construct grid models
% We use a sector model in given in the Eclipse input format (GRDECL). The
% model has five vertical layers in the Johansen formation and five shale
% layers above and one below in the Dunhil and Amundsen formations. The
% shale layers are removed and we construct the 2D VE grid of the top
% surface, assuming that the major fault is sealing, and identify all outer
% boundaries that are open to flow. Store grid and rock structures to file
% to avoid time-consuming processing.
[G, rock, ~, Gt, rock2D, bcIxVE] = makeJohansenVEgrid();
%% Set time and fluid parameters
gravity on
T = 1000*year();
stopInject = 50*year();
dT = 2*year();
dTplot = 1*dT;
% Fluid data at p = 300 bar
muw = 0.30860; rhow = 975.86; sw = 0.1;
muc = 0.056641; rhoc = 686.54; srco2 = 0.2;
kwm = [0.2142 0.85];
fluidVE = initVEFluidHForm(Gt, 'mu' , [muc muw] .* centi*poise, ...
'rho', [rhoc rhow] .* kilogram/meter^3, ...
'sr', srco2, 'sw', sw, 'kwm', kwm);
%% Set well and boundary conditions
% We use one well placed in the center of the model, perforated in layer 6.
% Injection rate is 1.4e4 m^3/day of supercritical CO2. Hydrostatic
% boundary conditions are specified on all outer boundaries that are not in
% contact with the shales; the latter are assumed to be no-flow boundaries.
% Set well in 3D model
wellIx = [54, 54, 6, 6];
rate = 1.4e4*meter^3/day;
W = verticalWell([], G, rock, wellIx(1), wellIx(2), wellIx(3):wellIx(4),...
'Type', 'rate', 'Val', rate, 'Radius', 0.1, 'comp_i', [1,0], ...
'name', 'I', 'InnerProduct', 'ip_simple');
% Well and BC in 2D model
WVE = convertwellsVE(W, G, Gt, rock2D);
bcVE = addBC([], bcIxVE, 'pressure', Gt.faces.z(bcIxVE)*rhow*norm(gravity));
bcVE = rmfield(bcVE,'sat');
bcVE.h = zeros(size(bcVE.face));
%% Prepare simulations
% Compute inner products and instantiate solution structure
% For the transport simulation, the default choice is to use a
% C-accelerated solver, but if this does not work, we use the standard
% solver from MRST.
SVE = computeMimeticIPVE(Gt, rock2D, 'Innerproduct','ip_simple');
preComp = initTransportVE(Gt, rock2D);
sol = initResSolVE(Gt, 0, 0);
sol.wellSol = initWellSol(W, 300*barsa());
sol.s = height2Sat(sol, Gt, fluidVE);
% select transport solver
try
mtransportVE;
cpp_accel = true;
catch me
d = fileparts(mfilename('fullpath'));
disp('mex-file for C++ acceleration not found');
disp(['See ', fullfile(mrstPath('co2lab'),'ve','VEmex','README'),...
' for building instructions']);
disp('Using matlab ve-transport');
cpp_accel = false;
end
%% Prepare plotting
% We will make a composite plot that consists of several parts: a 3D plot
% of the plume, a pie chart of trapped versus free volumes, a plane view of
% the plume from above, and two cross-sections in the x/y directions
% through the well
opts = {'slice', wellIx, 'Saxis', [0 1-fluidVE.res_water], 'maxH', 100, ...
'Wadd', 500, 'view', [-85 70], 'wireH', true, 'wireS', true};
plotPanelVE(G, Gt, W, sol, 0.0, zeros(1,4), opts{:});
%% Main loop
% Run the simulation using a sequential splitting with pressure and
% transport computed in separate steps. The transport solver is formulated
% with the height of the CO2 plume as the primary unknown and the relative
% height (or saturation) must therefore be reconstructed.
t = 0;
totVol = 0.0;
fprintf(1,'\nSimulating %d years on injection',convertTo(stopInject,year));
fprintf(1,' and %d years of migration\n', convertTo(T-stopInject,year));
fprintf(1,'Time: %4d years', convertTo(t,year));
tic;
while t<T
% Advance solution: compute pressure and then transport
sol = solveIncompFlowVE(sol, Gt, SVE, rock, fluidVE, ...
'bc', bcVE, 'wells', WVE);
if cpp_accel
[sol.h, sol.h_max] = mtransportVE(sol, Gt, dT, rock, ...
fluidVE, 'bc', bcVE, 'wells', WVE, ...
'gravity', norm(gravity), 'verbose', false);
else
sol = explicitTransportVE(sol, Gt, dT, rock, fluidVE, ...
'bc', bcVE, 'wells', WVE, ...
'preComp', preComp);
end
% Reconstruct 'saturation' defined as s=h/H, where h is the height of
% the CO2 plume and H is the total height of the formation
sol.s = height2Sat(sol, Gt, fluidVE);
assert( max(sol.s(:,1))<1+eps && min(sol.s(:,1))>-eps );
t = t + dT;
% Compute total injected, trapped and free volumes of CO2
if ~isempty(WVE)
totVol = totVol + WVE.val*dT;
end
vol = volumesVE(Gt, sol, rock2D, fluidVE);
% Check if we are to stop injecting
if t>= stopInject
WVE = []; bcVE = []; dT = 5*year(); dTplot = dT;
end
% Plotting
fprintf(1,'\b\b\b\b\b\b\b\b\b\b%4d years', convertTo(t,year));
if mod(t,dTplot)~= 0 && t<T
continue
else
plotPanelVE(G, Gt, W, sol, t, [vol totVol], opts{:});
drawnow
end
% break
end
fprintf(1,'\n\n');
% delete C++ simulator
if cpp_accel, mtransportVE(); end
etime = toc;
disp(['Elapsed simulation time: ', num2str(etime), ' seconds.']);