comfortFoam.C

Go to the documentation of this file.

#/*---------------------------------------------------------------------------*\
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
   \\  /    A nd           | Copyright held by original author
   \\/     M anipulation  |
   -------------------------------------------------------------------------------
   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 2 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, write to the Free Software Foundation,
   Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
   Application
   comfortFoam
 
   Description
   Das Comfort-Tool berechnet folgende Kennwerte:
   - Operative Raumlufttemperatur (TOp)
   - DR  (Draugth Rating oder Draft Risk), Wertebereich: 0 bis 100% (DR)
   - PMV (Predicted Percentage of Dissatisfied), Wertebereich: 0 bis 100% (PMV)
   - PPD (Predicted Mean Vote), Wertebereich: -3 (kalt) bis 3 (warm) (PPD)
   - Mittlere Verweildauer der Luft - Mean Age of Air (AoA)
   - Lüftungseffektivität (AE)
   - Berechnet die relative Raumluftfeuchtigkeit (Humidity), Wertebereich: 0 bis 100%
 
   内容描述
        舒适性工具需要计算以下参数：
            -1.操作室空气温度（TOp）
            -2.DR（危险等级或风险草案），值范围：0到100％（DR）
            -3.PMV（预测平均投票），取值范围：-3（冷）至3（暖）（PPD）
            -4.PPD（预测的不满意百分比），值范围：0到100％（PMV）
            -4.APMV（自适应PMV），值范围：0到100％（PMV）
            -5.平均空气年龄（AoA,AirAge）
            -6.通风效率（AE,air effecitve）
            -7.计算房间的相对湿度（湿度），值范围：0到100％
 
Grundlage:      "DIN EN ISO 7730"
Autor:      "Dipl.-Ing. Thomas Tian"
Version:        "1.5"
 
Source files:
createFields.H
 
\*---------------------------------------------------------------------------*/
 
#include "fvCFD.H"
#include "incompressible/singlePhaseTransportModel/singlePhaseTransportModel.H"
#include "wallFvPatch.H"
 
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
//计算平均风速
Foam::vector Uaverage(const Foam::fvMesh& mesh_)
{
    //获取体向量U
    const volVectorField& U = mesh_.lookupObject<volVectorField>("U");
    vector midU(0,0,0);
    scalar cellsum(0);
    scalar volume(0);
 
    //所有网格loop
    forAll (mesh_.cells(), cellI)
    {
        //sum 累加所有单个cell的速度*体积
        midU += U[cellI] * mesh_.V()[cellI];
        //累加每个cel的体积
        volume += mesh_.V()[cellI];
    };
 
    //将累加速度/累加体积，然后取标量获取大小值，获取速度大小.也就是使用提平均来计算风速平均值
    Info<< "Average velocity = "<< mag(midU/volume) << " m/s" << nl;
    //但是这里返回的居然还是平均风速的三个矢量
    return midU/(volume);
}
 
Foam::scalar Taverage(const Foam::fvMesh& mesh_)
{
    const volScalarField& T = mesh_.lookupObject<volScalarField>("T");
    scalar midT(0);
    scalar cellsum(0);
    scalar volume(0);
 
    forAll (mesh_.cells(), cellI)
    {
        //if (T[cellI]>283)
       // {
        midT += T[cellI] * mesh_.V()[cellI];
        volume += mesh_.V()[cellI];
       // }
 
    };
    Info<< "Average Temperature = "<< (midT/volume)-273.15 << " °C" << nl;
 
    return (midT/volume);
}
 
Foam::scalar radiationTemperature(const Foam::fvMesh& mesh_, const Foam::fvPatchList& Patches_)
{
    // Berechnet die Strahlungstemperatur der Umgebung
    // radiationTemperatur = Sum(Mittlere Wandoberflächentemperaturen / Anzahl Flächen)
    // Rückgabe: Temperatur in Kelvin
    //获取环境温度，是个标量
    const volScalarField& T = mesh_.lookupObject<volScalarField>("T");
    // 边界温度累计值，面积值
    scalar PatchSumTemp(0), area(0);
    //计数器
    int counter = 0;
    //循环所有的边界,每一个小的patch
    forAll (Patches_, patchI)
    {
        //当前patch的索引，是个整型数
        const label curPatch = Patches_[patchI].index();
        //如果类型是wall的patch
        //if (isType<wallFvPatch>( Patches_[patchI] ))
        if (isType<wallFvPatch>( Patches_[patchI] ))
        {
            //mesh_.magSf() 返回一个面face的面积标量大小
            //boundaryField()
            area = gSum(mesh_.magSf().boundaryField()[curPatch]);
            //如果面积不等于0
            if (area != 0)
            {
                PatchSumTemp +=
                    gSum
                    (
                     mesh_.magSf().boundaryField()[curPatch]
                     * T.boundaryField()[curPatch]
                    ) / area;
 
                counter++;
            }
        }
    }
    //return (T+PatchSumTemp / counter)/2 - 273.15;
    return PatchSumTemp / counter - 273.15;
}
 
