Laminar flow tutorials

Mesh generation

OpenFOAM always operates in a 3 dimensional Cartesian coordinate system and all geometries are generated in 3 dimensions. OpenFOAM solves the case in 3 dimensions by default but can be instructed to solve in 2 dimensions by specifying a ‘special’ empty boundary condition on boundaries normal to the (3rd) dimension for which no solution is required. In 2 dimensions, the directions x and y are always chosen even if the mesh is given in the x-z plane.

The sedimentation domain consists of a rectangle of side length L = 0.006 m in the x–y plane and of height h=0.06 m in the vertical direction z. The side length is chosen to ensure a unity cell aspect ratio. A uniform mesh of 1 by 1 by 100 cells will be used. The block structure is shown in Figure 1.

The mesh generator supplied with OpenFOAM, blockMesh, generates meshes from a description specified in an input dictionary, blockMeshDict located in the constant/polyMesh directory for a given case. The blockMeshDict entries for this case are as follows:

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  1.7.1                                 |
|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    object      blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

scale 0.06;

vertices
(
    (-0.005 0   0.005) //0
    ( 0.005 0   0.005) //1
    ( 0.005 1.  0.005) //2
    (-0.005 1.  0.005) //3
    (-0.005 0  -0.005) //4
    ( 0.005 0  -0.005) //5
    ( 0.005 1. -0.005) //6
    (-0.005 1. -0.005) //7
);

blocks
(
    hex (0 1 5 4 3 2 6 7) (1 1 200) simpleGrading (1 1 1)
);

edges
(
);

boundary
(
    inlet
    {
     type   cyclic;
     neighbourPatch outlet;
     faces ((0 3 7 4));
    }
    outlet
    {
     type   cyclic;
     neighbourPatch  inlet;
     faces  ((1 5 6 2));
    }
    top
    {
     type wall;
        faces ((7 6 2 3));
    }
    bottom
    {
     type wall;
     faces ((0 4 5 1));
    }
    frontAndBackPlanes
    {
    type empty;
    faces (
           (0 1 2 3)
           (4 7 6 5)
          );
    }
);

mergePatchPairs
(
);

// ************************************************************************* //

The file first contains header information in the form of a banner (lines 1-7), then file information contained in a FoamFile sub-dictionary, delimited by curly braces ({…}).

The file specifies coordinates of the block vertices; then, it defines the blocks (here, only 1) from the vertex labels and the number of cells within it; and finally, it defines the boundary patches. The user is encouraged to consult the official documentation here to understand the meaning of the entries in the blockMeshDict file.

The mesh is generated by running the command blockMesh. This is done simply by typing in the terminal:

blockMesh

The running status of blockMesh is reported in the terminal window. Any mistakes in the blockMeshDict file are picked up by blockMesh and the resulting error message directs the user to the line in the file where the problem occurred. There should be no error messages at this stage.

Boundary and initial conditions

Once the mesh generation is complete, the user can look at this initial fields set up for this case. The case is set up to start at time t = 0 s, so the initial field data is stored in a 0_org sub-directory of the 1DSedim directory. The 0_org sub-directory contains several files: alpha.a, U.a, U.b, p_rbgh, muI, pa, Theta, one for each of the solid phase volume fraction (alpha.a), velocities of solid and fluid phases (U.a and U.b), the reduced pressure (p_rbgh) fields whose initial values and boundary conditions must be set. Let us examine file U.a:

dimensions      [0 1 -1 0 0 0 0];

internalField   uniform (0 0 0);

boundaryField
{
    inlet
    {
        type            cyclic;
    }
    outlet
    {
        type            cyclic;
    }
    top
    {
        type            zeroGradient;
    }
    bottom
    {
        type            fixedValue;
        value           uniform (0 0 0);
    }
    frontAndBackPlanes
    {
        type            empty;
    }
}

There are 3 principal entries in field data files:

-dimensions: specifies the dimensions of the field, here velocity, i.e. m s ;

-internalField: the internal field data which are set to be uniform;

-boundaryField: the boundary field data that includes boundary conditions and data for all the boundary patches.

