LiquidEvaporationBoil.C

Go to the documentation of this file.

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     | Website:  https://openfoam.org
    \\  /    A nd           | Copyright (C) 2011-2022 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 "LiquidEvaporationBoil.H"
#include "specie.H"
#include "mathematicalConstants.H"

using namespace Foam::constant::mathematical;

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

template<class CloudType>
Foam::tmp<Foam::scalarField> Foam::LiquidEvaporationBoil<CloudType>::calcXc
(
    const label celli
) const
{
    scalarField Xc(this->owner().composition().carrier().Y().size());

    forAll(Xc, i)
    {
        Xc[i] =
            this->owner().composition().carrier().Y()[i][celli]
           /this->owner().composition().carrier().Wi(i);
    }

    return Xc/sum(Xc);
}


template<class CloudType>
Foam::scalar Foam::LiquidEvaporationBoil<CloudType>::Sh
(
    const scalar Re,
    const scalar Sc
) const
{
    return 2.0 + 0.6*Foam::sqrt(Re)*cbrt(Sc);
}


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

template<class CloudType>
Foam::LiquidEvaporationBoil<CloudType>::LiquidEvaporationBoil
(
    const dictionary& dict,
    CloudType& owner
)
:
    PhaseChangeModel<CloudType>(dict, owner, typeName),
    liquids_(owner.thermo().liquids()),
    activeLiquids_(this->coeffDict().lookup("activeLiquids")),
    liqToCarrierMap_(activeLiquids_.size(), -1),
    liqToLiqMap_(activeLiquids_.size(), -1)
{
    if (activeLiquids_.size() == 0)
    {
        WarningInFunction
            << "Evaporation model selected, but no active liquids defined"
            << nl << endl;
    }
    else
    {
        Info<< "Participating liquid species:" << endl;

        // Determine mapping between liquid and carrier phase species
        forAll(activeLiquids_, i)
        {
            Info<< "    " << activeLiquids_[i] << endl;
            liqToCarrierMap_[i] =
                owner.composition().carrierId(activeLiquids_[i]);
        }

        // Determine mapping between model active liquids and global liquids
        const label idLiquid = owner.composition().idLiquid();
        forAll(activeLiquids_, i)
        {
            liqToLiqMap_[i] =
                owner.composition().localId(idLiquid, activeLiquids_[i]);
        }
    }
}


template<class CloudType>
Foam::LiquidEvaporationBoil<CloudType>::LiquidEvaporationBoil
(
    const LiquidEvaporationBoil<CloudType>& pcm
)
:
    PhaseChangeModel<CloudType>(pcm),
    liquids_(pcm.owner().thermo().liquids()),
    activeLiquids_(pcm.activeLiquids_),
    liqToCarrierMap_(pcm.liqToCarrierMap_),
    liqToLiqMap_(pcm.liqToLiqMap_)
{}


// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //

template<class CloudType>
Foam::LiquidEvaporationBoil<CloudType>::~LiquidEvaporationBoil()
{}


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

