Help required in code for custom tracer concentration profile input

[piyush mohanty]

Hello,

I intend to show tracer flow with my own data for a time dependent tracer concentration profile. I curve fitted my data in a Fourier series form but I was unable to get the result by inputting a time dependent form of c in the loop. Please help in this regard

%% Define and discretize the domain

[nx,ny,nz] = deal(50, 1, 1);

[Lx,Ly,Lz] = deal(500*ft, 1*ft, 1*ft);

G = cartGrid([nx, ny, nz], [Lx, Ly, Lz]);

G = computeGeometry(G);

plotGrid(G); view(3); axis tight equal

a0=5.5802; % Constants for curve fitted fourier data

a1=4.1532;

b1=-1.9385;

a2=-0.2452;

b2=-1.3736;

a3=0.3142;

b3=-0.1856;

a4=-0.2297;

b4=-0.2484;

a5=-0.0054;

b5=0.2615;

a6=-0.0096;

b6=0.1223;

w=1.6465;

%% Define rock model

rock = makeRock(G, 30*milli*darcy, 0.3);

cr = 5e-6/psia;

Pref = 3000*psia;

Dref = 0.5*ft;

PVref = poreVolume(G, rock);

PVfunc = @(P) PVref.*exp(cr*(P - Pref));

P = linspace(1500*psia,Pref,50);

plot(P/psia, PVref(1)*exp(cr*(P-Pref))/ft^3,'LineWidth',2);

xlabel('Pressure [psia]')

ylabel('Single Gridblock Pore Volume [ft^3]')

%% Define model for an slightly compressible fluid

mu = 1*centi*poise;

cFluid = 3e-6/psia;

rhoRef = 62*pound/ft^3; % reservoir fluid density at the reference conditions

rhoFunc = @(P) rhoRef*exp(cFluid*(P - Pref));

plot(P/psia,rhoFunc(P)/(pound/ft^3),'LineWidth',2);

xlabel('Pressure [psia]')

ylabel('Fluid Density [lbm/ft^3]')

%% Specify boundary conditions

nc = G.cells.num;

W = addWell([], G, rock, 1, 'Name', 'Injection Well' , 'Radius', 0.5*ft);

W = addWell(W , G, rock, nc, 'Name', 'Production Well', 'Radius', 0.5*ft);

injCell = W(1).cells;

prodCell = W(2).cells;

PIprod = W(2).WI; % well productivity indices

dzProd = W(2).dZ; % perforated cell depth relative to bottom-hole

% Internal boundary conditions

qInj = 0.1*stb/day;

BHPprod = Pref;

influxConc = 1*mol/kilogram; % mol of solvent in 1 kilogram of solution

gravity reset on, g = norm(gravity);

Pperf = BHPprod + rhoFunc(BHPprod)*g.*dzProd; % wellbore pressures at perforated points

%% Set initial pressure and concentration

PG = rhoFunc(Pref)*g;

Pinit = Pref + PG.*(G.cells.centroids(:,3) - Dref);

Cinit = 0*mol/kilogram*ones(nc,1);

%% Compute transmissibilities

N = double(G.faces.neighbors); % Map: face -> cell

intInx = all(N ~= 0, 2); % Interior faces

N = N(intInx, :); % Interior neighbors

hT = computeTrans(G, rock); % Half-transmissibilities

cf = G.cells.faces(:,1); % Map: cell -> face

nf = G.faces.num; % Total number of faces

T = 1 ./ accumarray(cf, 1 ./ hT, [nf, 1]); % Harmonic average

T = T(intInx); % Restricted to interior

%% Define discrete operators to desretize the mathematical equations

n = size(N,1);

C = sparse( [(1:n)'; (1:n)'], N, ones(n,1)*[-1 1], n, G.cells.num);

grad = @(x)C*x;

div = @(x)-C'*x;

avg = @(x) 0.5 * (x(N(:,1)) + x(N(:,2)));

upwind = @(x,flag) x(N(:,1)).*double(flag)+x(N(:,2)).*double(~flag);

%% Initialize for solution loop

[P_AD, C_AD] = initVariablesADI(Pinit, Cinit);

[pIx, cIx] = deal(1:nc, (nc+1):(2*nc));

numSteps = 500; % number of time-steps

totTime = 270*day; % total simulation time

dt = totTime / numSteps; % constant time step

