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:

In [8]:
% Go to the FVTool folder and initialize the toolbox
cd('/home/ali/MyPackages/FVTool/')
FVToolStartUp()
Out[8]:
AGMG 3.x linear solver is NOT available.
warning: PVTtoolbox could not be found; Please run PVTinitialize.m manually
FiniteVolumeToolbox has started successfully.

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.

In [9]:
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');
Out[9]:

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.

In [10]:
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')
Out[10]:

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.

In [11]:
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;
Out[11]:

Define the nonlinear solver specifications

We don't use a nonlinear solver. We write one. The accuracies are defined below:

In [12]:
% 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
Out[12]:

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.

In [13]:
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;
Out[13]:

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.

Convection-diffusion in single-phase flow in porous media

Fingering and channeling

For single-phase flow in porous media, the mobility of a the fluid is defined by $k/\mu$, i.e., permeability (which is a property of the porous medium) over viscosity (a property of the fluid). Fingering and channeling terms refer to a phenomenon, where a fluid moves faster than the fluid in its surrounding zones due to the higher mobility in that particular zone of the porous medium. If this higher mobility is due to a higher permeability of the zone, the phenomenon is called channeling, and if the higher mobility is due to a lower viscosity, it is called fingering. Unlike channeling, viscous fingering can happen in an ideal homogeneous field. To initialize the fingers in numerical simulations, we often perturb slightly the permeability field. The other day, my supervisor ask me to investigate if this perturbation can indeed cause channeling in my simulations. One way to make a distinction between these two phenomenon is a tracer test, in which we inject a single phase fluid with a known concentration of a tracer into a heterogeneous field, saturated with the same fluid with zero concentration of the tracer. The concentration of the tracer does not change the properties of the fluid. We gradually increase the heterogeneity and observe its effect on the propagation of the concentration profile. We want to be certain that the small perturbation that we use to initialize the fingers does not lead to channeling phenomenon.

Single-phase flow in porous media

Flow of a Newtonian fluid in porous media is governed by the Darcy's law, which is basically an empirical relation between the permeability of the porous medium ($k$), the viscosity of the fluid ($\mu$), the superficial velocity of the fluid ($u$), and the pressure gradient ($\nabla p$): $$\mathbf{u}=-\frac{k}{\mu}\nabla p$$ For the flow of an incompressible single-phase fluid, the continuity equation reads $$\nabla.\mathbf{u}=\nabla.(-\frac{k}{\mu}\nabla p)=0$$ Finally, the convection-diffusion equation in porous media is defined by $$\varphi \frac{\partial c}{\partial t}+\nabla.(\mathbf{u}c-\varphi D \nabla c) =0,$$ where $\varphi$ is the porosity, $c$ is the concentration of the tracer or solute in mole per unit volume of the fluid, and $D$ is the diffusion coefficient. The initial condition is zero concentration, and the boundaries are Dirichlet on the left border, and Neumann on the right border.
single phase domain

Let me implement it here.

In [12]:
using JFVM, PyPlot;

Find the pressure profile

First, we solve the diffusion equation to find the pressure profile. Then, we will use the pressure gradient to find the velocity vector. The velocity vector will be part of the convection-diffusion equation.

In [9]:
Nx=50
Ny=50
Lx=1.0 # [m]
Ly=1.0 # [m]
# physical values
D_val=1.0e-9 # [m^2/s]
mu_val=1e-3 # [Pa.s]
poros=0.2
perm_val=1.0e-12 # [m^2]
clx=0.05
cly=0.05
V_dp=0.02
perm= permfieldlogrndg(Nx,Ny,perm_val,V_dp,clx,cly)
# physical system
u_inj=1.0/(3600*24) # [m/s]
c_init=0.0
c_inj=1.0
p_out=100e5 # [Pa]
# create mesh and assign values to the domain
m= createMesh2D(Nx, Ny, Lx, Ly)
k=createCellVariable(m, perm)
phi=createCellVariable(m,poros)
D=createCellVariable(m, D_val)
# Define the boundaries
BCp = createBC(m) # Neumann BC for pressure
BCc = createBC(m) # Neumann BC for concentration
# change the right boandary to constant pressure (Dirichlet)
BCp.right.a[:]=0.0
BCp.right.b[:]=1.0
BCp.right.c[:]=p_out
# left boundary
BCp.left.a[:]=-perm_val/mu_val
BCp.left.c[:]=u_inj
# change the left boundary to constant concentration (Dirichlet)
BCc.left.a[:]=0.0
BCc.left.b[:]=1.0
BCc.left.c[:]=1.0
labda_face=harmonicMean(k/mu_val)
Mdiffp=diffusionTerm(labda_face)
(Mbcp, RHSp) = boundaryConditionTerm(BCp)
Mp= Mdiffp+Mbcp
p_val=solveLinearPDE(m, Mp, RHSp)
figure(figsize=(5,4))
visualizeCells(p_val)
title("Pressure profile [Pa]")
colorbar();

Solve the convection-diffusion equation

Now we have the velocity profile, that we will use to define and solve the convection-diffusion equation.

In [10]:
# estimate a practical time step
n_loop=50
dt=Lx/u_inj/(5*100) # [s]
# find the velocity vector
u=-labda_face.*gradientTerm(p_val) # [m/s]
# find the matrices of coefficients
D_face=harmonicMean(phi.*D)
Mdiffc=diffusionTerm(D_face)
Mconvc=convectionUpwindTerm(u)
(Mbcc, RHSbcc) = boundaryConditionTerm(BCc)
# initialize
c_old = createCellVariable(m, c_init, BCc)
c_val = copyCell(c_old)
# start the loop
for i=1:n_loop
    (Mtransc, RHStransc) = transientTerm(c_old, dt, phi)
    Mc=-Mdiffc+Mconvc+Mbcc+Mtransc
    RHSc=RHSbcc+RHStransc
    c_val=solveLinearPDE(m, Mc, RHSc)
    c_old=copyCell(c_val)