For this case, boundary conditions consist of cyclic boundaries at the inlet and outlet (left and right boundaries) and walls at the top and bottom boundaries, the top boundary is zeroGradient whereas the bottom boundary is a no-slip boundary conditions i.e zero velocity at the wall. The frontAndBack patch represents the front and back planes of the 2D case and therefore must be set as empty.

The user can, similarly, examine the fluid phase velocity field in the 0_org/U.b file. The solid phase volume fraction has zeroGradient boundary conditions at the top and bottom boundaries whereas the pressure has fixedValue at the top boundary (equivalent to setting the atmospheric pressure) and a fixedFluxPressure at the bottom boundary. We shall explain below what this boundary condition means.

Physical properties

The physical properties for the case are stored in dictionaries whose names are given the suffix …Properties, located in the Dictionaries directory tree. For a SedFoam case, there are 8 property files: transportProperties, interfacialProperties, ppProperties, granularRheologyProperties, kineticTheoryProperties, RASProperties twophaseRASProperties

In the transport property file, the physical properties of the solid and fluid phases are set, such as density, kinematic viscosity and particle diameter.

// * * * * * * * * * * * * sediment properties * * * * * * * * * * * * //
phasea
{
    rho             rho [ 1 -3 0 0 0 ] 1050;
    nu              nu [ 0 2 -1 0 0 ] 1e-6;
    d               d [ 0 1 0 0 0 0 0 ] 290e-6;
    hExp            hExp [ 0 0  0 0 0 0 0 ] 3.15;  // hindrance exponent for drag: beta^(-hExp) (2.65 by default)
}
// * * * * * * * * * * * * fluid properties * * * * * * * * * * * * //
phaseb
{
    rho             rho [ 1 -3 0 0 0 ] 950;
    nu              nu [ 0 2 -1 0 0 ] 2.105e-05;
    d               d [ 0 1 0 0 0 0 0 ] 290e-6;
}
//*********************************************************************** //
transportModel  Newtonian;

nu              nu [ 0 2 -1 0 0 0 0 ] 2.105e-05;

nuMax           nuMax [0 2 -1 0 0 0 0] 1e0;      // viscosity limiter for the Frictional model (required for stability)

alphaSmall      alphaSmall [ 0 0 0 0 0 0 0 ] 1e-6;  // minimum volume fraction (phase a) for division by alpha

The hindrance function exponent of the Richardson-Zaki law can be tuned using the parameter hExp. This is important for this first test case which is in the laminar regime for the drag. The transportModel has to be prescribed when using turbulence model, it is present here for compatibility with the other test cases. The nuMax value is a viscosity limiter for the all the viscosities, it allows the user to control the numerical stability of the solver and will be discussed in the second test case. It is especially important as some shear stress terms are discretized explicitly and the frictional viscosity can reach huge values. The alphaSmall value is used to regularise the division by alpha.a in the source code, the default value is set to .

The dictionary interfacialProperties allows to specify the drag model:

dragModela      GidaspowSchillerNaumann;

dragModelb      GidaspowSchillerNaumann;

dragPhase       a;

The model chosen here is GidaspowSchilerNaumann that corresponds to the drag force around a spherical particle corrected by a hindrance function ( ).

The dictionary ppProperties contains the information concerning the permanent contact contribution to the solid phase pressure. In SedFoam, the model of Johnson and Jackson (1987) is used:

ppModel          JohnsonJackson;

alphaMax         alphaMax [ 0 0 0 0 0 0 0 ] 0.635;

alphaMinFriction alphaMinFriction [ 0 0 0 0 0 0 0 ] 0.57;

Fr               Fr [ 1 -1 -2 0 0 0 0 ] 5e-2;

eta0             eta0 [ 0 0 0 0 0 0 0 ] 3;

eta1             eta1 [ 0 0 0 0 0 0 0 ] 5;

packingLimiter  no;

And the dictionary granularRheologyProperties controls the granular rheology part:

granularRheology   on;

granularDilatancy  off;

granularCohesion off;

alphaMaxG       alphaMaxG [ 0 0 0 0 0 0 0 ] 0.6;

mus             mus  [ 0 0 0 0 0 0 0 ] 0.24;

mu2             mu2  [ 0 0 0 0 0 0 0 ] 0.39;

I0              I0  [ 0 0 0 0 0 0 0 ] 0.01;

Bphi            Bphi [ 0 0 0 0 0 0 0 ] 1.;

n               n [ 0 0 0 0 0 0 0 ] 1.;

Dsmall          Dsmall [ 0 0 -1 0 0 0 0 ] 1e-4;

relaxPa         relaxPa [ 0 0 0 0 0 0 0 ] 1;

FrictionModel   Coulomb;

PPressureModel  none;

FluidViscosityModel Einstein;

On the first line the switch granularRheology allows to turn on and off this module. The keywords FrictionModel, PPressureModel and FluidViscosityModel allow to choose the friction model, the shear induced particle pressure model and the effective viscosity model (see Input Parameters for details).

For the present case the granularRheology is set to Coulomb and the FluidViscosityModel is used for the effective viscosity model and an Einstein model is used. All the details concerning the different models can be found in the Geophysical Model Development paper on SedFoam.

The remaining dictionaries are not used for this tutorial and will not be detailed here. Their content will be explained in the next tutorials when they are used.

Control

Input data related to the control of time and reading and writing of the solution data are read from the controlDict dictionary. The user should view this file as a case control file, it is located in the system directory.

The start/stop times and the time step for the run must be set. In this tutorial we wish to start the run at time t = 0 which means that OpenFOAM needs to read field data from a directory named 0. Therefore we set the startFrom keyword to startTime and then specify the startTime keyword to be 0.

As the simulation progresses we wish to write results at certain intervals of time that we can later view with a post-processing package. The writeControl keyword presents several options for setting the time at which the results are written; here we select the adjustableRunTime option which specifies that results are written every writeInterval seconds of simulated time, adjusting the time steps to coincide with the writeInterval if necessary and the time step is automatically adjusted using adjustTimeStep = on. The time step controlled using different criteria on Courant numbers based on the maximum velocity, the relative velocity ( ) and the interfacial Courant number. The maximum values of the Courant numbers (maxCo and maxAlphaCo) are set in the controlDict and the maximum value of the time step is controlled by the parameter maxDeltaT.

OpenFOAM creates a new directory named after the current time, e.g. 20 s, on each occasion that it writes a set of data. For this case, the entries in the controlDict are shown below:

application     sedFoam;

startFrom       latestTime;

startTime       0;

stopAt          endTime;

endTime         1800;

deltaT          1e-5;

writeControl    adjustableRunTime;

writeInterval   20;

purgeWrite      0;

writeFormat     ascii;

writePrecision  6;

timeFormat      general;

timePrecision   6;

runTimeModifiable on;

adjustTimeStep  on;

maxCo           0.3;

maxAlphaCo      0.3;

maxDeltaT       2e-1;

Discretisation and linear-solver settings

The numerical schemes and linear solvers are set in the system/fvSchemes and system/fvSolution files respectively.

Concerning the numerical schemes, the time discretisation is an Euler implicit scheme, the divergence terms are of Gauss limitedLinear type (2nd order centered schemes with a Sweby limiter) and the laplacian scheme is Gauss linear corrected for all terms. The fvScheme file is shown below:

ddtSchemes
{
    default         Euler implicit;
}

gradSchemes
{
        default         Gauss linear;
}

divSchemes
{
    default         none;
// alphaEqn
    div(phi,alpha)  Gauss upwind;//limitedLinear01 1;
    div(phir,alpha) Gauss upwind;//limitedLinear01 1;
// UEqn
    div(phi.a,U.a)    Gauss limitedLinearV 1;
    div(phi.b,U.b)    Gauss limitedLinearV 1;
    div(phiRa,Ua)   Gauss limitedLinear 1;
    div(phiRb,Ub)   Gauss limitedLinear 1;
// pEqn
    div(alpha,nu)   Gauss linear;
// alphaPlastic
    div(phia,alphaPlastic)   Gauss limitedLinear01 1;
// pa
    div(phia,pa_new_value)   Gauss limitedLinear 1;
}

laplacianSchemes
{
    default         none;
// UEqn
    laplacian(nuEffa,U.a) Gauss linear corrected;
    laplacian(nuEffb,U.b) Gauss linear corrected;
    laplacian(nuFra,U.a)  Gauss linear corrected;
// pEqn
    laplacian((rho*(1|A(U))),p_rbgh) Gauss linear corrected;
}

