Merge branch 'master' of github.com-OpenFOAM:OpenFOAM/OpenFOAM-dev

bin/foamGet

@@ -3,7 +3,7 @@
 # =========                 |
 # \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
 #  \\    /   O peration     | Website:  https://openfoam.org
-#   \\  /    A nd           | Copyright (C) 2018-2020 OpenFOAM Foundation
+#   \\  /    A nd           | Copyright (C) 2018-2021 OpenFOAM Foundation
 #    \\/     M anipulation  |
 #------------------------------------------------------------------------------
 # License
@@ -77,7 +77,7 @@ findFiles () {
     _prefix="$2"
     _ext="$3"
 
-    if [ -z "$_ext" -o "$_ext" = "ANY" ] ; then
+    if [ -z "$_ext" ] || [ "$_ext" = "ANY" ] ; then
         findFilesInDirs "$_dirs" "$_prefix"
     else
         findFilesInDirs "$_dirs" "${_prefix}.${_ext}"
@@ -92,7 +92,7 @@ pruneFiles () {
              -maxdepth 1 -mindepth 1 \
              -type d -name "[[:digit:]]*" \
              -o -name "dev" 2> /dev/null | \
-             grep -v $WM_PROJECT_VERSION | \
+             grep -v "$WM_PROJECT_VERSION" | \
              awk -F "/" '{print $NF}')"
     [ "$_vers" ] && \
         for _v in $_vers
@@ -141,10 +141,37 @@ listArgs () {
 cpFile () {
     _file="$1"
     _dir="$2"
+
+    [ -d "$_dir" ] || return 0
+
     echo "Copying $_file to $_dir"
     cp "$_file" "$_dir"
 }
 
+setTgt () {
+    _prefix="$1"
+    _file="$2"
+    _tgt="$3"
+
+    ! [ "$_tgt" ] && \
+        echo "$_file" | grep -Eq "(fvModel|fvConstraint)" && \
+        _tgt="constant"
+
+    ! [ "$_tgt" ] && \
+        case "$prefix" in
+            All*) _tgt="." ;;
+            *Properties|*Cloud) _tgt="constant" ;;
+            s) _tgt="0" ;;
+            *) _tgt="system" ;;
+        esac
+
+    echo "$_tgt"
+
+    [ -d "$_tgt" ] && return 0
+    echo "target directory does not exist: '$_tgt'" >&2
+    return 0
+}
+
 setFile () {
     _files="$1"
     _n="$2"
@@ -207,25 +234,8 @@ done
 [ $# -eq 1 ] || error "missing argument: no file name/prefix <file> supplied"
 prefix="$1"
 
-[ "$tgt" ] || \
-    case "$prefix" in
-       All*)
-          tgt="."
-          ;;
-       *Properties|*Cloud)
-          tgt="constant"
-          ;;
-       s)
-          tgt="0"
-          ;;
-       *)
-          tgt="system"
-          ;;
-    esac
-
 [ -s "system/controlDict" ] || \
     echo "Warning: cannot find OpenFOAM case directory (no system/controlDict file)"
-[ -d "$tgt" ] || error "target directory does not exist: '$tgt'"
 
 files="$(findFiles "$searchDirs" "$prefix" "$ext")"
 
@@ -236,7 +246,9 @@ files="$(pruneFiles $files)"
     error "no file $prefix found $(noFilesMessage "$ext")"
 
 nFiles="$(nArgs "$files")"
-[ "$nFiles" -eq 1 ] && cpFile "$files" "$tgt" && exit 0
+[ "$nFiles" -eq 1 ] && \
+    tgt="$(setTgt "$prefix" "$files" "$tgt")" && \
+    cpFile "$files" "$tgt" && exit 0
 
 echo "Multiple files with \"$prefix\" prefix found:"
 suggest="$(listArgs "$files")"
@@ -246,15 +258,16 @@ echo "$files" | grep -q "annotated/" && \
 printf "\n%s" "Enter file number (1-$nFiles) to obtain description "
 [ -n "$suggest" ] && printf "%s" "(suggest $suggest) "
 printf "%s" ": "
-read nFile
+read -r nFile
 