end
In [11]:
figure(figsize=(11,4))
subplot(1,2,1)
visualizeCells(c_old)
title("concentration profile")
colorbar();
subplot(1,2,2)
visualizeCells(k)
colorbar()
title("perm field");

Conclusion

By changing the Dykstra-Parsons coefficient to a higher value, we can observe channeling phenomenon in porous media by looking at the concentration profile. A small value of Dykstra-Parsons, around 0.01, which makes a +/- 3% perturbation in the permeability field does not lead to a channeling phenomenon, and is enough to initialize the fingers. Beautiful figures can be made by assigning a value of 0.7 or 0.8 to the Dyktra-Parsons coefficient (see below)!
tracer profile

Update (May 24, 2015)

Find the matlab code for the above problem here

Choosing the right cell averaging method for linearized PDEs

Choosing the right averaging method

In the previous post, I discussed a linearization procedure that could be used to solve a nonlinear diffusion equation. We saw that the linearization converted the nonlinear diffusion equation to a PDE with linear diffusion and convection terms. Here, I'm going to solve that equation again, but this time on a nonuniform mesh. My objective is to see what happens if we choose an averaging term that is not consistent with the chosen method for the discretization of the convection term.

Another look at the equations

The equation that we solved was $$\frac{\partial\phi}{\partial t}+\nabla . (-D(\phi)\nabla \phi) =0$$ where we can linearize the diffusion term as $$D\left(\phi\right)\nabla \phi=D\left(\phi_{0}\right)\nabla \phi+\nabla \phi_{0}\left(\frac{\partial D}{\partial \phi}\right)_{\phi_{0}}\left(\phi-\phi_{0}\right).$$ Let's put it back in the original diffusion equation: $$\frac{\partial\phi}{\partial t}+\nabla . \left(-D\left(\phi_{0}\right)\nabla \phi-\nabla \phi_{0}\left(\frac{\partial D}{\partial \phi}\right)_{\phi_{0}}\left(\phi-\phi_{0}\right)\right) =0$$ Of course this is an approximation. The idea is that we gues a $\phi_0$, then solve the equation for $\phi$, and replace it in $\phi_0$ until the absolute or relative difference between them is less than a chosen small value. I write the above equation in a form that is actually used in a finite volume scheme with all the averaging on the cell faces. It reads $$\frac{\partial\phi}{\partial t}+\nabla . \left(-D\left(\phi_{0}\right)\nabla \phi\right)+\nabla.\left(-\nabla \phi_{0}\left(\frac{\partial D}{\partial \phi}\right)_{\phi_{0,face}}\phi\right) =\nabla.\left(-\nabla \phi_{0}\left(\frac{\partial D}{\partial \phi}\right)_{\phi_{0,face}}\phi_{0,face}\right)$$ Here's the trick: on the left hand side, we have, from left to right, a transient term, a diffusion term, a convection term, and on the right hand side, a divergence term. We can use any type of discretization for the convective term, as long as the same averaging term to calculate the face value on the right hand side, that would be:

Discretization Averaging
convectionTerm linearMean
convectionUpwindTerm upwindMean
convectionTvdTerm tvdMean

In this example, I'm going to show what will happen if you use an averaging technique that is not consistent with your chosen discretization method for the convective term.

Implementation

I'm repeating the nonlinear diffusion model of the previous post, this time on a nonuniform grid. Two functions are created: the direct substitution, which is assumed to give the right solution, and a linearized equation, with the possibility of using central difference or upwind convection discretization and linear, upwind, and arithmetic averaging methods.

Definition of functions in JFVM

In [1]:
using JFVM, PyPlot
INFO: Loading help data...

Linearized PDE (Newton method)

In [2]:
function diffusion_newton(conv_method, avg_method)
    L= 1.0 # domain length
    Nx= 100
    dx_min=L/Nx
        x=[0.0, dx_min]
    while x[end]<L
        push!(x, x[end]+1.05(x[end]-x[end-1]))
    end
    x[end]=L
    m= createMesh1D(x) # create a nonuniform 1D mesh
    D0= 1.0 # diffusion coefficient constant
    # Define the diffusion coefficientand its derivative
    D(phi)=D0*(1.0+phi.^2)
    dD(phi)=2.0*D0*phi
    # create boundary condition
    BC = createBC(m)
    BC.left.a[:]=0.0
    BC.left.b[:]=1.0
    BC.left.c[:]=5.0
    BC.right.a[:]=0.0
    BC.right.b[:]=1.0
    BC.right.c[:]=0.0
    (Mbc, RHSbc)= boundaryConditionTerm(BC)
    # define the initial condition
    phi0= 0.0 # initial value
    phi_old= createCellVariable(m, phi0, BC) # create initial cells
    phi_val= copyCell(phi_old)
    # initialize the solution for upwind case
    phi_face= linearMean(phi_val)
    u= -faceEval(dD, phi_face).*gradientTerm(phi_val)
    # define the time steps
    dt= 0.01*L*L/D0 # a proper time step for diffusion process
    for i=1:10
        err=100
        (Mt, RHSt) = transientTerm(phi_old, dt)
        while err>1e-10
            # calculate the diffusion coefficient
            Dface= harmonicMean(cellEval(D,phi_val))
            # calculate the face value of phi_0
            if avg_method==1
                phi_face= linearMean(phi_val)
            elseif avg_method==2
                phi_face= arithmeticMean(phi_val)
            elseif avg_method==3
                phi_face= upwindMean(phi_val, u)
            end
            # calculate the velocity for convection term
            u= -faceEval(dD, phi_face).*gradientTerm(phi_val)
            # diffusion term
            Mdif= diffusionTerm(Dface)
            # convection term
            if conv_method==1
                Mconv= convectionTerm(u)
            elseif conv_method==2
                Mconv= convectionUpwindTerm(u)
            end
            # divergence term on the RHS
            RHSlin= divergenceTerm(u.*phi_face)
            # matrix of coefficients
            M= Mt-Mdif+Mconv+Mbc
            # RHS vector
            RHS= RHSbc+RHSt+RHSlin
            # call the linear solver
            phi_new= solveLinearPDE(m, M, RHS)
            # calculate the error
            err= maximum(abs(phi_new.value-phi_val.value))
            # assign the new phi to the old phi
            phi_val=copyCell(phi_new)
        end
        phi_old= copyCell(phi_val)
        #visualizeCells(phi_val)
    end
    return phi_val
