View Issue Details

IDProjectCategoryView StatusLast Update
0000938OpenFOAMBugpublic2013-08-07 12:34
Reporterdkxls Assigned Touser2 
PrioritynormalSeveritymajorReproducibilityalways
Status resolvedResolutionfixed 
PlatformOpenSUSE 12.2OSLinuxOS Versionx86_64
Summary0000938: Wrong boiling temperature from Foam::liquidProperties::pvInvert(scalar p)
DescriptionThe boiling temperature can only be defined between the triple point and the critical point.

The pvInvert function in OpenFOAM does not obey these basic thermodynamic properties and hence returns wrong values if the pressure is higher than Pcrit or lower than Ptrip.

The code contains additionally minor flaws besides the critical issues (e.g. the unnecessarily defined variable 'pBoil'). Also the description in the .H file is wrong and should be updated.
Steps To ReproduceCall the function Foam::liquidProperties::pvInvert(scalar p) with values higher than the Pcrit or lower than Ptrip.
A utility that writes out the liquid properties for a mixture and the components is attached (for OpenFOAM version 2.2.x).
Additional InformationpvInvert is used in Lagrangian solver and surface film models, which should produce wrong results due to this issue:

lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/LiquidEvaporationBoil/LiquidEvaporationBoil.C

regionModels/surfaceFilmModels/submodels/thermo/phaseChangeModel/standardPhaseChange/standardPhaseChange.C
Tagsboiling, Evaporation, liquidProperties

Activities

dkxls

2013-08-01 21:01

reporter  

dkxls

2013-08-05 11:41

reporter   ~0002365

Last edited: 2013-08-05 11:56

Please find attached a reimplementation of the 'pvInvert' function that fixes this bug (I can also provide a patch-file if preferred).

I would also suggest to rename the function 'pvInvert' to 'Tboil' and the variables/functions 'Tb'/'Tb_' to 'TboilNorm'/'TboilNorm_'.
This is not implemented in the attached file, but would clarify the source code significantly!

dkxls

2013-08-05 11:42

reporter  

liquidProperties.C (9,780 bytes)   
/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
     \\/     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 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 "liquidProperties.H"
#include "HashTable.H"
#include "Switch.H"

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

namespace Foam
{
    defineTypeNameAndDebug(liquidProperties, 0);
    defineRunTimeSelectionTable(liquidProperties,);
    defineRunTimeSelectionTable(liquidProperties, Istream);
    defineRunTimeSelectionTable(liquidProperties, dictionary);
}

// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //

Foam::liquidProperties::liquidProperties
(
    scalar W,
    scalar Tc,
    scalar Pc,
    scalar Vc,
    scalar Zc,
    scalar Tt,
    scalar Pt,
    scalar Tb,
    scalar dipm,
    scalar omega,
    scalar delta
)
:
    W_(W),
    Tc_(Tc),
    Pc_(Pc),
    Vc_(Vc),
    Zc_(Zc),
    Tt_(Tt),
    Pt_(Pt),
    Tb_(Tb),
    dipm_(dipm),
    omega_(omega),
    delta_(delta)
{}


Foam::liquidProperties::liquidProperties(Istream& is)
:
    W_(readScalar(is)),
    Tc_(readScalar(is)),
    Pc_(readScalar(is)),
    Vc_(readScalar(is)),
    Zc_(readScalar(is)),
    Tt_(readScalar(is)),
    Pt_(readScalar(is)),
    Tb_(readScalar(is)),
    dipm_(readScalar(is)),
    omega_(readScalar(is)),
    delta_(readScalar(is))
{}


Foam::liquidProperties::liquidProperties(const dictionary& dict)
:
    W_(readScalar(dict.lookup("W"))),
    Tc_(readScalar(dict.lookup("Tc"))),
    Pc_(readScalar(dict.lookup("Pc"))),
    Vc_(readScalar(dict.lookup("Vc"))),
    Zc_(readScalar(dict.lookup("Zc"))),
    Tt_(readScalar(dict.lookup("Tt"))),
    Pt_(readScalar(dict.lookup("Pt"))),
    Tb_(readScalar(dict.lookup("Tb"))),
    dipm_(readScalar(dict.lookup("dipm"))),
    omega_(readScalar(dict.lookup("omega"))),
    delta_(readScalar(dict.lookup("delta")))
{}


