version 2.1.0.0 (1.95 MB) by
Ehsan

A simple yet general purpose FVM solver for transient convection diffusion PDE

A simple Finite volume tool

This code is the result of the efforts of a chemical/petroleum engineer to develop a simple tool to solve the general form of convection-diffusion equation:

α∂ϕ/∂t+∇.(uϕ)+∇.(−D∇ϕ)+βϕ=γ

on simple uniform/nonuniform mesh over 1D, 1D axisymmetric (radial), 2D, 2D axisymmetric (cylindrical), and 3D domains.

The code accepts Dirichlet, Neumann, and Robin boundary conditions (which can be achieved by changing a, b, and c in the following equation) on a whole or part of a boundary:

a∇ϕ.n+bϕ=c.

It also accepts periodic boundary conditions.

The main purpose of this code is to serve as a handy tool for those who try to play with mathematical models, solve the model numerically in 1D, compare it to analytical solutions, and extend their numerical code to 2D and 3D with the minimum number of modifications in the 1D code.

The discretizaion schemes include

* central difference

* upwind scheme for convective terms

* TVD schemes for convective terms with many flux limiters

To get started, go to the `Test` folder and run the test scripts.

A few calculus functions (divergence, gradient, etc) and averaging techniques (arithmetic average, harmonic average, etc) are available, which can be helpful specially for solving nonlinear or coupled equations or implementing explicit schemes.

I have used the code to solve coupled nonlinear systems of PDE. You can find some of them in the Examples/advanced folder.

There are a few functions in the 'PhysicalProperties' folder for the calculation of the physical properties of fluids. Some of them are not mine, which is specified inside the file.

I'll try to update the documents regularly, in the github repository. Please give me your feedback/questions by writing a comment in my weblog: <http://fvt.simulkade.com/>

Special thanks: I vastly benefited from the ideas behind Fipy <http://www.ctcms.nist.gov/fipy/>, a python-based finite volume solver.

To start the solver, download and extract the zip archive, open and run 'FVToolStartUp' function.

To see the code in action, copy and paste the following in your Matlab command window:

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; % Dirichlet for the 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(D); % calculate harmonic average of the diffusion coef on the cell faces

Mdiff = diffusionTerm(D_face); % matrix of coefficients for the diffusion term

[Mbc, RHSbc] = boundaryCondition(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(c); % visualize the results

You can find some animated results of this code in my youtube channel:

https://www.youtube.com/user/processsimulation/videos

Eftekhari, A.A. et al. (2015). FVTool: a finite volume toolbox for Matlab. Zenodo. http://doi.org/10.5281/zenodo.32745

Created with
R2014a

Compatible with any release

**Inspired by:**
IAPWS_IF97

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!Create scripts with code, output, and formatted text in a single executable document.

Meteb MejbelThanks @Ehsan for your responses. I've noticed something in 1-d spherical diffusion. Whenever I decrease diffusion coefficient, I get higher concentration at the boundary. For example, using Dirichlet boundary condition (a=0, b=1, c=2), as I decrease the diffusion coefficient from 0.5 to 0.005, the concentration at the first cell increases from 2.2 to 3.8. Do you think it is related to the ghost cell implementation or somewhere else?

Ehsan@Meteb yes it is possible. just define your initial condition as a cell by providing your array of values to the createCellVariable function. Remember to pass an array that is the same size as your mesh.

Meteb MejbelHi Ehsan, if I have several initial concentrations along x-domain, is it possible to apply a vector of initial concentrations and solve for the diffusion at the next time step?

Ehsan@Subhash Thanks for the nice words. Unfortunately, odd shape domains are not supported; only the simplest shapes in each coordinate :-)

subhash inavoluThank you @Ehsan for your wonderful contribution. I have been looking for something like this for a very longtime. Hats off, I have seen some of the examples I cannot stop appreciating your work I have one question.I would like to join shapes. Lets say I want to create an 'L' Shape. is it possible to do that using these codes? If so can you provide an example? I couldnt find any example like that.

EhsanExactly, @Meteb

Meteb MejbelThank you @Ehsan for your response. Removing transient terms, does that mean at line 37 and 38, the terms M and RHS should be equal to