interpolationSchemes
{
    default  linear;
}

snGradSchemes
{
    default         corrected;
}

fluxRequired
{
    default         no;
    p_rbgh               ;
}

In the system/fvSolution file, the user can specify the linear solver he/she wants to use for each equation and the tolerance for the resolution. The user can also control the number of corrector steps and non-orthogonal correctors (useless for cartesian grids) for the PIMPLE algorithm. In the present case, a single PIMPLE corrector is enough the residual of the pressure equation is shown to converge well. The fvSolution file is reproduced below:

solvers
{
    p_rbgh
    {
        solver          GAMG;
        tolerance       1e-9;
        relTol          0.0001;
        smoother        DIC;
        nPreSweeps      0;
        nPostSweeps     2;
        nFinestSweeps   2;
        cacheAgglomeration true;
        nCellsInCoarsestLevel 10;
        agglomerator    faceAreaPair;
        mergeLevels     1;
    }

    p_rbghFinal
    {
        $p_rbgh;
        tolerance       1e-9;
        relTol          0;
    }
    "(alpha.a|U.a|U.b|pa_new_value|alphaPlastic)"
    {
        solver          PBiCG;
        preconditioner  DILU;
        tolerance       1e-9;
        relTol          0;
    }
    "(alpha.aFinal|U.aFinal|U.bFinal|paFinal|alphaPlasticFinal)"
    {
        solver          PBiCG;
        preconditioner  DILU;
        tolerance       1e-9;
        relTol          0;
    }
}

PIMPLE
{
    momentumPredictor 0;
    nOuterCorrectors  1;
    nCorrectors     1;
    nNonOrthogonalCorrectors 0;
    correctAlpha    0;
    nAlphaCorr      1;
    pRefCell        0;
    pRefValue       0;
}


relaxationFactors
{
    fields
    {
        p_rbgh               1.;
        p_rbghFinal          1;
    }
    equations
    {
        "U.a|U.b"             1.;
        "(U.a|U.b)Final"      1;
    }

}

Using paraview

