SPE10 Quarter Five Spot Example
Contents
Document example dependciesrequire impes spe10 deckformat Define a oil/gas systemThe fluid data is contained in the file 'simpleOilGas.txt'. We read and process this using the deckformat module. current_dir = fileparts(mfilename('fullpath')); f = fullfile(current_dir, 'simpleOilGas.txt'); deck = readEclipseDeck(f); % Once the file has been read, we can use it to create the desired fluid. fluid = initEclipseFluid(deck); Set up grid and get petrophysical properties.We select three layers to layers = 25:28; cartDims = [60, 220, numel(layers)]; physDims = cartDims.*[20,10,2]*ft; rock = SPE10_rock(layers); rock.perm = convertFrom(rock.perm, milli*darcy); gravity off % Define a Cartesian Grid based on the SPE layers selected. G = cartGrid(cartDims, physDims); G = computeGeometry(G); % Cells with zero porosity are masked away to the minimum value to get % sensible results. is_pos = rock.poro>0; rock.poro(~is_pos) = min(rock.poro(is_pos)); Show grid and petrophysical propertiesSince the variations are large, use a logarithmic scale for permeability perm = log10(rock.perm); c = {'x', 'y', 'z'}; clf; for i = 1:3 subplot(2,2,i); plotCellData(G, perm(:,i)); colorbar(); axis equal tight; title(sprintf('Permeability in %s direction', c{i})); end % Plot porosity subplot(2,2,4); plotCellData(G, rock.poro); colorbar(); axis equal tight; title('Porosity') Set up four producers at reservoir pressure and an injector at 500 bar% Initialize empty well structure W = []; % Add an injector in the middle of the domain based on the midpoints of the % cartesian indices midpoint = round(cartDims)./2; % Use vertical well, with completion in all layers. To complete in a % different set of layers, enter the indices of the layers in completion completion = []; W = verticalWell(W, G, rock, midpoint(1), midpoint(2), completion, ... 'Type', 'bhp',... % Bottom hole pressure well 'Val', 500*barsa,... % Driving pressure 'Radius', .125*meter, ... % As per SPE10 Case B 'InnerProduct', 'ip_tpf',... % Use TPFA inner product 'Comp_i', [0, 1]); % Injects only gas % We want to place an producer in each corner of the domain, so we iterate % over the corners and add a 200 bar well at each corner. h_ind = [1, G.cartDims(1)]; v_ind = [1, G.cartDims(2)]; for i = 1:2 for j = 1:2 W = verticalWell(W, G, rock, h_ind(i), v_ind(j), [], ... 'Type', 'bhp',... 'Val', 200*barsa, ... 'Radius', .125*meter, ... 'InnerProduct', 'ip_tpf', ... 'Comp_i', [1, 0]); end end Compute one time values for the IMPES solverCompute TPFA transmissibilities Trans = computeTrans(G, rock);
% Compute pore volume
PV = poreVolume(G, rock);
Initialize a oil filled reservoir at 200 bar.Initially, it contains no gas x = initResSolComp(G, W, fluid, 200*barsa, [1, 0]); Do actual simulationsSimulate 100 days of production dT = 1*day; Nt = 100; f1 = figure(1); T = 0; tic() for kk = 1:Nt, x = impesTPFA(x, G, Trans, fluid, dT, PV, 'wells', W); T = T + dT; % We skip drawing the values for the cells containing wells, since they % are always outliers. notPerf = ~ismember(1:G.cells.num, vertcat(W.cells)); set(0, 'CurrentFigure', f1); clf; % Plot the pressure. subplot(1,3,1); plotCellData(G, (x.pressure), notPerf); title('Pressure') axis equal tight; % Plot the gas saturation subplot(1,3,2); plotCellData(G, log10(x.s(:,2)), notPerf); title('Saturation (Gas)') xlabel(sprintf('T = %1.0f days', convertTo(T, day))); axis equal tight; % Plot the x component of the permeability subplot(1,3,3); plotCellData(G, log10(rock.perm(:,1))); title('Permeability (x-direction)') axis equal tight; drawnow end toc() Elapsed time is 168.679231 seconds. |