-Mdiff+Mbc and RHSbc, respectively?

Ehsan@Meteb. Thanks for the comment. I'm not sure if I understand the part about conserving initial concentration on the surface. However, you can always solve the steady state equation. Simply remove the transient term (and the for loop) from the transient tutorial. You can get the concentration gradient using the function gradientTerm(c) where c is a cellvalue object that contains the concentration profile.

Meteb Mejbelsuper helpful tool. Thank you so much for that. I have a question @Ehsan, I am working on diffusion equation in 1d sphere. Is there anyway to simulate the concentration gradient within the sphere domain ( r or x) without dealing with transient? I just want to observe the distribution within the domain if I have conserved initial concentration on the surface. I am using diffusiontutorial_spherical.m

Andrea Bottacin BusolinEhsan@cmaurer, it is amazing that it works with CasADI. Could you open an issue on github about fixing combineBC functions? Also, could you add your CasADI example as a pull request in the FVTool repository?

cmaurerGreat Tool, I could also use it with Casadi (including small changes)

Therefore I use the combineBC function (similar to the use of Matlabs ode solver). I wondered if you could include the missing combineBC3D Function?

EhsanThanks, Burak. I just fixed it. BTW, you can always suggest these fixes on github, or even better change it yourself with a pull request. This way, you contribution will be registered under your name and you get the credit :-)

Burak DurkutThere is an issue in the "FVTool_functions_uniform_test.m" which is in the tests folder. The 78th line the "for" has not "end".

Ali MahmoodiLei Li@Ehsan, I want to create mesh for polyshape, could this be done in this tool? secondly, I would like to use symmetry BC and convection BC, is there any example for this kind of simulation? In the tutorial, it gives simple geometry, like a rectangle, but when the geometry becomes complex and have two domains and materials, how to deal with that in this tool box?

洪川 幸@Ehsan hello I want to simulate the acid wormholes in fractured media, how should I add cracks? Many literatures use the Monte Carlo method, but I don't know how to implement it in MATLAB

Ehsan@Ali, unfortunately, FVTool does not support triangular mesh.

Ali Hammouche@Ehsan Hello I want to solve a Poisson equation by FVM (-div(e*grad(A))=curl(M).M has tow components (Mx My).

By using triangulare mesh i dont khow how to dicretize term source by Fvm with tri mesh ?? I used (pdemesh,)

Mohammad Rahmani@Ehsan

Many thanks for sharing this! I believe this code has a lot of potential specially for solving reacting flow like chemical reactors.

EhsanI have solved several problems with nonlinear source term for my work (reactive transport in porous media) that I have not shared yet. You can find one example here: https://github.com/behzaadh/FVTool/tree/master/Examples/External/InjectionHClCoreFloodProblem

This is a relatively advanced case although I don't think the source term is linearized. I'll add a simple example later.

Mohammad RahmaniIs there any example to demonstrate such solution?

Sometimes this is not quite easy! The reason I asked for f(phi) instead of phi to let solver can solve nonlinear source term.

for example in chemical reaction, most of the time you have very nonlinear terms.

Ehsan@Mohammad Rahmani: you can do that easily by writing the Taylor expansion of your f(\phi) function :-)

Mohammad RahmaniGood job. if instead of \beta*\phi you can use f(\phi) then your code can be used for reaction-diffusion-convection cases!

Keep going on!

Good luck

mali MirBehzad HosseinzadehTamour Zubair@All

anyone please help me i am quite new in matlab. In very start, i am facing the following error in 1D meshing. Thanks in advance

Error:

Undefined function or variable 'MeshStructure'.

Error in createMesh1D (line 42)

MS=MeshStructure(1, [Nx,1], CellSize, CellLocation,

FaceLocation, [1], [1]);

Fernando FernandesHello! I'm really new in the simulation world. Your scripts are fantastic! I need to simulate a salt rock flow towards a petroleum well, How could I use your scripts for this?? I don't know to open the files, I think I need open the file "meshgeneration" together to the main file. Could u send a step by step???

Thank U

Amitava BiswasEhsan@Vahid,

Currently, internal boundaries are not possible. The Jacobian matrix is made manually in FVTool. It is not too hard to make them. I have one example here for a nonlinear diffusion coefficient: http://fvt.simulkade.com/posts/2015-04-06-solving-nonlinear-pdes-with-fvm.html

If you open issues on github page, I see them and react faster: https://github.com/simulkade/FVTool/issues

vahid mossHi. Is there any chance to apply your code for a geometry with internal boundaries. For instance, assume a large scale rectangular geometry and take a small rectangular out. See also, https://www.mathworks.com/help/matlab/ref/polyshape.regions.html . Another example can be found in : Fig. 4.3 (the page 89) of the book " The Finite Volume Method in Computational Fluid Dynamics An Advanced Introduction with OpenFOAM and Matlab".

vahid mossHi, Alireza, I am new with FVM but I have found you code wonderful. Here is my first question. I assume, a set of algebraic equation in FVM should be solved similar to FDM. In FDM, a nonlinear source term, e.g. nonlinear source in Poisson's equation, will produce nonlinear algebraic equations which must be solved, e.g. by Newton method and building a Jacobian matrix. Does your code include a function or class where the jacobian matrix being built? If yes, where it is located in the package?

Ehsan@Kannapiran:

Download the file, extract it, and run FVToolStartUP

Kannapiran ArasuHow to use or install this ? Thannks.

Junaid AhmadMITHRILHi,

Could you tell me how to define a non-uniform velocity field (say a Couette flow) in a 2D model, please?

Thank you.

Ann AnnmarziehSAHi adelAvaneesh NarlaI tried running the diffusiontutorial_spherical code, but I commented out line 20, and uncommented line 19 (switching BCs from left to right), and the computation no longer performs as expected. As I am just making a symmetric change in the BCs, is there a reason it doesn't work?

Thanks!

EhsanI mostly used "An Introduction to Computational Fluid Dynamics: The Finite Volume Method Approach" by Versteeg and Malalasekera, except for the boundary conditions that are my own design.

dinesh aravinthHi sir,

Can you specify which book you followed for writing this code?

Du Wangcatalina medinaHi!

I'm trying to use the FVM toolbox you created for MATLAB to calculate the temperature of a fluid being injected in an oil well that is closed on the bottom. The fluid goes down trough the space between the production tube and the casing of the well and then it goes up through the production tube. I need to calculate the heat transfer between the rocks and the fluid going down and also between the fluid going downwards and the fluid going upwards. There's a paper that explains the problem very well and it also has the needed equations. It's called "geothermal energy production utilizing abandoned oil and gas wells" by Xianbiao Bu. Do you think you'r code could do the job? I'm having trouble using it for this purpose. Any help will be very welcome!

Best regards,

Catalina.

Ali Akbar EftekhariHi Erdum,

The nonuniform mesh is created using the coordinates of the faces of the cells. You should use the cell centers coordinates:

```

x=[0 0.12 0.13 0.2 0.4 0.41 0.56 0.8 1]; %non-uniform grid

m = createMesh1D(x);

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

X = m.cellcenters.x;

D_val = sin(X)+2; % value of the diffusion coefficient

D = createCellVariable(m, D_val); % assign the diffusion coefficient to the cells

D_face = harmonicMean(D); % calculate harmonic average of the diffusion coef on the cell faces

Mdiff = diffusionTerm(D_face); % matrix of coefficients for the diffusion term

[Mbc, RHSbc] = boundaryCondition(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(c); % visualize the results

```

Erdem Uguzhello,

I am trying to use the toolbox for non-uniform grid with nonconstant diffusion coefficient. I can do the non-uniform grid (i.e. x).

I replaced D_val with sin(x)+2 and I get D values are 0 therefore result as NAN. How can I do this? I will extend this to 2D and 3D of course.

x=[0 0.12 0.13 0.2 0.4 0.41 0.56 0.8 1]; %non-uniform grid

>> m = createMesh1D(x);

>> 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 = sin(x)+2; % value of the diffusion coefficient

D = createCellVariable(m, D_val); % assign the diffusion coefficient to the cells