The reader is referred to the OpenFOAM tutorial for the lid driven cavity for the basic manipulation on paraview such as viewing loading the data or viewing the mesh (http://cfd.direct/openfoam/user-guide/cavity/).

Using fluidfoam

The postprocessing using python requires the installation of the package fluidfoam available at git clone https://github.com/fluiddyn/fluidfoam. If the package is installed properly with all the dependencies, you just have to run the python script plot_tuto1DSedim.py located in the folder tutorials/Py:

python plot_tuto1DSedim.py

and you should see the following figures:

Mesh generation

The numerical domain consists of a rectangle of side length L = 0.0065 m in the x-y plane and of height h=0.065 m in the vertical direction z. The side length is chosen to ensure a unity cell aspect ratio. A uniform mesh of 1 by 1 by 200 cells will be used. The block structure is the same as for the Sedimentation test case shown in Figure 1.

The mesh is generated by running blockMesh within the case directory:

blockMesh

Boundary and initial conditions

The case is set up to start at time t = 0 s, again the initial field data are stored in a 0_org sub-directory of the case directory. The 0_org sub-directory contains several files: alpha.a, alphaMinFriction, alphaPlastic, U.a, U.b, p_rbgh, muI, pa, pa_new, Theta, one for each of the solid phase volume fraction (alpha.a), velocities of solid and fluid phases (U.a and U.b), the reduced pressure (p_rbgh) fields whose initial values and boundary conditions must be set. Compared with the sedimentation test case, no-slip boundary conditions are imposed for both phase velocity at the top and bottom boundaries.

Physical properties

The transportProperties file for the 1DBedload case is shown below:

// * * * * * * * * * * * * sediment properties * * * * * * * * * * * * //
phasea
{
    rho             rho [ 1 -3 0 0 0 ] 1190;
    nu              nu [ 0 2 -1 0 0 ] 1e-6;
    d               d [ 0 1 0 0 0 0 0 ] 2.e-3;
}

phaseb
{
    rho             rho [ 1 -3 0 0 0 ] 1070;
    nu              nu [ 0 2 -1 0 0 ] 2.52e-4;
    d               d [ 0 1 0 0 0 0 0 ] 2.e-3;
}
//*********************************************************************** //
transportModel  Newtonian;

nu              nu [ 0 2 -1 0 0 0 0 ] 2.52e-4;

nuMax           nuMax [0 2 -1 0 0 0 0] 1e1;      // viscosity limiter for the Frictional model (required for stability)

alphaSmall      alphaSmall [ 0 0 0 0 0 0 0 ] 1e-6;  // minimum volume fraction (phase a) for division by alpha.a

The physical properties of the particles (diameter and density) and the fluid (density and kinematic viscosity) have been set to the experimental values from Aussillous et al. (2013).

In the dictionary interfacialProperties file the drag model from GidaspowShcillerNaumann has been kept. There is little influence of the drag model on this test case as the relative velocity between the fluid and the particles is negligible.

Compared with the sedimentation test case, the setting of granularRheologyProperties is critical. The file is shown below:

granularRheology   on;

granularDilatancy  off;

granularCohesion off;

alphaMaxG       alphaMaxG [ 0 0 0 0 0 0 0 ] 0.6;

mus             mus  [ 0 0 0 0 0 0 0 ] 0.32;

mu2             mu2  [ 0 0 0 0 0 0 0 ] 0.94;

I0              I0  [ 0 0 0 0 0 0 0 ] 0.0077;

Bphi            Bphi [ 0 0 0 0 0 0 0 ] 1.0;

n               n [ 0 0 0 0 0 0 0 ] 2.5;

Dsmall          Dsmall [ 0 0 -1 0 0 0 0 ] 1e-4;

relaxPa         relaxPa [ 0 0 0 0 0 0 0 ] 1;

FrictionModel  Coulomb;
//FrictionModel  MuIv;

PPressureModel  none;

FluidViscosityModel  Einstein;

To compare the numerical results with the analytical solution, the Coulomb Frictional model and the Einstein FluidViscosity model are used. No PPressureModel is specified meaning that the shear induced contribution to the solid phase pressure is not included i.e. the solid phase volume fraction profile must be sharp at the particle bed interface. The value of the rheological constants are chosen following Aussillous et al. (2013). The regularisation parameter Dsmall is taken as 1e-4, this is shown to be enough for this case but a sensitivity analysis could be performed in case of a too strong creeping flow in the particle bed.

Discretisation and linear-solver settings

The numerical schemes and the solvers are the same as in the 1DSedim case.

Control

The system/controlDict file is modified compared with the sedimentation test case. The endTime is fixed to 350 s with an output every 50s with an adjustable time step limited to 5e-3 s.

application     sedFoam;

startFrom       latestTime;

startTime       0;

stopAt          endTime;

endTime         350;

deltaT          0.005;

writeControl    adjustableRunTime;

writeInterval   50.;

purgeWrite      0;

writeFormat     ascii;

writePrecision  6;

timeFormat      general;

timePrecision   6;

runTimeModifiable on;

adjustTimeStep  yes;

maxCo           0.1;

maxAlphaCo      0.1;

maxDeltaT       5e-03;

Mesh generation

Because the avalanche is assumed to be in a stationary state, the simulation can be easily reduced to a one-dimensional problem. The numerical domain consists of a rectangle of side length L = 0.01 made dimensionless by the domain size in the x-y plane and of height h=1 in the vertical direction z.

A uniform mesh of 1 by 1 by 100 cells will be used (see =1 and =100 in figure 7). The block structure is the same as for the Sedimentation test case shown in figure 1.

The mesh is generated by running blockMesh within the case directory:

blockMesh

Boundary and initial conditions

The case is set up to start at time t = 0 s, again the initial field data are stored in a 0_org sub-directory of the case directory. The 0_org sub-directory contains several files: alpha.a, U.a, U.b, p_rbgh, muI, pa, Theta, one for each of the solid phase volume fraction (alpha.a), velocities of solid and fluid phases (U.a and U.b), the reduced pressure (p_rbgh) fields whose initial values and boundary conditions must be set. No-slip boundary conditions are imposed for both phases at the top and bottom boundaries (zero velocity at the solid boundaries). Additionally, cyclic boundaries are considered in the inlet and outlet (upstream and downstream boundaries in figure 7.). Cyclic boundaries are imposed to connect two equal meshes. The frontAndBack patch represents the front and back planes of the 2D case and therefore must be set as empty.

Physical properties

The following variables are considered:

The transportProperties file for the 1DAvalancheMuI case is shown below:

phasea
{
    rho             rho [ 1 -3 0 0 0 ] 1;
    nu              nu [ 0 2 -1 0 0 ] 1e-6;
    d               d [ 0 1 0 0 0 0 0 ] 0.02;
}

phaseb
{
    rho             rho [ 1 -3 0 0 0 ] 1e-3;
    nu              nu [ 0 2 -1 0 0 ] 1e-1;
    d               d [ 0 1 0 0 0 0 0 ] 1;
}
//*********************************************************************** //
transportModel  Newtonian;

nu              nu [ 0 2 -1 0 0 0 0 ] 1e-1;

// Diffusivity for mass conservation resolution (avoid num instab around shocks)
alphaDiffusion  alphaDiffusion [0 2 -1 0 0 0 0] 0e0;

nuMax           nuMax [0 2 -1 0 0 0 0] 1e2;      // viscosity limiter for the Frictional model (required for stability)

alphaSmall      alphaSmall [ 0 0 0 0 0 0 0 ] 1e-6;  // minimum volume fraction (phase a) for division by alpha

Dense granular flow rheology is used as shown in granularRheologyProperties:

granularRheology   on;

granularDilatancy  off;

granularCohesion off;

alphaMaxG       alphaMaxG [ 0 0 0 0 0 0 0 ] 0.6;

mus             mus  [ 0 0 0 0 0 0 0 ] 0.38;

mu2             mu2  [ 0 0 0 0 0 0 0 ] 0.64;

I0              I0  [ 0 0 0 0 0 0 0 ] 0.3;

Bphi            Bphi [ 0 0 0 0 0 0 0 ] 0.31;

n               n [ 0 0 0 0 0 0 0 ] 2.5;

BulkFactor      BulkFactor [ 0 0 0 0 0 0 0 ] 0e0;

PaMin           PaMin [ 0 0 0 0 0 0 0 ] 1e-6;

relaxPa         relaxPa [ 0 0 0 0 0 0 0 ] 1;

tau_inv_min     tau_inv_min [ 0 0 -1 0 0 0 0 ] 1e-3;

FrictionModel   MuI;

PPressureModel  MuI;

FluidViscosityModel  none;

Additionally, FluidViscosityModel model is set to none. Even though no FluidViscosityModel is specified, the shear induced contribution associated to the fluid phase is still present. However, it has insignificant influence on the results due to the low fluid viscosity. On the contrary, the analytical solution assumes that the shear induced contribution to the fluid phase is null.

Finally, the file constant/g is modified to take into account the and -axis component of gravity with magnitude . In this work the plane is inclined to , however, other angles may be considered.

dimensions      [0 1 -2 0 0 0 0];

// tilt angle = 26 deg
value           (0.43837 -0.89879 0);

Physical properties

The following variables are considered to reproduce the experimental results of Pailha et al. (2008):

The transportProperties file for the 1DWetAvalancheMuIv case is shown below:

phasea
{
    rho             rho [ 1 -3 0 0 0 ] 2500;
    nu              nu [ 0 2 -1 0 0 ] 1e-6;
    d               d [ 0 1 0 0 0 0 0 ] 160.e-6;
}

phaseb
{
    rho             rho [ 1 -3 0 0 0 ] 1041.;
    nu              nu [ 0 2 -1 0 0 ] 9.221902017291067e-05;
    d               d [ 0 1 0 0 0 0 0 ] 2.e-3;
}
//*********************************************************************** //
transportModel  Newtonian;

nu              nu [ 0 2 -1 0 0 0 0 ] 9.221902017291067e-05;

// Diffusivity for mass conservation resolution (avoid num instab around shocks)
alphaDiffusion  alphaDiffusion [0 2 -1 0 0 0 0] 0e-8;

nuMax           nuMax [0 2 -1 0 0 0 0] 1e0;      // viscosity limiter for the Frictional model (required for stability)

alphaSmall      alphaSmall [ 0 0 0 0 0 0 0 ] 1e-6;  // minimum volume fraction (phase a) for division by alpha