end
Out[2]:
diffusion_newton (generic function with 1 method)

Direct substitution

In [3]:
function diffusion_direct()
    L= 1.0 # domain length
    Nx= 100
    dx_min=L/Nx
        x=[0.0, dx_min]
    while x[end]<L
        push!(x, x[end]+1.05(x[end]-x[end-1]))
    end
    x[end]=L
    m= createMesh1D(x) # create a nonuniform 1D mesh
    D0= 1.0 # diffusion coefficient constant
    # Define the diffusion coefficientand its derivative
    D(phi)=D0*(1.0+phi.^2)
    dD(phi)=2.0*D0*phi
    # create boundary condition
    BC = createBC(m)
    BC.left.a[:]=0.0
    BC.left.b[:]=1.0
    BC.left.c[:]=5.0
    BC.right.a[:]=0.0
    BC.right.b[:]=1.0
    BC.right.c[:]=0.0
    (Mbc, RHSbc)= boundaryConditionTerm(BC)
    # define the initial condition
    phi0= 0.0 # initial value
    phi_old= createCellVariable(m, phi0, BC) # create initial cells
    phi_val= copyCell(phi_old)
    # define the time steps
    dt= 0.01*L*L/D0 # a proper time step for diffusion process
    for i=1:10
        err=100
        (Mt, RHSt) = transientTerm(phi_old, dt)
        while err>1e-10
            # calculate the diffusion coefficient
            Dface= harmonicMean(cellEval(D,phi_val))
            # diffusion term
            Mdif= diffusionTerm(Dface)
            # matrix of coefficients
            M= Mt-Mdif+Mbc
            # RHS vector
            RHS= RHSbc+RHSt
            # call the linear solver
            phi_new= solveLinearPDE(m, M, RHS)
            # calculate the error
            err= maximum(abs(phi_new.value-phi_val.value))
            # assign the new phi to the old phi
            phi_val=copyCell(phi_new)
        end
        phi_old= copyCell(phi_val)
        #visualizeCells(phi_val)
    end
    return phi_val
end
Out[3]:
diffusion_direct (generic function with 1 method)

Case 1: upwind convection with arithmetic average

In [8]:
phi_ds=diffusion_direct()
phi_n=diffusion_newton(2,2)
plot(phi_ds.domain.cellcenters.x, phi_ds.value[2:end-1], "y",
phi_n.domain.cellcenters.x, phi_n.value[2:end-1], "--r") 
legend(["direct subs", "Newton"])
title("Upwind/Arithmetic")
xlabel("x")
ylabel(L"\phi", fontsize=14);

As you can see, the converged solution of the Newton method, which uses inconsistent convection/averaging terms does not match the solution of the direct substitution solver.

Case 2: central convection with arithmetic average

In [9]:
phi_ds=diffusion_direct()
phi_n=diffusion_newton(1,2)
plot(phi_ds.domain.cellcenters.x, phi_ds.value[2:end-1], "y",
phi_n.domain.cellcenters.x, phi_n.value[2:end-1], "--r") 
legend(["direct subs", "Newton"])
title("Central/arithmetic")
xlabel("x")
ylabel(L"\phi", fontsize=14);

Again, choosing the inconsistent methods results in an erroneous solution (for the Newton method).

Case 3: central convection with linear average

In [10]:
phi_ds=diffusion_direct()
phi_n=diffusion_newton(1,1)
plot(phi_ds.domain.cellcenters.x, phi_ds.value[2:end-1], "y",
phi_n.domain.cellcenters.x, phi_n.value[2:end-1], "--r") 
legend(["direct subs", "Newton"])
title("Central/linear")
xlabel("x")
ylabel(L"\phi", fontsize=14);

Beautiful! The results match perfectly.

Case 4: upwind convection with upwind average

In [11]:
phi_ds=diffusion_direct()
phi_n=diffusion_newton(2,3)
plot(phi_ds.domain.cellcenters.x, phi_ds.value[2:end-1], "y",
phi_n.domain.cellcenters.x, phi_n.value[2:end-1], "--r") 
legend(["direct subs", "Newton"])
title("Upwind/upwind")
xlabel("x")
ylabel(L"\phi", fontsize=14);

Again, we get the correct result by choosing a right pair of discretization/averaging methods. Depending on the Peclet number (convective to diffusive rate ratio), using the upwind method gives a numerically more stable solution.

Solving nonlinear PDEs with FVM

There are many cases where the transfer coefficient in the convection-diffusion equations is a function of the independent variable. This makes the PDE nonlinear. While working on my thesis, I realized that it is not that difficult to use the Taylor expansion to linearize the nonlinear terms in a PDE. Here, I'm going to give a simple example. In another post, I will explain how the same idea can be used to numerically solve the Buckley-Leverett equation.

The PDE