template<class CloudType>
void Foam::LiquidEvaporationBoil<CloudType>::calculate
(
    const typename CloudType::parcelType& p,
    const typename CloudType::parcelType::trackingData& td,
    const scalar dt,
    const scalar Re,
    const scalar Pr,
    const scalar d,
    const scalar nu,
    const scalar T,
    const scalar Ts,
    const scalar pc,
    const scalar Tc,
    const scalarField& X,
    scalarField& dMassPC
) const
{
    // immediately evaporate mass that has reached critical condition
    if ((liquids_.Tc(X) - T) < small)
    {
        if (debug)
        {
            WarningInFunction
                << "Parcel reached critical conditions: "
                << "evaporating all available mass" << endl;
        }

        forAll(activeLiquids_, i)
        {
            const label lid = liqToLiqMap_[i];
            dMassPC[lid] = great;
        }

        return;
    }

    // droplet surface pressure assumed to surface vapour pressure
    scalar ps = liquids_.pv(pc, Ts, X);

    // vapour density at droplet surface [kg/m^3]
    scalar rhos = ps*liquids_.W(X)/(RR*Ts);

    // construct carrier phase species volume fractions for cell, celli
    const scalarField XcMix(calcXc(p.cell()));

    // carrier thermo properties
    scalar Hsc = 0.0;
    scalar Hc = 0.0;
    scalar Cpc = 0.0;
    scalar kappac = 0.0;
    forAll(this->owner().composition().carrier().Y(), i)
    {
        scalar Yc = this->owner().composition().carrier().Y()[i][p.cell()];
        Hc += Yc*this->owner().composition().carrier().Ha(i, pc, Tc);
        Hsc += Yc*this->owner().composition().carrier().Ha(i, ps, Ts);
        Cpc += Yc*this->owner().composition().carrier().Cp(i, ps, Ts);
        kappac += Yc*this->owner().composition().carrier().kappa(i, ps, Ts);
    }

    // calculate mass transfer of each specie in liquid
    forAll(activeLiquids_, i)
    {
        const label gid = liqToCarrierMap_[i];
        const label lid = liqToLiqMap_[i];

        // boiling temperature at cell pressure for liquid species lid [K]
        const scalar TBoil = liquids_.properties()[lid].pvInvert(pc);

        // limit droplet temperature to boiling/critical temperature
        const scalar Td = min(T, 0.999*TBoil);

        // saturation pressure for liquid species lid [Pa]
        const scalar pSat = liquids_.properties()[lid].pv(pc, Td);

        // carrier phase concentration
        const scalar Xc = XcMix[gid];


        if (Xc*pc > pSat)
        {
            // saturated vapour - no phase change
        }
        else
        {
            // vapour diffusivity [m^2/s]
            const scalar Dab = liquids_.properties()[lid].D(ps, Ts);

            // Schmidt number
            const scalar Sc = nu/(Dab + rootVSmall);

            // Sherwood number
            const scalar Sh = this->Sh(Re, Sc);


            if (pSat > 0.999*pc)
            {
                // boiling

                const scalar deltaT = max(T - TBoil, 0.5);

                // vapour heat of formation
                const scalar hv = liquids_.properties()[lid].hl(pc, Td);

                // empirical heat transfer coefficient W/m2/K
                scalar alphaS = 0.0;
                if (deltaT < 5.0)
                {
                    alphaS = 760.0*pow(deltaT, 0.26);
                }
                else if (deltaT < 25.0)
                {
                    alphaS = 27.0*pow(deltaT, 2.33);
                }
                else
                {
                    alphaS = 13800.0*pow(deltaT, 0.39);
                }

                // flash-boil vaporisation rate
                const scalar Gf = alphaS*deltaT*pi*sqr(d)/hv;

                // model constants
                // NOTE: using Sherwood number instead of Nusselt number
                const scalar A = (Hc - Hsc)/hv;
                const scalar B = pi*kappac/Cpc*d*Sh;

                scalar G = 0.0;
                if (A > 0.0)
                {
                    // heat transfer from the surroundings contributes
                    // to the vaporisation process
                    scalar Gr = 1e-5;

                    for (label i=0; i<50; i++)
                    {
                        scalar GrDash = Gr;

                        G = B/(1.0 + Gr)*log(1.0 + A*(1.0 + Gr));
                        Gr = Gf/G;

                        if (mag(Gr - GrDash)/GrDash < 1e-3)
                        {
                            break;
                        }
                    }
                }

                dMassPC[lid] += (G + Gf)*dt;
            }
            else
            {
                // evaporation

                // surface molar fraction - Raoult's Law
                const scalar Xs = X[lid]*pSat/pc;

                // molar ratio
                const scalar Xr = (Xs - Xc)/max(small, 1.0 - Xs);

                if (Xr > 0)
                {
                    // mass transfer [kg]
                    dMassPC[lid] += pi*d*Sh*Dab*rhos*log(1.0 + Xr)*dt;
                }
            }
        }
    }
}


template<class CloudType>
Foam::scalar Foam::LiquidEvaporationBoil<CloudType>::dh
(
    const label idc,
    const label idl,
    const scalar p,
    const scalar T
) const
{
    scalar dh = 0;

    scalar TDash = T;
    if (liquids_.properties()[idl].pv(p, T) >= 0.999*p)
    {
        TDash = liquids_.properties()[idl].pvInvert(p);
    }

    typedef PhaseChangeModel<CloudType> parent;
    switch (parent::enthalpyTransfer_)
    {
        case (parent::etLatentHeat):
        {
            dh = liquids_.properties()[idl].hl(p, TDash);
            break;
        }
        case (parent::etEnthalpyDifference):
        {
            scalar hc = this->owner().composition().carrier().Ha(idc, p, TDash);
            scalar hp = liquids_.properties()[idl].Ha(p, TDash);

            dh = hc - hp;
            break;
        }
        default:
        {
            FatalErrorInFunction
                << "Unknown enthalpyTransfer type" << abort(FatalError);
        }
    }

    return dh;
}


template<class CloudType>
Foam::scalar Foam::LiquidEvaporationBoil<CloudType>::Tvap
(
    const scalarField& X
) const
{
    return liquids_.Tpt(X);
}


template<class CloudType>
Foam::scalar Foam::LiquidEvaporationBoil<CloudType>::TMax
(
    const scalar p,
    const scalarField& X
) const
{
    return liquids_.pvInvert(p, X);
}


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