Modular Lagrangian in OpenFOAM

Contents

The following contains an overview of a typical case including a field-based Lagrangian cloud. The syntax and layout on disk are explained, as are many of the settings. Comparison with the particle-based Lagrangian cloud is made to aid transition for users familiar with the old system and syntax.

1 Case Layout

The following shows the directory layout for a basic case with a cloud named myCloud:

myCase
├── 0
│   ├── epsilon
│   ├── k
│   ├── Lagrangian
│   │   └── myCloud
│   │       ├── d
│   │       ├── number
│   │       └── U
│   ├── nut
│   ├── p
│   └── U
├── Allclean
├── Allrun
├── constant
│   ├── fvModels
│   ├── g
│   ├── Lagrangian
│   │   └── myCloud
│   │       ├── LagrangianModels
│   │       └── physicalProperties
│   ├── momentumTransport
│   └── physicalProperties
└── system
    ├── blockMeshDict
    ├── controlDict
    ├── decomposeParDict
    ├── functions
    ├── fvSchemes
    ├── fvSolution
    └── Lagrangian
        └── myCloud
            ├── LagrangianSchemes
            └── LagrangianSolution

The settings related to the cloud are all within subdirectories named Lagrangian/myCloud. The particle-based Lagrangian layout is similar, except that it uses an un-capitalised lagrangian outer directory name. This is necessary (in both cases) for it to be possible to have multiple clouds in operation in a single case.

The new Lagrangian files mirror those of finite-volume. Fields are defined in 0/Lagrangian/myCloud, models and physical properties in constant/Lagrangian/myCloud, and schemes and solution parameters in system/Lagrangian/myCloud. These sections will be explained in greater depth in the sections below.

1.1 FvModel and FunctionObject Clouds

A cloud is constructed and solved by either an fvModel or a functionObject.

An fvModel cloud is “two-way” coupled; i.e., the drag, heat transfer, mass transfer, etc…, applies to both the Lagrangian particles and the carrier phase. A functionObject cloud, by contrast, is only “one-way” coupled. The source terms that the Lagrangian particles experience are not applied to the carrier phase. The carrier is therefore unaffected by a functionObject cloud. One-way coupled clouds might be used when the particle loading is very low, or when the cloud is being used only for post-processing purposes.

The evolution sequence of the cloud and its carrier are also slightly different between the fvModel or a functionObject clouds. An fvModel cloud will be solved for at the start of the time-step before the carrier phase equations have been solved, whilst a functionObject cloud will be solved at the end. The principle is that an fvModel cloud is primarily used to provide sources to an Eulerian solution, so Eulerian is thought to be more “dependent” on Lagrangian than the other way around. The Lagrangian evolution and the resulting Eulerian sources should therefore be computed before the carrier equations are solved. A functionObject cloud, on the other hand, does not affect the carrier, so it makes sense for it to solve after the carrier so that it is using the most up-to-date carrier properties possible.

An fvModel can also be corrected by the outer PIMPLE iteration, whilst a functionObject cloud cannot. Outer-iteration is new to the field-based Lagrangian library; the particle-based library does not support it.

1.2 Cloud Types

There are currently seven different clouds. These control what fundamental properties are solved and are therefore analogous to the solvers in the Eulerian finite volume system. The clouds are:

• kinematicParticle: A cloud with spherical, constant density, particles in which the carrier interaction is characterised by a cloud/carrier density ratio. Neither the cloud nor the carrier needs the absolute value of the density to be defined. The carrier must therefore also have constant density. For use with incompressible fluid solvers.

• kinematicParcel: As kinematicParticle, but each Lagrangian element is a “parcel” which represents multiple particles

• dynamicParticle: A cloud with spherical, variable density, particles, but no thermodynamic modelling. Can be used with many fluid solvers, but likely to be of most use with incompressible multiphase.

• dynamicParcel: As dynamicParticle, but but each Lagrangian element is a “parcel” which represents multiple particles

• particle: A cloud with spherical particles with thermodynamic modelling. For use with compressible single-phase and multiphase solvers.

• parcel: As particle, but but each Lagrangian element is a “parcel” which represents multiple particles

• tracer: A cloud of massless and volume-less points that follows the carrier flow. For post-processing purposes.