I'm going to try a simple one here. I say try because I'm actually making it up right now. The equation reads $$\frac{\partial\phi}{\partial t}+\nabla . (-D\nabla \phi) =0$$ where $$D=D_0(1+\phi^2).$$ It means that the diffusion coefficient is increased when the concentration of the solute increases.
To linearize this equation, I assume that $\phi$ and $\nabla \phi$ are independent variables. Using Taylor expansion around a point $\phi_0$, we can obtain the following general form: $$f\left(\phi\right)\nabla \phi=f\left(\phi_{0}\right)\nabla \phi+\nabla \phi_{0}\left(\frac{\partial f}{\partial \phi}\right)_{\phi_{0}}\left(\phi-\phi_{0}\right).$$ In our case, this equation simplifies to $$D_0(1+\phi^2)\nabla \phi=D_0(1+\phi_0^2)\nabla \phi+2D_0\phi_0\nabla \phi_0(\phi-\phi_0).$$ Putting back everythin together, we obtain this linearized form for our PDE $$\frac{\partial\phi}{\partial t}+\nabla.(-D_{0}(1+\phi_{0}^{2})\nabla\phi)+\nabla.(-2D_{0}\phi_{0}\nabla\phi_{0}\phi)=\nabla.(-2D_{0}\phi_{0}\nabla\phi_{0}\phi_{0}).$$ Don't freak out! This equation is actually quite simple. By linearizing, we have added a linear convection term to our nonlinear diffusion equation. This equation is still an approximation of the real PDE. We have to solve the linear equation for $\phi$ by initializing $\phi_0$. Then, we assign the new value of $\phi$ to $\phi_0$ until it converges to a solution. The question that still remains is that "is this really a good idea to make a diffusion equation more complicated by adding a convection term to it?".

The solution procedure

Solving this equation is done in two loops. The outer loop is for the time steps, and the inner loop is for the nonlinear equation. We will assemple an equation which looks like this: $$\text{Transient}(\Delta t,\phi)-\text{Diffusion}(D(\phi_0),\phi)+\text{Convection}(\mathbf{u}(\phi_0),\phi)=\text{divergence}(D(\phi_0), \nabla\phi_0)$$

Implementation in JFVM

Here, I'll implement the procedure in JFVM. The same procedure can be done in FVtool.

In [1]:
using JFVM, PyPlot
INFO: Loading help data...
In [16]:
function diffusion_newton()
L= 1.0 # domain length
Nx= 100 # number of cells
m= createMesh1D(Nx, L) # create a 1D mesh
D0= 1.0 # diffusion coefficient constant
# Define the diffusion coefficientand its derivative
D(phi)=D0*(1.0+phi.^2)
dD(phi)=2.0*D0*phi
# create boundary condition
BC = createBC(m)
BC.left.a[:]=0.0
BC.left.b[:]=1.0
BC.left.c[:]=5.0
BC.right.a[:]=0.0
BC.right.b[:]=1.0
BC.right.c[:]=0.0
(Mbc, RHSbc)= boundaryConditionTerm(BC)
# define the initial condition
phi0= 0.0 # initial value
phi_old= createCellVariable(m, phi0, BC) # create initial cells
phi_val= copyCell(phi_old)
# define the time steps
dt= 0.001*L*L/D0 # a proper time step for diffusion process
for i=1:10
    err=100
    (Mt, RHSt) = transientTerm(phi_old, dt)
    while err>1e-10
        # calculate the diffusion coefficient
        Dface= harmonicMean(cellEval(D,phi_val))
        # calculate the face value of phi_0
        phi_face= linearMean(phi_val)
        # calculate the velocity for convection term
        u= faceEval(dD, phi_face).*gradientTerm(phi_val)
        # diffusion term
        Mdif= diffusionTerm(Dface)
        # convection term
        Mconv= convectionTerm(u)
        # divergence term on the RHS
        RHSlin= divergenceTerm(u.*phi_face)
        # matrix of coefficients
        M= Mt-Mdif-Mconv+Mbc
        # RHS vector
        RHS= RHSbc+RHSt-RHSlin
        # call the linear solver
        phi_new= solveLinearPDE(m, M, RHS)
        # calculate the error
        err= maximum(abs(phi_new.value-phi_val.value))
        # assign the new phi to the old phi
        phi_val=copyCell(phi_new)
    end
    phi_old= copyCell(phi_val)
    visualizeCells(phi_val)
end
    return phi_val
end
Out[16]:
diffusion_newton (generic function with 1 method)
In [33]:
@time phi_n=diffusion_newton()
elapsed time: 0.065438686 seconds (15208324 bytes allocated)
Out[33]:
CellValue{Float64}(MeshStructure(1,[100],CellSize{Float64}([0.01,0.01,0.01,0.01,0.01,0.01,0.01,0.01,0.01,0.01  …  0.01,0.01,0.01,0.01,0.01,0.01,0.01,0.01,0.01,0.01],[0.0],[0.0]),CellLocation{Float64}([0.005,0.015,0.025,0.035,0.045,0.055,0.065,0.075,0.085,0.095  …  0.905,0.915,0.925,0.935,0.945,0.955,0.965,0.975,0.985,0.995],[0.0],[0.0]),FaceLocation{Float64}([0.0,0.01,0.02,0.03,0.04,0.05,0.06,0.07,0.08,0.09  …  0.91,0.92,0.93,0.94,0.95,0.96,0.97,0.98,0.99,1.0],[0.0],[0.0]),[1],[1]),[5.0196,4.9804,4.94059,4.90018,4.85914,4.81745,4.77511,4.7321,4.68839,4.64398  …  0.000323843,0.000250452,0.000192389,0.000146146,0.000108911,7.84098e-5,5.27679e-5,3.04028e-5,9.92833e-6,-9.92833e-6])

The solution looks quite good. The diffusion coefficient is higher at higher concentrations, and it can be clearly observed in the above concentration profile (if you compare it with a diffusion problem with constant diffusion coefficient, in your head).

Some questions