tol = 1e-5; % Newton tolerance

maxIter = 10; % max number of Newton-Raphson iterations

sol = repmat(struct('time',[],'pressure',[], 'C', []),[numSteps+1,1]);

sol(1) = struct('time', 0, 'pressure', value(P_AD)/psia, 'C', value(C_AD));

%% Main loop

t = 0; step = 0;

hwb = waitbar(t,'Simulation ...');

while t < totTime

t = t + dt;

d = (t-177.7)/(55.2);

step = step + 1;

fprintf('\nTime step %d: Time %.2f -> %.2f days\n', ...

step, convertTo(t - dt, day), convertTo(t, day));

c = a0 + a1*cos(d*w) + b1*sin(d*w) + a2*cos(2*d*w) + b2*sin(2*d*w) + a3*cos(3*d*w) + b3*sin(3*d*w) + a4*cos(4*d*w) + b4*sin(4*d*w) + a5*cos(5*d*w) + b5*sin(5*d*w) + a6*cos(6*d*w) + b6*sin(6*d*w);

P0 = value(P_AD);

C0 = value(c/1000000);

% Newton loop

resNorm = 1e+99;

iterCounter = 0;

while (resNorm > tol) && (iterCounter <= maxIter)

%% Define pressure equation

[rho_AD, rho0] = deal(rhoFunc(P_AD), rhoFunc(P0));

[PV_AD , PV0] = deal(PVfunc(P_AD), PVfunc(P0));

gradz = grad(G.cells.centroids(:,3));

v = -(T/mu).*(grad(P_AD) - avg(rho_AD)*g.*gradz);

pressEq = (1/dt)*(PV_AD.*rho_AD - PV0.*rho0) + div(avg(rho_AD).*v);

% Adding boundary conditions

qProd = -PIprod.*(1/mu).*(P_AD(prodCell) - Pperf);

% The BC of the production well can be set in 2 different ways, either using the BHP directly or using the equivalent production flowrate.

pressEq(prodCell) = pressEq(prodCell) - qProd*rho_AD(prodCell); % Method 1

% pressEq(prodCell) = P_AD(prodCell) - BHPprod; % Method 2

pressEq(injCell) = pressEq(injCell) - qInj*rhoRef;

%% Define the mass transport equation of the desired species

% water upstream-index

upcw = v > 0;

accumTerm = (1/dt).*(PV_AD.*rho_AD.*C_AD - PV0.*rho0.*C0);

% Flux term without dispersion

% fluxTerm = div(upwind(C_AD.*rho_AD, upcw).*v);

% Flux term with dispersion

D = 10^(-9)*meter^2/second*ones;

fluxTerm = div(upwind(C_AD.*rho_AD, upcw).*v - avg(rock.poro.*rho_AD).*D.*grad(C_AD));

transEq = accumTerm + fluxTerm;

% Add boundary conditions

% transEq(injCell) = C_AD(injCell) - influxConc; % Method 1

transEq(injCell) = transEq(injCell) - qInj*rhoRef*influxConc; % Method 2

transEq(prodCell) = transEq(prodCell) - qProd*rho_AD(prodCell)*C_AD(prodCell);

%% Collect and concatenate all equations (i.e., assemble and linearize system)

eqs = {pressEq, transEq};

eq = cat(eqs{:});

J = eq.jac{1}; % Jacobian

res = eq.val; % residual

increment = -(J \ res); % Newton update

%% Update variables

P_AD.val = P_AD.val + increment(pIx);

C_AD.val = C_AD.val + increment(cIx);

resNorm = norm(res);

iterCounter = iterCounter + 1;

fprintf(' Iteration %3d: Res = %.4e\n', iterCounter, resNorm);

end

if iterCounter > maxIter

error('Newton solver did not converge')

else % Store solution

sol(step+1) = struct('time', t, 'pressure', value(P_AD)/psia, 'C', value(C_AD));

waitbar(t/totTime,hwb);

end

end

close(hwb);

%% For a 1-dimensional domain

clf

for i = 1:numel(sol)

plot(1:nc, sol(i).C, 'LineWidth', 1)

title(['Tracer Concentration Profile (', 'Time Step: ', num2str(i),')'])

xlim([1,nc])

ylim([Cinit(1), influxConc])

xlabel('Gridblock Index')

ylabel('Tracer Concentration [mole/kilogram]')

pause(0.25)

end