Consider a two-phase oil-water problem. Solve the two-phase pressure equation

where v is the Darcy velocity (total velocity) and lambda is the mobility, which depends on the water saturation S.

The saturation equation (conservation of the water phase) is given as:

where phi is the rock porosity, f is the Buckley-Leverett fractional flow function, and q_w is the water source term.

We will use the corner-point geometry of a synthetic reservoir model from the SAIGUP study that has faults and inactive cells and hetereogeneous anisotropic permeability.

We use MATLAB®'s built-in mldivide ("backslash") linear solver software to resolve all systems of linear equations that arise in the various discretzations. In many cases, a multigrid solver such as Notay's AGMG package may be prefereble. We refer the reader to Notay's home page at http://homepages.ulb.ac.be/~ynotay/AGMG/ for additional details on this package.

linsolve = @mldivide;

The model can be downloaded from the MRST page.

grdecl = fullfile(ROOTDIR, 'examples', 'data', 'saigup', 'SAIGUP.GRDECL'); if ~exist(grdecl, 'file'), error('Model data is not available.') end

The model data is provided as an ECLIPSE input file that can be read using the readGRDECL function.

grdecl = readGRDECL(grdecl);

MRST uses the strict SI conventions in all of its internal caluclations. The SAIGUP model, however, is provided using the ECLIPSE 'METRIC' conventions (permeabilities in mD and so on). We use the functions getUnitSystem and convertInputUnits to assist in converting the input data to MRST's internal unit conventions.

usys = getUnitSystem('METRIC'); grdecl = convertInputUnits(grdecl, usys);

We generate a space-filling geometry using the processGRDECL function and then compute a few geometric primitives (cell volumes, centroids, etc.) by means of the computeGeometry function.

G = processGRDECL (grdecl); G = computeGeometry(G);

The media (rock) properties can be extracted by means of the grdecl2Rock function.

rock = grdecl2Rock(grdecl, G.cells.indexMap);

The input data of the permeability in the SAIGUP realisation is an anisotropic tensor with zero vertical permeability in a number of cells. We work around this issue by (arbitrarily) assigning the minimum positive vertical (cross-layer) permeability to the grid blocks that have zero cross-layer permeability.

is_pos = rock.perm(:, 3) > 0; rock.perm(~is_pos, 3) = min(rock.perm(is_pos, 3));

Plot the logarithm of the permeability in the z-direction

clf, plotCellData(G,log10(rock.perm(:,3)),'EdgeColor','k','EdgeAlpha',0.1); axis off, view(-110,18); h=colorbar('horiz'); zoom(1.5) cs = [0.1 1 10 100 500 2000]; caxis(log10([min(cs) max(cs)]*milli*darcy)); set(h, 'XTick', log10(cs*milli*darcy), 'XTickLabel', num2str(round(cs)'));

For the two-phase fluid model, we use values:

gravity off fluid = initSimpleFluid('mu' , [ 1, 5]*centi*poise , ... 'rho', [1000, 700]*kilogram/meter^3, ... 'n' , [ 2, 2]);

The reservoir is produced using a set of production wells controlled by bottom-hole pressure and rate-controlled injectors. Wells are described using a Peacemann model, giving an extra set of equations that need to be assembled. For simplicity, all wells are assumed to be vertical and are assigned using the logical (i,j) subindex.

% Set vertical injectors, completed in the lowest 12 layers. nz = G.cartDims(3); I = [ 9, 8, 25, 25, 10]; J = [14, 35, 35, 95, 75]; R = [ 4, 4, 4, 4, 4, 4]*125*meter^3/day; nIW = 1:numel(I); W = []; for i = 1 : numel(I), W = verticalWell(W, G, rock, I(i), J(i), nz-11:nz, 'Type', 'rate', ... 'Val', R(i), 'Radius', 0.1, 'Comp_i', [1, 0], ... 'name', ['I$_{', int2str(i), '}$']); end % Set vertical producers, completed in the upper 14 layers I = [17, 12, 25, 33, 7]; J = [23, 51, 51, 95, 94]; nPW = (1:numel(I))+max(nIW); for i = 1 : numel(I), W = verticalWell(W, G, rock, I(i), J(i), 1:14, 'Type', 'bhp', ... 'Val', 300*barsa(), 'Radius', 0.1, ... 'name', ['P$_{', int2str(i), '}$'], 'Comp_i', [0, 1]); end % Plot grid outline and the wells clf subplot('position',[0.02 0.02 0.96 0.96]); plotGrid(G,'FaceColor','none','EdgeAlpha',0.1); axis tight off, view(-120,50) plotWell(G,W,'height',200); plotGrid(G, vertcat(W(nIW).cells), 'FaceColor', 'b'); plotGrid(G, vertcat(W(nPW).cells), 'FaceColor', 'r');

Initialize solution structures and assemble linear hybrid system from input grid, rock properties, and well structure.

S = computeMimeticIP(G, rock, 'Verbose', true); rSol = initState(G, W, 0, [0, 1]);

Using inner product: 'ip_simple'. Computing cell inner products ... Elapsed time is 9.530061 seconds. Assembling global inner product matrix ... Elapsed time is 0.087311 seconds.

Solve linear system construced from S and W to obtain solution for flow and pressure in the reservoir and the wells.

rSol = solveIncompFlow(rSol, G, S, fluid, 'wells', W, ... 'LinSolve', linsolve); clf plotCellData(G, convertTo(rSol.pressure(1:G.cells.num), barsa), ... 'EdgeColor', 'k', 'EdgeAlpha', 0.1); title('Initial pressure'), colorbar; plotWell(G, W, 'height', 200, 'color', 'k'); axis tight off; view(-80,36); caxis([300 400]) zoom(1.3)

In the main loop, we alternate between solving the transport and the flow equations. The transport equation is solved using the standard implicit single-point upwind scheme with a simple Newton-Raphson nonlinear solver.

T = 8*year(); dT = T/4; dTplot = 2*year(); pv = poreVolume(G,rock); % Prepare plotting of saturations clf plotGrid(G, 'FaceColor', 'none', 'EdgeAlpha', 0.1); plotWell(G, W, 'height', 200, 'color', 'c'); axis off, view(-80,30), colormap(flipud(jet)) colorbar; hs = []; ha=[]; zoom(1.3); % Start the main loop t = 0; plotNo = 1; while t < T,

rSol = implicitTransport(rSol, G, dT, rock, fluid, 'wells', W); % Check for inconsistent saturations assert(max(rSol.s(:,1)) < 1+eps && min(rSol.s(:,1)) > -eps); % Update solution of pressure equation. rSol = solveIncompFlow(rSol, G, S, fluid, 'wells', W); % Increase time and continue if we do not want to plot saturations t = t + dT; if ( t < plotNo*dTplot && t continue, end

Plot saturation

delete([hs, ha]) hs = plotCellData(G, rSol.s(:,1), find(rSol.s(:,1) > 0.01)); ha = annotation('textbox', [0.1714 0.8214 0.5000 0.1000], 'LineStyle', 'none', ... 'String', ['Water saturation at ', ... num2str(convertTo(t,year)), ' years']); view(-80, 36), drawnow, caxis([0 1]) plotNo = plotNo+1;

end

Published with MATLAB® 7.9

Published February 23, 2011

By clicking the "Agree"-button you are agreeing to our use of cookies. Find out more in our privacy policy (in Norwegian only).