This example aims to show the complete workflow for creating, running and analyzing a simulation model. Unlike the other examples, we will create all features of the model manually to get a self-contained script without any input files required.

The model we will setup is a two phase oil/water model with some compressibility and multiple wells. We begin by setting up the grid and rock structures. The grid is created by "makeModel3", which creates a structured model with intersecting faults.

We also define a layered permeability structure, with 300, 100 and 500 darcy in the lower, middle and top layers respectively.

Note that this example shows a simple conceptual model designed to show the workflow rather than a problem representing a realistic scenario in terms of well locations and fluid physics.

mrstModule add ad-core ad-blackoil ad-props mrst-gui % Define grid grdecl = makeModel3([50, 50, 5], [1000, 1000, 5]*meter); G = processGRDECL(grdecl); G = computeGeometry(G); % Set up permeability based on K-indices [I, J, K] = gridLogicalIndices(G); top = K < G.cartDims(3)/3; lower = K > 2*G.cartDims(3)/3; middle = ~(lower | top); px = ones(G.cells.num, 1); px(lower) = 300*milli*darcy; px(middle) = 100*milli*darcy; px(top) = 500*milli*darcy; % Introduce anisotropy by setting K_x = 10*K_z. perm = [px, px, 0.1*px]; rock = makeRock(G, perm, 0.3);

We will insert three producers, operating at fixed bottom hole pressure and perforated throughout all layers of the model along with a single injector injecting one pore volume over 10 years (the total simulation length). We also set up a schedule consisting of 5 small control steps initially, followed by 25 larger steps. We keep the well controls fixed throughout the simulation.

simTime = 10*year; nstep = 25; refine = 5; pv = poreVolume(G, rock); injRate = 1*sum(pv)/simTime; % Place wells offset = 10; W = []; % Producers W = verticalWell(W, G, rock, offset, offset, [],... 'Name', 'P1', 'comp_i', [1 0], 'Val', 250*barsa, 'Type', 'bhp'); W = verticalWell(W, G, rock, G.cartDims(1) - offset, G.cartDims(2) - offset, [],... 'Name', 'P2', 'comp_i', [1 0], 'Val', 250*barsa, 'Type', 'bhp'); W = verticalWell(W, G, rock, offset, G.cartDims(2) - offset, [], ... 'Name', 'P3', 'comp_i', [1 0], 'Val', 250*barsa, 'Type', 'bhp'); % Injectors W = verticalWell(W, G, rock, floor(G.cartDims(1)/2) + 4, floor(G.cartDims(1)/2) - 3, 1,... 'Name', 'Injector', 'comp_i', [1 0], 'Val', injRate, 'Type', 'rate'); % Compute the timesteps startSteps = repmat((simTime/(nstep + 1))/refine, refine, 1); restSteps = repmat(simTime/(nstep + 1), nstep, 1); timesteps = [startSteps; restSteps]; % Set up the schedule containing both the wells and the timesteps schedule = simpleSchedule(timesteps, 'Wells', W); % Plot horizontal permeability and wells clf; plotCellData(G, rock.perm(:, 1)/(milli*darcy)) plotWell(G, W) axis tight view(50, 50) colorbar title('K_x (mD)')

We set up a two phase simulation model based on automatic differentiation. We define a three phase fluid (with properties given for oil, water and gas) but the model will only use the oil and water values. Once we have defined the fluid, we further modify the oil compressibility by changing the bO-factor to have some constant compressibility using an anonymous function of pressure (using the standard Matlab builtin exp) to get a function on the form . Since the fluid model is a struct containing function handles, it is simple to modify the fluid to use alternate functions.

We then pass the fundamental structures (grid, rock and fluid) onto the two-phase oil/water model constructor.

fluid = initSimpleADIFluid('mu', [1, 5, 0]*centi*poise, ... 'rho', [1000, 700, 0]*kilogram/meter^3, ... 'n', [2, 2, 0]); c = 0.001/barsa; p_ref = 300*barsa; fluid.bO = @(p) exp((p - p_ref)*c); clf; p0 = (100:10:500)*barsa; plot(p0/barsa, fluid.bO(p0)) xlabel('Pressure (bar)') ylabel('Ratio') title('Reference to reservoir density ratio for the oil phase (bO)') gravity reset on model = TwoPhaseOilWaterModel(G, rock, fluid);

Once we have a model, we need to set up a initial state. We set up a very simple initial state. MRST uses water, oil, gas ordering internally, so in this case we have water in the first column and oil in the second for the saturations. We let the bottom part of the reservoir be completely water filled, and the top completely oil filled.

sW = ones(G.cells.num, 1); sW(G.cells.centroids(:, 3) < 5) = 0; sat = [sW, 1 - sW]; g = model.gravity(3); % Compute initial pressure p_res = p_ref + g*G.cells.centroids(:, 3).*[sW.*model.fluid.rhoWS + (1 - sW).*model.fluid.rhoOS]; state0 = initResSol(G, p_res, sat);

Once we have defined the schedule (dynamic controls and time), model (mathematical description of how to advance the solution) and the initial solution, we can simulate the problem.

[wellSols, states] = simulateScheduleAD(state0, model, schedule);

