Modelling coupled miscible fluid flow and transport processes

% Domain size along X, Y & Z directions

dom = RecDomain([365.76,670.56,0.6096]);

Dx = dom.Lx; Dy = dom.Ly; Dz = dom.Lz;

Grid.Nx = 60; Grid.Ny = 220; Grid.Nz = 1;

Nx = Grid.Nx; Ny = Grid.Ny; Nz = Grid.Nz;

Grid.hx = (Dx/Nx);

Grid.hy = (Dy/Ny);

Grid.hz = (Dz/Nz);

Grid.N = Grid.Nx*Grid.Ny*Grid.Nz;

N = Grid.N;

Grid.V = (Dx/Nx)*(Dy/Ny)*(Dz/Nz); % Cell volumes

Grid.K = 1.4e-11.*ones(3,Nx,Ny,Nz); % Unit-Darcy permeability

Grid.compr = 4.4e-10.*ones(Grid.Nx,Grid.Ny,Grid.Nz); % compressibility

load '..\data\spe_layer1_phi.mat' phi_layer;

phi = phi_layer; clear phi_layer;

phi = max (phi(:), 1e-3); % to avoid zero porosity cells

Grid.por = reshape(phi',Nx,Ny,Nz);

load '..\data\spe_layer1_perm.mat' perm_layer;

Grid.K = perm_layer*9.869232667160130e-13;

Grid.K = reshape(Grid.K',3,Nx,Ny,Nz);

G = cartGrid(0:Dx/Nx:Dx, 0:Dy/Ny:Dy);

figure;

subplot(1,2,1);

PlotCellData(G, phi, 'EdgeColor','none');

%colormap(jet(16));

axis tight equal; colorbar, title('Porosity field');

subplot(1,2,2);

PlotCellData(G, reshape(log10(Grid.K(1,:,:,:)),N,1), 'EdgeColor','none');

clear phi;

% Cell-centered injection/production flow rates

flow_rate = 100/3600;

Qw = zeros(N,1);

Qw([1 N]) = [+flow_rate -flow_rate];

% fluid properties units: density (Kg/m^3), initial viscosity (Pa.s)

water = Fluid('water',[1000, 1e-3]);

P = 100e5.*ones(Nx,Ny,Nz);

day = 3600*24; % seconds/day

T = 5*day; % 5 days

nt = 30; % number of time steps

dt = T/nt; % time step

opt.tol = 1e-9;

opt.maxiter = 50;

D = 1e-9.*ones(3,Nx,Ny,Nz); % uniform diffusion coefficient through the domain

C1 = zeros(N,1); C1(1) = 1; % unit input concentration at the injector & zero otherwise

Grid.sat = ones(Nx,Ny,Nz); % unit saturation everywhere

tic; % start timing the time marching loop

% call TPFA flow solver for single-phase flow (we assume steady state flow here)

[~,V] = Pressure(Grid,water,Qw);

for t=1:nt

fprintf('\nSolving solute transport problem. Time = %f days\n', t*dt/day);

% call the fully implicit first-order upwind convective transport for a half time step

C_bar = ImplicitConcentration(Grid,C1,V,Qw,dt/2,opt);

% solve implicitly the diffusion problem for this time step

[C_bar,~] = Diffusion(Grid,D,C_bar,dt);

% do advection step for the remaining half time step

C1 = ImplicitConcentration(Grid,C_bar,V,Qw,dt/2,opt);

% note: C_bar is the intermediate concentration after the first advection step

% and the diffusive step

figure;

% plot solute concentration for this first simulation

contourf(linspace((Dx/Nx)/2,Dx-(Dx/Nx)/2,Nx),...

linspace((Dy/Ny)/2,Dy-(Dy/Ny)/2,Ny),...

reshape(C1,Nx,Ny)','LevelList',0:0.1:1);

colormap(jet(16));

visc = @(c) 1e-3*(1+10.*c);

C2 = zeros(N,1); C2(1) = 1;

tic; % start timing the time marching loop

for t=1:nt

fprintf('\nUpdating the flow.\n');

% viscosity recalculation

water.viscosity = visc(C2);

% call TPFA flow solver for single-phase flow

[P,V] = Pressure(Grid,water,Qw);

fprintf('\nSolving solute transport problem. Time = %f days\n', t*dt/day);

% call the fully implicit first-order upwind convective transport for a half time step

C_bar = ImplicitConcentration(Grid,C2,V,Qw,dt/2,opt);

% solve implicitly the diffusion problem for this time step

[C_bar,~] = Diffusion(Grid,D,C_bar,dt);

% do advection step for the remaining half time step

C2 = ImplicitConcentration(Grid,C_bar,V,Qw,dt/2,opt);

figure;

% plot solute concentration

contourf(linspace((Dx/Nx)/2,Dx-(Dx/Nx)/2,Nx),...

linspace((Dy/Ny)/2,Dy-(Dy/Ny)/2,Ny),...

reshape(C2,Nx,Ny)','LevelList',0:0.1:1);

colormap(jet(16));

visc = @(c) 1e-3*(1-10.*c);

C3 = zeros(N,1); C3(1) = 1;

tic; % start timing the time marching loop

for t=1:nt

fprintf('\nUpdating the flow.\n');

% viscosity recalculation

water.viscosity = visc(C3);

% call TPFA flow solver for single-phase flow

[P,V] = Pressure(Grid,water,Qw);

fprintf('\nSolving solute transport problem. Time = %f days\n', t*dt/day);

% call the fully implicit first-order upwind convective transport for a half time step

C_bar = ImplicitConcentration(Grid,C3,V,Qw,dt/2,opt);

% solve implicitly the diffusion problem for this time step

[C_bar,~] = Diffusion(Grid,D,C_bar,dt);

% do advection step for the remaining half time step

C3 = ImplicitConcentration(Grid,C_bar,V,Qw,dt/2,opt);

figure;

% plot solute concentration

contourf(linspace((Dx/Nx)/2,Dx-(Dx/Nx)/2,Nx),...

linspace((Dy/Ny)/2,Dy-(Dy/Ny)/2,Ny),...

reshape(C3,Nx,Ny)','LevelList',0:0.1:1);

colormap(jet(16));