Foam::liquidProperties::liquidProperties(const liquidProperties& liq)
:
    W_(liq.W_),
    Tc_(liq.Tc_),
    Pc_(liq.Pc_),
    Vc_(liq.Vc_),
    Zc_(liq.Zc_),
    Tt_(liq.Tt_),
    Pt_(liq.Pt_),
    Tb_(liq.Tb_),
    dipm_(liq.dipm_),
    omega_(liq.omega_),
    delta_(liq.delta_)
{}


// * * * * * * * * * * * * * * * * * Selectors * * * * * * * * * * * * * * * //

Foam::autoPtr<Foam::liquidProperties> Foam::liquidProperties::New(Istream& is)
{
    if (debug)
    {
        Info<< "liquidProperties::New(Istream&): "
            << "constructing liquidProperties" << endl;
    }

    const word liquidPropertiesType(is);
    const word coeffs(is);

    if (coeffs == "defaultCoeffs")
    {
        ConstructorTable::iterator cstrIter =
            ConstructorTablePtr_->find(liquidPropertiesType);

        if (cstrIter == ConstructorTablePtr_->end())
        {
            FatalErrorIn("liquidProperties::New(Istream&)")
                << "Unknown liquidProperties type "
                << liquidPropertiesType << nl << nl
                << "Valid liquidProperties types are:" << nl
                << ConstructorTablePtr_->sortedToc()
                << abort(FatalError);
        }

        return autoPtr<liquidProperties>(cstrIter()());
    }
    else if (coeffs == "coeffs")
    {
        return autoPtr<liquidProperties>(new liquidProperties(is));
    }
    else
    {
        FatalErrorIn("liquidProperties::New(Istream&)")
            << "liquidProperties type " << liquidPropertiesType
            << ", option " << coeffs << " given"
            << ", should be coeffs or defaultCoeffs"
            << abort(FatalError);

        return autoPtr<liquidProperties>(NULL);
    }
}


Foam::autoPtr<Foam::liquidProperties> Foam::liquidProperties::New
(
    const dictionary& dict
)
{
    if (debug)
    {
        Info<< "liquidProperties::New(const dictionary&):"
            << "constructing liquidProperties" << endl;
    }

    const word& liquidPropertiesTypeName = dict.dictName();

    const Switch defaultCoeffs(dict.lookup("defaultCoeffs"));

    if (defaultCoeffs)
    {
        ConstructorTable::iterator cstrIter =
            ConstructorTablePtr_->find(liquidPropertiesTypeName);

        if (cstrIter == ConstructorTablePtr_->end())
        {
            FatalErrorIn
            (
                "liquidProperties::New(const dictionary&, const word&)"
            )   << "Unknown liquidProperties type "
                << liquidPropertiesTypeName << nl << nl
                << "Valid liquidProperties types are:" << nl
                << ConstructorTablePtr_->sortedToc()
                << abort(FatalError);
        }

        return autoPtr<liquidProperties>(cstrIter()());
    }
    else
    {
        return autoPtr<liquidProperties>
        (
            new liquidProperties
            (
                dict.subDict(liquidPropertiesTypeName + "Coeffs")
            )
        );
    }
}


// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //

Foam::scalar Foam::liquidProperties::rho(scalar p, scalar T) const
{
    notImplemented
    (
        "Foam::scalar Foam::liquidProperties::rho(scalar, scalar) const"
    );
    return 0.0;
}


Foam::scalar Foam::liquidProperties::pv(scalar p, scalar T) const
{
    notImplemented
    (
        "Foam::scalar Foam::liquidProperties::pv(scalar, scalar) const"
    );
    return 0.0;
}


Foam::scalar Foam::liquidProperties::hl(scalar p, scalar T) const
{
    notImplemented
    (
        "Foam::scalar Foam::liquidProperties::hl(scalar, scalar) const"
    );
    return 0.0;
}


Foam::scalar Foam::liquidProperties::Cp(scalar p, scalar T) const
{
    notImplemented
    (
        "Foam::scalar Foam::liquidProperties::Cp(scalar, scalar) const"
    );
    return 0.0;
}


Foam::scalar Foam::liquidProperties::h(scalar p, scalar T) const
{
    notImplemented
    (
        "Foam::scalar Foam::liquidProperties::h(scalar, scalar) const"
    );
    return 0.0;
}


Foam::scalar Foam::liquidProperties::Cpg(scalar p, scalar T) const
{
    notImplemented
    (
        "Foam::scalar Foam::liquidProperties::Cpg(scalar, scalar) const"
    );
    return 0.0;
}


