# Get from

### Highlights from A Simple Finite Volume Solver for Matlab

• BLanalytical(muw, muo)
BLANALYTICAL Analytical solution of Buckley-Leverett equation
• BLanalyticalCorey(muw, muo)
BLANALYTICAL Analytical solution of Buckley-Leverett equation
• BLanalyticalCoreyFoam(mug...
BLANALYTICAL Analytical solution of Buckley-Leverett equation
• BuckleyLeverett1D(MeshStr...
This function solves the continuity equation for an incompressible fluid
• BuckleyLeverett2D(MeshStr...
This function solves the continuity equation for an incompressible fluid
• BuckleyLeverett2D_foam(Me...
This function solves the continuity equation for an incompressible fluid
• BuckleyLeverett2D_pc(Mesh...
This function solves the continuity equation for an incompressible fluid
• BuckleyLeverett3D(MeshStr...
This function solves the continuity equation for an incompressible fluid
• CO2GasDensity(T, p)
This function uses the Span and Wagner equation of state for the pure CO2
• CouetteFlow( varargin )
UNTITLED Summary of this function goes here
• DuanLiquidDensity(T, p, x...
This function uses the model of Duan, Hu, Li, and Mao to calculate the
• DuanMollerWeare(Tr,Pr)
This is the EOS of Duan-Moller-Weare for the CO2-Methane-water system
• DuanSun(T, p, m_salt)
This function solves the Duan-Sun model to calculate the solubility of
• FVToolStartUp()
• GhiaSolution
GhiaSolution returns a benchmark solution for lid-driven cavity flow
• IAPWS_IF97(fun,in1,in2)
IAPWS_IF97(fun,in1,in2)
• JacobsenStewart_N2(p_pa, T)
JACOBSENSTEWART_N2 calculates the molar density of N2 as a function of
• MomentumEq(veloRelax,m, p...
MOMENTUMEQ Summary of this function goes here
• N2viscosity(p_pa, T)
this function calculates the viscosity of N2 as a function of temperature
• P_given_TVN(x)
• PoiseuilleFlow( varargin )
POISOILLEFLOW Summary of this function goes here
• RhieChow( U_cellVar,V_cel...
RHIECHOW Summary of this function goes here
• StefanProblemAnalyticalSo...
StefanProblemAnalyticalSolution returns the analytical solution (i.e.
• T_given_hP(x,h,P)
• V_given_TPN(x, phase)
Calculates volume for given temperature, pressure and mole number
• ViscosityCO2Water(T, p, x...
This function calculates the viscosity of the pure water, brine, and
• ViscosityCO2brine(T, p, x...
This function calculates the viscosity of the pure water, brine, and
• and(p,q)
this function compares the x, y, and z values of the structures that I use in
• arithmeticMean(phi)
This function gets the value of the field variable phi defined
• arithmeticMean1D(MeshStru...
This function gets the value of the field variable phi defined
• arithmeticMean2D(phi)
This function gets the value of the field variable phi defined
• arithmeticMean3D(MeshStru...
This function gets the value of the field variable phi defined
• atan(p)
UMINUS: this function calculates the atan of a x, y, and z values of
• boundaryCondition(BC)
creates the matrix of coefficient and RHS vector
• boundaryCondition1D(BC)
It creates the matrix of coefficient based on the BC structure provided
• boundaryCondition2D(BC)
It creates the matrix of coefficient based on the BC structureprovided
• boundaryCondition3D(BC)
It creates the matrix of coefficient based on the BC structureprovided
• boundaryConditionCylindri...
It creates the matrix of coefficient based on the BC structureprovided
It creates the matrix of coefficient based on the BC structureprovided
• cellBoundary(phi, BC)
This function calculates the value of the boundary cells and add them
• cellBoundary1D(phi, BC)
This function calculates the value of the boundary cells and add them
• cellBoundary2D(phi, BC)
function phiBC = cellBoundary2D(MeshStructure, BC, phi)
• cellBoundary3D(phi, BC)
function phiBC = cellBoundary2D(MeshStructure, BC, phi)
• cellBoundaryCylindrical3D...
function phiBC = cellBoundary2D(MeshStructure, BC, phi)
function phiBC = cellBoundary2D(MeshStructure, BC, phi)
• co2_sw(x)
co2_sw.m
• combineBC(MeshStructure, ...
COMBINEBC This function combines the boundary condition equations with the
• combineBC1D(MeshStructure...
COMBINEBC This function combines the boundary condition equations with the
• combineBC2D(MeshStructure...
COMBINEBC This function combines the boundary condition equations with the
• constantSourceTerm(phi)
RHS vector for an explicit source term
• constantSourceTerm1D(phi)
RHS vector for a Explicit source term
• constantSourceTerm2D(phi)
RHS vector for a Explicit source term
• constantSourceTerm3D(phi)
RHS vector for a Explicit source term
• convectionTerm(u)
This function uses the central difference scheme to discretize a 2D
• convectionTerm1D(u)
This function uses the central difference scheme to discretize a 1D
• convectionTerm2D(u)
This function uses the central difference scheme to discretize a 2D
• convectionTerm3D(u)
This function uses the central difference scheme to discretize a 2D
• convectionTermCylindrical...
This function uses the central difference scheme to discretize a 1D
• convectionTermCylindrical...
This function uses the central difference scheme to discretize a 2D
• convectionTermCylindrical...
This function uses the central difference scheme to discretize a 2D
This function uses the central difference scheme to discretize a 2D
• convectionTvdTerm(u, phi,...
This function uses the TVD scheme to discretize a
• convectionTvdTerm1D(u, ph...
This function uses the upwind scheme to discretize a 1D
• convectionTvdTerm2D(u, ph...
This function uses the upwind scheme to discretize a 1D
• convectionTvdTerm3D(u, ph...
This function uses the TVD scheme to discretize a 3D
• convectionTvdTermCylindri...
This function uses the upwind scheme to discretize a 1D
• convectionTvdTermCylindri...
This function uses the upwind scheme to discretize a 1D
• convectionTvdTermCylindri...
This function uses the TVD scheme to discretize a 3D
This function uses the upwind scheme to discretize a 1D
• convectionUpwindTerm(u)
This function uses the upwind scheme to discretize a 2D
• convectionUpwindTerm1D(u)
This function uses the upwind scheme to discretize a 1D
• convectionUpwindTerm2D(u)
This function uses the upwind scheme to discretize a 2D
• convectionUpwindTerm3D(u)
This function uses the upwind scheme to discretize a 2D
• convectionUpwindTermCylin...
This function uses the upwind scheme to discretize a 1D
• convectionUpwindTermCylin...
This function uses the upwind scheme to discretize a 2D
• convectionUpwindTermCylin...
This function uses the upwind scheme to discretize a 2D
This function uses the upwind scheme to discretize a 2D
• createBC(meshvar)
function BC = createBC(meshvar)
• createBC1D(meshvar)
function BC = createBC1D(meshvar)
• createBC2D(meshvar)
function BC = createBC2D(meshvar)
• createBC3D(meshvar)
Creates a boundary condition structure from a mesh structure
• createCellVariable(meshva...
cellvar = createCellVariable(MeshStructure, cellval)
• createCellVector(meshvar,...
this function creates a vector field based on the geometry and mesh
• createFaceVariable(meshva...
this function creates a face variable based on the geometry and mesh
• createMesh1D(varargin)
MeshStructure = createMesh1D(Nx, Width)
• createMesh2D(varargin)
• createMesh3D(varargin)
MeshStructure = buildMesh3D(Nx, Ny, Nz, Width, Height, Depth)
• createMeshCylindrical1D(v...
MeshStructure = createMeshCylindrical1D(Nr, Lr)
• createMeshCylindrical2D(v...
MeshStructure = buildMeshCylindrical2D(Nr, Ny, Lr, Ly)
• createMeshCylindrical3D(v...
MeshStructure = createMesh3D(Nr, Ntheta, Nz, Radius, theta, height)
• createMeshCylindrical3D(v...
MeshStructure = createMesh3D(Nr, Ntheta, Nz, Radius, theta, height)
MeshStructure = createMeshRadial2D(Nr, Ntetta, Lr, Tetta)
• ddtTerm(MeshStructure, dt...
function ddt = ddtTerm(MeshStructure, dt, phi)
• ddtTerm1D(MeshStructure, ...
function ddt = ddtTerm1D(MeshStructure, h, dt, phi)
• ddtTerm2D(MeshStructure, ...
function ddt = ddtTerm2D(MeshStructure, dt, phi)
• ddtTerm3D(MeshStructure, ...
Matrix of coefficients and the RHS vector for a transient term
• diffusionTerm(D)
This function uses the central difference scheme to discretize a 2D
• diffusionTerm1D(D)
This function uses the central difference scheme to discretize a 1D
• diffusionTerm2D(D)
This function uses the central difference scheme to discretize a 2D
• diffusionTerm3D(D)
This function uses the central difference scheme to discretize a 2D
• diffusionTermCylindrical1...
This function uses the central difference scheme to discretize a 1D
• diffusionTermCylindrical2...
This function uses the central difference scheme to discretize a 2D
• diffusionTermCylindrical3...
This function uses the central difference scheme to discretize a 2D
This function uses the central difference scheme to discretize a 2D
• divergenceTerm(F)
This function calculates the divergence of a field using its face
• divergenceTerm1D(F)
This function calculates the divergence of a field using its face
• divergenceTerm2D(F)
This function calculates the divergence of a field using its face
• divergenceTerm3D(F)
This function calculates the divergence of a field using its face
• divergenceTermCylindrical...
This function calculates the divergence of a field using its face
• divergenceTermCylindrical...
This function calculates the divergence of a field using its face
• divergenceTermCylindrical...
This function calculates the divergence of a field using its face
This function calculates the divergence of a field using its face
• domainInt(MeshStructure, ...
INTEGRATE integrate the value of phi over the domain defined by the mesh
• domainInt1D(MeshStructure...
INTEGRATE integrate the value of phi over the domain defined by the mesh
• domainInt2D(MeshStructure...
INTEGRATE integrate the value of phi over the domain defined by the mesh
• domainInt3D(MeshStructure...
INTEGRATE integrate the value of phi over the domain defined by the mesh
• domainIntCylindrical1D(Me...
INTEGRATE integrate the value of phi over the domain defined by the mesh
• domainIntCylindrical2D(Me...
INTEGRATE integrate the value of phi over the domain defined by the mesh
• domainIntCylindrical3D(Me...
INTEGRATE integrate the value of phi over the domain defined by the mesh
INTEGRATE integrate the value of phi over the domain defined by the mesh
• eq(p,q)
this function compares the equality condition of x, y, and z values of the structures that I use in
• excludeGhostRHS(MS, RHS)
this function cuts out the RHS values related to ghost cells and returns
• field2d(Nx,Ny,k_avg,V_dp,...
Gaussian filter
• field3d(Nx,Ny,Nz,k_avg,V_...
Gaussian filter
• fluxLimiter(varargin)
This function returns a function handle to a flux limiter of user's
• funceval(f, varargin)
FUNCEVAL calculates the value of the function f at the face variable phi
• ge(p,q)
this function compares the greater than or equality condition of x, y, and z values of the structures that I use in
• geometricMean(phi)
This function gets the value of the field variable phi defined
• geometricMean1D(MeshStruc...
This function gets the value of the field variable phi defined
• geometricMean2D(MeshStruc...
This function gets the value of the field variable phi defined
• geometricMean3D(MeshStruc...
This function gets the value of the field variable phi defined
this function calculates the gradient of a variable in x direction, in
this function calculates the gradient of a variable in x direction
this function calculates the gradient of a variable in x direction
this function calculates the gradient of a variable in x and y direction
this function calculates the gradient of a variable in x and y direction
this function calculates the gradient of a variable in x and y direction
this function calculates the gradient of a variable in x and y direction
• gt(p,q)
this function compares the greater than condition of x, y, and z values of the structures that I use in
• h_given_TVN(x)
Calculates spesific entalphy from given temperature, volume and
• harmonicMean(phi)
This function gets the value of the field variable phi defined
• harmonicMean1D(MeshStruct...
This function gets the value of the field variable phi defined
• harmonicMean2D(MeshStruct...
This function gets the value of the field variable phi defined
• harmonicMean3D(MeshStruct...
This function gets the value of the field variable phi defined
• internalCells(phi)
this function extracts the internal cells from variable phi, which
• le(p,q)
this function compares the less than or equal to condition of x, y, and z values of the structures that I use in
• linearMean(phi)
This function gets the value of the field variable phi defined
• linearSourceTerm(k)
Matrix of coefficients for a linear source term in the form of k \phi
• linearSourceTerm1D(k)
Matrix of coefficients for a linear source term in the form of k \phi
• linearSourceTerm2D(k)
Matrix of coefficients for a linear source term in the form of k \phi
• linearSourceTerm3D(k)
Matrix of coefficients for a linear source term in the form of k \phi
• lt(p,q)
this function compares the less than condition of x, y, and z values of the structures that I use in
This function masks the specified cells by giving them a constant value.
• minus(p,q)
MINUS: this function subtracts the x, y, and z values of the structures that I use in
• mpower(p,q)
• mrdivide(p,q)
RDIVIDE this function divides the x, y, and z values of the structures that I use in
• mtimes(p,q)
TIMES this function multiplies the x, y, and z values of the structures that I use in
• ne(p,q)
this function compares the inequality condition of x, y, and z values of the structures that I use in
• normCellVector(cellvec)
this function calculates the second norm of a cell vector field
• not(p)
UMINUS: this function applies a unary minus to the x, y, and z values of the structures that I use in
• or(p,q)
this function compares the x, y, and z values of the structures that I use in
• phaseChangeEnthalpyMethod...
phaseChangeEnthalpyMethodExample() demonstrates the 2D phase change using
• phaseChangeEnthalpyMethod...
phaseChangeEnthalpyMethodExample() demonstrates the 2D phase change using
• plus(p,q)
this function adds the x, y, and z values of the structures that I use in
• power(p,q)
• rdivide(p,q)
RDIVIDE this function divides the x, y, and z values of the structures that I use in
• reshapeCell(MS, phi)
this function reshapes a vetorized cell variable to its domain shape
• solveExplicitPDE(MeshStru...
SolveExplicitPDE solves for the new value of variable \phi in an explicit
• solvePDE(MS, M, RHS, vara...
SOLVEPDE solves the linear system M x \phi = RHS and returns the value of
• times(p,q)
TIMES this function multiplies the x, y, and z values of the structures that I use in
• tracer_tvd(m, V_dp)
physical values
• transientTerm(phi_old, dt...
• transientTerm1D(MeshStruc...
Matrix of coefficients and the RHS vector for a transient term
• transientTerm2D(MeshStruc...
Matrix of coefficients and the RHS vector for a transient term
• transientTerm3D(MeshStruc...
Matrix of coefficients and the RHS vector for a transient term
• tvdMean(phi, u, FL)
This function gets the value of the field variable phi defined
• tvdMean1D(phi, u, FL)
This function gets the value of the field variable phi defined
• tvdMean2D(phi, u, FL)
This function gets the value of the field variable phi defined
• tvdMean3D(phi, u, FL)
This function gets the value of the field variable phi defined
• uminus(p)
UMINUS: this function applies a unary minus to the x, y, and z values of the structures that I use in
• uplus(p)
UMINUS: this function applies a unary minus to the x, y, and z values of the structures that I use in
• upwindMean(phi, u)
This function gets the value of the field variable phi defined
• upwindMean1D(phi, u)
This function gets the value of the field variable phi defined
• upwindMean2D(phi, u)
This function gets the value of the field variable phi defined
• upwindMean3D(phi, u)
This function gets the value of the field variable phi defined
• visualizeCellVectors(phi_...
VISUALIZECELLS plots the values of cell variable phi.value
• visualizeCellVectors2D(ph...
VISUALIZECELLS plots the values of cell variable phi
• visualizeCellVectors3D(ph...
VISUALIZECELLS plots the values of cell variable phi
• visualizeCellVectorsCylin...
VISUALIZECELLS plots the values of cell variable phi
VISUALIZECELLS plots the values of cell variable phi
• visualizeCells(phi)
VISUALIZECELLS plots the values of cell variable phi.value
• visualizeCells1D(phi)
VISUALIZECELLS plots the values of cell variable phi
• visualizeCells2D(phi)
VISUALIZECELLS plots the values of cell variable phi
• visualizeCells3D(phi)
VISUALIZECELLS plots the values of cell variable phi
• visualizeCellsCylindrical...
VISUALIZECELLS plots the values of cell variable phi
VISUALIZECELLS plots the values of cell variable phi
• visualizeMesh(MS)
VISUALIZECELLS plots the values of cell variable phi.value
• visualizeMesh2D(m)
• BoundaryCondition
• CellVariable
• CellVector
• FaceVariable
• MeshStructure
• diffusiontutorial.m
Transient diffusion equation
• FVTdemo.m
Apparently it must be here!
• wave_equation_1D.m
• tracer_hetrogeneous.m
• poisson_example.m
Poisson equation
• heatconductionfin.m
Transient diffusion equation
• diffusiontutorialExplicit.m
a tutorial adapted from the fipy diffusion 1D example
• diffusiontutorial3D.m
Transient diffusion equation
• diffusiontutorial2D.m
Transient diffusion equation
• FVTool_functions_nonunifo...
FVTool test script
• diffusionODEtutorial.m
a tutorial adapted from the fipy diffusion 1D example
• diffusionNonlinearTutorial.m
Nonlinear diffusion equation: a tutorial
• diffconv_variable_diff.m
Transient diffusion equation
• diffconvCyl3D.m
Transient diffusion equation
• diffconv3D.m
Transient diffusion equation
• convectionTVDexample.m
another simple 1D convection example with periodic BC
• convectionTVDburger.m
BURGER's equation; there are still a few strange outputs when I use TVD
• convectionODEexample.m
another simple convection example using method of lines, still a bit
• convectiondiffusiontutori...
a tutorial adapted from the fipy convection diffusion 1D example
• convectioncoupledExample.m
1D convection coupled example with source terms
Transient diffusion equation
• BL2DcoupledTVD.m
Coupled nonlinear PDE's
• BL2DcoupledCentral.m
Coupled nonlinear PDE's
• BL2Dcoupled.m
Coupled nonlinear PDE's
• FVTool_functions_uniform_...
FVTool test script
• View all files
4.66667
4.7 | 3 ratings Rate this file 208 Downloads (last 30 days) File Size: 577 KB File ID: #46637 Version: 2.0

# A Simple Finite Volume Solver for Matlab

by

### Ehsan (view profile)

• 1 file
• 4.66667

16 May 2014 (Updated )

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

File Information
Description

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:

Acknowledgements

Iapws If97 Functional Form With No Slip inspired this file.

Required Products MATLAB
MATLAB release MATLAB 8.3 (R2014a)
29 Oct 2015 Peng Cao

### Peng Cao (view profile)

08 Sep 2015 Nikhil Kumar

### Nikhil Kumar (view profile)

18 Jul 2014 Martin

### Martin (view profile)

02 Aug 2015 2.0

Support for nonuniform grid types
New objects with backward incompatible APIs

05 Aug 2015 2.0

Edited the example above to work with the new API

18 Aug 2015 1.1

18 Aug 2015 1.2

18 Aug 2015 1.3

updated descriptions

18 Aug 2015 1.4

added support for 2D radial (r, theta) and 3D cylindrical (r, thetta, z)

28 Jan 2016 2.0

showdemo function is not available. I will update it later.