The convection term that appears after linearization of the problem has a velocity term, which reads $$\mathbf{u}(\phi_0)=−2.0D_0\phi_0\nabla\phi_0$$ I have written a function, called `gradientTerm` which gives you the value of the gradient on the cell faces. This gradient is multiplied by a term that needs to be averaged on each cell face. Here, I have calculated the arithmetic average of $\phi_0$ using a function called suprise suprise linearMean. This is of course one way of doing it. The other ways are geometricMean, harmonicMean, upwindMean, tvdMean, and aritheticMean. Here, we will discuss them.

Averaging techniques for the face values in FVM

The following figure shows two neigbor cells with different length and values.
neighbor cells FVM The averaging terms that are used in JFVM are as follows:

Linear mean

It is equivalent to a linear interpolation or the lever rule. It reads $$\phi_{face}=\frac{\Delta x_1\phi_2+\Delta x_2\phi_1}{\Delta x_1+\Delta x_2}$$

Arithmetic mean

The value of each cell is weighted by its length: $$\phi_{face}=\frac{\Delta x_1\phi_1+\Delta x_2\phi_2}{\Delta x_1+\Delta x_2}$$

Geometric mean

It reads $$\phi_{face}=\exp\left(\frac{\Delta x_1\log(\phi_1)+\Delta x_2\log(\phi_2)}{\Delta x_1+\Delta x_2}\right)$$

Harmonic mean

It reads $$\phi_{face}=\frac{\Delta x_1+\Delta x_2}{\Delta x_1/\phi_1+\Delta x_2/\phi_2}$$ Note that if the value of one of the cells is zero, the geometric and harmonic means are also equal to zero.
These methods have different applications for different ptoblems. We will discuss them in future posts.

Direct substitution

The diffusion problem that we solved above can also be solved using a direct substitution method. Here, we don't use the derivatives and no linearization is required. The code is written here:

In [18]:
function diffusion_direct()
L= 1.0 # domain length
Nx= 100 # number of cells
m= createMesh1D(Nx, L) # create a 1D mesh
D0= 1.0 # diffusion coefficient constant
# Define the diffusion coefficientand its derivative
D(phi)=D0*(1.0+phi.^2)
dD(phi)=2.0*D0*phi
# create boundary condition
BC = createBC(m)
BC.left.a[:]=0.0
BC.left.b[:]=1.0
BC.left.c[:]=5.0
BC.right.a[:]=0.0
BC.right.b[:]=1.0
BC.right.c[:]=0.0
(Mbc, RHSbc)= boundaryConditionTerm(BC)
# define the initial condition
phi0= 0.0 # initial value
phi_old= createCellVariable(m, phi0, BC) # create initial cells
phi_val= copyCell(phi_old)
# define the time steps
dt= 0.001*L*L/D0 # a proper time step for diffusion process
for i=1:10
    err=100
    (Mt, RHSt) = transientTerm(phi_old, dt)
    while err>1e-10
        # calculate the diffusion coefficient
        Dface= harmonicMean(cellEval(D,phi_val))
        # diffusion term
        Mdif= diffusionTerm(Dface)
        # matrix of coefficients
        M= Mt-Mdif+Mbc
        # RHS vector
        RHS= RHSbc+RHSt
        # call the linear solver
        phi_new= solveLinearPDE(m, M, RHS)
        # calculate the error
        err= maximum(abs(phi_new.value-phi_val.value))
        # assign the new phi to the old phi
        phi_val=copyCell(phi_new)
    end
    phi_old= copyCell(phi_val)
    visualizeCells(phi_val)
end
    return phi_val
end
Out[18]:
diffusion_direct (generic function with 1 method)
In [34]:
@time phi_ds=diffusion_direct()
elapsed time: 0.115909663 seconds (22458788 bytes allocated, 26.14% gc time)
Out[34]:
CellValue{Float64}(MeshStructure(1,[100],CellSize{Float64}([0.01,0.01,0.01,0.01,0.01,0.01,0.01,0.01,0.01,0.01  …  0.01,0.01,0.01,0.01,0.01,0.01,0.01,0.01,0.01,0.01],[0.0],[0.0]),CellLocation{Float64}([0.005,0.015,0.025,0.035,0.045,0.055,0.065,0.075,0.085,0.095  …  0.905,0.915,0.925,0.935,0.945,0.955,0.965,0.975,0.985,0.995],[0.0],[0.0]),FaceLocation{Float64}([0.0,0.01,0.02,0.03,0.04,0.05,0.06,0.07,0.08,0.09  …  0.91,0.92,0.93,0.94,0.95,0.96,0.97,0.98,0.99,1.0],[0.0],[0.0]),[1],[1]),[5.0196,4.9804,4.94059,4.90018,4.85914,4.81745,4.77511,4.7321,4.68839,4.64398  …  0.000323843,0.000250452,0.000192389,0.000146146,0.000108911,7.84098e-5,5.27679e-5,3.04028e-5,9.92833e-6,-9.92833e-6])

Show and compare the results

Here, I'm plotting the results of the Newton and the direct substitution methods for comparison.

In [40]:
plot(phi_ds.domain.cellcenters.x, phi_ds.value[2:end-1], "y",
phi_n.domain.cellcenters.x, phi_n.value[2:end-1], "--r") 
legend(["direct subs", "Newton"])
Out[40]:
PyObject <matplotlib.legend.Legend object at 0x7ffa453f94d0>

Conclusion

The final results shows that we can linearize a nonlinear PDE, relatively easily and solve the equations using a linear PDE solver. No difference is observed between the final results of the direct substitution method (with no linearization) and the Newton method (with a linearized PDE). We can also observe that, as expected, the Newton method is relatively faster.

Update

The same problem implemented in Matlab:

2D transient diffusion equation; numerical FVM solution

Important Update: the codes in this post will not work with the new version of FVTool. Download the old version of FVTool here.

A time dependent diffusion problem

