ENH: BilgerMixtureFraction: add new reactionThermo FO

Signed-off-by: Sergio Ferraris <s.ferraris@opencfd.co.uk>
Signed-off-by: Kutalmis Bercin <kutalmis.bercin@esi-group.com>

CONTRIBUTORS.md

@@ -50,3 +50,7 @@ It is likely incomplete...
 - Norbert Weber
 - Henry Weller
 - Niklas Wikstrom
+- Thorsten Zirwes
+
+
+<!----------------------------------------------------------------------------->

src/thermophysicalModels/chemistryModel/Make/files

@@ -7,5 +7,6 @@ chemistryModel/TDACChemistryModel/tabulation/makeChemistryTabulationMethods.C
 chemistrySolver/chemistrySolver/makeChemistrySolvers.C
 
 functionObjects/specieReactionRates/specieReactionRates.C
+functionObjects/BilgerMixtureFraction/BilgerMixtureFraction.C
 
 LIB = $(FOAM_LIBBIN)/libchemistryModel

src/thermophysicalModels/chemistryModel/functionObjects/BilgerMixtureFraction/BilgerMixtureFraction.C

@@ -0,0 +1,409 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2020 Thorsten Zirwes
+    Copyright (C) 2020 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 "BilgerMixtureFraction.H"
+#include "basicThermo.H"
+#include "reactingMixture.H"
+#include "thermoPhysicsTypes.H"
+#include "scalarRange.H"
+#include "basicChemistryModel.H"
+#include "psiReactionThermo.H"
+#include "rhoReactionThermo.H"
+#include "BasicChemistryModel.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace functionObjects
+{
+    defineTypeNameAndDebug(BilgerMixtureFraction, 0);
+    addToRunTimeSelectionTable
+    (
+        functionObject,
+        BilgerMixtureFraction,
+        dictionary
+    );
+}
+}
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::functionObjects::BilgerMixtureFraction::calcBilgerMixtureFraction()
+{
+    if (!mesh_.foundObject<volScalarField>(resultName_, false))
+    {
+        auto tCo = tmp<volScalarField>::New
+        (
+            IOobject
+            (
+                resultName_,
+                mesh_.time().timeName(),
+                mesh_,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh_,
+            dimensionedScalar(dimless, Zero)
+        );
+        mesh_.objectRegistry::store(tCo.ptr());
+    }
+
+    auto& f_Bilger = mesh_.lookupObjectRef<volScalarField>(resultName_);
+
+    auto& Y = thermo_.Y();
+
+    f_Bilger = -o2RequiredOx_;
+    forAll(Y, i)
+    {
+        f_Bilger +=
+            Y[i]
+           *(nAtomsC_[i] + nAtomsS_[i] + 0.25*nAtomsH_[i] - 0.5*nAtomsO_[i])
+           /thermo_.W(i);
+    }
+    f_Bilger /= o2RequiredFuelOx_;
+    f_Bilger.clip
+    (
+        dimensionedScalar(dimless, 0),
+        dimensionedScalar(dimless, 1)
+    );
+}
+
+
+bool Foam::functionObjects::BilgerMixtureFraction::readComposition
+(
+    const dictionary& subDict,
+    scalarField& comp
+)
+{
+    auto& Y = thermo_.Y();
+    const speciesTable& speciesTab = thermo_.species();
+
+    // Read mass fractions of all species for the oxidiser or fuel
+    forAll(Y, i)
+    {
+        comp[i] =
+            subDict.getCheckOrDefault<scalar>
+            (
+                speciesTab[i],
+                0,
+                scalarRange::ge0()
+            );
+    }
+
+    if (sum(comp) < SMALL)
+    {
+        FatalIOErrorInFunction(subDict)
+            << "No composition is given" << nl
+            << "Valid species are:" << nl
+            << speciesTab
+            << exit(FatalIOError);
+
+        return false;
+    }
+
+    const word fractionBasisType
+    (
+        subDict.getOrDefault<word>("fractionBasis", "mass")
+    );
+
+    if (fractionBasisType == "mass")
+    {
+        // Normalize fractionBasis to the unity
+        comp /= sum(comp);
+    }
+    else if (fractionBasisType == "mole")
+    {
+        // In case the fractionBasis is given in mole fractions,
+        // convert from mole fractions to normalized mass fractions
+        scalar W(0);
+        forAll(comp, i)
+        {
+            comp[i] *= thermo_.W(i);
+            W += comp[i];
+        }
+        comp /= W;
+    }
+    else
+    {
+        FatalIOErrorInFunction(subDict)
+            << "The given fractionBasis type is invalid" << nl
+            << "Valid fractionBasis types are" << nl
+            << "  \"mass\" (default)" << nl
+            << "  \"mole\""
+            << exit(FatalIOError);
+
+        return false;
+    }
+
+    return true;
+}
+
+
+Foam::scalar Foam::functionObjects::BilgerMixtureFraction::o2Required
+(
+    const scalarField& comp
+) const
+{
+    scalar o2req(0);
+    forAll(thermo_.Y(), i)
+    {
+        o2req +=
+            comp[i]/thermo_.W(i)*(nAtomsC_[i] + nAtomsS_[i] + 0.25*nAtomsH_[i]);
+    }
+
+    return o2req;
+}
+
+
+Foam::scalar Foam::functionObjects::BilgerMixtureFraction::o2Present
+(
+    const scalarField& comp
+) const
+{
+    scalar o2pres(0);
+    forAll(thermo_.Y(), i)
+    {
+        o2pres += comp[i]/thermo_.W(i)*nAtomsO_[i];
+    }
+
+    return 0.5*o2pres;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::functionObjects::BilgerMixtureFraction::BilgerMixtureFraction
+(
+    const word& name,
+    const Time& runTime,
+    const dictionary& dict
+)
+:
+    fvMeshFunctionObject(name, runTime, dict),
+    phaseName_(dict.getOrDefault<word>("phase", word::null)),
+    resultName_
+    (
+        dict.getOrDefault<word>
+        (
+            "result",
+            IOobject::groupName("f_Bilger", phaseName_)
+        )
+    ),
+    thermo_
+    (
+        mesh_.lookupObject<basicSpecieMixture>
+        (
+            IOobject::groupName(basicThermo::dictName, phaseName_)
+        )
+    ),
+    nSpecies_(thermo_.Y().size()),
+    o2RequiredOx_(0),
+    o2RequiredFuelOx_(0),
+    nAtomsC_(nSpecies_, 0),
+    nAtomsS_(nSpecies_, 0),
+    nAtomsH_(nSpecies_, 0),
+    nAtomsO_(nSpecies_, 0),
+    Yoxidiser_(nSpecies_, 0),
+    Yfuel_(nSpecies_, 0)
+{
+    read(dict);
+
+    calcBilgerMixtureFraction();
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+bool Foam::functionObjects::BilgerMixtureFraction::read(const dictionary& dict)
+{
+    if (!fvMeshFunctionObject::read(dict))
+    {
+        return false;
+    }
+
+    Info<< nl << type() << " " << name() << ":" << nl;
+
+    phaseName_ = dict.getOrDefault<word>("phase", word::null);
+    resultName_ =
+        dict.getOrDefault<word>
+        (
+            "result",
+            IOobject::groupName("f_Bilger", phaseName_)
+        );
+
+    nSpecies_ = thermo_.Y().size();
+
+    if (nSpecies_ == 0)
+    {
+        FatalErrorInFunction
+            << "Number of input species is zero"
+            << exit(FatalError);
+    }
+
+    nAtomsC_.resize(nSpecies_, 0);
+    nAtomsS_.resize(nSpecies_, 0);
+    nAtomsH_.resize(nSpecies_, 0);
+    nAtomsO_.resize(nSpecies_, 0);
+
+    auto& Y = thermo_.Y();
+    const speciesTable& speciesTab = thermo_.species();
+
+    typedef BasicChemistryModel<psiReactionThermo> psiChemistryModelType;
+    typedef BasicChemistryModel<rhoReactionThermo> rhoChemistryModelType;
+
+    const auto* psiChemPtr =
+        mesh_.cfindObject<psiChemistryModelType>("chemistryProperties");
+
+    const auto* rhoChemPtr =
+        mesh_.cfindObject<rhoChemistryModelType>("chemistryProperties");
+
+    autoPtr<HashTable<List<specieElement>>> speciesCompPtr;
+
+    if (psiChemPtr)
+    {
+        speciesCompPtr.reset((*psiChemPtr).thermo().specieComposition());
+    }
+    else if (rhoChemPtr)
+    {
+        speciesCompPtr.reset((*rhoChemPtr).thermo().specieComposition());
+    }
+    else
+    {
+        FatalErrorInFunction
+            << "BasicChemistryModel not found"
+            << exit(FatalError);
+    }
+
+    forAll(Y, i)
+    {
+        const List<specieElement>& curSpecieComposition =
+            (speciesCompPtr.ref())[speciesTab[i]];
+
+        forAll(curSpecieComposition, j)
+        {
+            const word& e = curSpecieComposition[j].name();
+            const label nAtoms  = curSpecieComposition[j].nAtoms();
+
+            if (e == "C")
+            {
+                nAtomsC_[i] = nAtoms;
+            }
+            else if (e == "S")
+            {
+                nAtomsS_[i] = nAtoms;
+            }
+            else if (e == "H")
+            {
+                nAtomsH_[i] = nAtoms;
+            }
+            else if (e == "O")
+            {
+                nAtomsO_[i] = nAtoms;
+            }
+        }
+    }
+
+    if (sum(nAtomsO_) == 0)
+    {
+        FatalErrorInFunction
+            << "No specie contains oxygen"
+            << exit(FatalError);
+    }
+
+    Yoxidiser_.resize(nSpecies_, 0);
+    Yfuel_.resize(nSpecies_, 0);
+
+    if
+    (
+        !readComposition(dict.subDict("oxidiser"), Yoxidiser_)
+     || !readComposition(dict.subDict("fuel"), Yfuel_)
+    )
+    {
+        return false;
+    }
+
+    o2RequiredOx_ = o2Required(Yoxidiser_) - o2Present(Yoxidiser_);
+
+    if (o2RequiredOx_ > 0)
+    {
+        FatalErrorInFunction
+            << "Oxidiser composition contains not enough oxygen" << endl
+            << "Mixed up fuel and oxidiser compositions?"
+            << exit(FatalError);
+    }
+
+    const scalar o2RequiredFuel = o2Required(Yfuel_) - o2Present(Yfuel_);
+
+    if (o2RequiredFuel < 0)
+    {
+        FatalErrorInFunction
+            << "Fuel composition contains too much oxygen" << endl
+            << "Mixed up fuel and oxidiser compositions?"
+            << exit(FatalError);
+    }
+
+    o2RequiredFuelOx_ = o2RequiredFuel - o2RequiredOx_;
+
+    if (mag(o2RequiredFuelOx_) < SMALL)
+    {
+        FatalErrorInFunction
+            << "Fuel and oxidiser have the same composition"
+            << exit(FatalError);
+    }
+
+    return true;
+}
+
+
+bool Foam::functionObjects::BilgerMixtureFraction::execute()
+{
+    calcBilgerMixtureFraction();
+
+    return true;
+}
+
+
+bool Foam::functionObjects::BilgerMixtureFraction::clear()
+{
+    return clearObject(resultName_);
+}
+
+
+bool Foam::functionObjects::BilgerMixtureFraction::write()
+{
+    Log << type() << " " << name() << " write:" << nl
+        << "    writing field " << resultName_ << endl;
+
+    return writeObject(resultName_);
+}
+
+
+// ************************************************************************* //

src/thermophysicalModels/chemistryModel/functionObjects/BilgerMixtureFraction/BilgerMixtureFraction.H

@@ -0,0 +1,301 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2020 Thorsten Zirwes
+    Copyright (C) 2020 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::functionObjects::BilgerMixtureFraction
+
+Group
+    grpThermophysicalFunctionObjects
+
+Description
+    Calculates the Bilger mixture-fraction field (i.e.
+    \f$f_{\mathrm{Bilger}}\f$) based on the elemental composition
+    of the mixture. Elements \c C, \c H, \c S and \c O are considered.
+
+    \f$f_{\mathrm{Bilger}}\f$ is the mass mixing ratio of fuel and oxidiser
+    in [kg fuel / kg mixture], and is invariant to the reaction progress of
+    the mixture (e.g. the same for unburnt and burnt mixtures).
+
+    \f$f_{\mathrm{Bilger}}\f$ equals to the unity
+    for pure fuel and to zero for pure oxidiser.
+
+    The Bilger mixture-fraction field is computed based on the following:
+
+    \f[
+        f_{\mathrm{Bilger}} =
+            \frac{\beta - \beta_{\mathrm{ox}}}{\beta_{\mathrm{fuel}}
+          - \beta_{\mathrm{ox}}}
+    \f]
+
+    with
+
+    \f[
+        \beta =
+            2\frac{Y_C}{W_C}
+          + 2\frac{Y_S}{W_S}
+          + \frac{1}{2}\frac{Y_H}{W_H}
+          - \frac{Y_O}{W_O}
+    \f]
+
+    where
+    \vartable
+      \beta               | Coupling function                        [kmol/kg]
+      Y_e                 | Elemental mass fraction of element, e
+      W_e                 | Atomic weight of element, e
+      {.}_{\mathrm{ox}}   | Subscript to denote oxidiser composition
+      {.}_{\mathrm{fuel}} | Subscript to denote fuel composition
+    \endvartable
+
+    Operands:
+    \table
+      Operand        | Type           | Location
+      input          | -              | -
+      output file    | -              | -
+      output field   | volScalarField | $FOAM_CASE/\<time\>/\<outField\>
+    \endtable
+
+    References:
+    \verbatim
+        Original method:
+            Bilger, R. W. (1979).
+            Turbulent jet diffusion flames.
+            Energy and Combustion Science, p. 109-131.
+            DOI:10.1016/B978-0-08-024780-9.50011-3
+
+        Implementation:
+            Zirwes, T., Zhang, F., Habisreuther, P., Hansinger, M.,
+            Bockhorn, H., Pfitzner, M., & Trimis, D. (2019).
+            Quasi-DNS dataset of a piloted flame
+            with inhomogeneous inlet conditions.
+            Flow, Turbulence and Combustion, vol 104, p. 997-1027.
+            DOI:10.1007/s10494-019-00081-5
+    \endverbatim
+
+Usage
+    Minimal example by using \c system/controlDict.functions:
+    \verbatim
+    BilgerMixtureFraction1
+    {
+        // Mandatory entries (unmodifiable)
+        type                 BilgerMixtureFraction;
+        libs                 (fieldFunctionObjects);
+
+        // Mandatory entries (runtime modifiable)
+        fuel
+        {
+            // Optional entries (runtime modifiable)
+            fractionBasis    mass;
+
+            // Conditional mandatory entries (runtime modifiable)
+            CH4              1;
+        }
+
+        oxidiser
+        {
+            // Optional entries (runtime modifiable)
+            fractionBasis    mole;
+
+            // Conditional mandatory entries (runtime modifiable)
+            O2               0.23;
+            N2               0.77;
+        }
+
+        // Optional entries (runtime modifiable)
+        phase                <phaseName>;
+        result               <resultName>;
+
+        // Optional (inherited) entries
+        ...
+    }
+    \endverbatim
+
+    where the entries mean:
+    \table
+      Property | Description                         | Type | Reqd | Dflt
+      type     | Type name: BilgerMixtureFraction    | word | yes  | -
+      libs     | Library name: fieldFunctionObjects  | word | yes  | -
+      fuel     | Dictionary for fuel composition     | dict | yes  | -
+      oxidiser | Dictionary for oxidiser composition | dict | yes  | -
+      phase    | Name of phase (e.g. "gas")          | word | no   | ""
+      result   | Name of resulting field             | word | no   | f_Bilger
+      fractionBasis | Species-fraction interpretation method | word | no | mass
+    \endtable
+
+    Options for the \c fractionBasis entry:
+    \verbatim
+      mass   | Interpret species fractions as mass fractions
+      mole   | Interpret species fractions as mole fractions
+    \endverbatim
+
+    The inherited entries are elaborated in:
+     - \link functionObject.H \endlink
+
+    Usage by the \c postProcess utility is not available.
+
+Note
+  - The mole or mass fractions are automatically normalized to the unity.
+
+See also
+  - Foam::functionObject
+  - Foam::functionObjects::fvMeshFunctionObject
+
+SourceFiles
+    BilgerMixtureFraction.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef BilgerMixtureFraction_H
+#define BilgerMixtureFraction_H
+
+#include "fvMeshFunctionObject.H"
+#include "specieElement.H"
+#include "basicSpecieMixture.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace functionObjects
+{
+/*---------------------------------------------------------------------------*\
+                    Class BilgerMixtureFraction Declaration
+\*---------------------------------------------------------------------------*/
+
+class BilgerMixtureFraction
+:
+    public fvMeshFunctionObject
+{
+    // Private Data
+
+        //- Name of the phase (e.g. "gas" for multiphase applications)
+        word phaseName_;
+
+        //- Name of the resulting mixture-fraction field
+        word resultName_;
+
+        //- Reference to thermo object
+        const basicSpecieMixture& thermo_;
+
+        //- Number of species
+        label nSpecies_;
+
+        // Amount of oxygen required to fully oxidise the oxidiser
+        scalar o2RequiredOx_;
+
+        // Amount of oxygen required to oxidise the fuel minus the oxidiser
+        scalar o2RequiredFuelOx_;
+
+        //- Number of carbon atoms for each species
+        labelField nAtomsC_;
+
+        //- Number of sulphur atoms for each species
+        labelField nAtomsS_;
+
+        //- Number of hydrogen atoms for each species
+        labelField nAtomsH_;
+
+        //- Number of oxygen atoms for each species
+        labelField nAtomsO_;
+
+        //- Mass fractions of species in the oxidiser
+        scalarField Yoxidiser_;
+
+        //- Mass fractions of species in the fuel
+        scalarField Yfuel_;
+
+
+    // Private Member Functions
+
+        //- Calculate the Bilger mixture-fraction
+        void calcBilgerMixtureFraction();
+
+        //- Read composition of fuel and oxidiser from subdictionary
+        bool readComposition
+        (
+            const dictionary& subDict,
+            scalarField& comp
+        );
+
+        //- Compute amount of oxygen required to oxidise a mixture
+        scalar o2Present(const scalarField&) const;
+
+        //- Compute amount of oxygen present in a mixture
+        scalar o2Required(const scalarField&) const;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("BilgerMixtureFraction");
+
+
+    // Constructors
+
+        //- Construct from Time and dictionary
+        BilgerMixtureFraction
+        (
+            const word& name,
+            const Time& runTime,
+            const dictionary& dict
+        );
+
+        //- No copy construct
+        BilgerMixtureFraction(const BilgerMixtureFraction&) = delete;
+
+        //- No copy assignment
+        void operator=(const BilgerMixtureFraction&) = delete;
+
+
+    //- Destructor
+    virtual ~BilgerMixtureFraction() = default;
+
+
+    // Member Functions
+
+        //- Read the BilgerMixtureFraction data
+        virtual bool read(const dictionary&);
+
+        //- Calculate the Bilger mixture-fraction field
+        virtual bool execute();
+
+        //- Clear the Bilger mixture-fraction field from registry
+        virtual bool clear();
+
+        //- Write Bilger mixture-fraction field
+        virtual bool write();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace functionObjects
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //