INT: age/comfort: add new field function objects

src/functionObjects/field/Make/files

@@ -128,4 +128,7 @@ DMD/DMDModels/DMDModel/DMDModel.C
 DMD/DMDModels/DMDModel/DMDModelNew.C
 DMD/DMDModels/derived/STDMD/STDMD.C
 
+age/age.C
+comfort/comfort.C
+
 LIB = $(FOAM_LIBBIN)/libfieldFunctionObjects

src/functionObjects/field/age/age.C

@@ -0,0 +1,296 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2018-2021 OpenFOAM Foundation
+    Copyright (C) 2021 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 "age.H"
+#include "fvmDiv.H"
+#include "fvmLaplacian.H"
+#include "fvOptions.H"
+#include "turbulentTransportModel.H"
+#include "turbulentFluidThermoModel.H"
+#include "turbulenceModel.H"
+#include "inletOutletFvPatchField.H"
+#include "wallFvPatch.H"
+#include "zeroGradientFvPatchField.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace functionObjects
+{
+    defineTypeNameAndDebug(age, 0);
+    addToRunTimeSelectionTable(functionObject, age, dictionary);
+}
+}
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+Foam::wordList Foam::functionObjects::age::patchTypes() const
+{
+    wordList result
+    (
+        mesh_.boundary().size(),
+        inletOutletFvPatchField<scalar>::typeName
+    );
+
+    forAll(mesh_.boundary(), patchi)
+    {
+        if (isA<wallFvPatch>(mesh_.boundary()[patchi]))
+        {
+            result[patchi] = zeroGradientFvPatchField<scalar>::typeName;
+        }
+    }
+
+    return result;
+}
+
+
+bool Foam::functionObjects::age::converged
+(
+    const int nCorr,
+    const scalar initialResidual
+) const
+{
+    if (initialResidual < tolerance_)
+    {
+        Info<< "Field " << typeName
+            << " converged in " << nCorr << " correctors"
+            << nl << endl;
+
+        return true;
+    }
+
+    return false;
+}
+
+
+template<class GeoField>
+Foam::autoPtr<GeoField>
+Foam::functionObjects::age::newField
+(
+    const word& baseName,
+    const wordList patches
+) const
+{
+    return autoPtr<GeoField>::New
+    (
+        IOobject
+        (
+            scopedName(baseName),
+            time_.timeName(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh_,
+        dimensioned<typename GeoField::value_type>(dimTime, Zero),
+        patches
+    );
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::functionObjects::age::age
+(
+    const word& name,
+    const Time& runTime,
+    const dictionary& dict
+)
+:
+    fvMeshFunctionObject(name, runTime, dict)
+{
+    read(dict);
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+bool Foam::functionObjects::age::read(const dictionary& dict)
+{
+    if (fvMeshFunctionObject::read(dict))
+    {
+        phiName_ = dict.getOrDefault<word>("phi", "phi");
+        rhoName_ = dict.getOrDefault<word>("rho", "rho");
+        schemesField_ = dict.getOrDefault<word>("schemesField", typeName);
+        tolerance_ = dict.getOrDefault<scalar>("tolerance", 1e-5);
+        nCorr_ = dict.getOrDefault<int>("nCorr", 5);
+        diffusion_ = dict.getOrDefault<bool>("diffusion", false);
+
+        return true;
+    }
+
+    return false;
+}
+
+
+bool Foam::functionObjects::age::execute()
+{
+    auto tage = tmp<volScalarField>::New
+    (
+        IOobject
+        (
+            typeName,
+            mesh_.time().timeName(),
+            mesh_,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE,
+            false
+        ),
+        mesh_,
+        dimensionedScalar(dimTime, 0),
+        patchTypes()
+    );
+    volScalarField& age = tage.ref();
+
+    const word divScheme("div(phi," + schemesField_ + ")");
+
+    // Set under-relaxation coeff
+    scalar relaxCoeff = 0;
+    if (mesh_.relaxEquation(schemesField_))
+    {
+        relaxCoeff = mesh_.equationRelaxationFactor(schemesField_);
+    }
+
+    Foam::fv::options& fvOptions(Foam::fv::options::New(mesh_));
+
+
+    // This only works because the null constructed inletValue for an
+    // inletOutletFvPatchField is zero. If we needed any other value we would
+    // have to loop over the inletOutlet patches and explicitly set the
+    // inletValues. We would need to change the interface of inletOutlet in
+    // order to do this.
+
+    const auto& phi = mesh_.lookupObject<surfaceScalarField>(phiName_);
+
+    if (phi.dimensions() == dimMass/dimTime)
+    {
+        const auto& rho = mesh_.lookupObject<volScalarField>(rhoName_);
+
+        tmp<volScalarField> tmuEff;
+        word laplacianScheme;
+
+        if (diffusion_)
+        {
+            tmuEff =
+                mesh_.lookupObject<compressible::turbulenceModel>
+                (
+                    turbulenceModel::propertiesName
+                ).muEff();
+
+            laplacianScheme =
+                "laplacian(" + tmuEff().name() + ',' + schemesField_ + ")";
+        }
+
+        for (int i = 0; i <= nCorr_; ++i)
+        {
+            fvScalarMatrix ageEqn
+            (
+                fvm::div(phi, age, divScheme) == rho //+ fvOptions(rho, age)
+            );
+
+            if (diffusion_)
+            {
+                ageEqn -= fvm::laplacian(tmuEff(), age, laplacianScheme);
+            }
+
+            ageEqn.relax(relaxCoeff);
+
+            fvOptions.constrain(ageEqn);
+
+            if (converged(i, ageEqn.solve().initialResidual()))
+            {
+                break;
+            };
+
+            fvOptions.correct(age);
+        }
+    }
+    else
+    {
+        tmp<volScalarField> tnuEff;
+        word laplacianScheme;
+
+        if (diffusion_)
+        {
+            tnuEff =
+                mesh_.lookupObject<incompressible::turbulenceModel>
+                (
+                    turbulenceModel::propertiesName
+                ).nuEff();
+
+            laplacianScheme =
+                "laplacian(" + tnuEff().name() + ',' + schemesField_ + ")";
+        }
+
+        for (int i = 0; i <= nCorr_; ++i)
+        {
+            fvScalarMatrix ageEqn
+            (
+                fvm::div(phi, age, divScheme)
+             == dimensionedScalar(1) + fvOptions(age)
+            );
+
+            if (diffusion_)
+            {
+                ageEqn -= fvm::laplacian(tnuEff(), age, laplacianScheme);
+            }
+
+            ageEqn.relax(relaxCoeff);
+
+            fvOptions.constrain(ageEqn);
+
+            if (converged(i, ageEqn.solve().initialResidual()))
+            {
+                break;
+            }
+
+            fvOptions.correct(age);
+        }
+    }
+
+    Info<< "Min/max age:"
+        << min(age).value() << ' '
+        << max(age).value()
+        << endl;
+
+    // Workaround
+    word fieldName = typeName;
+
+    return store(fieldName, tage);
+}
+
+
+bool Foam::functionObjects::age::write()
+{
+    return writeObject(typeName);
+}
+
+
+// ************************************************************************* //

src/functionObjects/field/age/age.H

@@ -0,0 +1,196 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2018-2021 OpenFOAM Foundation
+    Copyright (C) 2021 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::age
+
+Description
+    Calculates and writes out the time taken for a particle to travel from an
+    inlet to the location. Solves the following equation when incompressible:
+
+    \f[
+        \div (\phi t) = 1
+    \f]
+
+    where
+    \vartable
+        t    | Age [s]
+        \phi | Volumetric flux [m^3/s]
+    \endvartable
+
+    Boundary conditions are generated automatically as \c zeroGradient
+    on all walls and \c inletOutlet everywhere else.
+
+Usage
+    Minimal example by using \c system/controlDict.functions:
+    \verbatim
+    age1
+    {
+        // Mandatory entries
+        type            age;
+        libs            (fieldFunctionObjects);
+
+        // Optional entries
+        phi             <word>;
+        rho             <word>;
+        schemesField    <word>;
+        tolerance       <scalar>;
+        nCorr           <int>;
+        diffusion       <bool>;
+
+        // Inherited entries
+        ...
+    }
+    \endverbatim
+
+    where the entries mean:
+    \table
+      Property  | Description                        | Type | Reqd | Deflt
+      type      | Type name: age                     | word | yes  | -
+      libs      | Library name: fieldFunctionObjects | word | yes  | -
+      phi       | Name of flux field                 | word | no   | phi
+      rho       | Name of density field              | word | no   | rho
+      schemesField | Name of the field from which schemes are taken <!--
+                -->                                  | word | no   | age
+      tolerance | Solver residual control            | scalar | no | 1e-5
+      nCorr     | Maximum number of correctors       | int  | no   | 5
+      diffusion | Flag to turn on/off the diffusion term | bool  | no | false
+    \endtable
+
+    The inherited entries are elaborated in:
+     - \link functionObject.H \endlink
+
+See also
+  - Foam::functionObjects::comfort
+
+SourceFiles
+    age.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef functionObjects_age_H
+#define functionObjects_age_H
+
+#include "fvMeshFunctionObject.H"
+#include "volFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace functionObjects
+{
+
+/*---------------------------------------------------------------------------*\
+                             Class age Declaration
+\*---------------------------------------------------------------------------*/
+
+class age
+:
+    public fvMeshFunctionObject
+{
+    // Private Data
+
+        //- The name of the flux field
+        word phiName_;
+
+        //- The name of the density field
+        word rhoName_;
+
+        //- Name of field from which schemes are taken
+        word schemesField_;
+
+        //- Convergence tolerance
+        scalar tolerance_;
+
+        //- Number of corrections
+        int nCorr_;
+
+        //- Flag to turn on/off the diffusion term
+        bool diffusion_;
+
+
+    // Private Member Functions
+
+        //- The list of patch types for the age field
+        wordList patchTypes() const;
+
+        //- Return true if convergence is reached
+        bool converged(const int nCorr, const scalar initialResidual) const;
+
+        //- Create and allocate a new zero geometric field
+        template<class GeoField>
+        autoPtr<GeoField> newField
+        (
+            const word& baseName,
+            const wordList
+        ) const;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("age");
+
+
+    // Constructors
+
+        //- Construct from Time and dictionary
+        age
+        (
+            const word& name,
+            const Time& runTime,
+            const dictionary& dict
+        );
+
+
+    //- Destructor
+    virtual ~age() = default;
+
+
+    // Member Functions
+
+        //- Read the data
+        virtual bool read(const dictionary&);
+
+        //- Execute
+        virtual bool execute();
+
+        //- Write
+        virtual bool write();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace functionObjects
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/functionObjects/field/comfort/comfort.C

@@ -0,0 +1,478 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2019 OpenFOAM Foundation
+    Copyright (C) 2021 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 "comfort.H"
+#include "wallFvPatch.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace functionObjects
+{
+    defineTypeNameAndDebug(comfort, 0);
+    addToRunTimeSelectionTable(functionObject, comfort, dictionary);
+}
+}
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+Foam::tmp<Foam::volScalarField> Foam::functionObjects::comfort::magU() const
+{
+    tmp<volScalarField> tmagU = mag(lookupObject<volVectorField>("U"));
+    volScalarField& magU = tmagU.ref();
+
+    // Flag to use the averaged velocity field in the domain.
+    // Consistent with EN ISO 7730 but does not make physical sense
+    if (meanVelocity_)
+    {
+        magU = magU.weightedAverage(mesh_.V());
+    }
+
+    return tmagU;
+}
+
+
+Foam::dimensionedScalar Foam::functionObjects::comfort::Trad() const
+{
+    dimensionedScalar Trad(Trad_);
+
+    // The mean radiation is calculated by the mean wall temperatures
+    // which are summed and divided by the area | only walls are taken into
+    // account. This approach might be correct for a squared room but will
+    // defintely be inconsistent for complex room geometries. The norm does
+    // not provide any information about the calculation of this quantity.
+    if (!TradSet_)
+    {
+        const volScalarField::Boundary& TBf =
+            lookupObject<volScalarField>("T").boundaryField();
+
+        scalar areaIntegral = 0;
+        scalar TareaIntegral = 0;
+
+        forAll(TBf, patchi)
+        {
+            const fvPatchScalarField& pT = TBf[patchi];
+            const fvPatch& pTBf = TBf[patchi].patch();
+            const scalarField& pSf = pTBf.magSf();
+
+            if (isType<wallFvPatch>(pTBf))
+            {
+                areaIntegral += gSum(pSf);
+                TareaIntegral += gSum(pSf*pT);
+            }
+        }
+
+        Trad.value() = TareaIntegral/areaIntegral;
+    }
+
+    // Bounds based on EN ISO 7730
+    if ((Trad.value() < 283.15) || (Trad.value() > 313.15))
+    {
+        WarningInFunction
+            << "The calculated mean wall radiation temperature is out of the\n"
+            << "bounds specified in EN ISO 7730:2005\n"
+            << "Valid range is 10 degC < T < 40 degC\n"
+            << "The actual value is: " << Trad - 273.15 << nl << endl;
+    }
+
+    return Trad;
+}
+
+
+Foam::tmp<Foam::volScalarField> Foam::functionObjects::comfort::pSat() const
+{
+    static const dimensionedScalar kPaToPa(dimPressure, 1000);
+    static const dimensionedScalar A(dimless, 16.6563);
+    static const dimensionedScalar B(dimTemperature, 4030.183);
+    static const dimensionedScalar C(dimTemperature, -38.15);
+
+    tmp<volScalarField> tpSat = volScalarField::New("pSat", mesh_, pSat_);
+
+    // Calculate the saturation pressure if no user input is given
+    if (pSat_.value() == 0)
+    {
+        const auto& T = lookupObject<volScalarField>("T");
+
+        // Equation based on ISO 7730:2006
+        tpSat = kPaToPa*exp(A - B/(T + C));
+    }
+
+    return tpSat;
+}
+
+
+Foam::tmp<Foam::volScalarField> Foam::functionObjects::comfort::Tcloth
+(
+    volScalarField& hc,
+    const dimensionedScalar& metabolicRateSI,
+    const dimensionedScalar& extWorkSI,
+    const volScalarField& T,
+    const dimensionedScalar& Trad
+)
+{
+    const dimensionedScalar factor1(dimTemperature, 308.85);
+
+    const dimensionedScalar factor2
+    (
+        dimTemperature/metabolicRateSI.dimensions(),
+        0.028
+    );
+
+    const dimensionedScalar factor3
+    (
+        dimMass/pow3(dimTime)/pow4(dimTemperature),
+        3.96e-8
+    );
+
+    // Heat transfer coefficient based on forced convection [W/m^2/K]
+    const volScalarField hcForced
+    (
+        dimensionedScalar(hc.dimensions()/sqrt(dimVelocity), 12.1)
+       *sqrt(magU())
+    );
+
+    // Tcl [K] (surface cloth temperature)
+    tmp<volScalarField> tTcl
+    (
+        volScalarField::New
+        (
+            "Tcl",
+            T.mesh(),
+            dimTemperature
+        )
+    );
+    volScalarField& Tcl = tTcl.ref();
+
+    // Initial guess
+    Tcl = T;
+
+    label i = 0;
+
+    Tcl.storePrevIter();
+
+    // Same temperatures as for the radiation
+    const dimensionedScalar Tmin(dimTemperature, 283.15);
+    const dimensionedScalar Tmax(dimTemperature, 313.15);
+
+    // Iterative solving of equation (2)
+    do
+    {
+        Tcl = (Tcl + Tcl.prevIter())/2;
+        Tcl.storePrevIter();
+
+        // Heat transfer coefficient based on natural convection
+        volScalarField hcNatural
+        (
+            dimensionedScalar(hc.dimensions()/pow025(dimTemperature), 2.38)
+           *pow025(mag(Tcl - T))
+        );
+
+        // Set heat transfer coefficient based on equation (3)
+        hc =
+            pos(hcForced - hcNatural)*hcForced
+          + neg0(hcForced - hcNatural)*hcNatural;
+
+        // Calculate surface temperature based on equation (2)
+        Tcl =
+            factor1
+          - factor2*(metabolicRateSI - extWorkSI)
+          - Icl_*factor3*fcl_*(pow4(Tcl) - pow4(Trad))
+          - Icl_*fcl_*hc*(Tcl - T);
+
+        // Make sure that Tcl is in some physical limit (same range as we used
+        // for the radiative estimation - based on ISO EN 7730:2005)
+        Tcl.clip(Tmin, Tmax);
+
+    } while (!converged(Tcl) && i++ < maxClothIter_);
+
+    if (i == maxClothIter_)
+    {
+        WarningInFunction
+            << "The surface cloth temperature did not converge within " << i
+            << " iterations" << nl;
+    }
+
+    return tTcl;
+}
+
+
+bool Foam::functionObjects::comfort::converged
+(
+    const volScalarField& phi
+) const
+{
+    return
+        max(mag(phi.primitiveField() - phi.prevIter().primitiveField()))
+      < tolerance_;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::functionObjects::comfort::comfort
+(
+    const word& name,
+    const Time& runTime,
+    const dictionary& dict
+)
+:
+    fvMeshFunctionObject(name, runTime, dict),
+    clothing_("clothing", dimless, 0),
+    metabolicRate_("metabolicRate", dimMass/pow3(dimTime), 0.8),
+    extWork_("extWork", dimMass/pow3(dimTime), 0),
+    Trad_("Trad", dimTemperature, 0),
+    relHumidity_("relHumidity", dimless, 0.5),
+    pSat_("pSat", dimPressure, 0),
+    Icl_("Icl", pow3(dimTime)*dimTemperature/dimMass, 0),
+    fcl_("fcl", dimless, 0),
+    tolerance_(1e-4),
+    maxClothIter_(100),
+    TradSet_(false),
+    meanVelocity_(false)
+{
+    read(dict);
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+bool Foam::functionObjects::comfort::read(const dictionary& dict)
+{
+    if (fvMeshFunctionObject::read(dict))
+    {
+        clothing_.readIfPresent(dict);
+        metabolicRate_.readIfPresent(dict);
+        extWork_.readIfPresent(dict);
+        pSat_.readIfPresent(dict);
+        tolerance_ = dict.getOrDefault("tolerance", 1e-4);
+        maxClothIter_ = dict.getOrDefault("maxClothIter", 100);
+        meanVelocity_ = dict.getOrDefault<bool>("meanVelocity", false);
+
+        // Read relative humidity if provided and convert from % to fraction
+        if (dict.found(relHumidity_.name()))
+        {
+            relHumidity_.read(dict);
+            relHumidity_ /= 100;
+        }
+
+        // Read radiation temperature if provided
+        if (dict.found(Trad_.name()))
+        {
+            TradSet_ = true;
+            Trad_.read(dict);
+        }
+
+        Icl_ = dimensionedScalar(Icl_.dimensions(), 0.155)*clothing_;
+
+        fcl_.value() =
+            Icl_.value() <= 0.078
+          ? 1.0 + 1.290*Icl_.value()
+          : 1.05 + 0.645*Icl_.value();
+
+        return true;
+    }
+
+    return false;
+}
+
+
+bool Foam::functionObjects::comfort::execute()
+{
+    // Assign and build fields
+    const dimensionedScalar Trad(this->Trad());
+    const volScalarField pSat(this->pSat());
+
+    const dimensionedScalar metabolicRateSI(58.15*metabolicRate_);
+    const dimensionedScalar extWorkSI(58.15*extWork_);
+
+    const auto& T = lookupObject<volScalarField>("T");
+
+    // Heat transfer coefficient [W/m^2/K]
+    // This field is updated in Tcloth()
+    volScalarField hc
+    (
+        IOobject
+        (
+            "hc",
+            mesh_.time().timeName(),
+            mesh_
+        ),
+        mesh_,
+        dimensionedScalar(dimMass/pow3(dimTime)/dimTemperature, 0)
+    );
+
+    // Calculate the surface temperature of the cloth by an iterative
+    // process using equation (2) from DIN EN ISO 7730 [degC]
+    const volScalarField Tcloth
+    (
+        this->Tcloth
+        (
+            hc,
+            metabolicRateSI,
+            extWorkSI,
+            T,
+            Trad
+        )
+    );
+
+    // Calculate the PMV quantity
+    const dimensionedScalar factor1(pow3(dimTime)/dimMass, 0.303);
+    const dimensionedScalar factor2
+    (
+        dimless/metabolicRateSI.dimensions(),
+        -0.036
+    );
+    const dimensionedScalar factor3(factor1.dimensions(), 0.028);
+    const dimensionedScalar factor4(dimLength/dimTime, 3.05e-3);
+    const dimensionedScalar factor5(dimPressure, 5733);
+    const dimensionedScalar factor6(dimTime/dimLength, 6.99);
+    const dimensionedScalar factor8(metabolicRateSI.dimensions(), 58.15);
+    const dimensionedScalar factor9(dimless/dimPressure, 1.7e-5);
+    const dimensionedScalar factor10(dimPressure, 5867);
+    const dimensionedScalar factor11(dimless/dimTemperature, 0.0014);
+    const dimensionedScalar factor12(dimTemperature, 307.15);
+    const dimensionedScalar factor13
+    (
+        dimMass/pow3(dimTime)/pow4(dimTemperature),
+        3.96e-8
+    );
+
+    const scalar factor7
+    (
+        // Special treatment of Term4
+        // if metaRate - extWork < factor8, set to zero
+        (metabolicRateSI - extWorkSI).value() < factor8.value() ? 0 : 0.42
+    );
+
+    Info<< "Calculating the predicted mean vote (PMV)" << endl;
+
+    // Equation (1)
+    tmp<volScalarField> PMV =
+        (
+            // Term1: Thermal sensation transfer coefficient
+            (factor1*exp(factor2*metabolicRateSI) + factor3)
+           *(
+                (metabolicRateSI - extWorkSI)
+
+                // Term2: Heat loss difference through skin
+              - factor4
+               *(
+                    factor5
+                  - factor6*(metabolicRateSI - extWorkSI)
+                  - pSat*relHumidity_
+                )
+
+                // Term3: Heat loss through sweating
+              - factor7*(metabolicRateSI - extWorkSI - factor8)
+
+                // Term4: Heat loss through latent respiration
+              - factor9*metabolicRateSI*(factor10 - pSat*relHumidity_)
+
+                // Term5: Heat loss through dry respiration
+              - factor11*metabolicRateSI*(factor12 - T)
+
+                // Term6: Heat loss through radiation
+              - factor13*fcl_*(pow4(Tcloth) - pow4(Trad))
+
+                // Term7: Heat loss through convection
+              - fcl_*hc*(Tcloth - T)
+            )
+        );
+
+    Info<< "Calculating the predicted percentage of dissatisfaction (PPD)"
+        << endl;
+
+    // Equation (5)
+    tmp<volScalarField> PPD =
+        100 - 95*exp(-0.03353*pow4(PMV()) - 0.21790*sqr(PMV()));
+
+    Info<< "Calculating the draught rating (DR)\n";
+
+    const dimensionedScalar Umin(dimVelocity, 0.05);
+    const dimensionedScalar Umax(dimVelocity, 0.5);
+    const dimensionedScalar pre(dimless, 0.37);
+    const dimensionedScalar C1(dimVelocity, 3.14);
+
+    // Limit the velocity field to the values given in EN ISO 7733
+    volScalarField Umag(mag(lookupObject<volVectorField>("U")));
+    Umag.clip(Umin, Umax);
+
+    // Calculate the turbulent intensity if turbulent kinetic energy field k
+    // exists
+    volScalarField TI
+    (
+        IOobject
+        (
+            "TI",
+            mesh_.time().timeName(),
+            mesh_
+        ),
+        mesh_,
+        dimensionedScalar(dimless, 0)
+    );
+
+    if (foundObject<volScalarField>("k"))
+    {
+        const auto& k = lookupObject<volScalarField>("k");
+        TI = sqrt(2/3*k)/Umag;
+    }
+
+    // For unit correctness
+    const dimensionedScalar correctUnit
+    (
+        dimensionSet(0, -1.62, 1.62, -1, 0, 0, 0),
+        1
+    );
+
+    // Equation (6)
+    tmp<volScalarField> DR =
+        correctUnit*(factor12 - T)*pow(Umag - Umin, 0.62)*(pre*Umag*TI + C1);
+
+    // Workaround
+    word fieldNamePMV = "PMV";
+    word fieldNamePPD = "PPD";
+    word fieldNameDR = "DR";
+
+    return
+        store(fieldNamePMV, PMV)
+     && store(fieldNamePPD, PPD)
+     && store(fieldNameDR, DR);
+}
+
+bool Foam::functionObjects::comfort::write()
+{
+    return
+        writeObject("PMV")
+     && writeObject("PPD")
+     && writeObject("DR");
+}
+
+
+// ************************************************************************* //

src/functionObjects/field/comfort/comfort.H

@@ -0,0 +1,246 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2019-2021 OpenFOAM Foundation
+    Copyright (C) 2021 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::comfort
+
+Description
+    Calculates the thermal comfort quantities predicted mean vote (PMV),
+    predicted percentage of dissatisfaction (PPD) and the draught rate (DR)
+    based on DIN ISO EN 7730:2005.
+
+    The draught rate is defined for velocities between 0 m/s and 0.5 m/s. Values
+    larger than 0.5 m/s will be set to 0.5 m/s. Furthermore, the draught rate is
+    defined between 20 degC and 26 degC. A temperature limitation is not
+    implemented. The draught rate is mainly used for HVAC analysis in rooms.
+
+Usage
+    Minimal example by using \c system/controlDict.functions:
+    \verbatim
+    comfort1
+    {
+        // Mandatory entries
+        type            comfort;
+        libs            (fieldFunctionObjects);
+
+        // Optional entries
+        clothing        <scalar>;
+        metabolicRate   <scalar>;
+        extWork         <scalar>;
+        Trad            <scalar>;
+        relHumidity     <scalar>;
+        pSat            <scalar>;
+        tolerance       <scalar>;
+        maxClothIter    <int>;
+        meanVelocity    <bool>;
+
+        // Inherited entries
+        ...
+    }
+    \endverbatim
+
+    where the entries mean:
+    \table
+      Property  | Description                        | Type | Reqd | Deflt
+      type      | Type name: comfort                 | word | yes  | -
+      libs      | Library name: fieldFunctionObjects | word | yes  | -
+      clothing  | The insulation value of the cloth  | scalar | no | 0
+      metabolicRate | The metabolic rate             | scalar | no | 0.8
+      extWork   | The external work                  | scalar | no | 0
+      Trad      | Radiation temperature              | scalar | no | 0
+      relHumidity   | Relative humidity of the air   | scalar | no | 0.5
+      pSat      | Saturation pressure of water       | scalar | no | -1
+      tolerance | Residual control for the cloth temperature <!--
+                -->                                  | scalar | no | 1e-4
+      maxClothIter  | Maximum number of iterations   | int    | no | 100
+      meanVelocity  | Flag to use a constant mean velocity   <!--
+                -->             in the whole domain  | bool   | no | false
+    \endtable
+
+    The inherited entries are elaborated in:
+     - \link functionObject.H \endlink
+
+    \table
+        Predicted Mean Vote (PMV)   | evaluation
+        + 3                         | hot
+        + 2                         | warm
+        + 1                         | slightly warm
+        + 0                         | neutral
+        - 1                         | slightly cool
+        - 2                         | cool
+        - 3                         | cold
+    \endtable
+
+    \table
+        Draught rate based on 7730  | category
+        0 - 10                      | I - fine
+        10 - 20                     | II - okay
+        20 - 30                     | III - intermedian
+        > 30                        | bad - commonly too high
+    \endtable
+
+See also
+  - Foam::functionObjects::age
+
+SourceFiles
+    comfort.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef comfort_H
+#define comfort_H
+
+#include "fvMeshFunctionObject.H"
+#include "volFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace functionObjects
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class comfort Declaration
+\*---------------------------------------------------------------------------*/
+
+class comfort
+:
+    public fvMeshFunctionObject
+{
+    // Private Data
+
+        //- Clothing [-]
+        dimensionedScalar clothing_;
+
+        //- Metabolic rate [kg/s^3]
+        dimensionedScalar metabolicRate_;
+
+        //- External work [kg/s^3]
+        dimensionedScalar extWork_;
+
+        //- Mean radiation temperature [K]
+        dimensionedScalar Trad_;
+
+        //- Relative humidity [percentage]
+        dimensionedScalar relHumidity_;
+
+        //- Saturation pressure of water [Pa]
+        dimensionedScalar pSat_;
+
+        //- Thermal insulation of clothing [W/m^2/K]
+        dimensionedScalar Icl_;
+
+        //- Prefactor of cloth area [-]
+        dimensionedScalar fcl_;
+
+        //- Tolerance criteria for iterative process to find Tcl
+        scalar tolerance_;
+
+        //- Maximum number of correctors for cloth temperature
+        int maxClothIter_;
+
+        //- Flag to set to true if the radiation temperature is provided
+        bool TradSet_;
+
+        //- Flag to use volume weighted velocity field for caluclation
+        bool meanVelocity_;
+
+
+    // Private Member Functions
+
+        //- Calculate the magnitude of the velocity [m/s]
+        tmp<volScalarField> magU() const;
+
+        //- Calculate the radiation temperature in the domain using a simple
+        //- approach [K]
+        dimensionedScalar Trad() const;
+
+        //- Calculate the saturation pressure based on 7730:2005
+        //- Possible options: adding different calculation methods such as
+        //- the formulation based on Magnus or others [Pa]
+        tmp<volScalarField> pSat() const;
+
+        //- Calculate and return the surface temperature of the cloth [K]
+        //- and the heat transfer coefficient hc [W/m^2/K]
+        tmp<volScalarField> Tcloth
+        (
+            volScalarField& hc,
+            const dimensionedScalar& metabolicRateSI,
+            const dimensionedScalar& extWorkSI,
+            const volScalarField& TdegC,
+            const dimensionedScalar& Trad
+        );
+
+        //- Return true if the cloth temperature iteration has converged
+        bool converged(const volScalarField&) const;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("comfort");
+
+
+    // Constructors
+
+        //- Construct from Time and dictionary
+        comfort
+        (
+            const word& name,
+            const Time& runTime,
+            const dictionary& dict
+        );
+
+
+    //- Destructor
+    virtual ~comfort() = default;
+
+
+    // Member Functions
+
+        //- Read the data needed for the comfort calculation
+        virtual bool read(const dictionary&);
+
+        //- Calculate the predicted mean vote (PMV)
+        //- and predicted percentage dissatisfaction (PPD) fields
+        virtual bool execute();
+
+        //- Write the PPD and PMV fields
+        virtual bool write();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace functionObjects
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //