# 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:

# 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:

# 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{e}+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.