int main(int argc, char *argv[])
{
 
//#    argList::noParallel();
 
#   include "addTimeOptions.H"
#   include "setRootCase.H"
#   include "createTime.H"
 
    // Get times list
    instantList Times = runTime.times();
 
    // set startTime and endTime depending on -time and -latestTime options
#   include "checkTimeOptions.H"
 
    runTime.setTime(Times[startTime], startTime);
    Foam::instantList timeDirs = Foam::timeSelector::select0(runTime, args);
 
#   include "createMesh.H"
 
    forAll(timeDirs, timeI)
    {
        runTime.setTime(timeDirs[timeI], timeI);
 
        Info<< "Time = " << runTime.timeName() << endl;
 
#include "createFields.H"
 
        // Operative Raumluftemperatur Top
        volScalarField T(THeader, mesh);
        const fvPatchList& Patches = T.mesh().boundary();
 
        scalar STemp(20);
 
        if ((w.headerOk()==1) )
        {
            RHswitch = true;
        }
 
        if (QrHeader.headerOk()==1 )
        {
            volScalarField Qr(THeader, mesh);
            const fvPatchList& Patches_ = Qr.mesh().boundary();
            scalar PatchSumTemp(0);
 
            forAll (Patches_, patchI)
            {
                const label curPatch = Patches_[patchI].index();
 
                if (isType<wallFvPatch>( Patches_[patchI] ))
                {
                    PatchSumTemp +=
                        gSum
                        (
                         mesh.magSf().boundaryField()[curPatch]
                         * Qr.boundaryField()[curPatch]
                        );
                }
            }
            Info<< "Sum of all heat flows: "<< PatchSumTemp-273.15 << " oC" << endl;
        }
        else
        {
            STemp = radiationTemperature(mesh, Patches);
            Info << "function for Average Radiation Temperature: " << STemp << " oC" << endl;
        };
 
        if (AoAHeader.headerOk()==1 )
        {
 
            volScalarField AoA(AoAHeader, mesh);
            scalar cellsum(0);
            scalar midAoA(0);
 
            forAll (mesh.cells(), cellI)
            {
                midAoA += AoA[cellI];
                cellsum++;
            };
 
            scalar maxValue = max(AoA).value();
 
            Info << "Average age of air " << (midAoA/cellsum) << " s" << endl;
            Info << "Maximum age of air " << maxValue << " s" << endl;
 
 
            // Lüftungseffektivität AE = (Tn/2TM) x 100 %
            Info << "Room volume " << gSum( mesh.V() ) << " m3" << endl;
            Info << "Ventilation efficiency AE = " << ( gSum( mesh.V() ) / sumZuluft )  / (2 * ( maxValue / 3600) )  << " %" << endl;
        };
 
        volVectorField U(UHeader,mesh);
        vector wl = Uaverage(mesh);
        scalar Tr = Taverage(mesh);
        volScalarField p_rgh(p_rghHeader,mesh);
 
        scalar cellsum(0);
        scalar P1(0), P2(0), P3(0), P4(0), P5(0), XN(0), XF(0), HCN(0), HCF(0), HC(0), PA(0), FCL(0), EPS(0), ICL(0);
        scalar Tu(0), HL1(0), HL2(0), HL3(0), HL4(0), HL5(0), HL6(0), TCL(0), TS(0), TCLA(0), midPMV(0), midPPD(0),midAPMV(0), midDR(0), midTO(0), midRH(0);
        scalar p_w(0), p_ws(0.01);
 
        scalar volume(0);
 
        int N;
        double a_korr; // Korrekturfaktor a nach DIN EN 7730;
 
        // FANGER - Gleichungen (PMV)
 
        forAll (mesh.cells(), cellI)
        {
 
            // Zum Testen
            /*
               STemp= 21;
               Tr = 23.0+273.15;
               T[cellI] = Tr;
               U[cellI].x() = 0.3;
               U[cellI].y() = 0;
               U[cellI].z() = 0;
               */
 
            if (T[cellI] < 0)
                T[cellI] = 273;
 
            if (T[cellI] > 350)
                T[cellI] = 350;
 
 
            if (RHswitch)
            {
                // Berechne die Raumluftfeuchte
                p_w = ((101325 + p_rgh[cellI] ) * w[cellI])/(0.62198 + 0.37802 * w[cellI]);
 
                p_ws =  Foam::exp(-0.58002206*Foam::pow(10.0,4)*Foam::pow(T[cellI],-1)
                        +0.13914993*10
                        -0.48640239*Foam::pow(10.0,-1)*T[cellI]
                        +0.41764768*Foam::pow(10.0,-4)*Foam::pow(T[cellI],2)
                        -0.14452093*Foam::pow(10.0,-7)*Foam::pow(T[cellI],3)
                        +0.65459673*10*Foam::log(T[cellI]) );
 
                RH[cellI] = p_w/p_ws*100;
                PA = RH[cellI] * 10 * Foam::exp( 16.6563 - (4030.183 / (T[cellI] - 273.15 + 235) )); // water vapour pressure, Pa
                midRH += RH[cellI];
            } else {
                PA = RH1 * 10 * Foam::exp( 16.6563 - (4030.183 / (T[cellI] - 273.15 + 235) )); // water vapour pressure, Pa
                midRH += RH1;
            }
 
            // Berechne den Turbulenzgrad %
            // Tu = Wurzel ( 1/3 * (Ux'² + Uy'² + Uz'²) ) / Mittelwert U
            if ( mag(wl) >0 )
            {
                Tu =  (Foam::sqrt(scalar(1)/scalar(3) * (
                                Foam::pow( U[cellI].x() ,2) +
                                Foam::pow( U[cellI].y() ,2) +
                                Foam::pow( U[cellI].z() ,2)
                                )) / mag(wl)) * 100;
                // ToDo!!!!
                Tu = 40;
 
            }
            else
            {
                Tu = 0;
            };
 
            // Berechne DR-Wert
            // bei wl<0,05 m/s ist wl=0,05 m/s einzusetzen!
            if (mag(U[cellI]) >= 0.05)
            {
                DR[cellI] = ( 34 - (T[cellI] - 273.15) ) * ( Foam::pow( mag(U[cellI]) - 0.05 ,0.62 ) * ( (0.37 * mag(U[cellI]) * Tu) + 3.14) );
            }
            else
            {
                DR[cellI] = (34 - T[cellI] - 273.15) * Foam::pow( 0.05 ,0.6223) * ((0.37 * 0.05 * Tu) + 3.14);
            };
            if (DR[cellI] > 100)
            {
                DR[cellI] = 100;
            };
            if (DR[cellI] < 0)
            {
                DR[cellI] = 0;
            };
 
 
 
            ICL = 0.155 * clo;
 
            if (ICL < 0.078) // ok
            { FCL = 1 + 1.29 * ICL; } // ok
            else
            { FCL = 1.05 + 0.645 * ICL; }; // clothing area factor
 
            HCF = 12.1 * Foam::sqrt(mag( U[cellI] )); // heat transf. coeff. by forced convection
 
            TCLA = T[cellI] + (35.5 - (T[cellI] - 273.15)) / (3.5 * 6.45 * ( ICL + 0.1));
            XN = TCLA / 100; // ok
 
            P1 = ICL * FCL; // ok
            P2 = P1 * 3.96; // ok
            P3 = P1 * 100; // ok
            P4 = P1 * T[cellI]; // ok
            P5 = 308.7 - 0.028 * (met*58.15 - wme*58.15) + P2 * Foam::pow( (STemp + 273.0)/100, 4);
 
            // geändert Tian 06.10.2012
            // XF = XN;
            XF = TCLA / 50;
            EPS = 0.0015;
 
            N=0;
 
            do
            {
 
                N++;
                XF = (XF + XN)/2; // stop criteria by iteration
                HCF = 12.1 * Foam::sqrt(mag( U[cellI] ));
                HCN = 2.38 * Foam::pow( mag( 100 * XF - T[cellI] ), 0.25); // heat transf. coeff. by natural convection;
 
                if (HCF>HCN) { HC = HCF; } else { HC = HCN; }; // ok
                XN = (P5 + P4 * HC - P2 * Foam::pow(XF,4.0) ) / (100 + P3 * HC);
                if (N>150) { break; };
            } while (mag(XN-XF)>EPS);
 
 
            TCL = 100 * XN - 273; // surface temperature of clothing
            HL1 = 3.05 * 0.001 * (5733 - (6.99 * (met*58.15 - wme*58.15)) - PA); // heat loss diff. through skin
            if ( (met*58.15 - wme*58.15) > 58.15) { HL2 = 0.42 * ( (met*58.15 - wme*58.15) - 58.15); }
            else { HL2 = 0; }; // heat loss by sweating (comfort)
            HL3 = 1.7 * 0.00001 * met * 58.15 * (5867 - PA); // latent respiration heat loss
            HL4 = 0.0014 * met * 58.15 * (34 - (T[cellI] - 273.15) ); // dry respiration heat loss
            HL5 = 3.96 * FCL * (Foam::pow(XN,4) - Foam::pow( (((STemp+Troom)/2 + 273.0)/100),4)) ; // heat lose by radiation
            HL6 = FCL * HC * (TCL - (T[cellI] - 273.15)); // heat lose by convection
            TS = 0.303 * Foam::exp(-0.036 * met * 58.15) + 0.028; // thermal sensation trans coeff
 
            // Berechne PMV-Wert
            PMV[cellI] = TS * ( (met*58.15 - wme*58.15) - HL1 - HL2 - HL3 - HL4 - HL5 - HL6 );
 
//if (T[cellI]>290)
//{
//Info << " thermal sensation trans coedd == " << TS <<endl;
//Info << " wm == " << (met*58.15 - wme*58.15)  <<endl;
//Info << " TCL == " << TCL <<endl;
//Info << " TS == " << TS <<endl;
//Info << " XN == " << XN <<endl;
//Info << " Stemp == " << STemp <<endl;
//Info << " T == " << T[cellI] <<endl;
//Info << " HC == " << HC <<endl;
//Info << " FCL == " << FCL <<endl;
//Info << " HL1 == " << HL1 <<endl;
//Info << " HL2 == " << HL2 <<endl;
//Info << " HL3 == " << HL3 <<endl;
//Info << " HL4 == " << HL4 <<endl;
//Info << " HL5 == " << HL5 <<endl;
//Info << " HL6 == " << HL6 <<endl;
//Info << " PMV of this cell == " << PMV[cellI] <<endl;
//}
 
            // Berechne PPD-Wert
            PPD[cellI] = 100 - 95 * Foam::exp( -0.03353*pow(PMV[cellI],4) - 0.2179 * pow(PMV[cellI],2) );
 
         scalar new_lamda=0.24;
            if (PMV[cellI]>= 0)
            {
               new_lamda=0.24;
            }
            else
            {
               new_lamda=-0.5;
            }
            APMV[cellI] =PMV[cellI]/(1+new_lamda*PMV[cellI]);
 
            if (mag(U[cellI])<0.2) { a_korr=0.5; };
            if ( (mag(U[cellI])>=0.2) || (mag(U[cellI])<=0.6) ) { a_korr=0.6; };
            if (mag(U[cellI])>0.6) { a_korr=0.7; };
            // Berechne die operative Raumlufttemperatur
            TOp[cellI] = a_korr * T[cellI] + (1 - a_korr) * (STemp+273.15);
 
            midPMV += PMV[cellI];
            midPPD += PPD[cellI];
            midAPMV += APMV[cellI];
            midDR += DR[cellI];
            midTO += TOp[cellI];
 
            volume += mesh.V()[cellI];
            cellsum++;
 
 
        };
 
        DR.write();
        PMV.write();
        PPD.write();
        APMV.write();
        TOp.write();
        RH.write();
 
        Info<< "Mean Radiation temperature "<< STemp << " °C" << endl;
        Info<< "Water vapour pressure "<< PA << " Pa" <<nl;
        Info << "Average PMV-Value = " << (midPMV/cellsum) << nl;
        Info << "Average PPD-Value = " << (midPPD/cellsum) << " %" << nl;
        Info << "Average APMV-Value = " << (midAPMV/cellsum) << " %" << nl;
        Info << "Average DR-Value = " << (midDR/cellsum) << " %" << nl;
        Info<< "Average air velocity = "<< mag(wl) << " m/s" << nl;
        Info<< "Average room temperature = "<< mag(Tr)-273.15 << " °C" << nl;
        Info<< "Average relative room humidity = "<< mag(midRH/cellsum) << " %" << nl;
        Info << "Average operative temperature = " << (midTO/cellsum)-273.15 << " °C" << endl;
 
        if ( ((midPMV/cellsum > -0.2) || (midPMV/cellsum < 0.2)) && (midDR/cellsum < 10) && (midPPD/cellsum < 6))
        {
            Info << "Analysis: Category A" << endl;
        } else {
            if ( ((midPMV/cellsum > -0.5) || (midPMV/cellsum < 0.5)) && (midDR/cellsum < 20) && (midPPD/cellsum < 10))
            {
                Info << "Analysis: Category B" << endl;
            } else {
                if ( ((midPMV/cellsum > -0.7) || (midPMV/cellsum < 0.7)) && (midDR/cellsum < 30) && (midPPD/cellsum < 15))
                {
                    Info << "Analysis: Category C" << endl;
                } else
                {
                    Info << "Analysis: Category D" << endl;
                }}};
 
        // Info<< "Average Ux = "<< midUx / cellsum << endl;
        // Info<< "Average Uy = "<< midUy / cellsum << endl;
        // Info<< "Average Uz = "<< midUz / cellsum << endl;
    }
 
    Info<< "\nEnd\n" << endl;
 
    return(0);
}
 
 
// ************************************************************************* //