From dff4c3da05e04f548a2746a1fdbd961511537e66 Mon Sep 17 00:00:00 2001
From: Andrew Heather <>
Date: Tue, 29 Nov 2022 09:18:57 +0000
Subject: [PATCH 1/2] ENH: limitHeight: new faOption to limit film height

---
 src/faOptions/Make/files                      |   1 +
 .../corrections/limitHeight/limitHeight.C     | 152 +++++++++++++++++
 .../corrections/limitHeight/limitHeight.H     | 153 ++++++++++++++++++
 3 files changed, 306 insertions(+)
 create mode 100644 src/faOptions/corrections/limitHeight/limitHeight.C
 create mode 100644 src/faOptions/corrections/limitHeight/limitHeight.H

diff --git a/src/faOptions/Make/files b/src/faOptions/Make/files
index af26eff51b..43e1de42e4 100644
--- a/src/faOptions/Make/files
+++ b/src/faOptions/Make/files
@@ -17,5 +17,6 @@ $(derivedSources)/externalFileSource/externalFileSource.C
 corrections=corrections
 
 $(corrections)/limitVelocity/limitVelocity.C
+$(corrections)/limitHeight/limitHeight.C
 
 LIB = $(FOAM_LIBBIN)/libfaOptions
diff --git a/src/faOptions/corrections/limitHeight/limitHeight.C b/src/faOptions/corrections/limitHeight/limitHeight.C
new file mode 100644
index 0000000000..89e1787c16
--- /dev/null
+++ b/src/faOptions/corrections/limitHeight/limitHeight.C
@@ -0,0 +1,152 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2022 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "limitHeight.H"
+#include "areaFields.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace fa
+{
+    defineTypeNameAndDebug(limitHeight, 0);
+    addToRunTimeSelectionTable
+    (
+        option,
+        limitHeight,
+        dictionary
+    );
+}
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::fa::limitHeight::limitHeight
+(
+    const word& name,
+    const word& modelType,
+    const dictionary& dict,
+    const fvMesh& mesh
+)
+:
+    faceSetOption(name, modelType, dict, mesh),
+    hName_("h"),
+    max_(0)  // overwritten later
+{
+    read(dict);
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+bool Foam::fa::limitHeight::read(const dictionary& dict)
+{
+    if (!faceSetOption::read(dict))
+    {
+        return false;
+    }
+
+    coeffs_.readIfPresent("h", hName_);
+    coeffs_.readEntry("max", max_);
+
+    fieldNames_.resize(1, hName_);
+
+    applied_.resize(1, false);
+
+    return true;
+}
+
+
+void Foam::fa::limitHeight::correct(areaScalarField& h)
+{
+    // Count nTotFaces ourselves
+    // (maybe only applying on a subset)
+    label nFacesAbove = 0;
+    const label nTotFaces = returnReduce(faces_.size(), sumOp<label>());
+
+    scalarField& hif = h.primitiveFieldRef();
+
+    for (const label facei : faces_)
+    {
+        auto& hval = hif[facei];
+
+        if (hval > max_)
+        {
+            hval *= max_/max(hval, SMALL);
+            ++nFacesAbove;
+        }
+    }
+
+    // Handle boundaries in the case of 'all'
+    label nEdgesAbove = 0;
+
+    if (!faceSetOption::useSubMesh())
+    {
+        for (faPatchScalarField& hp : h.boundaryFieldRef())
+        {
+            if (!hp.fixesValue())
+            {
+                for (auto& hval : hp)
+                {
+                    if (hval > max_)
+                    {
+                        hval *= max_/max(hval, SMALL);
+                        ++nEdgesAbove;
+                    }
+                }
+            }
+        }
+    }
+
+    // Percent, max 2 decimal places
+    const auto percent = [](scalar num, label denom) -> scalar
+    {
+        return (denom ? 1e-2*round(1e4*num/denom) : 0);
+    };
+
+
+    reduce(nFacesAbove, sumOp<label>());
+    reduce(nEdgesAbove, sumOp<label>());
+
+    Info<< type() << ' ' << name_ << " Limited "
+        << nFacesAbove << " ("
+        << percent(nFacesAbove, nTotFaces)
+        << "%) of faces, with max limit " << max_ << endl;
+
+    if (nFacesAbove || nEdgesAbove)
+    {
+        // We've changed internal values so give
+        // boundary conditions opportunity to correct
+        h.correctBoundaryConditions();
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/src/faOptions/corrections/limitHeight/limitHeight.H b/src/faOptions/corrections/limitHeight/limitHeight.H
new file mode 100644
index 0000000000..d9872ff677
--- /dev/null
+++ b/src/faOptions/corrections/limitHeight/limitHeight.H
@@ -0,0 +1,153 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2022 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::fa::limitHeight
+
+Description
+    Limits the film height to a specified \c max value.
+
+    The \c limitHeight applies an explicit correction to limit the max film
+    height to the user-specified value. This will ensure that the max value
+    is adhered to, but may become unstable if the limiting instigates an
+    oscillatory/destabilising effect into neighbouring faces that then spreads.
+
+Usage
+    Minimal example by using \c constant/faOptions:
+    \verbatim
+    <faOption>
+    {
+        // Mandatory entries
+        type            limitHeight;
+        active          yes;
+        selectionMode   all;
+        max             <scalar>;
+
+        // Optional entries
+        h               <word>;
+
+        // Inherited entries
+        ...
+    }
+    \endverbatim
+
+    where the entries mean:
+    \table
+      Property  | Description                   | Type   | Reqd  | Deflt
+      type      | Type name: limitHeight        | word   | yes   | -
+      max       | Maximum height limit [m]      | scalar | yes   | -
+      h         | Name of operand height field  | word   | no    | h
+    \endtable
+
+    The inherited entries are elaborated in:
+      - \link faceSetOption.H \endlink
+
+SourceFiles
+    limitHeight.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef fa_limitHeight_H
+#define fa_limitHeight_H
+
+#include "faceSetOption.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace fa
+{
+
+/*---------------------------------------------------------------------------*\
+                        Class limitHeight Declaration
+\*---------------------------------------------------------------------------*/
+
+class limitHeight
+:
+    public faceSetOption
+{
+protected:
+
+    // Protected Data
+
+        //- Name of operand height field
+        word hName_;
+
+        //- Maximum height [m]
+        scalar max_;
+
+
+    // Protected Member Functions
+
+        //- No copy construct
+        limitHeight(const limitHeight&) = delete;
+
+        //- No copy assignment
+        void operator=(const limitHeight&) = delete;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("limitHeight");
+
+
+    // Constructors
+
+        //- Construct from components
+        limitHeight
+        (
+            const word& name,
+            const word& modelType,
+            const dictionary& dict,
+            const fvMesh& mesh
+        );
+
+
+    //- Destructor
+    virtual ~limitHeight() = default;
+
+
+    // Member Functions
+
+        //- Read dictionary
+        virtual bool read(const dictionary& dict);
+
+        //- Correct the height field
+        virtual void correct(areaScalarField& h);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace fv
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

From 59811983db969d03b15746e0922bac27d0b49c51 Mon Sep 17 00:00:00 2001
From: Andrew Heather <>
Date: Tue, 29 Nov 2022 10:45:10 +0000
Subject: [PATCH 2/2] ENH: BrunDrippingInjection: new injection model for
 liquid films

---
 src/regionFaModels/Make/files                 |   1 +
 .../BrunDrippingInjection.C                   | 147 ++++++++++++++
 .../BrunDrippingInjection.H                   | 185 ++++++++++++++++++
 3 files changed, 333 insertions(+)
 create mode 100644 src/regionFaModels/liquidFilm/subModels/kinematic/injectionModel/BrunDrippingInjection/BrunDrippingInjection.C
 create mode 100644 src/regionFaModels/liquidFilm/subModels/kinematic/injectionModel/BrunDrippingInjection/BrunDrippingInjection.H

diff --git a/src/regionFaModels/Make/files b/src/regionFaModels/Make/files
index a4f1d35563..4e761c227d 100644
--- a/src/regionFaModels/Make/files
+++ b/src/regionFaModels/Make/files
@@ -24,6 +24,7 @@ liquidFilm/subModels/kinematic/injectionModel/injectionModel/injectionModel.C
 liquidFilm/subModels/kinematic/injectionModel/injectionModel/injectionModelNew.C
 
 liquidFilm/subModels/kinematic/injectionModel/curvatureSeparation/curvatureSeparation.C
+liquidFilm/subModels/kinematic/injectionModel/BrunDrippingInjection/BrunDrippingInjection.C
 
 liquidFilm/subModels/kinematic/force/forceList/forceList.C
 liquidFilm/subModels/kinematic/force/force/force.C
diff --git a/src/regionFaModels/liquidFilm/subModels/kinematic/injectionModel/BrunDrippingInjection/BrunDrippingInjection.C b/src/regionFaModels/liquidFilm/subModels/kinematic/injectionModel/BrunDrippingInjection/BrunDrippingInjection.C
new file mode 100644
index 0000000000..d0e4417b0e
--- /dev/null
+++ b/src/regionFaModels/liquidFilm/subModels/kinematic/injectionModel/BrunDrippingInjection/BrunDrippingInjection.C
@@ -0,0 +1,147 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2016-2017 OpenFOAM Foundation
+    Copyright (C) 2020-2022 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "BrunDrippingInjection.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace regionModels
+{
+namespace areaSurfaceFilmModels
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(BrunDrippingInjection, 0);
+addToRunTimeSelectionTable
+(
+    injectionModel,
+    BrunDrippingInjection,
+    dictionary
+);
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+BrunDrippingInjection::BrunDrippingInjection
+(
+    liquidFilmBase& film,
+    const dictionary& dict
+)
+:
+    injectionModel(type(), film, dict),
+    ubarStar_
+    (
+        coeffDict_.getCheckOrDefault<scalar>
+        (
+            "ubarStar",
+            1.62208,
+            scalarMinMax::ge(SMALL)
+        )
+    ),
+    dCoeff_(coeffDict_.getOrDefault<scalar>("dCoeff", 3.3)),
+    deltaStable_(coeffDict_.getOrDefault<scalar>("deltaStable", 0)),
+    diameter_(film.regionMesh().nFaces(), -1.0)
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+void BrunDrippingInjection::correct
+(
+    scalarField& availableMass,
+    scalarField& massToInject,
+    scalarField& diameterToInject
+)
+{
+    // Calculate available dripping mass
+    tmp<areaScalarField> tsinAlpha = -film().gn()/mag(film().g());
+    const scalarField& sinAlpha = tsinAlpha();
+
+    const areaScalarField& delta = film().h();
+    const areaScalarField& rho = film().rho();
+    const areaScalarField& sigma = film().sigma();
+    const scalar magg = mag(film().g().value());
+
+    forAll(delta, facei)
+    {
+        bool dripping = false;
+
+        if (sinAlpha[facei] > SMALL && delta[facei] > deltaStable_)
+        {
+            const scalar rhoc = rho[facei];
+            const scalar lc = sqrt(sigma[facei]/(rhoc*magg));
+            const scalar deltaStable = max
+            (
+                3*lc*sqrt(1 - sqr(sinAlpha[facei]))
+               /(ubarStar_*sqrt(sinAlpha[facei])*sinAlpha[facei]),
+                deltaStable_
+            );
+
+            if (delta[facei] > deltaStable)
+            {
+                const scalar massDrip =
+                    availableMass[facei]*(delta[facei] - deltaStable);
+
+                if (massDrip > 0)
+                {
+                    const scalar diam = dCoeff_*lc;
+                    diameter_[facei] = diam;
+
+                    massToInject[facei] += massDrip;
+                    availableMass[facei] -= massDrip;
+
+                    diameterToInject[facei] = diam;
+                    addToInjectedMass(massDrip);
+
+                    dripping = true;
+                }
+            }
+        }
+
+        if (!dripping)
+        {
+            diameterToInject[facei] = 0;
+            massToInject[facei] = 0;
+        }
+    }
+
+    injectionModel::correct();
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace areaSurfaceFilmModels
+} // End namespace regionModels
+} // End namespace Foam
+
+// ************************************************************************* //
diff --git a/src/regionFaModels/liquidFilm/subModels/kinematic/injectionModel/BrunDrippingInjection/BrunDrippingInjection.H b/src/regionFaModels/liquidFilm/subModels/kinematic/injectionModel/BrunDrippingInjection/BrunDrippingInjection.H
new file mode 100644
index 0000000000..b94ac6b080
--- /dev/null
+++ b/src/regionFaModels/liquidFilm/subModels/kinematic/injectionModel/BrunDrippingInjection/BrunDrippingInjection.H
@@ -0,0 +1,185 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2016-2017 OpenFOAM Foundation
+    Copyright (C) 2022 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::regionModels::areaSurfaceFilmModels::BrunDrippingInjection
+
+Description
+    Film-dripping mass-transfer model.
+
+    If the film thickness exceeds the critical value needed to generate one or
+    more drops, the equivalent mass is removed from the film.  The critical film
+    thickness is calculated from the Rayleigh-Taylor stability analysis of film
+    flow on an inclined plane, proposed by (Brun et. al., 2015).
+
+    Reference:
+    \verbatim
+        Brun, P. T., Damiano, A., Rieu, P., Balestra, G., & Gallaire, F. (2015).
+        Rayleigh-Taylor instability under an inclined plane.
+        Physics of Fluids, 27(8), 084107.
+        DOI:10.1063/1.4927857
+    \endverbatim
+
+    The diameter of the drops formed are obtained from the local
+    capillary length multiplied by the \c dCoeff coefficient.
+
+    Reference:
+    \verbatim
+        Lefebvre, A. (1988).
+        Atomization and sprays
+        (Vol. 1040, No. 2756). CRC press.
+    \endverbatim
+
+Usage
+    Minimal example:
+    \verbatim
+    injectionModels
+    (
+        BrunDrippingInjection
+    );
+
+    BrunDrippingInjectionCoeffs
+    {
+        // Optional entries
+        ubarStar        <scalar>;
+        dCoeff          <scalar>;
+        deltaStable     <scalar>;
+    }
+    \endverbatim
+
+    where the entries mean:
+    \table
+      Property  | Description                      | Type   | Reqd  | Deflt
+      type      | Type name: BrunDrippingInjection | word   | yes   | -
+      ubarStar  | Critical non-dimensional interface velocity <!--
+                -->                                | scalar | yes   | 1.62208
+      dCoeff    | Coefficient relating the diameter of the drops <!--
+                --> formed to the capillary length | scalar | yes   | 3.3
+      deltaStable  | Stable film thickness         | scalar | yes   | 0
+    \endtable
+
+    The inherited entries are elaborated in:
+      - \link injectionModel.H \endlink
+
+SourceFiles
+    BrunDrippingInjection.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef areaSurfaceFilmModels_BrunDrippingInjection_H
+#define areaSurfaceFilmModels_BrunDrippingInjection_H
+
+#include "injectionModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace regionModels
+{
+namespace areaSurfaceFilmModels
+{
+
+/*---------------------------------------------------------------------------*\
+                    Class BrunDrippingInjection Declaration
+\*---------------------------------------------------------------------------*/
+
+class BrunDrippingInjection
+:
+    public injectionModel
+{
+    // Private Member Functions
+
+        //- No copy construct
+        BrunDrippingInjection(const BrunDrippingInjection&) = delete;
+
+        //- No copy assignment
+        void operator=(const BrunDrippingInjection&) = delete;
+
+
+protected:
+
+    // Protected Data
+
+        //- Critical non-dimensional interface velocity
+        //  Coefficient in the film angle stability function.
+        scalar ubarStar_;
+
+        //- Coefficient relating the diameter of the drops formed to
+        //- the capillary length.
+        scalar dCoeff_;
+
+        //- Stable film thickness - drips only formed if thickness
+        //- exceeds this threshold value
+        scalar deltaStable_;
+
+        //- Diameters of particles to inject into the dripping
+        scalarList diameter_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("BrunDrippingInjection");
+
+
+    // Constructors
+
+        //- Construct from surface film model
+        BrunDrippingInjection
+        (
+            liquidFilmBase& film,
+            const dictionary& dict
+        );
+
+
+    //- Destructor
+    virtual ~BrunDrippingInjection() = default;
+
+
+    // Member Functions
+
+        //- Correct
+        virtual void correct
+        (
+            scalarField& availableMass,
+            scalarField& massToInject,
+            scalarField& diameterToInject
+        );
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace areaSurfaceFilmModels
+} // End namespace regionModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //