Coupled nonlinear PDEs: two phase flow in porous media
In this post, I'm going to quickly explain how to solve a system of nonlinear equations. In the previous posts (here and here), I explained how to linearize a nonlinear PDE. Now, we are going to linearize two PDEs, which means we need to use the Taylor expansion again, which for a function with two independent variables reads: $$f\left(x,y\right)=f\left(x_{0},y_{0}\right)+\frac{\partial f}{\partial x}\left(x_{0},y_{0}\right)\left(x-x_{0}\right)+\frac{\partial f}{\partial y}\left(x_{0},y_{0}\right)\left(y-y_{0}\right)+h.o.t$$ We are going to use the above equation to linearize the following equations, which describe the two phase flow (oil, and water) in porous media. We ignore the gravity term and assume that the fluids and the porous medium are incompressible: $$\nabla.\left(-k\left(\frac{k_{rw}\left(S_{w}\right)}{\mu_{w}}+\frac{k_{ro}\left(S_{w}\right)}{\mu_{o}}\right)\nabla p\right)=0,$$ and $$\varphi\frac{\partial S_{w}}{\partial t}+\nabla.\left(-\frac{k_{rw}\left(S_{w}\right)k}{\mu_{w}}\nabla p\right)=0.$$ The linearized form, after some algebraic operations, read $$\nabla.\left(-k\left(\frac{k_{rw}\left(S_{w,0}\right)}{\mu_{w}}+\frac{k_{ro}\left(S_{w,0}\right)}{\mu_{o}}\right)\nabla p\right)+\\\nabla.\left(k\left(\frac{\partial k_{rw}\left(S_{w,0}\right)}{\mu_{w}\partial S_{w}}+\frac{\partial k_{ro}\left(S_{w,0}\right)}{\mu_{o}\partial S_{w}}\right)\nabla p_{0}S_{w}\right)=\\\nabla.\left(k\left(\frac{\partial k_{rw}\left(S_{w,0}\right)}{\mu_{w}\partial S_{w}}+\frac{\partial k_{ro}\left(S_{w,0}\right)}{\mu_{o}\partial S_{w}}\right)\nabla p_{0}S_{w,0}\right)$$ and $$\varphi\frac{\partial S_{w}}{\partial t}+\nabla.\left(-\frac{k_{rw}\left(S_{w,0}\right)k}{\mu_{w}}\nabla p\right)+\\\nabla.\left(-\frac{k}{\mu_{w}}\frac{\partial k_{rw}\left(S_{w,0}\right)}{\partial S_{w}}\nabla p_{0}S_{w}\right)=\\\nabla.\left(-\frac{k}{\mu_{w}}\frac{\partial k_{rw}\left(S_{w,0}\right)}{\partial S_{w}}\nabla p_{0}S_{w,0}\right)$$ We will discretize the above equations using the FVTool, which gives us a system of linear equations, i.e., $$\left[\begin{array}{cc} J_{p1}+M_{bc,p} & J_{S_{w}1}\\ J_{p2} & J_{S_{w}2}+M_{bc,S_{w}} \end{array}\right]\left[\begin{array}{c} \mathbf{p}\\ \mathbf{S_{w}} \end{array}\right]=\left[\begin{array}{c} RHS_{1}+RHS_{bc,p}\\ RHS_{2}+RHS_{bc,S_{w}} \end{array}\right]$$ Subscript 1 denotes the first equation (continuity) and 2 denotes the second equation (water saturation balance). Note how we add the boundary condition terms to the equations. I'm not going to explain it in more details, because I'm too tired and it's 1:00 o'clock, but I give you my word that this is the right way to do it! For each time step, the above equation needs to be solved in a loop until the new values of $S_w$ converges to the initial estimates $S_{w,0}$. Let me start the implementation:
% Go to the FVTool folder and initialize the toolbox
cd('/home/ali/MyPackages/FVTool/')
FVToolStartUp()
Rel-perm curves¶
First, we need to define the relative permeability functions, and their derivatives. I then plot the functions, just for the fun of it.
krw0 = 1;
kro0 = 1;
nw = 2;
no = 2;
krw = @(sw)(krw0*sw.^nw);
dkrwdsw = @(sw)(krw0*nw*sw.^(nw-1));
kro = @(sw)(kro0*(1-sw).^no);
dkrodsw = @(sw)(-kro0*no*(1-sw).^(no-1));
lw = geometricMean(k)/mu_water;
lo = geometricMean(k)/mu_oil;
Lw = @(sw)(krw(sw));
Lo = @(sw)(k/mu_oil*kro(sw));
dLwdsw = @(sw)(k/mu_water*dkrwdsw(sw));
dLodsw = @(sw)(k/mu_oil*dkrodsw(sw));
sw_0=linspace(0,1,100);
plot(sw_0, krw(sw_0), sw_0, kro(sw_0), '--');
xlabel('S_w')
ylabel('Rel-perms')
legend('water', 'oil');
Create the mesh¶
Now, we create a mesh and initialize all the variables over the mesh, except for pressure and saturation. I use a large value of Dykstra-Parsons coefficient to create a hetrogeneous permeability field. Let's visualize it as well.
Nx = 100; % number of cells in x direction
Ny = 70; % number of cells in y direction
W = 50; % [m] length of the domain in x direction
H = 30; % [m] length of the domain in y direction
m = createMesh2D(Nx, Ny, W, H); % creates a 2D mesh
%% define the physical parametrs
p0 = 1e5; % [bar] pressure
pin = 50e5; % [bar] injection pressure at the left boundary
sw0 = 0; % initial water saturation
swin = 1;
mu_oil = 10e-3; % [Pa.s] oil viscosity
mu_water = 1e-3; % [Pa.s] water viscosity
% reservoir
k0 = 2e-12; % [m^2] average reservoir permeability
phi0 = 0.2; % average porosity
clx=0.05;
cly=0.05;
V_dp=0.4; % Dykstra-Parsons coef.
perm_val= field2d(Nx,Ny,k0,V_dp,clx,cly);
k=createCellVariable(m, perm_val);
phi=createCellVariable(m, phi0);
visualizeCells(k);
title('Permeability field')
Define the boundary and initial conditions¶
We want to solve a problem, where the domain is initially saturated with oil (Sw=0). We inject water (Sw=1) from the left boundary, with an injection pressure of 50 bar. The right boundary is at a constant pressure of 1 bar.
BCp = createBC(m); % Neumann BC for pressure
BCs = createBC(m); % Neumann BC for saturation
% change the left and right boandary to constant pressure (Dirichlet)
BCp.left.a(:)=0; BCp.left.b(:)=1; BCp.left.c(:)=pin;
BCp.right.a(:)=0; BCp.right.b(:)=1; BCp.right.c(:)=p0;
% change the left boundary to constant saturation (Dirichlet)
BCs.left.a(:)=0; BCs.left.b(:)=1; BCs.left.c(:)=1;
%% define the variables
sw_old = createCellVariable(m, sw0, BCs);
p_old = createCellVariable(m, p0, BCp);
sw = sw_old;
p = p_old;
Define the nonlinear solver specifications¶
We don't use a nonlinear solver. We write one. The accuracies are defined below:
% define the time step and solver properties
dt = 1000; % [s] time step
t_end = 50*dt; % [s] final time
eps_p = 1e-5; % pressure accuracy
eps_sw = 1e-5; % saturation accuracy
uw = -gradientTerm(p_old); % an estimation of the water velocity
Main loop¶
Here we write the external time loop and the internal nonlinear solver loop. I've played with the time steps a bit to find a suitable value. You can adapt the time step inside the loop if you prefer.
t = 0;
while (t<t_end)
error_p = 1e5;
error_sw = 1e5;
% Implicit loop
while ((error_p>eps_p) || (error_sw>eps_sw))
% calculate parameters
pgrad = gradientTerm(p);
sw_face = upwindMean(sw, -pgrad); % average value of water saturation
labdao = lo.*funceval(kro, sw_face);
labdaw = lw.*funceval(krw, sw_face);
dlabdaodsw = lo.*funceval(dkrodsw, sw_face);
dlabdawdsw = lw.*funceval(dkrwdsw, sw_face);
labda = labdao+labdaw;
dlabdadsw = dlabdaodsw+dlabdawdsw;
% compute [Jacobian] matrices
Mdiffp1 = diffusionTerm(-labda);
Mdiffp2 = diffusionTerm(-labdaw);
Mconvsw1 = convectionUpwindTerm(-dlabdadsw.*pgrad);
Mconvsw2 = convectionUpwindTerm(-dlabdawdsw.*pgrad);
[Mtranssw2, RHStrans2] = transientTerm(sw_old, dt, phi);
% Compute RHS values
RHS1 = divergenceTerm(-dlabdadsw.*sw_face.*pgrad);
RHS2 = divergenceTerm(-dlabdawdsw.*sw_face.*pgrad);
% include boundary conditions
[Mbcp, RHSbcp] = boundaryCondition(BCp);
[Mbcsw, RHSbcsw] = boundaryCondition(BCs);
% Couple the equations; BC goes into the block on the main diagonal
M = [Mdiffp1+Mbcp Mconvsw1; Mdiffp2 Mconvsw2+Mtranssw2+Mbcsw];
RHS = [RHS1+RHSbcp; RHS2+RHStrans2+RHSbcsw];
% solve the linear system of equations
x = M\RHS;
% separate the variables from the solution
p_new = reshapeCell(m,full(x(1:(Nx+2)*(Ny+2))));
sw_new = reshapeCell(m,full(x((Nx+2)*(Ny+2)+1:end)));
% calculate error values
error_p = max(abs((p_new(:)-p.value(:))./p_new(:)));
error_sw = max(abs(sw_new(:)-sw.value(:)));
% assign new values of p and sw
p.value = p_new;
sw.value = sw_new;
end
t=t+dt;
p_old = p;
sw_old = sw;
end
visualizeCells(sw);shading interp;
This is a very simple example. There are a few functions that I have used here, like reshapeCell which should be explained later in the documents. All in all, we could show that with a few lines of code, we have solved a system of nonlinear coupled equations using a fully implicit scheme. Now, I'm going to watch the lunar eclipse! Here, you can find the full Matlab/Octave code for the above problem.
Exciting news¶
Kai Schüller has solved a phase change problem using the FVTool. You can find his implemetation here.