All these clouds have a velocity field, though tracer can initialise that field automatically and does not therefore need the user to provide it. All clouds except tracer also have a diameter field. The dynamic clouds also have a density field, and the particle and parcel clouds have a temperature field. The parcel clouds also have a number field. The number field controls how many identical physical particles are represented by each Lagrangian element, or “parcel”. Parcels are useful because they help manage the expense of simulating large numbers of physical particles, and they also provide the ability to simultaneously specify the number, size (distribution) and volume or mass flow rate of particles at an injection.

The case shown above has d, number and U fields, so it is set up to simulate an kinematicParcel cloud.

As the field-based Lagrangian library is developed, more clouds with additional fundamental properties may be added, such as species fractions and multicomponent thermodynamic modelling, and multiple phases (e.g., to represent solid particles covered in a liquid film).

Each cloud has a corresponding fvModel and a functionObject with the same name as the cloud, plus the suffix Cloud. So, to enable an kinematicParcel cloud as an fvModel, the following can be added in constant/fvModels:

myCloud
{
    type        kinematicParcelCloud;
    libs        ("libLagrangianCloud.so");
}

And the equivalent to create a functionObject cloud in system/functions is:

myCloud
{
    type        kinematicParcelCloud;
    libs        ("libLagrangianCloud.so");
    executeControl timeStep;
    writeControl writeTime;
}

Note that the cloud name is taken from the name of the dictionary, which is also the name of the fvModel or functionObject.

This system, where every cloud has its own fvModel/functionObject is somewhat different to the selection of the particle-based clouds. Those clouds have a single fvModel and functionObject for all cloud types, and then a sub-selection system for the specific cloud model. This additional level of indirection was considered confusing and not necessary in the field-based system.

1.3 Multiphase

When run in conjunction with a multiphase solver module there are options regarding which Eulerian phase is the carrier (i.e., the continuous phase which surrounds the particles) and which Eulerian phase (if any) is considered to be the same phase as the particles.

These options are controlled with carrierPhase and phase entries in the top-level fvModel or functionObject specification. In both cases, these settings are optional.

If the carrierPhase entry is provided, then that phase’s properties (velocity, viscosity, temperature, etc…) will be accessed and used when evaluating most interaction models (drag, lift, heat transfer, etc…). In some cases, to support VoF-like solvers, a fallback to the mixture property will be used if a phase-specific property is not found (e.g., U will be used if U.air is not available). If the carrierPhase entry is not provided, then properties relating to the full mixture of phases will be accessed instead.

If the phase entry is provided then the particles will be considered to be the same phase as the identified Eulerian phase. This allows properties to be looked up from the corresponding phase, rather than being re-defined, and it also allows models to be used which facilitate transfers between the Eulerian and Lagrangian representations. If the phase entry is not provided then the cloud will be considered to be a separate phase entirely.

So, for example, to create an fvModel cloud that represents water droplets within an air-water VoF simulation, the following model specification could be used:

droplets
{
    type        kinematicParticleCloud;
    libs        ("libLagrangianCloud.so");
    phase       water;
    carrierPhase air;
}

A solver module typically provides properties for each phase, and properties for the mixture of all phases. If there are more than two phases, it might be that the carrier is logically a subset of those (e.g., bubbles in an air-water-oil simulation). In that case, the solver will probably not provide information about the mixture of the subset (water-oil, in this example). Additional development would be needed to make that possible.

2 Lagrangian Fields

Lagrangian properties are stored in fields. As for finite-volume, the user needs to provide these fields at the start time to initialise the simulation. The syntax of these fields is very similar to that for finite-volume.

An example velocity field file (0/Lagrangian/myCloud/U), is shown below:

/*--------------------------------*- C++ -*----------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     | Website:  https://openfoam.org
    \\  /    A nd           | Version:  dev
     \\/     M anipulation  |
\*---------------------------------------------------------------------------*/
FoamFile
{
    format      ascii;
    class       LagrangianVectorField;
    location    "0/Lagrangian/myCloud";
    object      U;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

dimensions      [velocity];

internalField   uniform (0 0 0);

boundaryField
{
    #includeEtc "caseDicts/setConstraintTypes"

    walls
    {
        type            reboundVelocity;
        e               0.97;
        mu              0.09;
    }

    open
    {
        type            escapeVelocity;
    }
}


sources
{
    injection1
    {
        type            coneVelocity;
        Umean           (1 0 0);
        thetaInner      0 [deg];
        thetaOuter      5 [deg];
    }
}

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

The four top level entries are as follows:

• The dimensions entry states the dimensions are that of velocity.

• The internalField entry states that all particles are initially at rest. It is possible/probable that there are no particles at the start. In which case an empty list could be specified instead. E.g.:

       internalField   nonuniform 0();
  

• The boundaryField entry determines how the property is affected when a particle encounters a patch. In this case, on the walls patch the velocity is modified to represent a rebound, and on the open patch the velocity is not modified and the particle escapes.

• The sources entry defines what property values are applied to elements that are created by injection-like models and/or what values are added to by transfer-like models. In this case, it states that the injection1 model should initialise the velocities of injected elements so that they spray in a cone shape with a given angle and mean velocity.

2.1 Lagrangian Boundary Conditions

There are fifteen Lagrangian boundary conditions that can be used in the boundaryField section described above. Most of these (eleven) are the constraint conditions (cyclic, empty, processor, wedge, etc …). These can be specified without any additional entries, but, as with finite volume, they are also constructed in patch groups with names matching the boundary condition type. This means that they can be set in bulk using an #includeEtc "caseDicts/setConstraintTypes" entry as shown above.

There are three non-constraint velocity boundary conditions; reboundVelocity, escapeVelocity and stickVelocity. These govern how the velocity is affected by interaction with the boundary, and also indirectly whether or not the particle is retained in the simulation after said interaction. Of these, reboundVelocity requires the user to provide a coefficient of restitution and a coefficient for the fractional tangential momentum loss (see the field example above). The other two require no additional parameters.

The last condition is the calculated condition. This condition leaves whatever parameter it is applied to unchanged following interaction with the boundary. It, too, does not require the user to specify any additional parameters.

In most simple cases, velocity is the only field that is instantaneously affected by the boundary. In which case, all other non-constraint boundary conditions will be of calculated type. This can be concisely specified as follows:

boundaryField
{
    #includeEtc "caseDicts/setConstraintTypes"

    ".*"
    {
        type            calculated;
    }
}

Whilst velocity is, at the moment, the only field with a non-trivial boundary field specification, the system provides the opportunity to enable boundary interactions to affect other fields in future. For example, if it was desired for a boundary interaction to transfer away some heat, or break up a parcel (i.e., reduce diameter and increase number), then these requirements could be realised conveniently by implementing further boundary conditions.

2.2 Lagrangian Source Conditions

Boundary conditions for Lagrangian are less varied and are of lower complexity than for finite volume. Source conditions, conversely, are more complicated and there is a larger variety of options on offer. This is representative of the fact that most particles are introduced as the result of injection, rather than by flow rate through a boundary, as is the case with finite volume. Source conditions are used to define the property values of all injected elements.

There are thirty source conditions in total, at present, but about half of these are designed to be selected automatically by the code in order to initialise properties that the user does not need to define directly; for thermodynamic properties for example, or for initialising age or volume or similar.

There are fourteen conditions that the user is likely to interact with directly. A summary of these is as follows:

• uniformFixedValue sets a user-defined value, which may be a Function1 of time. For example, in the 0/Lagrangian/myCloud/d file, the following could be specified to set the diameters of all introduced particles to 1 mm.

       <modelName>
       {
           type            uniformFixedValue;
           uniformValue    1 [mm];
       }
  

• zero sets a hard-coded value of zero. This is also selected automatically by the code in a number of places. This has no additional entries.

• carrier interpolates the value from a corresponding finite-volume carrier field. The name of the field can be specified, but it will default to the name of the Lagrangian field if it is omitted:

       <modelName>
       {
           type            carrier;
           fieldc          U.air;
       }
  

• totalPressureVelocity sets the velocity based on the difference between a specified total pressure and the local interpolated carrier pressure, in a given direction:

       <modelName>
       {
           type            totalPressureVelocity;
           direction       (0 0 0);
           p0              1 [bar];
       }
  

We then have five velocity conditions that relate to injection a point or over a disc. These conditions all specify the inner and outer angles of the cone or fan in which the particles are injected. Some conditions also require a centreline velocity or velocity magnitude or the direction of the injection. Others can look that information up from corresponding model.

• coneVelocity

       <modelName>
       {
           type            coneVelocity;
           Ucentre         (0.1 0 0);
           thetaInner      15 [deg];
           thetaOuter      22.5 [deg];
       }
  

• totalPressureConeVelocity

       <modelName>
       {
           type            totalPressureConeVelocity;
           direction       (1 0 0);
           thetaInner      15 [deg];
           thetaOuter      22.5 [deg];
           p0              1 [bar];
       }
  

• coneDiskVelocity

       <modelName>
       {
           type            coneDiskVelocity;
           Umag            0.4 [m/s];
           thetaInner      5 [deg];
           thetaOuter      30 [deg];
       }
  

• flowRateConeDiskVelocity

       <modelName>
       {
           type            flowRateConeDiskVelocity;
           thetaInner      5 [deg];
           thetaOuter      30 [deg];
       }
  

• fanVelocity

       <modelName>
       {
           type            fanVelocity;
           Ucentre         (0.1 0 0);
           normal          (0 0 1);
           thetaInner      15 [deg];
           thetaOuter      22.5 [deg];
       }
  

Next we have conditions which facilitate initialising property values by randomly sampling user-defined distributions.

• distribution sets values for a property which are randomly sampled from a given distribution:

       <modelName>
       {
           type            distribution;
           distribution
           {
               type            normal;
               Q               0;
               min             200 [kg/m^3];
               max             1800 [kg/m^3];
               mu              1000 [kg/m^3];
               sigma           400 [];
           }
       }
  

Whilst the distribution condition could be used to set a distribution of diameter, in practice it is a common requirement to also require a specified mass or mass flow-rate of particles. Both these requirements can be achieved for parcel clouds by adjusting the number of particles per parcel to match the desired quantity of material. The following source conditions do this. The distributionDiameter condition has an identical specification to distribution, but additionally interfaces with a source condition for the number field; either flowRateNumber or totalNumber. These number conditions specify the quantity of material introduced.

• flowRateNumber

       <modelName>
       {
           type            flowRateNumber;
           uniformSize     volume; // number, surfaceArea, mass
           volumetricFlowRate 1 [ml/s]; // massFlowRate
       }
  

• totalNumber

       <modelName>
       {
           type            totalNumber;
           uniformSize     volume; // number, surfaceArea, mass
           mass            1 [g]; // volume
       }
  

Either condition can specify the amount of material in terms of either volume or mass. Additionally, to fully specify the condition, a uniformSize entry must be provided. This tells the condition what size is kept uniform across all the created parcels. It can take a value of either number, surfaceArea, volume or mass. If it is set to number then all the parcels will have the same number of particles, but if it is set to volume then they will all have the same volume, and the parcels representing larger particles will have a correspondingly smaller number. This uniformSize control can be thought of a means of grading the Lagrangian discretisation. If the physics of interest scales in proportion with the volume of the particles, then it would be wasteful to give them equal number and thereby devote a significant proportion of the discretisation to the small particles that do not have a significant affect the simulation.

• Finally, LaplacePressure is a condition which sets the pressure within the particles to that of the interpolated carrier pressure, plus a correction to account for the compressive effect of surface tension.

       <modelName>
       {
           type        LaplacePressure;
           sigma       0.07 [N/m];
       }
  

3 Lagrangian Physical Properties

Lagrangian physical properties are defined in a physicalProperties in the relevant Lagrangian subdirectory of constant (in the above, for example, the file is constant/Lagrangian/myCloud/physicalProperties). This is consistent with where physical properties are defined for Eulerian simulations.

The kinematic clouds require this file, as do the particle and parcel clouds. In the case of the kinematic clouds, this file only contains one setting; the density ratio between the particles and the carrier, rhoByRhoc. For example:

/*--------------------------------*- C++ -*----------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     | Website:  https://openfoam.org
    \\  /    A nd           | Version:  dev
     \\/     M anipulation  |
\*---------------------------------------------------------------------------*/
FoamFile
{
    format      ascii;
    class       dictionary;
    location    "constant/Lagrangian/myCloud";
    object      physicalProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

rhoByRhoc       833.3333333333334;

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

With particle and parcel clouds, this file is used to specify the thermodynamic model. For example, for a cloud of air bubbles, the following model might be used:

/*--------------------------------*- C++ -*----------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     | Website:  https://openfoam.org
    \\  /    A nd           | Version:  dev
     \\/     M anipulation  |
\*---------------------------------------------------------------------------*/
FoamFile
{
    format      ascii;
    class       dictionary;
    location    "constant/Lagrangian/myCloud";
    object      physicalProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

thermoType
{
    type            heRhoThermo;
    mixture         pureMixture;
    transport       const;
    thermo          eConst;
    equationOfState perfectGas;
    specie          specie;
    energy          sensibleInternalEnergy;
}

mixture
{
    specie
    {
        molWeight   28.9;
    }
    thermodynamics
    {
        Cv          718;

        hf          0;
    }
    transport
    {
        mu          1.84e-05;
        Pr          0.7;
    }
}

pressure
{
    type        LaplacePressure;
    sigma       0.07 [N/m];
}


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

The format of the thermodynamic specification is identical to that used in compressible finite-volume solvers. In many cases, settings can be shared and re-used between cases by using #include entries and/or other standard dictionary functions.

The one addition to this file, compared with the equivalent settings required for finite-volume, is the pressure entry. This defines how the internal pressure of the Lagrangian elements is constructed, given the carrier pressure and the particles’ properties. This is a standard source condition for the pressure, but it is specified here rather than in a pressure field file because the pressure field file does not otherwise need to be specified by the user, and because this source condition does not apply to any particular model. This condition is likely either to be a carrier condition, or LaplacePressure as in the example above.

4 Lagrangian Models

Lagrangian models are specified in a file called LagrangianModels in the relevant Lagrangian subdirectory of constant (in the above, for example, the file is constant/Lagrangian/myCloud/LagrangianModels). The models are specified in a simple flat list in a manner similar to that which is used for fvModels.

All models are specified in the LagrangianModels, and are intended to be specified here in the future. The LagrangianModel interface is designed to be general enough to facilitate injections, forces, heat and mass transfers and phase changes, turbulent dispersion, breakup, and more. This is in contrast to the particle-based system in which different interfaces and syntaxes and locations are used for different model categories. Creating a truly general interface to such models is challenging, but it is hoped that the resulting consistency of implementation and syntax provides substantial net benefit in the understanding both for users and those maintaining and developing the code.

4.1 Injection Models

The following is an example of a LagrangianModels file. It contains all the currently available injection models; that is, models that create Lagrangian elements.

/*--------------------------------*- C++ -*----------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     | Website:  https://openfoam.org
    \\  /    A nd           | Version:  dev
     \\/     M anipulation  |
\*---------------------------------------------------------------------------*/
FoamFile
{
    format      ascii;
    class       dictionary;
    location    "constant/Lagrangian/myCloud";
    object      LagrangianModels;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

// Continuously inject at a point
pointInjection
{
    type        pointInjection;
    point       (1 0 0) [cm];
    numberRate  1e5;
}

// Instantaneously inject at a list of locations
manualInjection
{
    type        manualInjection;
    file        "manualInjectionPositions";
    units       [mm];
    time        0.01;
}

// Instantaneously inject to fill a cell-zone to a given density
volumeInjection
{

    type        volumeInjection;
    cellZone    injection;
    numberDensity 2 [mm^-3];
    time        0.025;
}

// Continuously inject through a patch
patchInjection
{
    type        patchInjection;
    patch       top;
    numberRate  100000;
}

// Continuously inject across a disk
diskInjection
{
    type        diskInjection;
    centre      (2 0 0) [cm];
    axis        (1 0 0);
    innerDiameter 2 [mm];
    outerDiameter 2 [mm];
    numberRate  100000;
}

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

These models largely fall into two categories; instantaneous and continuous. Instantaneous injections create all of their elements at an instant, which is typically the start of simulation. Continuous injections create elements at a specified rate, often over the entire duration of the simulation. The controls determining the rate and location of continuously generated elements are usually Function1-s of time.

Note that these models only specify the location and number of elements created. The property values of those elements are controlled entirely using the sources entry in the relevant field file. This separation means the injection models are of limited scope and relatively modest complexity. The principle here is that a larger number of smaller, truly interoperable components results in a flexible system that is more functional than the sum of its parts.

4.2 Force Models

The following models, also specified in the LagrangianModels file, all represent forces that the particles experience. They apply these forces by creating source terms in the particles’ momentum equations, in much the same way as fvModel-based force models create source terms in a momentum equation of one of the finite volume fluid solvers. These models represent simple and intuitively understood dynamic effects and are largely self-explanatory given their names.

gravity
{
    type        gravity;
}

drag
{
    // Standard drag force for isolated spherical particles
    /*
    type        SchillerNaumannDrag;
    */

    // Drag force for densely packed spherical particles
    type        GidaspowErgunWenYuDrag;
    alphaMax    0.65;
}

lift
{
    type        SaffmanMeiLift;
}

pressureGradientForce
{
    type        pressureGradientForce;
}

virtualMass
{
    type        constantCoefficientVirtualMass;
    Cvm         0.5;
}

turbulentDispersion
{
    type        turbulentDispersion;

}

4.3 Thermal Models

The following model enables heat transfer between the particle and its carrier phase in proportion with the temperature difference between them, according to the correlation of Ranz and Marshall for spherical particles in a turbulent carrier flow. It creates a source term in the particles’ energy equations.

heatTransfer
{
    type        RanzMarshallHeatTransfer;
}

This is the only heat-transfer model available at present, but further models can be conveniently added as necessary. The implementation of the Ranz-Marshall model is concise and readable, and only contains about 15 lines of functional code.

The other thermal model is a pressure-work model. This model subtracts the pressure-volume work associated with particle expansion from the particle energy, and adds it to the carrier energy. This model is needed to make (e.g.,) bubbles reduce in temperature as they expand.

pressureWork
{
    type        pressureWork;
}

4.4 Mass Transfer Models

Mass transfer models enable the transfer of mass between a cloud and the surrounding continuous fluid. Two such models are currently available.

The first model transfers mass into the carrier phase with a constant surface volume or mass flux. This model provides a crude approximation of phase change which neglects any thermal effects. It can be used by clouds and in solvers that do not feature any thermodynamic modelling.

constantFluxCarrierTransfer
{
    type        constantFluxCarrierTransfer;
    volumeFlux  1 [ml/cm^2/s];
}

The second model represents the absorption of droplets or bubbles into a corresponding Eulerian phase as a result of collisions with that phase (and potentially other non-carrier Eulerian phases). This can be used to represent (e.g.) a spray of droplets hitting and being absorbed by a VoF interface.

collisionPhaseTransfer
{
    type         collisionPhaseTransfer;
}

These models can be used simultaneously and either can result in a particle being removed entirely. The system is constructed so that the mass or volume will all be accounted for and property equations are not rendered ill-posed in the event of particle removal, even in the presence of multiple active mass-transfer models.

4.5 Further Models …

The intended next steps are to add phase change models (i.e., thermally-coupled mass-transfer models), and breakup models, similar to those that exist in the particle-based Lagrangian library.

5 Lagrangian Schemes

Lagrangian schemes are specified in a LagrangianSchemes file, in the relevant Lagrangian subdirectory of system (in the above, for example, the file is system/Lagrangian/myCloud/LagrangianSchemes). Much like for fvSchemes, this file contains settings related to numerics. The way in which tracking, interpolation, time-integration, averaging, accumulation and application of sources is done is all controlled here.

An example of this file is given below:

/*--------------------------------*- C++ -*----------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     | Website:  https://openfoam.org
    \\  /    A nd           | Version:  dev
     \\/     M anipulation  |
\*---------------------------------------------------------------------------*/
FoamFile
{
    format      ascii;
    class       dictionary;
    location    "system/Lagrangian/myCloud";
    object      LagrangianSchemes;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

tracking    linear; // parabolic;

ddtSchemes
{
    default     Euler; // CrankNicolson;
}

SpSchemes
{
    default     none;
    Sp(D,U)     implicit;
    Sp(D,Uc)    explicit;
    Sp(L,U)     implicit;
    Sp(L,Uc)    explicit;
}

averagingSchemes
{
    default     none;
    alpha       cellPoint; // cell

}

interpolationSchemes
{
    default     cell;
    Uc          cellPoint;
}

accumulationSchemes
{
    default     none;
    v           cell;
    S(v*U)      cellPoint;
}

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

The top level entries are as follows:

• Tracking is the method with which particles are moved through the mesh. linear is a first order method in which particles are assumed to move through each cell with a constant velocity. parabolic is a second order method in which the velocity varies with a constant acceleration.

• Time schemes control how the properties of the particles are advanced in time. Euler is a first order implicit scheme. CrankNicolson is a second order scheme in which implicit and explicit contributions are equally weighted.

  CrankNicolson is required for parabolic tracking. The combination of parabolic tracking and CrankNicolson time-evolution is notable in that it conserves the mechanical energy of the particles.

• Source schemes control how the sources are applied to the particles and to the carrier. They can be either implicit or explicit.

  Note that the conservation of properties (momentum, heat, etc…) that are transferred between the particles and their carrier is only maintained if the phase which is solved first has implicit sources and the phase which is solved last has explicit sources.

  In the above example, Sp(D,U) means the (potentially) implicit drag source term as applied to the particles’ velocity equation, whilst Sp(D,Uc) is the corresponding term applied to the carrier’s velocity equation. The terms with L are similar, but relate to lift rather than drag.

• Averaging schemes control how particle properties are interpolated back to the particles. cell does unweighted averages into the cells to create a piecewise constant average. cellPoint averages at the cell centres and at the points using the coordinates within the tetrahedral decomposition to weight the contributions. The result is a continuous piecewise linear function for the average value, but this comes at a higher cost (notably, it involves parallel communication).

• Interpolation schemes control how carrier properties are interpolated to the particles. Typically either cell or cellPoint is chosen, but there are many other options. These options are the same as for the particle-based Lagrangian library. Note that cellPoint is almost universally used for the carrier velocity, as it is important not to create a velocity that has a component into a wall, as this can cause particles to hit the wall at an infinite frequency.

• Accumulation schemes control how Lagrangian data is summed in order to create a integral quantity on the Eulerian mesh. cell does a basic sum in which a particle contributes to its containing cell only. cellPoint, by contrast, adds both to the cell and to the nearby points, using the coordinates within the tetrahedral decomposition to weight the contributions. The result is more distributed and rather smoother, but it comes at a higher cost.

  In the above example, v is the volume (the accumulation of which is likely to be informing a particle volume-fraction calculation), and S(v*U) denotes the source term into the volume-weighted velocity equation.

6 Lagrangian Solution

Lagrangian solution parameters are specified in a LagrangianSolution file, in the relevant Lagrangian subdirectory of system (in the above, for example, the file is system/Lagrangian/myCloud/LagrangianSolution). This file contains controls relating to the solution algorithm. Currently all the settings here relate to controlling the numbers of iterations and sub-cycles.

At the moment, there are only four controls in this file, though it is thought that this is likely to increase as the system is developed further. An example file with all four controls is shown below:

/*--------------------------------*- C++ -*----------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     | Website:  https://openfoam.org
    \\  /    A nd           | Version:  dev
     \\/     M anipulation  |
\*---------------------------------------------------------------------------*/
FoamFile
{
    format      ascii;
    class       dictionary;
    location    "system/Lagrangian/myCloud";
    object      LagrangianSolution;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

outerCorrectors yes;

maxTimeStepFraction 0.3;

maxCellLengthScaleFraction 0.3;

nCorrectors     1;

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

The entries are as follows:

• The outerCorrectors control determines if the cloud is solved every time the outer PIMPLE loop is executed, or if it is just solved once at the start. This setting is only required by fvModel clouds and when multiple outer correctors are being used.

• The maxTimeStepFraction control specifies the maximum fraction of the time-step that particles may be tracked before their equations are solved and their properties updated. The lower this value, the more sub-steps the particles do.

• The maxCellLengthScaleFraction is similar to maxTimeStepFraction, but it specifies a maximum fraction of the local cell size that may be tracked between each property update.

• The nCorrectors setting controls the number of times the equations are solved and the properties updated for each tracking sub-step. This can improve the convergence and implicit-ness of the inner solution in a similar manner to inner PIMPLE iterations. It does not result in additional tracking.

Together, the two max*Fraction controls effectively specify the number of sub-steps that the particles do. In the parcel-based system, these controls were merged into a single number, which was referred to as maxCo; presumably meaning a maximum Courant number. In the field-based system it was concluded that this isn’t formally a Courant number, it would be better to name it more explicitly, and that it was only possible to do so if the two ways in which it was used were specified separately.

cellValueSourceCorrection

The particle-based system has a setting called cellValueSourceCorrection. This adds the exchanged momentum (heat, mass, …) of the current particle to the interpolated value of the velocity (temperature, density, …) so that subsequent uses of the interpolated value were more up to date. The idea (it is assumed) was that this reduces the splitting error associated with solving the particles and the carrier separately.

In the field-based system, this idea has been discarded. It might be of benefit if the particle sources are directly causing sub-time-step transient behaviour, but in many other cases (e.g., the carrier and the particles are at equilibrium) then it is likely to have a detrimental effect.

Instead, coupling can be improved by means of outer iteration, which is now possible for fvModel based clouds. In addition, carrier properties are now interpolated in time as well as in space; e.g.; if a particle is one quarter of the way through its evolution, then the carrier velocity it experiences is 75% of the old-time value and 25% of the new-time value. These two properties mean that, if necessary, it it possible to iterate the system to a point where the splitting error is formally second-order.

7 Lagrangian Function Objects

A number of function objects have been implemented that interact with the field-based Lagrangian system. These functions are “normal” function objects and are specified in system/functions in the same way as function objects that interact with the Eulerian system.

7.1 Non-Cloud Function Objects

Two function objects are “non-cloud” in that they don’t relate to the cloud solution procedure; they only operate on and post-process Lagrangian field data. These are LagrangianFieldValue, which performs reduction operations on Lagrangian values, and LagrangianDistribution, which outputs a plots of the distribution of Lagrangian values.

The following functions entry will, for example, sum the kinetic energy (KE) and gravitational potential energy (GPE) of all the particles in the cloud, and write out a plot of these sums against time:

sumEnergy
{
    type        LagrangianFieldValue;
    libs        ("libLagrangianFunctionObjects.so");
    Lagrangian  myCloud;
    operation   sum;
    fields      (KE GPE);
    weightField number;
}

And this functions entry will write a plot of the diameter distribution:

diameterDistribution
{
    type        LagrangianDistribution;
    libs        ("libLagrangianFunctionObjects.so");
    Lagrangian  myCloud;
    field       d;
    weightField number;
    nBins       20;
    setFormat   raw;
}

7.2 Cloud Function Objects

There are another sixteen “cloud” functions which do interact with the cloud solution procedure (i.e., they have to know when particles move, hit faces, etc…). They are still, however, specified in system/functions. For all but one of these functions, the specification is trivial; none of them require any settings other than the model type, the library and the name of the cloud. So, specific examples will not be given here. They can all be specified in the following general way.

<cloudFunctionName>
{
    type        <cloudFunctionType>;
    libs        ("libLagrangianCloudFunctionObjects.so");
    cloud       myCloud;
}

Where <cloudFunctionName> should be set to a unique name by the user to meaningfully identify the function, and <cloudFunctionType> is the type of the function itself; i.e., one of the following:

• cloudPosition writes a Lagrangian field of the positions of the particles. The positions are not written explicitly otherwise. The particle location is stored and written as addressing defining the tetrahedron containing the particle centre and the local (Barycentric) coordinates within that tetrahedron.

• cloudVolume writes a Lagrangian field of the volumes of the particles

• cloudSurfaceArea writes a Lagrangian field of the surface areas of the particles

• cloudMass writes a Lagrangian field of the masses of the particles

• cloudVolumeFraction writes an Eulerian vol-field of the volume fraction of the particles. I.e., the result of an accumulation scheme being applied to the particle volumes, and then divided by the cell volumes.

• cloudSurfaceAreaPerUnitVolume writes an Eulerian vol-field of the surface area of the particles per unit volume. I.e., the result of an accumulation scheme being applied to the particle surface areas, and then divided by the cell volumes.

• cloudLagrangianVolumeFraction writes a Lagrangian field of the particulate volume fraction. I.e., the result of an averaging scheme being applied to the particle volumes.

• cloudAge writes a Lagrangian field containing the “age” of the particles; i.e., the amount of physical time that has elapsed since they were injected

• cloudKineticEnergy writes a Lagrangian field of the kinetic energy of the particles

• cloudGravitationalPotentialEnergy writes a Lagrangian field of the gravitational potential energy of the particles

• cloudNumberFlux writes an Eulerian surface-field of the number flux of the particles

• cloudVolumeFlux writes an Eulerian surface-field of the volume flux of the particles

• cloudMassFlux writes an Eulerian surface-field of the mass flux of the particles

• cloudBoundaryCollisionNumberFlux writes an Eulerian boundary-field of the number flux of particles that rebound off the patch faces

• cloudBoundaryCollisionForce writes an Eulerian boundary-field of the force exerted on the patch faces

The last function is cloudSurfaceDistribution, which writes a plot of the distribution of the values of particles that pass through a face-zone, face-set or patch. It is similar to LagrangianDistribution, but restricted to a surface, and accumulated over time. This does require some specification, though most of this is identical to the LagrangianDistribution function described above. The only difference is that a patch entry must be specified, and (as this is a cloud function object) a cloud must be specified rather than Lagrangian. For example:

outletDiameterDistribution
{
    type        LagrangianDistribution;
    libs        ("libLagrangianFunctionObjects.so");
    cloud       myCloud;
    patch       outlet;
    field       d;
    weightField number;
    nBins       20;
    setFormat   raw;
}