In my last post, I promised to solve a time dependent conduction problem. Indeed I tried to keep my promise by writing a full post about single phase compressible flow in porous media. But in the middle of writing, I noticed that there are too many things that need to be clarified before we can jump into that problem. Therefore, I decided to hold my horses and only solve a simple time-dependent diffusion problem. I this post, you are going to learn how to define the initial conditions and use a for loop for time steps.

This blog's icon

If you look carefully at this open tab in your browser, you see a colorful icon that looks like a bad design for a gay flag. It is nothing but a square domain, which has been initially at a concentration c=1.0, and suddenly its boundaries are exposed to an environment at a zero concentration, c=0.0 and mass starts moving out of the domain only by diffusion mechanism. As simple as that!

We gaan beginnen

The equation for this problem reads

$$\frac{\partial c}{\partial t} +\nabla.(-D \nabla c) = 0$$

where D [m^2/s] is the diffusion coefficient and c [mol/m^3] is the concentration. The boundary conditions are all Dirichlet, i.e.,

$$ c=0 $$

The coding steps are as always in the following sequence:

  1. Geometry and mesh
  2. Boundary condition
  3. Initial condition
  4. Matrix of coefficients
  5. Linear solver

The geometry can be defined as

clc; clear; close all;
L = 0.1; % [m] length of the domain
H = 0.1; % [m] height of the domain
Nx = 50; % number of grids in x direction
Ny = 50; % number of grids in y direction
m = createMesh2D(Nx, Ny, L, H); % create the mesh

Then the boundary conditions are defined. Just as a reminder, the boundaries are defined by the following general relation:

$$ a \nabla \phi.\mathbf{n}+b \phi = c$$

The code for specifying the boundary condition and finding the matrix of coefficients for the boundary nodes is written as

bc = createBC(m);
bc.left.a(:)=0; bc.left.b(:)=1; bc.left.c(:)=0;
bc.right.a(:)=0; bc.right.b(:)=1; bc.right.c(:)=0;
bc.top.a(:)=0; bc.top.b(:)=1; bc.top.c(:)=0;
bc.bottom.a(:)=0; bc.bottom.b(:)=1; bc.bottom.c(:)=0;
[M_bc, RHS_bc] = boundaryCondition(m, bc);

Now is time to learn something new: defining initial conditions. We remember that the variables are build over the domain using the createCellVariable function. We used it before to assign a conductivity value to each cell in the domain. We don't need to assign transfer coefficients to the ghost cells. Now, if we need to define a cell variable which also has value on the ghost cell. The value on the ghost cells should be consistent with the boundary condition. Therefore, we send the boundary condition as an extra input the the createCellVariable function. The function will take care of the rest! One more thing. The initial condition must be assigned to a structure that is called Old. In fact, the transitionTerm that returns the matrix of coefficients for the transient term needs the initial conditon to be written this way.

D_val = 1e-5; % [m^2/s] diffusion coefficient
D = createCellVariable(m, D_val); % assign diff. coef. to each cell
D_face = geometricMean(m, D); % average of diff. coef. on cell faces
M_diff = diffusionTerm(m, D_face);
c_init = 1.0; % [mol/m^3] initial concentration
c.Old = createCellVariable(m, c_init, bc); % initial condition
alfa = createCellVariable(m, 1.0); % it will be required later

The next step is to define the time steps and the final time of the simulation. Normally, for a diffusion problem, we can choose a time step based on the length of the domain and the diffusion coefficient, i.e., a fraction of $L^2/D$. In Matlab, we can write

t_end = 0.2*L^2/D_val; % [s] final time
dt = t_end/100; % [s] time step 

And the final part. We have to solve the PDE literally step by step, in a loop. We start from the initial condition (t=0) and find the concentration profile for t=dt. Then we use the new concentration profile as the initial condition and solve the PDE for the next time step and so on. The Matlab code for this procedure reads

for t = 0:dt:t_end
        % This part must be inside the time loop:
    [M_trans, RHS_trans] = transientTerm(m, alfa, dt, c);
    c_new = solvePDE(m, M_bc+M_trans-M_diff, RHS_trans+RHS_bc);
    c.Old = c_new; % replace the old value with the new time step
    visualizeCells(m, c_new); title(['t= ' num2str(t) ' s']);
        drawnow;
end

The code will show you how the material escapes the square domain through the boundaries. We can do a bit of mass balance test on the whole thing, which I will explain later.
This is the code in one piece:

clc; clear; close all;
L = 0.1; % [m] length of the domain
H = 0.1; % [m] height of the domain
Nx = 50; % number of grids in x direction
Ny = 50; % number of grids in y direction
m = createMesh2D(Nx, Ny, L, H); % create the mesh
bc = createBC(m);
bc.left.a(:)=0; bc.left.b(:)=1; bc.left.c(:)=0;
bc.right.a(:)=0; bc.right.b(:)=1; bc.right.c(:)=0;
bc.top.a(:)=0; bc.top.b(:)=1; bc.top.c(:)=0;
bc.bottom.a(:)=0; bc.bottom.b(:)=1; bc.bottom.c(:)=0;
[M_bc, RHS_bc] = boundaryCondition(m, bc);
D_val = 1e-5; % [m^2/s] diffusion coefficient
D = createCellVariable(m, D_val); % assign diff. coef. to each cell
D_face = geometricMean(m, D); % average of diff. coef. on cell faces
M_diff = diffusionTerm(m, D_face);
c_init = 1.0; % [mol/m^3] initial concentration
c.Old = createCellVariable(m, c_init, bc); % initial condition
alfa = createCellVariable(m, 1.0); % it will be required later
t_end = 0.1*L^2/D_val; % [s] final time
dt = t_end/100; % [s] time step 
for t = 0:dt:t_end
    [M_trans, RHS_trans] = transientTerm(m, alfa, dt, c);
    c_new = solvePDE(m, M_bc+M_trans-M_diff, RHS_trans+RHS_bc);
    c.Old = c_new;
    visualizeCells(m, c_new); title(['t= ' num2str(t) ' s']);
        drawnow;