D_face = harmonicMean(D); % calculate harmonic average of the diffusion coef on the cell faces

Mdiff = diffusionTerm(D_face); % matrix of coefficients for the diffusion term

[Mbc, RHSbc] = boundaryCondition(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(c); % visualize the results

Mehul FadnavisMd Shujan AliImpressive and helpful

SfalahatiHi guys,

I need to write a code for CFD to solve the difference heat equation and conduct 6 cases simulations.

Equation: (

Thanh Hung VoDear Sir. Could Sir help me explain how can i use these code for enhanced oil recovery. Thank you so much

Yongjia ZhangMarco dos Santos BernardesAli Akbar EftekhariHi all,

Please ask your questions by opening an issue in the github repository (https://github.com/simulkade/FVTool) orby writing a comment in my blog (http://fvt.simulkade.com). I don't have access to this page anymore, since my previous email at TU Delft is not active.

I'm adding a few examples to the github page to answer your questions regarding the Poisson equation and the source terms.

In short, you can have a derivative boundary condition by changing the BC lines in the above code to:

BC.left.a(:) = 1; BC.left.b(:)=0; BC.left.c(:)=1; % Neumann for the left boundary

BC.right.a(:) = 1; BC.right.b(:)=0; BC.right.c(:)=0; % Neumann for the right boundary

A source term can be added by calling the function constantSourceTerm(q) where q is a cell variable. Don't forget to divide the source term by the cell volume.

Yue SUNJust want to say thx

Chu-Yao ChouRodward HewlinHi I would like to know how to add a source term to the diffusion equation using your solver.

mohammadrezahi, how I can solve 3d poison equation with derivative boundary conditions in pdetool in rectangular channel domain?

Ali Akbar EftekhariHi Shijie,

I don't visit this page regularly. You can always download the latest version from github:

https://github.com/simulkade/FVTool

Ali

Shijie liuHi Ali,

I don't know why I can't download the attached .zip file. And I have tried many times. Could you send the files to me? My email is shi.jieliu@163.com. Thanks a lot!

Saif ManjiYuri FeldmanHi Ali,

Can you please explain the following corrections in your upwind scheme for convection term:

% Also correct for the boundary cells (not the ghost cells)

% Left boundary:

APx(1,:) = APx(1,:)-uw_max(1,:)/(2*DXp(1)); AW(1,:) = AW(1,:)/2;

% Right boundary:

AE(end,:) = AE(end,:)/2; APx(end,:) = APx(end,:) + ue_min(end,:)/(2*DXp(end));

% Bottom boundary:

APy(:,1) = APy(:,1)-vs_max(:,1)/(2*DYp(1)); AS(:,1) = AS(:,1)/2;

% Top boundary:

AN(:,end) = AN(:,end)/2; APy(:,end) = APy(:,end) + vn_min(:,end)/(2*DYp(end));

I will also appreciate if you could refer me to the relevant literature.

Thank you,

Yuri

Ali Akbar EftekhariHi Everyone,

Please write your questions or comments preferably in the Github page of the code. My mathworks account is changed and I don't receive notifications anymore.

https://github.com/simulkade/FVTool

Ali Akbar EftekhariSorry Alireza for this late reply. I have a new mathworks account so I don't receive notifications anymore. Yes, you can define a heterogeneous field. For instance, in the example above, replace the mesh creation and D definition lines with:

m=createMesh2D(Nx, Nx, L, L);

D=createCellVariable(m, rand(Nx, Nx));

You can replace the rand(Nx, Nx) with any matrix of Nx x Nx size.

Antoniosuper

AGThanks for the magnificent work. Does your code support heterogeneous material properties as well? I am trying to solve a 2D transient heat equation on a domain that has different conductivities and heat capacities and I was hoping your framework could be of help.

Many thanks!

Ali Akbar EftekhariHi Hongwei, I'm glad you find the code useful. Please let me know if you like to add your application to the example folder. You are always welcome to send a pull request on github.

Hongwei GuoThanks a lot !!! Very professional and general code !!! I will try to apply this in electron transport problems !!!

Hongwei GuoGuoxi HEWilliam SchillPeng CaoNikhil KumarMartin