Solving timestep 01/30: -> 28 Days, 1 Hour, 3046 Seconds Solving timestep 02/30: 28 Days, 1 Hour, 3046 Seconds -> 56 Days, 3 Hours, 2492 Seconds Solving timestep 03/30: 56 Days, 3 Hours, 2492 Seconds -> 84 Days, 5 Hours, 1938 Seconds Solving timestep 04/30: 84 Days, 5 Hours, 1938 Seconds -> 112 Days, 7 Hours, 1384 Seconds Solving timestep 05/30: 112 Days, 7 Hours, 1384 Seconds -> 140 Days, 9 Hours, 830 Seconds ...

We launch a plotting tool for both the reservoir quantities (pressures and saturations, located in states) and the well solution (well rates and bottom hole pressures, located in wellSols).

clf plotToolbar(G, states) view(50, 50); % Wells simtime = cumsum(schedule.step.val); plotWellSols(wellSols, simtime, 'field', 'qOs');

The fine scale model has approximately 10000 cells. If we want a smaller model we can easily define an upscaled model. Here we set up a simple uniform partition of approximately 50 cells based on the IJK-indices of the grid.

mrstModule add coarsegrid cdims = [5, 5, 2]; p0 = partitionUI(G, cdims); close all plotCellData(G, mod(p0, 13), 'EdgeColor', 'none') axis tight off view(50, 50)

We see that several coarse blocks cross the fault lines. To get hexahedral coarse blocks, we create a grid where the faults act as barriers and apply the "processPartition" routine to split any coarse blocks intersected by faults.

Afterwards, we show the new partition and highlight blocks created due to the modification of the fault.

G_fault = makeInternalBoundary(G, find(G.faces.tag > 0)); p = processPartition(G_fault, p0); clf plotCellData(G, mod(p, 13), 'EdgeColor', 'none') plotGrid(G, p ~= p0, 'EdgeColor', 'r', 'LineWidth', 2, 'FaceColor', 'none') axis tight off view(50, 50) title('Partition after processing by faults');

We can now directly upscale the model, schedule and initial state automatically. By default, the upscaling routines uses the simplest possible options, i.e. harmonic averaging of permeabilities. It is possible to use more advanced options, but for the purpose of this example we will use the defaults.

Once we have a upscaled model, we can again simulate the new schedule and observe that the time taken is greatly reduced. For instance, on a Intel Core i7 desktop computer, the fine model takes little under a minute to run, while the upscaled model takes 4 seconds.

model_c = upscaleModelTPFA(model, p); G_c = model_c.G; rock_c = model_c.rock; schedule_c = upscaleSchedule(model_c, schedule); state0_c = upscaleState(model_c, model, state0); [wellSols_c, states_c] = simulateScheduleAD(state0_c, model_c, schedule_c);

We plot the coarse solutions and compare the well solutions. Note that the upscaling will result in only 70 cells, which is unlikely to give good results with only simple harmonic averaging of permeabilities.

clf plotToolbar(G_c, states_c) view(50, 50); plotWellSols({wellSols, wellSols_c}, simtime, 'DatasetNames', {'Fine scale', 'Upscaled'}, 'Field', 'qOs');

As an alternative to looking at well curves, we can also look at the flow diagnostics of the models. Flow diagnostics are simple routines based on time-of-flight and tracer equations, which aim to give a qualitative understanding of reservoir dynamics. Here, we take the end-of-simulation states as a snapshot for both the fine and coarse model and compute time-of-flight and well tracers.

mrstModule add diagnostics D = computeTOFandTracer(states{end}, G, rock, 'Wells', schedule.control.W); D_c = computeTOFandTracer(states_c{end}, G_c, rock_c, 'Wells', schedule_c.control.W);

We plot the sum of forward and backwards time of flight, showing which regions have large amounts of flow. Large values indicate that a cell has little flow going through it.

Since the values vary by several orders of magnitude, we take the transform before plotting. We also use the same color axis to ensure that the plots can be compared.

close all figure(1); clf plotCellData(G, log10(sum(D.tof, 2))); view(50, 50); title('Log of total travel time, fine model'); c = caxis(); figure(2); clf plotCellData(G_c, log10(sum(D_c.tof, 2))); view(50, 50); title('Log of total travel time, coarse model'); caxis(c)

We can also look at the tracer partitioning for the producers, showing the drainage regions for the different wells.

See the diagnostics module for more examples and more in-depth discussions of how flow diagnostics can be used.

figure(1); clf plotCellData(G, D.ppart); view(50, 50); title('Drainage regions, fine model'); figure(2); clf plotCellData(G_c, D_c.ppart); view(50, 50); title('Drainage regions, coarse model');

We can also examine the diagnostics interactively using the diagnostics viewer.

close all; interactiveDiagnostics(G, rock, schedule.control.W, 'state', states{end}, 'computeFlux', false, 'name', 'Fine model'); view(50, 50); interactiveDiagnostics(G_c, rock_c, schedule_c.control.W, 'state', states_c{end}, 'computeFlux', false, 'name', 'Coarse model'); view(50, 50);

Published May 1, 2015

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