Function equationCompositional.m
Residual equations functionsIn this function, the residual equations are assembled and evaluated. By using Automatic Differentiation (AD), we compute automatically the derivatives of the residual. Examples of use of automatic differentiation can be found in MRST book. The governing equations are the mass conservation equations for the 4 components: He, Ne, CO_{2} and water, and the equilibrium equations, or flash equations. The mass conservation equations are discretized in space using a Two Point Flux Approximation. The transmissibilities are computed before hand when the system is set up in the function setupSystem. The equations are solved implicitly in time. function eqs = equationCompositional(state0, state, dt, bc, system, varargin)
opt = struct('Verbose', mrstVerbose,... 'scaling', [],... 'history', [], ... 'iteration', 1, ... 'stepOptions', []); opt = merge_options(opt, varargin{:}); s = system.s; f = system.fluid; Temp = system.Temp; T = s.T; % Transmissibilities. nComp = system.nComp; G = system.G; The main variables are pressure, liquid saturation and total concentrations. They are converted to AD class objects by the function initVariablesADI. Since they are initialized as main variables, their Jacobian are identity matrices. By operator overloading, any expression which contains these variables (p, sL and C) will also be an AD object and have a derivative. In this way, the derivatives of the residuals are automatically obtained. p = state.pressure;
sL = state.sL;
C = state.C;
[p, C{:}, sL] = initVariablesADI(p, C{:}, sL);
p0 = state0.pressure;
sL0 = state0.sL;
C0 = state0.C;
% The previous values should be stored in state
[cg0, cl0] = computeComposition(C0, p0, sL0, system);
[Cw0, cwg0, cwl0] = computeWaterComp(p0, C0, cl0, sL0, system);
Compute the relative permeabilities. We also need the values at the boundaries where Dirichlet conditions hold. [krL, krG] = f.relPerm(sL); % Liquid, gas relative permeability bd = bc.dirichlet; % shortcut for dirichlet bc values. [bc_krL, bc_krG] = f.relPerm(bd.s(:, 2)); % Liquid and gas relative permeability as % the boundary Set up the gravity term. The molar mass are needed. g = norm(gravity); dz = s.dz; mmW = f.mmW; % f.mmW returns molar mass for water. mmC = f.mmC; % f.mmC returns a vector of molar mass for each component. Compute the mobilities mobL = krL./f.muL; mobG = krG./f.muG; bc_mobL = bc_krL./f.muL; bc_mobG = bc_krG./f.muG; Compute the upstream direction for each component. dpC = cell(nComp, 1); upC = cell(nComp, 1); for ic = 1:nComp dpC{ic} = s.p_grad(p)  g*(mmC(ic).*dz); upC{ic} = (double(dpC{ic})>=0); end dpW = s.p_grad(p)  g*(mmW.*dz); upW = (double(dpW)>=0); Components the composition of the phases: the gas cg and liquid cl concentrations. Note that cg and cl are AD objects and contain the Jacobian with respect to the main AD variables C, p and sL. [cg, cl] = computeComposition(C, p, sL, system); fluxC = cell(nComp, 1); Compute the residual of the mass conservation equation for each component (He, Ne and CO_{2}) for ic = 1 : nComp
The function s.faceConcentrations computes concentration on the faces given cell concentrations, using upwind directions. bc_val = bd.cg{ic}.*bc_mobG + bd.cl{ic}.*bc_mobL; fluxC{ic} = s.faceConcentrations(upC{ic}, cg{ic}.*mobG + cl{ic}.*mobL, bc_val); eqs{ic} = (s.pv/dt).*(C{ic}  C0{ic}) + s.div(fluxC{ic}.*T.*dpC{ic}); end
Compute the residual of the mass conservation equation for water. [Cw, cwg, cwl] = computeWaterComp(p, C, cl, sL, system); bc_val = bd.cw(:, 1).*bc_mobG + bd.cw(:, 2).*bc_mobL; fluxW = s.faceConcentrations(upW, cwg.*mobG + cwl.*mobL, bc_val); eqs{nComp + 1} = (s.pv/dt).*(Cw  Cw0) + s.div(fluxW.*T.*dpW); Compute the residual for the flash equation for the saturation. Note that, since the function equationCompositional is always called after the flash are equations are solved, this residual is always equal to zero in every cells. But, the residual variable is an AD objects which also contains the derivative. For this particular residual, we obtain in this way, given in an implicit form, the derivative of the liquid saturation sL with respect to the total concentrations C and the pressure p. vp = system.vp; eqs{nComp + 2} = (1  p/vp); R = system.R; k = system.k; for ic = 1:nComp eqs{nComp + 2} = eqs{nComp + 2} + 1/vp*R*Temp.*C{ic}./(1 + (R*Temp/k(ic)  1).*sL); end Check the status for each cell
st = getCellStatus(state, system); We do not handle the case only gas is present. assert(all(st ~= 2), 'pure gas phase case is not handled');
For the cells where only liquid is present (st == 1), the liquid saturation equation is simly sL = 1. is_st_one = (st == 1); if any(is_st_one) eqs{nComp + 2}(is_st_one) = sL(is_st_one)  1; end Add input fluxes. for ic = 1 : nComp eqs{ic}(bc.influx_cells) = eqs{ic}(bc.influx_cells)  bc.C_influx{ic}; end eqs{nComp + 1}(bc.influx_cells) = eqs{nComp + 1}(bc.influx_cells)  bc.water_influx; end function st = getCellStatus(state, system) This function computes the status in each cell
[R, k, Temp, vp, nComp] = deal(system.R, system.k, system.Temp, system.vp, system.nComp); p = state.pressure; Initiate status by assuming that gas and liquid are both present. st = 3*ones(numel(p), 1); Set pure liquid status, if pressure is above threshold. C = cell2mat({state.C{:}}); st(C*k + vp <= p) = 1; Set pure liquid status, if pressure is below threshold. omega_g = R*Temp*ones(nComp, 1); min_pressure = C*omega_g + vp; st(p < min_pressure) = 2; end
