Function equationOWPolymer.m

Set up the polymer equations

In this function, the modeling equations (in this case mass conservation of oil, water and polymer) are assembled. The function returns the residual vector (one value per cell) and the Jacobians with respect to the main variables (pressure, saturation and polymer concentration).

```function eqs = equationOWPolymer(state0, state, dt, bc, system, param)
```
```   sys = system; % short-cut
f = sys.fluid;
poro = sys.rock.poro;
```

Initialize then main variables as Automatic Differentiation (AD) variables

```   p    = state.pressure;
sW   = state.s(:,1);
c    = state.c;
cmax = state.cmax;
nc   = numel(p);
[p, sW, c] = initVariablesADI(p, sW, c);
```

Initialize previous time step variables.

```   p0  = state0.pressure;
sW0 = state0.s(:,1);
c0  = state0.c;
cmax0 = state0.cmax;
```

Compute formation volume factors

```   bW     = f.bW(p);
bO     = f.bO(p);
bc_bW  = f.bW(bc.dirichlet.pressure); % boundary values
bc_bO  = f.bO(bc.dirichlet.pressure); %
```

```   dpO = sys.grad(p, bc.dirichlet.pressure);
dpW = dpO; % no capillary pressure
```

Computation of faces formation volume factor, given as average.

```   fbW = sys.averageFaceValues(bW, bc_bW);
fbO = sys.averageFaceValues(bO, bc_bO);
```

Computation of the mobilities, in the cells. The mobilities at the dirichlet boundary have been already computed.

For an implicit solver, we compute the mobility using the current values for pressure, saturation and polymer concentration. For an explicit solver, we use the previous values. Note that this is the only modifications which are needed in the code to switch between an implicit and explicit transport solver.

```   if system.implicit_transport
[mobW, mobO, mobP] = computeMobilities(p, sW, c, cmax, f);
else
[mobW, mobO, mobP] = computeMobilities(p0, sW0, c0, cmax0, f);
end
```

Computation of upstream phases mobilities.

```   wupc = (double(dpW)>=0);
oupc = (double(dpO)>=0);
fmobW = sys.upwindFaceValues(wupc, mobW, bc.mobW);
fmobO = sys.upwindFaceValues(oupc, mobO, bc.mobO);
fmobP = sys.upwindFaceValues(wupc, mobP, bc.mobP);
```

Computation of fluxes

```   vW = fbW.*fmobW.*sys.T.*dpW;
vO = fbO.*fmobO.*sys.T.*dpO;
vP = fbW.*fmobP.*sys.T.*dpW;
```

Computation of accumulation terms for oil, water and polymer The accumulation term for polymer includes adsorption

```   accumO = (sys.pv/dt).*bO.*(1 - sW);
accumW = (sys.pv/dt).*bW.*sW;
accumP = (sys.pv/dt).*bW.*sW.*c + (sys.pv/dt).*(f.rhoR* ( (1-poro)./poro)).* ...
```

Accumulation term for previous time step.

```   accumO0 = 1/dt*bO.*sys.pv.*(1 - sW0);
accumW0 = 1/dt*bW.*sys.pv.*sW0;
accumP0 = (sys.pv/dt).*bW.*sW0.*c0 + (sys.pv/dt).*(f.rhoR* ( (1-poro)./poro)).* ...
```

Compute source term.

```   sInj = bc.injection.s;
cInj = bc.injection.c;
influxO = zeros(nc, 1);
influxW = zeros(nc, 1);
influxP = zeros(nc, 1);
influxO(bc.injection.cells) = (1 - sInj)*bc.injection.rate;
influxW(bc.injection.cells) = sInj*bc.injection.rate;
influxP(bc.injection.cells) = cInj*sInj*bc.injection.rate;
```

assemble the equations as the sum of accumulation terms, flux terms and influxes.

```   eqs{1} = accumO - accumO0 + sys.div(vO) - influxO;
eqs{2} = accumW - accumW0 + sys.div(vW) - influxW;
eqs{3} = accumP - accumP0 + sys.div(vP) - influxP;
```
```end
```

Published November 6, 2014