Get from Ico-github-logo

Highlights from
A Simple Finite Volume Solver for Matlab

  • Apparently it must be here!
  • ...
    This function uses the upwind scheme to discretize a 1D
  • ...
    This function uses the TVD scheme to discretize a 3D
  • ...
    This function calculates the divergence of a field using its face
  • ...
    This function uses the upwind scheme to discretize a 1D
  • 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
  • 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()
  • 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
  • N2viscosity(p_pa, T)
    this function calculates the viscosity of N2 as a function of temperature
  • P_given_TVN(x)
    Downloaded from the personal page of Professor Sigurd Skogestad
  • RandPermField(k_avrg, phi...
    generates permeability field on a m x n matrix, using the Dykstra-Parsons coefficient.
  • RandPermField3D(k_avrg, p...
    generates permeability field on a m x n x w matrix, using the Dykstra-Parsons coefficient.
  • T_given_hP(x,h,P)
    Downloaded from the personal page of Professor Sigurd Skogestad
  • 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(MeshStruct...
    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(MeshStru...
    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(MeshStr...
    creates the matrix of coefficient and RHS vector
  • boundaryCondition1D(MeshS...
    It creates the matrix of coefficient based on the BC structure provided
  • boundaryCondition2D(MeshS...
    It creates the matrix of coefficient based on the BC structureprovided
  • boundaryCondition3D(MeshS...
    It creates the matrix of coefficient based on the BC structureprovided
  • buildMesh1D(Nx, Width)
    MeshStructure = buildMesh1D(Nx, Width)
  • buildMesh2D(Nx, Ny, Lx, Ly)
  • buildMesh3D(Nx, Ny, Nz, L...
    MeshStructure = buildMesh3D(Nx, Ny, Nz, Width, Height, Depth)
  • buildMeshCylindrical1D(Nr...
    MeshStructure = buildMeshCylindrical1D(Nr, Lr)
  • buildMeshCylindrical2D(Nr...
    MeshStructure = buildMeshCylindrical2D(Nr, Ny, Lr, Ly)
  • capillaryPressureOptions()
    RELPERMOPIONS creates the options structure for the relperm functions
  • cellBoundary(MeshStructur...
    This function calculates the value of the boundary cells and add them
  • cellBoundary1D(MeshStruct...
    This function calculates the value of the boundary cells and add them
  • cellBoundary2D(MeshStruct...
    function phiBC = cellBoundary2D(MeshStructure, BC, phi)
  • cellBoundary3D(MeshStruct...
    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
  • convectionTerm(MeshStruct...
    This function uses the central difference scheme to discretize a 2D
  • convectionTerm1D(MeshStru...
    This function uses the central difference scheme to discretize a 1D
  • convectionTerm2D(MeshStru...
    This function uses the central difference scheme to discretize a 2D
  • convectionTerm3D(MeshStru...
    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
  • convectionTvdTerm(MeshStr...
    This function uses the TVD scheme to discretize a
  • convectionTvdTerm1D(MeshS...
    This function uses the upwind scheme to discretize a 1D
  • convectionTvdTermCylindri...
    This function uses the upwind scheme to discretize a 1D
  • convectionUpwindTerm(Mesh...
    This function uses the upwind scheme to discretize a 2D
  • convectionUpwindTerm1D(Me...
    This function uses the upwind scheme to discretize a 1D
  • convectionUpwindTerm2D(Me...
    This function uses the upwind scheme to discretize a 2D
  • convectionUpwindTerm3D(Me...
    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
  • createBC(MeshStructure)
    Creates a boundary condition structure from a mesh structure
  • createBC1D(MeshStructure)
    Creates a boundary condition structure from a mesh structure
  • createBC2D(MeshStructure)
    Creates a boundary condition structure from a mesh structure
  • createBC3D(MeshStructure)
    Creates a boundary condition structure from a mesh structure
  • createCellVariable(MeshSt...
    cellvar = createCellVariable(MeshStructure, cellval)
  • createCellVector(MeshStru...
    this function creates a vector field based on the geometry and mesh
  • createFaceVariable(MeshSt...
    this function creates a face variable based on the geometry and mesh
  • createMesh1D(Nx, Width)
    MeshStructure = buildMesh1D(Nx, Width)
  • createMesh2D(Nx, Ny, Lx, Ly)
  • createMesh3D(Nx, Ny, Nz, ...
    MeshStructure = buildMesh3D(Nx, Ny, Nz, Width, Height, Depth)
  • createMeshCylindrical1D(N...
    MeshStructure = buildMeshCylindrical1D(Nr, Lr)
  • createMeshCylindrical2D(N...
    MeshStructure = buildMeshCylindrical2D(Nr, Ny, Lr, Ly)
  • 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(MeshStructu...
    This function uses the central difference scheme to discretize a 2D
  • diffusionTerm1D(MeshStruc...
    This function uses the central difference scheme to discretize a 1D
  • diffusionTerm2D(MeshStruc...
    This function uses the central difference scheme to discretize a 2D
  • diffusionTerm3D(MeshStruc...
    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
  • divergenceTerm(MeshStruct...
    This function calculates the divergence of a field using its face
  • divergenceTerm1D(MeshStru...
    This function calculates the divergence of a field using its face
  • divergenceTerm2D(MeshStru...
    This function calculates the divergence of a field using its face
  • divergenceTerm3D(MeshStru...
    This function calculates the divergence of a field using its face
  • divergenceTermCylindrical...
    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
  • eq(p,q)
    this function compares the equality condition of x, y, and z values of the structures that I use in
  • excludeGhostRHS(MeshStruc...
    this function cuts out the RHS values related to ghost cells and returns
  • 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
  • gaswaterrelperm(s, varargin)
    GASRELPERM returns the relative permeability of the gas phase as a
  • 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(MeshStructu...
    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
  • gradientCellTerm(MeshStru...
    this function calculates the gradient of a variable in x direction, in
  • gradientCellTerm1D(MeshSt...
    this function calculates the gradient of a variable in x direction
  • gradientCellTerm2D(MeshSt...
    this function calculates the gradient of a variable in x and y direction
  • gradientCellTerm3D(MeshSt...
    this function calculates the gradient of a variable in x and y direction
  • gradientTerm(MeshStructur...
    this function calculates the gradient of a variable in x direction
  • gradientTerm1D(MeshStruct...
    this function calculates the gradient of a variable in x direction
  • gradientTerm2D(MeshStruct...
    this function calculates the gradient of a variable in x and y direction
  • gradientTerm3D(MeshStruct...
    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=h_given_TVN(x)
    Calculates spesific entalphy from given temperature, volume and
  • harmonicMean(MeshStructur...
    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(MeshStructu...
    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
  • lt(p,q)
    this function compares the less than condition of x, y, and z values of the structures that I use in
  • maskCells(meshstruct, M, ...
    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(MeshStruct...
    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
  • oilwaterrelperm(s, varargin)
    OILRELPERM returns the relative permeability of the oil phase as a
  • or(p,q)
    this function compares the x, y, and z values of the structures that I use in
  • 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
  • relpermOptions()
    RELPERMOPIONS creates the options structure for the relperm functions
  • reshapeCell(MeshStructure...
    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(MeshStructure, M...
    SOLVEPDE solves the linear system M x \phi = RHS and returns the value of
  • sourceExplicitTerm(MeshSt...
    RHS vector for an explicit source term
  • sourceExplicitTerm1D(Mesh...
    RHS vector for a Explicit source term
  • sourceExplicitTerm2D(Mesh...
    RHS vector for a Explicit source term
  • sourceExplicitTerm3D(Mesh...
    RHS vector for a Explicit source term
  • sourceTerm(MeshStructure, k)
    Matrix of coefficients for a linear source term in the form of k \phi
  • sourceTerm1D(MeshStructur...
    Matrix of coefficients for a linear source term in the form of k \phi
  • sourceTerm2D(MeshStructur...
    Matrix of coefficients for a linear source term in the form of k \phi
  • sourceTerm3D(MeshStructur...
    Matrix of coefficients for a linear source term in the form of k \phi
  • times(p,q)
    TIMES this function multiplies the x, y, and z values of the structures that I use in
  • transientCellVariable(phi...
    this function is indeed a very simple tool to create a transient cell
  • transientTerm(MeshStructu...
  • 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(MeshStructure, ph...
    This function gets the value of the field variable phi defined
  • tvdMean1D(MeshStructure, ...
    This function gets the value of the field variable phi defined
  • tvdMean2D(MeshStructure, ...
    This function gets the value of the field variable phi defined
  • tvdMean3D(MeshStructure,u...
    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(MeshStructure,...
    This function gets the value of the field variable phi defined
  • upwindMean1D(MeshStructur...
    This function gets the value of the field variable phi defined
  • upwindMean2D(MeshStructur...
    This function gets the value of the field variable phi defined
  • upwindMean3D(MeshStructur...
    This function gets the value of the field variable phi defined
  • visualizeCells(MeshStruct...
    VISUALIZECELLS plots the values of cell variable phi
  • visualizeCells1D(MeshStru...
    VISUALIZECELLS plots the values of cell variable phi
  • visualizeCells2D(MeshStru...
    VISUALIZECELLS plots the values of cell variable phi
  • visualizeCells3D(MeshStru...
    VISUALIZECELLS plots the values of cell variable phi
  • watergascapillarypressure...
    CAPILLARYPRESSURE gives the capillary pressure and its derivative with
  • watergasrelperm(s, varargin)
    LIQUIDRELPERM returns the relative permeability of the water phase as a
  • wateroilcapillarypressure...
    CAPILLARYPRESSURE gives the capillary pressure and its derivative with
  • wateroilrelperm(s, varargin)
    OILRELPERM returns the relative permeability of the water phase as a
  • BL2D_pc_tvd.m
    2D Buckley-Leverett solution including source and
  • BL2Dcoupled.m
    Coupled nonlinear PDE's
  • BL2DcoupledCentral.m
    Coupled nonlinear PDE's
  • BL2DcoupledTVD.m
    Coupled nonlinear PDE's
  • convectionODEexample.m
    another simple convection example using method of lines, still a bit
  • convectionTVDburger.m
    BURGER's equation; there are still a few strange outputs when I use TVD
  • convectionTVDexample.m
    another simple 1D convection example with periodic BC
  • convectiondiffusiontutori...
    a tutorial adapted from the fipy convection diffusion 1D example
  • diffusionODEtutorial.m
    a tutorial adapted from the fipy diffusion 1D example
  • diffusiontutorial.m
    Transient diffusion equation
  • diffusiontutorial2D.m
    Transient diffusion equation
  • diffusiontutorial3D.m
    Transient diffusion equation
  • diffusiontutorialExplicit.m
    a tutorial adapted from the fipy diffusion 1D example
  • poisson_example.m
    Poisson equation
  • View all files

5.0

5.0 | 1 rating Rate this file 154 Downloads (last 30 days) File Size: 377 KB File ID: #46637
image thumbnail

A Simple Finite Volume Solver for Matlab

by

 

16 May 2014 (Updated )

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

| Watch this File

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

A short html document is available to get you started. Simply type
showdemo FVTdemo
after downloading and extracting the code.

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(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 some animated results of this code in my youtube channel:
https://www.youtube.com/user/processsimulation/videos

Acknowledgements

Iapws If97 Functional Form With No Slip inspired this file.

Required Products MATLAB
MATLAB release MATLAB 8.3 (R2014a)
Tags for This File   Please login to tag files.
Please login to add a comment or rating.
Comments and Ratings (1)
18 Jul 2014 Martin  
Updates
17 May 2014

add youtube channel link to descriptions

20 May 2014

update my weblog address

24 May 2014

updated descriptions

Contact us