In the following we give a short introduction to the MaFloT functions and how to use the InputFile. For a more detailed description of Finite-Volume discretizations of the Flow and Transport equation and the realization in MaFLoT you can download the following script:


The equations that are solved in the MaFloT function are the mass conservation equation of the fluid using Darcy’s law and the mass conservation equation of the solvent that is transported with the fluid. We assume the fluid to be incompressible, whereas the fluid’s density and viscosity are linear functions of the dissolved substance. The equations are discretized by a standard Finite-Volume scheme.  As the two equations are coupled by the concentration, a sequentially implicit scheme is provided for the solution. For the convective terms in the transport equation a first order upwind scheme is used and the accumulation term is discretized by a backward Euler scheme using constant time steps. Both molecular diffusion and dispersion can be assigned. Note, that a higher order TVD scheme as well as an adaptive time step function exist and can be asked for by email.

A 2009 or newer MATLAB version is needed to run the function. The whole package includes the following functions:

MaFloT.m                                            Main function

    InputFile.m                                           Set-up the flow problem

    Initialize.m                                           Initializes hydraulic conductivities and gravity terms

    PresMat.m                                           Creates the linear pressure system

    Velocity.m                                           Calculates interface based velocities

    Transport.m                                          Sets-up and solves the transport problem

            UpMat.m                                            Creates the linear transport upwind system

            Diffusion.m                                          Constructs diffusive linear system

            Dispersion.m                                        Constructs dispersive linear system

    Output.m                                             Function to control the desired output

        DisplayVariable.m                                    Displays pressure and saturation profiles

The simulator can be started by entering the command “MaFloT” in the command window of Matlab (the functions need to be placed in the actual directory). MaFloT can be used by only editing the Inputfile that has the following self explanatory form and leaving all other functions untouched:

%% GRID PARAMETERS ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------%
global Nf len dx 
len     = [10 10];                                  % physical length of the domain in x and y direction [m]
Nf      = [100 100];  	                        % number of cells in x and y direction
dx      = len./Nf;                                  % cell length [m]
%% PERMEABILITY -----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------%
K       = ones(Nf(1),Nf(2))*5e-9;          % permeability field [m2]
%% Porosity -------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------%
global phi
phi     = ones(Nf(1),Nf(2))*0.44;	% porosity field
% INITIAL CONDITIONS------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------%
global s0 smax
s0 = zeros(Nf(1),Nf(2));                     % Initial saturation (normalized concentration) [-]
smax = 1;                                          % maximum concentration [kg/m3] for normalization
%% PHASE PROPERTIES------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------%
global viscosity density gravity
viscosity = [0.001 0.0001];                % viscosity [kg/m/s] [viscosity(s=0) viscosity(s=smax)]
density   = [1000 1000];                    % density [kg/m3] [density(s=0) density(s=smax)]
gravity   = 9.81;                                % gravity acceleration in y [m/s2]
%% SIMULATION PARAMETER FOR TRANSPORT -----------------------------------------------------------------------------------------------------------------------------%
global dt
timeSim  = 5000;                                % total simulation time [s]
dt          = 60;                                    % time step length [s]
tolpS      = 1.e-4;                                % saturation tolerance of pressure saturation loop  [-]
maxpS  = 50;                                     % maximum number of pressure saturation loops to converge                
%% BC FLUID --------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------%
global Fix ibcs                                
ibcs = zeros(2*sum(Nf),1);                 % type 0:Neumann(N); 1:Dirichlet(D)
Fix  = zeros(2*sum(Nf),1);                  % value N [m2/s] (inflow>0); D [Pa]        
%% BC SOLVENT -----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------%
global FixT
FixT = zeros(2*sum(Nf),1);
%% DIFFUSION AND DISPERSION ----------------------------------------------------------------------------------------------------------------------------------------------------------%
global Dif FlagDif alphal alphat
Dif       = 1e-9;                                % [m2/s] molecular diffusion 
ibcD      = zeros(2*sum(Nf),1);;       % 1 -> Diffusion on boundary cell
alphal   = 1e-2;                              % longitudinal dispersivity [m]
alphat   = 1e-3;                              % transversal dispersivity [m]
%% SOURCE TERMS ----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------%
global Q QT
Q       = zeros(Nf);                        %  source term; inflow positive [m3/s]
QT      = zeros(Nf);                        % normalized concentration of flux [-]

Units are always given in brackets. The form of the Input file respectively all variables shown here must be preserved to avoid errors. Boundary vectors, initial values, source terms, etc. always have to be initialized in the dimensions shown above. Density effects can be neglected by setting gravity = 0. “tolpS” is the user defined maximum residual for the coupling loop and “maxpS” the maximum number of iterations allowed per time step. Viscosity and density are defined as a two entry vector with the first value referring to concentration zero and the second to max concentration. For a detailed explanation on how to assign boundary conditions see the BOUNDARY page.

So far the pressure as well as the saturation profiles at the last time step are written in ascii format to the folder “results” and the corresponding graphics are displayed through the whole simulation but not saved. Additional variables as velocities can be saved. These manipulations can be done in the main file.