-[ -z "$nFile" -a -n "$suggest" ] && nFile="$suggest"
+[ -z "$nFile" ] && [ -n "$suggest" ] && nFile="$suggest"
 [ -z "$nFile" ] && \
     echo "Cannot specify nothing; re-run and enter a file number" && exit 1
 ! [ "$nFile" -eq "$nFile" ] 2>/dev/null && \
     echo "\"$nFile\" is not a number between 1 and $nFiles" && exit 1
-[ "$nFile" -lt 1 -o "$nFile" -gt "$nFiles" ] && \
+[ "$nFile" -lt 1 ] || [ "$nFile" -gt "$nFiles" ] && \
     echo "\"$nFile\" is not a number between 1 and $nFiles" && exit 1
 
 file="$(setFile "$files" "$nFile")"
+tgt="$(setTgt "$prefix" "$file" "$tgt")"
 cpFile "$file" "$tgt"

etc/caseDicts/general/fvConstraints/limitPressure

@@ -0,0 +1,29 @@
+/*--------------------------------*- 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;
+    object      limitPressure;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+limitp
+{
+    type       limitPressure;
+
+    // limits can be specified absolute
+    // min        200;
+    // max        500;
+
+    // ... or as multipliers on initial min/max values
+    minFactor  0.1;
+    maxFactor  1.5;
+}
+
+//************************************************************************** //

etc/caseDicts/general/fvConstraints/limitTemperature

@@ -0,0 +1,24 @@
+/*--------------------------------*- 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;
+    object      limitTemperature;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+limitT
+{
+    type       limitTemperature;
+    min        200;
+    max        6000;
+    selectionMode all;
+}
+
+//************************************************************************** //

etc/caseDicts/solvers/lagrangian/cloudProperties

@@ -13,6 +13,8 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+type            cloud;
+
 solution
 {
     transient       true;

etc/templates/closedVolumeRotating/constant/dynamicMeshDict

@@ -15,17 +15,14 @@ FoamFile
 
 dynamicFvMesh   dynamicMotionSolverFvMesh;
 
-solver solidBody;
+motionSolverLibs ( "libfvMotionSolvers.so" );
 
-solidBodyCoeffs
+motionSolver    solidBody;
-{
-    cellZone        rotatingZone;
 
-    solidBodyMotionFunction  rotatingMotion;
+cellZone        rotatingZone;
-    rotatingMotionCoeffs
+
-    {
+solidBodyMotionFunction  rotatingMotion;
-        #include "rotatingZoneProperties"
+
-    }
+#include        "rotatingZoneProperties"
-}
 
 // ************************************************************************* //

etc/templates/closedVolumeRotating/system/snappyHexMeshDict

@@ -87,8 +87,8 @@ castellatedMeshControls
 
 snapControls
 {
-//  explicitFeatureSnap    on;
+    explicitFeatureSnap    on;
-//  implicitFeatureSnap    off;
+    implicitFeatureSnap    off;
 }
 
 addLayersControls

etc/templates/compressibleInflowOutflow/system/fvSolution

@@ -42,8 +42,6 @@ SIMPLE
     }
 
     nNonOrthogonalCorrectors 0;
-    pMinFactor      0.1;
-    pMaxFactor      1.5;
 }
 
 relaxationFactors

etc/templates/inflowOutflowRotating/constant/dynamicMeshDict

@@ -15,17 +15,14 @@ FoamFile
 
 dynamicFvMesh   dynamicMotionSolverFvMesh;
 
-solver solidBody;
+motionSolverLibs ( "libfvMotionSolvers.so" );
 
-solidBodyCoeffs
+motionSolver    solidBody;
-{
-    cellZone        rotatingZone;
 
-    solidBodyMotionFunction  rotatingMotion;
+cellZone        rotatingZone;
-    rotatingMotionCoeffs
+
-    {
+solidBodyMotionFunction  rotatingMotion;
-        #include "rotatingZoneProperties"
+
-    }
+#include        "rotatingZoneProperties"
-}
 
 // ************************************************************************* //

etc/templates/inflowOutflowRotating/system/snappyHexMeshDict

@@ -116,8 +116,8 @@ castellatedMeshControls
 
 snapControls
 {
-//  explicitFeatureSnap    on;
+    explicitFeatureSnap    on;
-//  implicitFeatureSnap    off;
+    implicitFeatureSnap    off;
 }
 
 addLayersControls

etc/templates/singleFluidCHT/system/snappyHexMeshDict

@@ -107,8 +107,8 @@ castellatedMeshControls
 
 snapControls
 {
-//  explicitFeatureSnap    on;
+    explicitFeatureSnap    on;
-//  implicitFeatureSnap    off;
+    implicitFeatureSnap    off;
 }
 
 addLayersControls

etc/templates/singleFluidCHT/templates/materials/aluminium/radiationProperties

@@ -13,7 +13,7 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-radiationModel  opaqueSolid;
+radiationModel  none; // opaqueSolid;
 
 absorptionEmissionModel constant;
 

etc/templates/singleFluidCHT/templates/materials/copper/radiationProperties

@@ -13,7 +13,7 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-radiationModel  opaqueSolid;
+radiationModel  none; // opaqueSolid;
 
 absorptionEmissionModel constant;
 

src/OpenFOAM/db/error/IOerror.C

@@ -41,10 +41,6 @@ Foam::IOerror::IOerror(const string& title)
 {}
 
 
-Foam::IOerror::~IOerror() throw()
-{}
-
-
 Foam::OSstream& Foam::IOerror::operator()
 (
     const char* functionName,
@@ -186,7 +182,7 @@ void Foam::IOerror::exit(const int)
             IOerror errorException(*this);
 
             // Rewind the message buffer for the next error message
-            messageStreamPtr_->rewind();
+            messageStream_.rewind();
 
             throw errorException;
         }
@@ -231,7 +227,7 @@ void Foam::IOerror::abort()
             IOerror errorException(*this);
 
             // Rewind the message buffer for the next error message
-            messageStreamPtr_->rewind();
+            messageStream_.rewind();
 
             throw errorException;
         }

src/OpenFOAM/db/error/error.C

@@ -42,9 +42,9 @@ Foam::error::error(const string& title)
     sourceFileLineNumber_(0),
     abort_(env("FOAM_ABORT")),
     throwExceptions_(false),
-    messageStreamPtr_(new OStringStream())
+    messageStream_()
 {
-    if (!messageStreamPtr_->good())
+    if (!messageStream_.good())
     {
         Perr<< endl
             << "error::error(const string& title) : cannot open error stream"
@@ -54,12 +54,6 @@ Foam::error::error(const string& title)
 }
 
 
-Foam::error::~error() throw()
-{
-    delete messageStreamPtr_;
-}
-
-
 Foam::OSstream& Foam::error::operator()
 (
     const char* functionName,
@@ -93,7 +87,7 @@ Foam::OSstream& Foam::error::operator()
 
 Foam::OSstream& Foam::error::operator()()
 {
-    if (!messageStreamPtr_->good())
+    if (!messageStream_.good())
     {
         Perr<< endl
             << "error::operator OSstream&() : error stream has failed"
@@ -101,7 +95,7 @@ Foam::OSstream& Foam::error::operator()()
         abort();
     }
 
-    return *messageStreamPtr_;
+    return messageStream_;
 }
 
 
@@ -124,7 +118,7 @@ Foam::error::operator Foam::dictionary() const
 
 Foam::string Foam::error::message() const
 {
-    return messageStreamPtr_->str();
+    return messageStream_.str();
 }
 
 
@@ -155,7 +149,7 @@ void Foam::error::exit(const int errNo)
             error errorException(*this);
 
             // Rewind the message buffer for the next error message
-            messageStreamPtr_->rewind();
+            messageStream_.rewind();
 
             throw errorException;
         }
@@ -200,7 +194,7 @@ void Foam::error::abort()
             error errorException(*this);
 
             // Rewind the message buffer for the next error message
-            messageStreamPtr_->rewind();
+            messageStream_.rewind();
 
             throw errorException;
         }

src/OpenFOAM/db/error/error.H

@@ -48,6 +48,7 @@ SourceFiles
 #ifndef error_H
 #define error_H
 
+#include "OStringStream.H"
 #include "messageStream.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -81,7 +82,9 @@ protected:
         bool abort_;
 
         bool throwExceptions_;
-        OStringStream* messageStreamPtr_;
+
+        OStringStream messageStream_;
+
 
 public:
 
@@ -91,10 +94,6 @@ public:
         error(const string& title);
 
 
-    //- Destructor
-    virtual ~error() throw();
-
-
     // Member Functions
 
         string message() const;
@@ -205,10 +204,6 @@ public:
         IOerror(const string& title);
 
 
-    //- Destructor
-    virtual ~IOerror() throw();
-
-
     // Member Functions
 
         const string& ioFileName() const

src/functionObjects/lagrangian/particles/particles.H

@@ -25,8 +25,8 @@ Class
     Foam::functionObjects::particles
 
 Description
-    This functionObject tracks a uncoupled kinematic particle cloud in the
+    This functionObject tracks a particle cloud in the specified velocity field
-    specified velocity field of an incompressible flow (laminar, RANS or LES).
+    of an incompressible flow (laminar, RANS or LES).
 
     It may be used in conjunction with any transient single-phase incompressible
     flow solver such as pisoFoam or pimpleFoam and tracks the particles or

src/lagrangian/parcel/fvModels/clouds.C

@@ -57,16 +57,87 @@ Foam::fv::clouds::clouds
 )
 :
     fvModel(sourceName, modelType, dict, mesh),
-    carrierThermo_
+    g_
+    (
+        IOobject
+        (
+            "g",
+            mesh.time().constant(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::NO_WRITE
+        ),
+        dimensionedVector(dimAcceleration, Zero)
+    ),
+    carrierHasThermo_
+    (
+        mesh.foundObject<fluidThermo>(physicalProperties::typeName)
+    ),
+    tCarrierThermo_
+    (
+        carrierHasThermo_
+      ? tmpNrc<fluidThermo>
         (
             mesh.lookupObject<fluidThermo>(physicalProperties::typeName)
+        )
+      : tmpNrc<fluidThermo>(nullptr)
     ),
-    clouds_
+    tCarrierViscosity_
     (
-        mesh.lookupObject<volScalarField>("rho"),
+        carrierHasThermo_
-        mesh.lookupObject<volVectorField>("U"),
+      ? tmpNrc<viscosityModel>(nullptr)
-        mesh.lookupObject<uniformDimensionedVectorField>("g"),
+      : tmpNrc<viscosityModel>
-        carrierThermo_
+        (
+            mesh.lookupObject<viscosityModel>(physicalProperties::typeName)
+        )
+    ),
+    tRho_
+    (
+        carrierHasThermo_
+      ? tmp<volScalarField>(nullptr)
+      : tmp<volScalarField>
+        (
+            new volScalarField
+            (
+                IOobject
+                (
+                    "rho",
+                    mesh.time().timeName(),
+                    mesh
+                ),
+                mesh,
+                dimensionedScalar("rho", dimDensity, tCarrierViscosity_())
+            )
+        )
+    ),
+    tMu_
+    (
+        carrierHasThermo_
+      ? tmp<volScalarField>(nullptr)
+      : tmp<volScalarField>
+        (
+            new volScalarField("mu", tRho_()*tCarrierViscosity_().nu())
+        )
+    ),
+    rhoName_(dict.lookupOrDefault<word>("rho", "rho")),
+    UName_(dict.lookupOrDefault<word>("U", "U")),
+    cloudsPtr_
+    (
+        carrierHasThermo_
+      ? new parcelCloudList
+        (
+            mesh.lookupObject<volScalarField>(rhoName_),
+            mesh.lookupObject<volVectorField>(UName_),
+            g_,
+            tCarrierThermo_()
+        )
+      : new parcelCloudList
+        (
+            tRho_(),
+            mesh.lookupObject<volVectorField>(UName_),
+            tMu_(),
+            g_
+        )
     ),
     curTimeIndex_(-1)
 {}
@@ -76,12 +147,20 @@ Foam::fv::clouds::clouds
 
 Foam::wordList Foam::fv::clouds::addSupFields() const
 {
-    wordList fieldNames({"rho", "U", carrierThermo_.he().name()});
+    wordList fieldNames(1, UName_);
 
-    if (isA<basicSpecieMixture>(carrierThermo_))
+    if (carrierHasThermo_)
+    {
+        const fluidThermo& carrierThermo = tCarrierThermo_();
+
+        fieldNames.append(rhoName_);
+
+        fieldNames.append(carrierThermo.he().name());
+
+        if (isA<basicSpecieMixture>(carrierThermo))
         {
             const basicSpecieMixture& composition =
-            refCast<const basicSpecieMixture>(carrierThermo_);
+                refCast<const basicSpecieMixture>(carrierThermo);
 
             const PtrList<volScalarField>& Y = composition.Y();
 
@@ -93,6 +172,7 @@ Foam::wordList Foam::fv::clouds::addSupFields() const
                 }
             }
         }
+    }
 
     return fieldNames;
 }
@@ -105,7 +185,12 @@ void Foam::fv::clouds::correct()
         return;
     }
 
-    clouds_.evolve();
+    if (!carrierHasThermo_)
+    {
+        tMu_.ref() = tRho_()*tCarrierViscosity_().nu();
+    }
+
+    cloudsPtr_().evolve();
 
     curTimeIndex_ = mesh().time().timeIndex();
 }
@@ -122,16 +207,18 @@ void Foam::fv::clouds::addSup
         Info<< type() << ": applying source to " << eqn.psi().name() << endl;
     }
 
-    if (fieldName == "rho")
+    if (!carrierHasThermo_)
-    {
-        eqn += clouds_.Srho(eqn.psi());
-    }
-    else
     {
         FatalErrorInFunction
-            << "Support for field " << fieldName << " is not implemented"
+            << "Applying source to compressible equation when carrier thermo "
+            << "is not available"
             << exit(FatalError);
     }
+
+    if (fieldName == rhoName_)
+    {
+        eqn += cloudsPtr_().Srho(eqn.psi());
+    }
 }
 
 
@@ -147,35 +234,65 @@ void Foam::fv::clouds::addSup
         Info<< type() << ": applying source to " << eqn.psi().name() << endl;
     }
 
-    if (fieldName == "rho")
+    if (!carrierHasThermo_)
     {
-        eqn += clouds_.Srho(eqn.psi());
+        FatalErrorInFunction
+            << "Applying source to compressible equation when carrier thermo "
+            << "is not available"
+            << exit(FatalError);
     }
-    else if (fieldName == carrierThermo_.he().name())
+
+    const fluidThermo& carrierThermo = tCarrierThermo_();
+
+    if (fieldName == rhoName_)
     {
-        eqn += clouds_.Sh(eqn.psi());
+        eqn += cloudsPtr_().Srho(eqn.psi());
+    }
+    else if (fieldName == carrierThermo.he().name())
+    {
+        eqn += cloudsPtr_().Sh(eqn.psi());
     }
     else if
     (
-        isA<basicSpecieMixture>(carrierThermo_)
+        isA<basicSpecieMixture>(carrierThermo)
-     && refCast<const basicSpecieMixture>(carrierThermo_).contains
+     && refCast<const basicSpecieMixture>(carrierThermo).contains
         (
             eqn.psi().name()
         )
     )
     {
-        eqn += clouds_.SYi
+        eqn += cloudsPtr_().SYi
         (
-            refCast<const basicSpecieMixture>(carrierThermo_).index(eqn.psi()),
+            refCast<const basicSpecieMixture>(carrierThermo).index(eqn.psi()),
             eqn.psi()
         );
     }
-    else
+}
+
+
+void Foam::fv::clouds::addSup
+(
+    fvMatrix<vector>& eqn,
+    const word& fieldName
+) const
+{
+    if (debug)
+    {
+        Info<< type() << ": applying source to " << eqn.psi().name() << endl;
+    }
+
+    if (carrierHasThermo_)
     {
         FatalErrorInFunction
-            << "Support for field " << fieldName << " is not implemented"
+            << "Applying source to incompressible equation when carrier thermo "
+            << "is available"
             << exit(FatalError);
     }
+
+    if (fieldName == UName_)
+    {
+        eqn += cloudsPtr_().SU(eqn.psi())/tRho_();
+    }
 }
 
 
@@ -191,23 +308,25 @@ void Foam::fv::clouds::addSup
         Info<< type() << ": applying source to " << eqn.psi().name() << endl;
     }
 
-    if (fieldName == "U")
+    if (!carrierHasThermo_)
-    {
-        eqn += clouds_.SU(eqn.psi());
-    }
-    else
     {
         FatalErrorInFunction
-            << "Support for field " << fieldName << " is not implemented"
+            << "Applying source to compressible equation when carrier thermo "
+            << "is not available"
             << exit(FatalError);
     }
+
+    if (fieldName == UName_)
+    {
+        eqn += cloudsPtr_().SU(eqn.psi());
+    }
 }
 
 
 void Foam::fv::clouds::preUpdateMesh()
 {
     // Store the particle positions
-    clouds_.storeGlobalPositions();
+    cloudsPtr_().storeGlobalPositions();
 }
 
 

src/lagrangian/parcel/fvModels/clouds.H

@@ -25,17 +25,46 @@ Class
     Foam::fv::clouds
 
 Description
-    Lagrangian clouds fvModel
+    This fvModel adds any number of Lagrangian clouds to any single-phase
+    solver. The particles are tracked through, and exchange sources with, the
+    Eulerian flow field.
+
+    As well as the fvModel controls, properties must be specified for each
+    cloud. For a single cloud, these should be provided in the
+    constant/cloudProperties file. For multiple clouds, the list of cloud names
+    must first be provided in the constant/clouds file. Then, each named cloud
+    has its properties specified in its constant/<cloudName>Properties file.
+
+    The application of sources to the Eulerian fields is controlled by the
+    solution/coupled switch in each cloud's properties file. If set to "true"
+    then the Eulerian phase will have forces, and heat and mass sources applied
+    to it by the Lagrangian phase. If set to "false" then these will be omitted,
+    and the Lagrangian phase will not affect the Eulerian phase.
+
+    If this model is used with an incompressible solver, then the density of
+    the Eulerian phase must be specified in the constant/physicalProperties
+    dictionary.
+
+    Gravity will be read from the constant/g file if present, otherwise it will
+    default to zero.
 
 Usage
     Example usage:
     \verbatim
     clouds
     {
+        libs        ("liblagrangianParcel.so");
         type        clouds;
     }
     \endverbatim
 
+    \table
+        Property | Description                     | Required   | Default value
+        type     | Type name: clouds               | yes        |
+        rho      | Name of the density field       | no         | rho
+        U        | Name of the velocity field      | no         | U
+    \endtable
+
 SourceFiles
     clouds.C
 
@@ -46,6 +75,7 @@ SourceFiles
 
 #include "fvModel.H"
 #include "fluidThermo.H"
+#include "viscosityModel.H"
 #include "uniformDimensionedFields.H"
 #include "parcelCloudList.H"
 
@@ -66,11 +96,35 @@ class clouds
 {
     // Private Data
 
-        //- Reference to the carrier phase thermo
+        //- Optional acceleration due to gravity
-        const fluidThermo& carrierThermo_;
+        const uniformDimensionedVectorField g_;
+
+        //- Flag to indicate whether a thermo model is available for the
+        //  carrier
+        const bool carrierHasThermo_;
+
+        //- Reference to the carrier phase thermo. Valid only if the carrier
+        //  has thermo.
+        const tmpNrc<fluidThermo> tCarrierThermo_;
+
+        //- Reference to the carrier viscosity model. Valid only if the carrier
+        //  does not have thermo.
+        const tmpNrc<viscosityModel> tCarrierViscosity_;
+
+        //- Density field. Valid only if the carrier does not have thermo.
+        const tmp<volScalarField> tRho_;
+
+        //- Viscosity field. Valid only if the carrier does not have thermo.
+        tmp<volScalarField> tMu_;
+
+        //- Name of the density field
+        const word rhoName_;
+
+        //- Name of the velocity field
+        const word UName_;
 
         //- The Lagrangian cloud list
-        mutable parcelCloudList clouds_;
+        mutable autoPtr<parcelCloudList> cloudsPtr_;
 
         //- Current time index (used for updating)
         mutable label curTimeIndex_;
@@ -127,7 +181,7 @@ public:
                 const word& fieldName
             ) const;
 
-            //- Add source to pressure or enthalpy equation
+            //- Add source to enthalpy or species equation
             virtual void addSup
             (
                 const volScalarField& rho,
@@ -135,7 +189,14 @@ public:
                 const word& fieldName
             ) const;
 
-            //- Add source to momentum equation
+            //- Add source to incompressible momentum equation
+            virtual void addSup
+            (
+                fvMatrix<vector>& eqn,
+                const word& fieldName
+            ) const;
+
+            //- Add source to compressible momentum equation
             virtual void addSup
             (
                 const volScalarField& rho,

src/lagrangian/parcel/parcels/Templates/MomentumParcel/MomentumParcel.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2020 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2021 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -313,7 +313,7 @@ bool Foam::MomentumParcel<ParcelType>::move
         // maxCo times the total value.
         scalar f = 1 - p.stepFraction();
         f = min(f, maxCo);
-        f = min(f, maxCo*l/max(small*l, mag(s)));
+        f = min(f, maxCo/min(max(mag(s)/l, rootSmall), rootGreat));
         if (p.moving())
         {
             // Track to the next face
@@ -421,7 +421,13 @@ void Foam::MomentumParcel<ParcelType>::hitWallPatch
     trackingData&
 )
 {
-    // wall interactions are handled by the generic hitPatch method
+    const polyPatch& pp = this->mesh().boundaryMesh()[this->patch()];
+
+    FatalErrorInFunction
+        << "Particle " << this->origId() << " hit " << pp.type() << " patch "
+        << pp.name() << " at " << this->position()
+        << " but no interaction model was specified for this patch"
+        << exit(FatalError);
 }
 
 

tutorials/incompressible/pimpleFoam/RAS/TJunction/constant/cloudProperties

@@ -0,0 +1,91 @@
+/*--------------------------------*- 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";
+    object      cloudProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+type        cloud;
+
+solution
+{
+    coupled         false;
+    transient       yes;
+    cellValueSourceCorrection off;
+    maxCo           0.3;
+
+    interpolationSchemes
+    {
+        rho             cell;
+        U               cellPoint;
+        mu              cell;
+    }
+
+    integrationSchemes
+    {
+        U               Euler;
+    }
+}
+
+constantProperties
+{
+    rho0            964;
+}
+
+subModels
+{
+    particleForces
+    {
+        sphereDrag;
+    }
+
+    injectionModels
+    {
+        model1
+        {
+            type            patchInjection;
+            massTotal       0;
+            SOI             0;
+            duration        5;
+            parcelBasisType fixed;
+            nParticle       1;
+            parcelsPerSecond 1000;
+            patchName       inlet;
+            U0              (1 0 0);
+            flowRateProfile constant 1;
+            sizeDistribution
+            {
+                type        fixedValue;
+                fixedValueDistribution
+                {
+                    value       1e-4;
+                }
+            }
+        }
+    }
+
+    dispersionModel none;
+
+    patchInteractionModel standardWallInteraction;
+
+    standardWallInteractionCoeffs
+    {
+        type rebound;
+    }
+
+    surfaceFilmModel none;
+
+    stochasticCollisionModel none;
+}
+
+
+// ************************************************************************* //

tutorials/incompressible/pimpleFoam/RAS/TJunction/constant/fvModels

@@ -0,0 +1,23 @@
+/*--------------------------------*- 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";
+    object      fvModels;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+clouds
+{
+    type    clouds;
+    libs    ("liblagrangianParcel.so");
+}
+
+// ************************************************************************* //

tutorials/incompressible/pimpleFoam/RAS/TJunction/constant/physicalProperties

@@ -18,4 +18,6 @@ viscosityModel  constant;
 
 nu              [0 2 -1 0 0 0 0] 1e-05;
 
+rho             [1 -3 0 0 0 0 0] 1.2;
+
 // ************************************************************************* //