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 start by reading the model from a file in the Eclipse formate (GRDECL). In Real-Field Model: I, we examined the model in more detail, and showed that the grid has two components, of which we will only use the first one.

grdecl = readGRDECL(fullfile(ROOTDIR, 'examples', 'grids', 'GSmodel.grdecl')); G = processGRDECL(grdecl, 'Verbose', false); clear grdecl; G = computeGeometry(G(1), 'Verbose', false);

The permeability is lognormal and isotropic within nine distinct layers and is generated using our simplified 'geostatistics' function and then transformed to lay in the interval 200-2000 mD. For the permeability-porosity relationship we use the simple relationship that phi~0.25*K^0.11, porosities in the interval 0.25-0.32. For the two-phase fluid model, we use values:

gravity off K = logNormLayers(G.cartDims, rand(9,1), 'sigma', 2); K = K(G.cells.indexMap); K = 200 + (K-min(K))/(max(K)-min(K))*1800; rock.perm = K*milli*darcy; rock.poro = 0.25*(K/200).^0.1; clear K; fluid = initSimpleFluid('mu', [1 5]); clf, plotCellData(G,log10(rock.perm),'EdgeColor','k'); axis off, view(15,60), h=colorbar('horiz'); cs = [200 400 700 1000 1500 2000]; caxis(log10([min(cs) max(cs)]*milli*darcy)); set(h, 'XTick', log10(cs*milli*darcy), 'XTickLabel', num2str(round(cs)')); zoom(2.5), title('Log_{10} of x-permeability [mD]');

Initialize solution structure with reservoir pressure equal 350 bar. Assemble linear hybrid system from input grid and rock properties.

The reservoir is produced using a set 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, 26, 8, 25, 35, 10]; J = [14, 14, 35, 35, 68, 75]; R = [ 4, 4, 4, 4, 4, 4]*1000*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,0], ... 'name', ['I$_{', int2str(i), '}$']); end % Set vertical producers, completed in the upper 14 layers I = [17, 12, 25, 35, 15]; 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), '}$']); 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, zoom(1.1), view(-5,58) 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 = assembleMimeticSystem(G, rock, 'Verbose', true); W = assembleWellSystem(G, W); rSol = initResSol(G, 350*barsa, 0.0); wSol = initWellSol(W, 300*barsa());

Using inner product: 'ip_simple'. Computing component matrices C and D ... Elapsed time is 0.018762 seconds. Computing cell inner products ... Elapsed time is 6.298584 seconds. Assembling global inner product matrix ... Elapsed time is 0.097781 seconds.

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

[rSol, wSol] = solveIncompFlow(rSol, wSol, G, S, fluid, 'wells', W); clf plotCellData(G, convertTo(rSol.cellPressure, barsa), 'EdgeColor','k'); title('Initial pressure'), colorbar('horiz') plotWell(G,W,'height',200,'color','w'); axis tight off; view(20,80); zoom(2)

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 = 30*year(); dT = T/12; dTplot = 5*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(30,50), colormap(flipud(jet)) colorbar('horiz'); hs = []; ha=[]; zoom(2.5) % Start the main loop t = 0; plotNo = 1; while t < T,

rSol = implicitTransport(rSol, wSol, G, dT, rock, fluid, 'wells', W); % Check for inconsistent saturations assert(max(rSol.s) < 1+eps && min(rSol.s) > -eps); % Update solution of pressure equation. [rSol, wSol] = solveIncompFlow(rSol, wSol, 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

Plot saturation

delete([hs, ha]) hs = plotCellData(G, rSol.s,find(rSol.s>0.01)); ha = annotation('textbox',[0.6 0.2 0.5 0.1], 'LineStyle','none', ... 'String', ['Water saturation at ',num2str(convertTo(t,year)),' years']); view(30, 50+7*(plotNo-1)), drawnow plotNo = plotNo+1;

end

Published March 28, 2009

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