Foam::scalar Foam::liquidProperties::mu(scalar p, scalar T) const
{
    notImplemented
    (
        "Foam::scalar Foam::liquidProperties::mu(scalar, scalar) const"
    );
    return 0.0;
}


Foam::scalar Foam::liquidProperties::mug(scalar p, scalar T) const
{
    notImplemented
    (
        "Foam::scalar Foam::liquidProperties::mug(scalar, scalar) const"
    );
    return 0.0;
}


Foam::scalar Foam::liquidProperties::K(scalar p, scalar T) const
{
    notImplemented
    (
        "Foam::scalar Foam::liquidProperties::K(scalar, scalar) const"
    );
    return 0.0;
}


Foam::scalar Foam::liquidProperties::Kg(scalar p, scalar T) const
{
    notImplemented
    (
        "Foam::scalar Foam::liquidProperties::Kg(scalar, scalar) const"
    );
    return 0.0;
}


Foam::scalar Foam::liquidProperties::sigma(scalar p, scalar T) const
{
    notImplemented
    (
        "Foam::scalar Foam::liquidProperties::sigms(scalar, scalar) const"
    );
    return 0.0;
}


Foam::scalar Foam::liquidProperties::D(scalar p, scalar T) const
{
    notImplemented
    (
        "Foam::scalar Foam::liquidProperties::D(scalar, scalar) const"
    );
    return 0.0;
}


Foam::scalar Foam::liquidProperties::D(scalar p, scalar T, scalar Wb) const
{
    notImplemented
    (
        "Foam::scalar Foam::liquidProperties::D(scalar, scalar) const"
    );
    return 0.0;
}


Foam::scalar Foam::liquidProperties::pvInvert(scalar p) const
{
    // Check for critical and solid phase conditions
    if (p >= Pc_)
    {
        return Tc_;
    }
    else if (p < Pt_)
    {
        if (debug)
        {
            WarningIn
            (
                "Foam::scalar Foam::liquidProperties::pvInvert(scalar p) const"
            )   << "Pressure below triple point pressure: "
                << "p = "<<p<<" < Pt = "<<Pt_ <<  nl << endl;
        }
        return -1;
    }

    // Set initial upper and lower bounds
    scalar Thi = Tc_;
    scalar Tlo = Tt_;

    // Boiling temperature at normal conditions seems a good initial guess
    scalar T = Tb_;

    while ((Thi - Tlo) > 1.0e-4)
    {
        if ((pv(p, T) - p) <= 0.0)
        {
            Tlo = T;
        }
        else
        {
            Thi = T;
        }
        T = (Thi + Tlo)*0.5;
    }

    return T;
}


void Foam::liquidProperties::writeData(Ostream& os) const
{

    os  << W_ << token::SPACE
        << Tc_ << token::SPACE
        << Pc_ << token::SPACE
        << Vc_ << token::SPACE
        << Zc_ << token::SPACE
        << Tt_ << token::SPACE
        << Pt_ << token::SPACE
        << Tb_ << token::SPACE
        << dipm_ << token::SPACE
        << omega_<< token::SPACE
        << delta_;
}


// ************************************************************************* //
liquidProperties.C (9,780 bytes)   

user2

2013-08-06 11:36

  ~0002369

Thanks for the report and suggested changes - resolved by commit 21c0050

Issue History

Date Modified Username Field Change
2013-08-01 21:01 dkxls New Issue
2013-08-01 21:01 dkxls File Added: liqMixProperties.tar.gz
2013-08-05 11:41 dkxls Note Added: 0002365
2013-08-05 11:42 dkxls File Added: liquidProperties.C
2013-08-05 11:56 dkxls Note Edited: 0002365
2013-08-06 11:35 user2 Priority immediate => normal
2013-08-06 11:35 user2 Severity crash => major
2013-08-06 11:36 user2 Note Added: 0002369
2013-08-06 11:36 user2 Status new => resolved
2013-08-06 11:36 user2 Fixed in Version => 2.2.x
2013-08-06 11:36 user2 Resolution open => fixed
2013-08-06 11:36 user2 Assigned To => user2
2013-08-07 12:34 dkxls Tag Attached: Evaporation
2013-08-07 12:34 dkxls Tag Attached: boiling
2013-08-07 12:34 dkxls Tag Attached: liquidProperties