end

One of the frames in your final result should look like this: 2d diffusion fvm

Conduction and diffusion: a brief tutorial

Important Update: the codes in this post will not work with the new version of FVTool. Download the old version of FVTool here.

Warning

This post is not edited. You may find horrible English mistakes.

Diffusion/conduction equation; mass and heat transfer

As a chemical engineer mass transfer is my favorite topic, which is sort of strange given the fact that heat transfer is way easier to feel and understand. It has to be due to the very busy schedule of my heat transfer Professor, who was involved in a project at the time and did not spend enough time on getting ready for his teaching duties, unlike my mass transfer Professor. Back to business...

Temperature gradient is the driving force behind the conductive heat transfer, and the gradient of chemical potential is the driving force behind the diffusive mass transfer. At low concentrations, the gradient of the chemical potential can be replaced by the concentration gradient. Here I'm going to show you how to solve a conservation equation when the only flux term is the diffusive heat or conductive mass (doesn't sound as funny as I expected). The equation reads

$$ \nabla . (-\lambda \nabla \phi)=q, $$

where for mass transfer, $\lambda$ denotes the diffusivity (m^2/s) and $\phi$ denotes the concentration (mol/m^3), and q is a mass source term (mol/(m^3.s)) and for heat transfer, $\lambda$ denotes the conductivity (J/(m.K.s)) and $\phi$ denotes the temperature (K), and q is a heat source term (J/(m^3.s)). There are other physical phenomena that can be described by the above relation, e.g., Poisson equation or flow in porous media, in which $\lambda$ denotes the total mobility (permeability divided by the viscosity for single phase flow), and $\phi$ denotes pressure.

Boundary conditions

I think Dirichlet (constant value) and Neumann (constant flux) boundary conditions are quite clear in terms of their physical meaning. However, there are two important situations that lead to a Robin boundary condition. In heat transfer, when a boundary is gaining/losing heat to a medium with a constant temperature of $$ T_{\infty} $$ by convection mechanism with a heat transfer coefficient h (J/(m^2.K.s)), the energy balance equation at the boundary reads

$$ -\lambda (\mathbf{n}.\nabla T) = h(T-T_{\infty}), $$

which can be rearranged to the following form that can be formulated in the FVMtool:

$$\frac{\lambda}{h} (\mathbf{n}.\nabla T) + T = T_{\infty}$$

For the mass transfer, if mass is produced or consumed at a boundary by a first order reaction rate, the equation reads

$$ -\lambda (\mathbf{n}.\nabla c) = k_0 c, $$

where $k_0$ is the rate constant of the reaction that happens on the boundary. For a chapter of my thesis, I encountered this sort of boundaries and it was the main reason that I rewrote the implementation of boundary condition from scratch. I will talk about the implementation of the boundary conditions later in a separate post.

A (sort of) real problem

Fins are used to increase the heat transfer area and thus the rate of heat transfer. Here we are going to model a fin in 2D. The base of the rectangular fin is attached to a surface with a constant temperature of 100 degree Celsius. The fin is made of Aluminum, with a thermal conductivity of 237 W/(m.K), with a thickness of 0.1 cm and a length of 10 cm. The fin is exposed to a air at a constant temperature of 25 degree Celsius, and the heat transfer coefficient is 10 W/(m^2.K). This number may be off so you can estimate it using one of these correlations.
Let us jump into our Matlab FVTool, and solve this problem numerically using the finite volume method.

Solution procedure

First, we start by defining the domain and creating the mesh structure:

clc; clear; % clean the command prompt, clean the memory
L = 0.1; % 10 cm length
W = 0.01; % 1 cm thickness
Nx = 50; % number of cell in x direction
Ny = 20; % number of cells in the y direction
m = createMesh2D(Nx, Ny, L, W); % creates a 2D Cartesian grid

Now, we define the the transfer coefficients. We assume that the thermal conductivity is constant on the whole domain. However, we need to assign this constant value to each individual cell. This is done by the function createCellVariable. After creating this thermal conductivity field, we have to calculate the average values on the interface between cells or the faces. Again, in this special case, the average values are equal to a constant. But in other cases when the conductivity varies with space, the averaging becomes more important. Three averaging techniques are available in the FVTool, viz. arithmetic, geometric, and harmonic. The output of these averaging functions is always a face variable. Let's see that in action:

T_inf = 25+273.15; % [K] ambient temperature
T_base = 100+273.15; % [K] temperature at the base of the fin
k_val = 237; % W/(m.K) thermal conductivity
h_val = 10; % W/(m^2.K) heat transfer coefficient
k = createCellVariable(m, k_val); % assign thermal cond. value to all cells
k_face = geometricMean(m, k); % geometric average of the thermal conductivity values on the cell faces

Now, we can define the boundary conditions. Here, we have one Dirichlet (constant temperature) and three Robin boundaries. The general boundary condition is defined as

$$ a (\mathbf{n}.\nabla \phi)+b\phi = c. $$

a, b, and c must be defined in the program. Here, we have $a=\lambda/h$, $b=1$, and $c=T_{\infty}$. Don't forget to include the sign of the normal vector for the bottom boundary. The normal vector is in the opposite direction of the y axis and therefore its sign must be included in the temperature gradient term. I will talk about it in more details later.

BC = createBC(m); % creates a BC structure for the domain m; all Neumann boundaries
BC.left.a(:)=0; BC.left.b(:)=1; BC.left.c(:)=T_base; % convert the left boundary to constant temperature
BC.right.a(:)=k_val/h_val; BC.right.b(:)=1; BC.right.c(:)=T_inf; % right boundary to Robin
BC.top.a(:)=k_val/h_val; BC.top.b(:)=1; BC.top.c(:)=T_inf; % top boundary to Robin
BC.bottom.a(:)=-k_val/h_val; BC.bottom.b(:)=1; BC.bottom.c(:)=T_inf; % bottom boundary to Robin

Now, we have a domain with fully specified transfer coefficients and defined boundary. Next and final step is to find the matrix of coefficients for the conduction term and the boundary conditions, and solve the linear system of discretized linear equations:

M_cond = diffusionTerm(m, k_face); % matrix of coefficients for the heat diffusion
[M_bc, RHS_bc] = boundaryCondition(m, BC); % matrix of coefficients and RHS vector for the boundary conditions
T = solvePDE(m, M_cond+M_bc, RHS_bc); % solve the linear system of discretized PDE
visualizeCells(m, T); % visualize the results

Little fun

Please find the source code for this tutorial here. To have some numerical fun, try to change the problem from 2D to 3D by activating the appropriate line in the Matlab source code. You can also increase the heat transfer coefficient and see its effect on the temperature profile. In the next post, I will explain how to convert this example to a unsteady-state simulation.
This is a snapshot of the 3D result:

3d conduction in a fin

A simple finite volume toolbox for Matlab

An introduction to a Matlab FVM (toy) toolbox

For some reasons, I had to solve a few PDE's including single/multi phase flow in porous media, heat transfer in saturated porous media, multi-component mass transfer, and so on. My job never included the development of numerical method. In fact, I was supposed to come up with simple and useful models for the physical system that I was/am studying. Solving those models could be done in PDE solver of my choice. At the end of the day, I chose to write my own codes in Matlab.

Why [in general]?

Let me count my reasons here. In case you don't like to listen to me bragging about it, you can skip to the middle of this post.

Black-boxes

These numerical solvers, in particular COMSOL always was like a black-box to me. I could not understand its error messages, which was frustrating, and even when it generated nice results (which most of the times happens), I was really unable to have a physical feeling for it. I will talk about the physical feeling later, and don't expect a post about love-making.

Boundary conditions

For some reason, choosing and handling boundary conditions is not explained in books, or papers, or lectures. The same situation applies to the documentation of PDE solvers. Id it really that difficult for you guys to bring us a general Robin boundary condition?

Geometry

Honestly, one of the most attractive features of PDE solvers is the graphical pre-processor (with all sort of different mesh generation techniques) and CAD modules. For me, with the simple rectangular or cylindrical geometries of my experimental set-up, there was no need of a fancy pre-processor.

Learning/Teaching

I'll be doing a lot of teaching soon, so I needed to learn numerical methods. What's better than learning by doing? Now, I can share my experiences and my coode with my students.

Why [in particular]?

I needed a mass conservative scheme (e.g., finite volume method), which is implemented in an understandable language (yes, I know. C++ is quite beautiful and elegant and understandable even for a kid with the right genes, but I prefer Matlab), with some flexibility for specifying boundary conditions and changing the physics. I actually found a code. It's called FiPy. But I was already comfortable with Matlab and, don't tell anyone, I couldn't understand Python and NumPy. So having FiPy's syntax in mind I decided to write my own code in Matlab. I think it is in good enough shape to be shared with other FVM/Matlab users.

Really! Why?

For most of our important PDE's in chemical and petroleum engineering (did I mention that I'm a chemical/petroleum engineer?) we have analytical solutions. Also, most of the experiments that we do in the lab can be modeled with a one dimensional PDE. However, so many curious things happen when we go to a two- or three-dimensional domain. I wanted to have something, like FiPy, to make me able to solve a 1D equation, verify it by comparing it to my analytical solution or experimental data, and switch it to 2D and 3D domains without too many modifications. I have it now.

What do we solve?

We solve this general form of transient convection-diffusion equation:

$$ \alpha\frac{\partial\phi}{\partial t}+\nabla.\left(\mathbf{u}\phi\right)+\nabla.\left(-D\nabla\phi\right)+\beta\phi=\gamma$$

with the following general (Robin) boundary condition:

$$a\nabla\phi.\mathbf{n}+b\phi=c.$$

All of the coefficients can be defined explicitly for each control volume or on the surface of a control volume.

Where to find/How to use 'the code'?

You can download the code from this github repository (click on the download zip button): FVTool Or alternatively, if you are on linux (and hopefully you are), use the command

git clone https://github.com/simulkade/FVTool

How to start

Start Matlab (or Octave), go to the FVTool folder, and type

FVToolStartUp 

You must see a few messages and finally you should see FiniteVolumeToolbox has started successfully. in Matlab command prompt. With the following command, you can see a short document that introduces you to the code:

showdemo FVTdemo

If you want to jump into it, you can run the following script to solve a diffusion equation with Dirichlet boundary conditions:

Important Update: this code does not work with the new version of FVTool. Download the old version here.

clc; clear;
L = 50;  % domain length
Nx = 20; % number of cells
m = createMesh3D(Nx,Nx,Nx, L,L,L);
BC = createBC(m); % all Neumann boundary condition structure
BC.left.a(:) = 0; BC.left.b(:)=1; BC.left.c(:)=1; % Dirichlet for the left boundary
BC.right.a(:) = 0; BC.right.b(:)=1; BC.right.c(:)=0; % right boundary
D_val = 1; % value of the diffusion coefficient
D = createCellVariable(m, D_val); % assign the diffusion coefficient to the cells
D_face = harmonicMean(m, D); % calculate harmonic average of the diffusion coef on the cell faces
Mdiff = diffusionTerm(m, D_face); % matrix of coefficients for the diffusion term
[Mbc, RHSbc] = boundaryCondition(m, BC); % matix of coefficients and RHS vector for the BC
M = Mdiff + Mbc; % matrix of cefficients for the PDE
c = solvePDE(m,M, RHSbc); % send M and RHS to the solver
visualizeCells(m, c); % visualize the results

You can find more examples here.