View Issue Details

Jump to NotesJump to History
IDProjectCategoryView StatusDate SubmittedLast Update
0000942OpenFOAMBugpublic2013-08-05 09:472013-08-09 13:12
Reporterdkxls Assigned Touser2 
PrioritynormalSeveritymajorReproducibilityalways
Status resolvedResolutionfixed 
PlatformOpenSUSE 12.2OSLinuxOS Versionx86_64
Summary0000942: [SprayParcel]: No limiter for max. parcel temperature to critical temperature
DescriptionThe current implementation of SprayParcel (and hence also sprayFoam or sprayEngineFoam) does not limit the maximum parcel temperature to the critical temperature.
This leads for most spray (combustion) simulation to wrong results and/or crashes.

This bug is related to bugs #938 and #939, but requires additional changes in the PhaseChange classes.
Although this bug is relevant to liquids under boiling conditions, the most sever effects will be experienced in non-boiling conditions, where (compressed) liquid fuel reaches super-critical conditions (P > Pcrit).
Steps To ReproduceThe "Spray A" test case, specified by the Engine Combustion Network (ECN), resembles typical engine spray conditions and leads to reproducible crashes.
http://www.sandia.gov/ecn/cvdata/targetCondition/sprayA.php
Additional Informationhttp://www.openfoam.org/mantisbt/view.php?id=938
http://www.openfoam.org/mantisbt/view.php?id=939

A correct implementation is given in the dieselSpray classes:
https://github.com/OpenFOAM/OpenFOAM-2.0.x/blob/master/src/lagrangian/dieselSpray/parcel/parcel.C
https://github.com/OpenFOAM/OpenFOAM-2.0.x/blob/master/src/lagrangian/dieselSpray/parcel/setRelaxationTimes.C
TagsEvaporation, Lagrangian, spray, sprayEngineFoam, sprayFoam

Activities

user2

2013-08-06 11:37

  ~0002370

Thanks for the report - fixed by commit f3bdac5

dkxls

2013-08-06 17:36

reporter   ~0002374

Thanks for taking care of this. I still think this is not yet done correctly. The current implementation assumes that a parcel consists of independent liquid components instead of a liquid mixture!
I don't think this is what we want, or at least I think it's not correct for liquid fuel sprays.

A correct way IMO of doing it, was implemented in the dieselSpray class (see the two links earlier). The important part is that the we consider mixing quantities for max temperature considerations.
Here the steps:
1) Limit the temperature to the critical temperature for the mixture (do this in any case)
2) Check if the mixture is at boiling conditions; if true limit to the boiling temperature for the mixture

So instead of having a 'TMax' in 'liquidMixtureProperties' we need to have 'pvInvert(p,X)' for the mixture quantities. I implemented and tested this in the attached files (liquidMixtureProperties.*). I also implemented a pseudo triple point temperature Tpt(X), implementation is done the same way as Tpc(X). I needed Tpt(X) for pvInvert(p,X), but it can also be used to fix bug #945. Again this is tested and works!

Also attached is a new SprayParcel.C with the required changes. However, I didn't test the SprayParcel code at all, but idea should be clear! (I also cleaned up the code a bit, would be nice if that could be included as well, which would also fix bug #946).

I think that having a function 'TMax' in the PhaseChange model is not only unnecessary, but also confusion. 'TMax' should be removed entirely from the PhaseChange models (as from liquidMixtureProperties) and 'TMax(p)' in ReactingParcel.C should be replaced by td.cloud().constProps().TMax(). Of course, given td.cloud().constProps().TMax() is set correctly as described above.

dkxls

2013-08-06 17:36

reporter  

liquidMixtureProperties.C (9,912 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 "liquidMixtureProperties.H"
#include "dictionary.H"
#include "specie.H"

// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //

const Foam::scalar Foam::liquidMixtureProperties::TrMax = 0.999;


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

Foam::liquidMixtureProperties::liquidMixtureProperties
(
    const dictionary& dict
)
:
    components_(),
    properties_()
{
    components_ = dict.toc();
    properties_.setSize(components_.size());

    forAll(components_, i)
    {
        properties_.set
        (
            i,
            liquidProperties::New(dict.subDict(components_[i]))
        );
    }
}


Foam::liquidMixtureProperties::liquidMixtureProperties
(
    const liquidMixtureProperties& lm
)
:
    components_(lm.components_),
    properties_(lm.properties_.size())
{
    forAll(properties_, i)
    {
        properties_.set(i, lm.properties_(i)->clone());
    }
}


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

Foam::autoPtr<Foam::liquidMixtureProperties> Foam::liquidMixtureProperties::New
(
    const dictionary& thermophysicalProperties
)
{
    return autoPtr<liquidMixtureProperties>
    (
        new liquidMixtureProperties(thermophysicalProperties)
    );
}


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

Foam::scalar Foam::liquidMixtureProperties::Tc(const scalarField& x) const
{
    scalar vTc = 0.0;
    scalar vc = 0.0;

    forAll(properties_, i)
    {
        scalar x1 = x[i]*properties_[i].Vc();
        vc += x1;
        vTc += x1*properties_[i].Tc();
    }

    return vTc/vc;
}


Foam::scalar Foam::liquidMixtureProperties::Tpt(const scalarField& x) const
{
    scalar Tpt = 0.0;
    forAll(properties_, i)
    {
        Tpt += x[i]*properties_[i].Tt();
    }

    return Tpt;
}


Foam::scalar Foam::liquidMixtureProperties::pvInvert
(
    const scalar p,
    const scalarField& x
) const
{
    // Set upper and lower bounds and initial guess
    scalar Thi = Tc(x);
    scalar Tlo = Tpt(x);

    // Check for critical and solid phase conditions
    if (p >= pv(p, Thi, x))
    {
        return Thi;
    }
    else if (p < pv(p, Tlo, x))
    {
        WarningIn
        (
            "Foam::scalar Foam::liquidMixtureProperties::pvInvert"
            "("
            "    const scalar,"
            "    const scalarField&"
            ") const"
        )   << "Pressure below triple point pressure: "
            << "p = " << p << " < Pt = " << pv(p, Tlo, x) <<  nl << endl;
        return -1;
    }

    // Set initial guess
    scalar T = (Thi + Tlo)*0.5;

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

        T = (Thi + Tlo)*0.5;
    }

    return T;
}


Foam::scalar Foam::liquidMixtureProperties::Tpc(const scalarField& x) const
{
    scalar Tpc = 0.0;
    forAll(properties_, i)
    {
        Tpc += x[i]*properties_[i].Tc();
    }

    return Tpc;
}


Foam::scalar Foam::liquidMixtureProperties::Ppc(const scalarField& x) const
{
    scalar Vc = 0.0;
    scalar Zc = 0.0;
    forAll(properties_, i)
    {
        Vc += x[i]*properties_[i].Vc();
        Zc += x[i]*properties_[i].Zc();
    }

    return specie::RR*Zc*Tpc(x)/Vc;
}


Foam::scalar Foam::liquidMixtureProperties::omega(const scalarField& x) const
{
    scalar omega = 0.0;
    forAll(properties_, i)
    {
        omega += x[i]*properties_[i].omega();
    }

    return omega;
}


Foam::scalarField Foam::liquidMixtureProperties::Xs
(
    const scalar p,
    const scalar Tg,
    const scalar Tl,
    const scalarField& xg,
    const scalarField& xl
) const
{
    scalarField xs(xl.size(), 0.0);

    // Raoult's Law
    forAll(xs, i)
    {
        scalar Ti = min(TrMax*properties_[i].Tc(), Tl);
        xs[i] = properties_[i].pv(p, Ti)*xl[i]/p;
    }
    return xs;
}


Foam::scalar Foam::liquidMixtureProperties::W(const scalarField& x) const
{
    scalar W = 0.0;
    forAll(properties_, i)
    {
        W += x[i]*properties_[i].W();
    }

    return W;
}


Foam::scalarField Foam::liquidMixtureProperties::Y(const scalarField& X) const
{
    scalarField Y(X/W(X));

    forAll(Y, i)
    {
        Y[i] *= properties_[i].W();
    }

    return Y;
}


Foam::scalarField Foam::liquidMixtureProperties::X(const scalarField& Y) const
{
    scalarField X(Y.size());
    scalar Winv = 0.0;
    forAll(X, i)
    {
        Winv += Y[i]/properties_[i].W();
        X[i] = Y[i]/properties_[i].W();
    }

    tmp<scalarField> tfld = X/Winv;
    return tfld();
}


Foam::scalar Foam::liquidMixtureProperties::rho
(
    const scalar p,
    const scalar T,
    const scalarField& x
) const
{
    scalar v = 0.0;

    forAll(properties_, i)
    {
        if (x[i] > SMALL)
        {
            scalar Ti = min(TrMax*properties_[i].Tc(), T);
            scalar rho = SMALL + properties_[i].rho(p, Ti);
            v += x[i]*properties_[i].W()/rho;
        }
    }

    return W(x)/v;
}


Foam::scalar Foam::liquidMixtureProperties::pv
(
    const scalar p,
    const scalar T,
    const scalarField& x
) const
{
    scalar pv = 0.0;

    forAll(properties_, i)
    {
        if (x[i] > SMALL)
        {
            scalar Ti = min(TrMax*properties_[i].Tc(), T);
            pv += x[i]*properties_[i].pv(p, Ti)*properties_[i].W();
        }
    }

    return pv/W(x);
}


Foam::scalar Foam::liquidMixtureProperties::hl
(
    const scalar p,
    const scalar T,
    const scalarField& x
) const
{
    scalar hl = 0.0;

    forAll(properties_, i)
    {
        if (x[i] > SMALL)
        {
            scalar Ti = min(TrMax*properties_[i].Tc(), T);
            hl += x[i]*properties_[i].hl(p, Ti)*properties_[i].W();
        }
    }

    return hl/W(x);
}


Foam::scalar Foam::liquidMixtureProperties::Cp
(
    const scalar p,
    const scalar T,
    const scalarField& x
) const
{
    scalar Cp = 0.0;

    forAll(properties_, i)
    {
        if (x[i] > SMALL)
        {
            scalar Ti = min(TrMax*properties_[i].Tc(), T);
            Cp += x[i]*properties_[i].Cp(p, Ti)*properties_[i].W();
        }
    }

    return Cp/W(x);
}


Foam::scalar Foam::liquidMixtureProperties::sigma
(
    const scalar p,
    const scalar T,
    const scalarField& x
) const
{
    // sigma is based on surface mole fractions
    // which is estimated from Raoult's Law
    scalar sigma = 0.0;
    scalarField Xs(x.size(), 0.0);
    scalar XsSum = 0.0;
    forAll(properties_, i)
    {
        scalar Ti = min(TrMax*properties_[i].Tc(), T);
        scalar Pvs = properties_[i].pv(p, Ti);
        scalar xs = x[i]*Pvs/p;
        XsSum += xs;
        Xs[i] = xs;
    }

    forAll(properties_, i)
    {
        if (Xs[i] > SMALL)
        {
            scalar Ti = min(TrMax*properties_[i].Tc(), T);
            sigma += (Xs[i]/XsSum)*properties_[i].sigma(p, Ti);
        }
    }

    return sigma;
}


Foam::scalar Foam::liquidMixtureProperties::mu
(
    const scalar p,
    const scalar T,
    const scalarField& x
) const
{
    scalar mu = 0.0;

    forAll(properties_, i)
    {
        if (x[i] > SMALL)
        {
            scalar Ti = min(TrMax*properties_[i].Tc(), T);
            mu += x[i]*log(properties_[i].mu(p, Ti));
        }
    }

    return exp(mu);
}


Foam::scalar Foam::liquidMixtureProperties::K
(
    const scalar p,
    const scalar T,
    const scalarField& x
) const
{
    // calculate superficial volume fractions phii
    scalarField phii(x.size(), 0.0);
    scalar pSum = 0.0;

    forAll(properties_, i)
    {
        scalar Ti = min(TrMax*properties_[i].Tc(), T);

        scalar Vi = properties_[i].W()/properties_[i].rho(p, Ti);
        phii[i] = x[i]*Vi;
        pSum += phii[i];
    }

    forAll(phii, i)
    {
        phii[i] /= pSum;
    }

    scalar K = 0.0;

    forAll(properties_, i)
    {
        scalar Ti = min(TrMax*properties_[i].Tc(), T);

        forAll(properties_, j)
        {
            scalar Tj = min(TrMax*properties_[j].Tc(), T);

            scalar Kij =
                2.0
               /(
                    1.0/properties_[i].K(p, Ti)
                  + 1.0/properties_[j].K(p, Tj)
                );
            K += phii[i]*phii[j]*Kij;
        }
    }

    return K;
}


Foam::scalar Foam::liquidMixtureProperties::D
(
    const scalar p,
    const scalar T,
    const scalarField& x
) const
{
    // Blanc's law
    scalar Dinv = 0.0;

    forAll(properties_, i)
    {
        if (x[i] > SMALL)
        {
            scalar Ti = min(TrMax*properties_[i].Tc(), T);
            Dinv += x[i]/properties_[i].D(p, Ti);
        }
    }

    return 1.0/Dinv;
}


// ************************************************************************* //
liquidMixtureProperties.C (9,912 bytes)   

dkxls

2013-08-06 17:36

reporter  

liquidMixtureProperties.H (6,887 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/>.

Class
    Foam::liquidMixtureProperties

Description
    A mixture of liquids

    An example of a two component liquid mixture:
    \verbatim
        <parentDictionary>
        {
            H2O
            {
                defaultCoeffs   yes;     // employ default coefficients
            }
            C7H16
            {
                defaultCoeffs   no;
                C7H16Coeffs
                {
                    ... user defined properties for C7H16
                }
            }
        }
    \endverbatim

SourceFiles
    liquidMixtureProperties.C

SeeAlso
    Foam::liquidProperties

\*---------------------------------------------------------------------------*/

#ifndef liquidMixtureProperties_H
#define liquidMixtureProperties_H

#include "word.H"
#include "scalarField.H"
#include "PtrList.H"
#include "liquidProperties.H"
#include "autoPtr.H"

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

namespace Foam
{

/*---------------------------------------------------------------------------*\
                Class liquidMixtureProperties Declaration
\*---------------------------------------------------------------------------*/

class liquidMixtureProperties
{
    // Private data

        //- Maximum reduced temperature
        static const scalar TrMax;

        //- The names of the liquids
        List<word> components_;

        //- The liquid properties
        PtrList<liquidProperties> properties_;


public:


    // Constructors

        //- Construct from dictionary
        liquidMixtureProperties(const dictionary& dict);

        //- Construct copy
        liquidMixtureProperties(const liquidMixtureProperties& lm);

        //- Construct and return a clone
        virtual autoPtr<liquidMixtureProperties> clone() const
        {
            return autoPtr<liquidMixtureProperties>
            (
                new liquidMixtureProperties(*this)
            );
        }


    //- Destructor
    virtual ~liquidMixtureProperties()
    {}


    // Selectors

        //- Select construct from dictionary
        static autoPtr<liquidMixtureProperties> New(const dictionary&);


    // Member Functions

        //- Return the liquid names
        inline const List<word>& components() const
        {
            return components_;
        }

        //- Return the liquid properties
        inline const PtrList<liquidProperties>& properties() const
        {
            return properties_;
        }

        //- Return the number of liquids in the mixture
        inline label size() const
        {
            return components_.size();
        }

        //- Calculate the critical temperature of mixture
        scalar Tc(const scalarField& x) const;

        //- Invert the vapour pressure relationship to retrieve the boiling 
        //  temperature of the mixture as a function of pressure
        scalar pvInvert(const scalar p, const scalarField& x) const;

        //- Return pseudocritical temperature according to Kay's rule
        scalar Tpc(const scalarField& x) const;

        //- Return pseudocritical pressure (modified Prausnitz and Gunn)
        scalar Ppc(const scalarField& x) const;

        //- Return pseudo triple point temperature (mole averaged formulation)
        scalar Tpt(const scalarField& x) const;

        //- Return mixture accentric factor
        scalar omega(const scalarField& x) const;

        //- Return the surface molar fractions
        scalarField Xs
        (
            const scalar p,
            const scalar Tg,
            const scalar Tl,
            const scalarField& xg,
            const scalarField& xl
        ) const;


        //- Calculate the mean molecular weight [kg/kmol]
        // from mole fractions
        scalar W(const scalarField& x) const;

        //- Returns the mass fractions, given mole fractions
        scalarField Y(const scalarField& X) const;

        //- Returns the mole fractions, given mass fractions
        scalarField X(const scalarField& Y) const;

        //- Calculate the mixture density [kg/m^3]
        scalar rho
        (
            const scalar p,
            const scalar T,
            const scalarField& x
        ) const;

        //- Calculate the mixture vapour pressure [Pa]
        scalar pv
        (
            const scalar p,
            const scalar T,
            const scalarField& x
        ) const;

        //- Calculate the mixture latent heat [J/kg]
        scalar hl
        (
            const scalar p,
            const scalar T,
            const scalarField& x
        ) const;

        //- Calculate the mixture heat capacity [J/(kg K)]
        scalar Cp
        (
            const scalar p,
            const scalar T,
            const scalarField& x
        ) const;

        //- Estimate mixture surface tension [N/m]
        scalar sigma
        (
            const scalar p,
            const scalar T,
            const scalarField& x
        ) const;

        //- Calculate the mixture viscosity [Pa s]
        scalar mu
        (
            const scalar p,
            const scalar T,
            const scalarField& x
        ) const;

        //- Estimate thermal conductivity  [W/(m K)]
        // Li's method, Eq. 10-12.27 - 10.12-19
        scalar K
        (
            const scalar p,
            const scalar T,
            const scalarField& x
        ) const;

        scalar D
        (
            const scalar p,
            const scalar T,
            const scalarField& x
        ) const;
};


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

} // End namespace Foam

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

#endif

// ************************************************************************* //
liquidMixtureProperties.H (6,887 bytes)   

dkxls

2013-08-06 17:37

reporter  

SprayParcel.C (13,830 bytes)    
/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright (C) 2011-2013 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 "SprayParcel.H"
#include "CompositionModel.H"
#include "AtomizationModel.H"

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

template<class ParcelType>
template<class TrackData>
void Foam::SprayParcel<ParcelType>::setCellValues
(
    TrackData& td,
    const scalar dt,
    const label cellI
)
{
    ParcelType::setCellValues(td, dt, cellI);
}


template<class ParcelType>
template<class TrackData>
void Foam::SprayParcel<ParcelType>::cellValueSourceCorrection
(
    TrackData& td,
    const scalar dt,
    const label cellI
)
{
    ParcelType::cellValueSourceCorrection(td, dt, cellI);
}


template<class ParcelType>
template<class TrackData>
void Foam::SprayParcel<ParcelType>::calc
(
    TrackData& td,
    const scalar dt,
    const label cellI
)
{
    typedef typename TrackData::cloudType::reactingCloudType reactingCloudType;
    const CompositionModel<reactingCloudType>& composition =
        td.cloud().composition();

    const bool coupled = td.cloud().solution().coupled();

    // check if parcel belongs to liquid core
    if (liquidCore() > 0.5)
    {
        // liquid core parcels should not interact with the gas
        td.cloud().solution().coupled() = false;
    }

    // get old mixture composition
    const scalarField& Y0(this->Y());
    scalarField X0(composition.liquids().X(Y0));

    // check if we have critical or boiling conditions
    scalar TMax = composition.liquids().Tc(X0);
    const scalar T0 = this->T();
    const scalar pc0 = this->pc_;
    if (composition.liquids().pv(pc0, T0, X0) >= pc0*0.999)
    {
        // set TMax to boiling temperature
        TMax = composition.liquids().pvInvert(pc0,X0);
    }

    // set the maximum temperature limit
    td.cloud().constProps().TMax() = TMax;

    this->Cp() = composition.liquids().Cp(pc0, T0, X0);
    scalar rho0 = composition.liquids().rho(pc0, T0, X0);
    this->rho() = rho0;

    ParcelType::calc(td, dt, cellI);

    if (td.keepParticle)
    {
        // update Cp, diameter and density due to change in temperature
        // and/or composition
        scalar T1 = this->T();
        const scalarField& Y1(this->Y());
        scalarField X1(composition.liquids().X(Y1));

        this->Cp() = composition.liquids().Cp(this->pc_, T1, X1);

        scalar rho1 = composition.liquids().rho(this->pc_, T1, X1);
        this->rho() = rho1;
        scalar d1 = this->d()*cbrt(rho0/rho1);
        this->d() = d1;

        if (liquidCore() > 0.5)
        {
            calcAtomization(td, dt, cellI);

            // preserve the total mass/volume by increasing the number of
            // particles in parcels due to breakup
            scalar d2 = this->d();
            this->nParticle() *= pow3(d1/d2);
        }
        else
        {
            calcBreakup(td, dt, cellI);
        }
    }

    // restore coupled
    td.cloud().solution().coupled() = coupled;
}


template<class ParcelType>
template<class TrackData>
void Foam::SprayParcel<ParcelType>::calcAtomization
(
    TrackData& td,
    const scalar dt,
    const label cellI
)
{
    typedef typename TrackData::cloudType::reactingCloudType reactingCloudType;
    const CompositionModel<reactingCloudType>& composition =
        td.cloud().composition();

    typedef typename TrackData::cloudType::sprayCloudType sprayCloudType;
    const AtomizationModel<sprayCloudType>& atomization =
        td.cloud().atomization();


    // cell state info is updated in ReactingParcel calc
    const scalarField& Y(this->Y());
    scalarField X(composition.liquids().X(Y));

    scalar rho = composition.liquids().rho(this->pc(), this->T(), X);
    scalar mu = composition.liquids().mu(this->pc(), this->T(), X);
    scalar sigma = composition.liquids().sigma(this->pc(), this->T(), X);

    // Average molecular weight of carrier mix - assumes perfect gas
    scalar Wc = this->rhoc_*specie::RR*this->Tc()/this->pc();
    scalar R = specie::RR/Wc;
    scalar Tav = atomization.Taverage(this->T(), this->Tc());

    // calculate average gas density based on average temperature
    scalar rhoAv = this->pc()/(R*Tav);

    scalar soi = td.cloud().injectors().timeStart();
    scalar currentTime = td.cloud().db().time().value();
    const vector& pos = this->position();
    const vector& injectionPos = this->position0();

    // disregard the continous phase when calculating the relative velocity
    // (in line with the deactivated coupled assumption)
    scalar Urel = mag(this->U());

    scalar t0 = max(0.0, currentTime - this->age() - soi);
    scalar t1 = min(t0 + dt, td.cloud().injectors().timeEnd() - soi);

    // this should be the vol flow rate from when the parcel was injected
    scalar volFlowRate = td.cloud().injectors().volumeToInject(t0, t1)/dt;

    scalar chi = 0.0;
    if (atomization.calcChi())
    {
        chi = this->chi(td, X);
    }

    atomization.update
    (
        dt,
        this->d(),
        this->liquidCore(),
        this->tc(),
        rho,
        mu,
        sigma,
        volFlowRate,
        rhoAv,
        Urel,
        pos,
        injectionPos,
        td.cloud().pAmbient(),
        chi,
        td.cloud().rndGen()
    );
}


template<class ParcelType>
template<class TrackData>
void Foam::SprayParcel<ParcelType>::calcBreakup
(
    TrackData& td,
    const scalar dt,
    const label cellI
)
{
    typedef typename TrackData::cloudType::reactingCloudType reactingCloudType;
    const CompositionModel<reactingCloudType>& composition =
        td.cloud().composition();

    typedef typename TrackData::cloudType cloudType;
    typedef typename cloudType::parcelType parcelType;
    typedef typename cloudType::forceType forceType;

    const parcelType& p = static_cast<const parcelType&>(*this);
    const forceType& forces = td.cloud().forces();

    if (td.cloud().breakup().solveOscillationEq())
    {
        solveTABEq(td, dt);
    }

    // cell state info is updated in ReactingParcel calc
    const scalarField& Y(this->Y());
    scalarField X(composition.liquids().X(Y));

    scalar rho = composition.liquids().rho(this->pc(), this->T(), X);
    scalar mu = composition.liquids().mu(this->pc(), this->T(), X);
    scalar sigma = composition.liquids().sigma(this->pc(), this->T(), X);

    // Average molecular weight of carrier mix - assumes perfect gas
    scalar Wc = this->rhoc()*specie::RR*this->Tc()/this->pc();
    scalar R = specie::RR/Wc;
    scalar Tav = td.cloud().atomization().Taverage(this->T(), this->Tc());

    // calculate average gas density based on average temperature
    scalar rhoAv = this->pc()/(R*Tav);
    scalar muAv = this->muc();
    vector Urel = this->U() - this->Uc();
    scalar Urmag = mag(Urel);
    scalar Re = this->Re(this->U(), this->d(), rhoAv, muAv);

    const scalar mass = p.mass();
    const forceSuSp Fcp = forces.calcCoupled(p, dt, mass, Re, muAv);
    const forceSuSp Fncp = forces.calcNonCoupled(p, dt, mass, Re, muAv);
    scalar tMom = mass/(Fcp.Sp() + Fncp.Sp());

    const vector g = td.cloud().g().value();

    scalar massChild = 0.0;
    scalar dChild = 0.0;
    if
    (
        td.cloud().breakup().update
        (
            dt,
            g,
            this->d(),
            this->tc(),
            this->ms(),
            this->nParticle(),
            this->KHindex(),
            this->y(),
            this->yDot(),
            this->d0(),
            rho,
            mu,
            sigma,
            this->U(),
            rhoAv,
            muAv,
            Urel,
            Urmag,
            tMom,
            dChild,
            massChild
        )
    )
    {
        scalar Re = rhoAv*Urmag*dChild/muAv;
        this->mass0() -= massChild;

        // Add child parcel as copy of parent
        SprayParcel<ParcelType>* child = new SprayParcel<ParcelType>(*this);
        child->mass0() = massChild;
        child->d() = dChild;
        child->nParticle() = massChild/rho*this->volume(dChild);

        const forceSuSp Fcp =
            forces.calcCoupled(*child, dt, massChild, Re, muAv);
        const forceSuSp Fncp =
            forces.calcNonCoupled(*child, dt, massChild, Re, muAv);

        child->liquidCore() = 0.0;
        child->KHindex() = 1.0;
        child->y() = td.cloud().breakup().y0();
        child->yDot() = td.cloud().breakup().yDot0();
        child->tc() = -GREAT;
        child->ms() = 0.0;
        child->injector() = this->injector();
        child->tMom() = massChild/(Fcp.Sp() + Fncp.Sp());
        child->user() = 0.0;
        child->setCellValues(td, dt, cellI);

        td.cloud().addParticle(child);
    }
}


template<class ParcelType>
template<class TrackData>
Foam::scalar Foam::SprayParcel<ParcelType>::chi
(
    TrackData& td,
    const scalarField& X
) const
{
    // modifications to take account of the flash boiling on primary break-up

    typedef typename TrackData::cloudType::reactingCloudType reactingCloudType;
    const CompositionModel<reactingCloudType>& composition =
        td.cloud().composition();

    scalar chi = 0.0;
    scalar T0 = this->T();
    scalar p0 = this->pc();
    scalar pAmb = td.cloud().pAmbient();

    scalar pv = composition.liquids().pv(p0, T0, X);

    forAll(composition.liquids(), i)
    {
        if (pv >= 0.999*pAmb)
        {
            // liquid is boiling - calc boiling temperature

            const liquidProperties& liq = composition.liquids().properties()[i];
            scalar TBoil = liq.pvInvert(p0);

            scalar hl = liq.hl(pAmb, TBoil);
            scalar iTp = liq.h(pAmb, T0) - liq.rho(pAmb, T0);
            scalar iTb = liq.h(pAmb, TBoil) - pAmb/liq.rho(pAmb, TBoil);

            chi += X[i]*(iTp - iTb)/hl;
        }
    }

    chi = min(1.0, max(chi, 0.0));

    return chi;
}


template<class ParcelType>
template<class TrackData>
void Foam::SprayParcel<ParcelType>::solveTABEq
(
    TrackData& td,
    const scalar dt
)
{
    typedef typename TrackData::cloudType::reactingCloudType reactingCloudType;
    const CompositionModel<reactingCloudType>& composition =
        td.cloud().composition();

    const scalar& TABCmu = td.cloud().breakup().TABCmu();
    const scalar& TABWeCrit = td.cloud().breakup().TABWeCrit();
    const scalar& TABComega = td.cloud().breakup().TABComega();

    scalar r = 0.5*this->d();
    scalar r2 = r*r;
    scalar r3 = r*r2;

    const scalarField& Y(this->Y());
    scalarField X(composition.liquids().X(Y));

    scalar rho = composition.liquids().rho(this->pc(), this->T(), X);
    scalar mu = composition.liquids().mu(this->pc(), this->T(), X);
    scalar sigma = composition.liquids().sigma(this->pc(), this->T(), X);

    // inverse of characteristic viscous damping time
    scalar rtd = 0.5*TABCmu*mu/(rho*r2);

    // oscillation frequency (squared)
    scalar omega2 = TABComega*sigma/(rho*r3) - rtd*rtd;

    if (omega2 > 0)
    {
        scalar omega = sqrt(omega2);
        scalar rhoc = this->rhoc();
        scalar Wetmp = this->We(this->U(), r, rhoc, sigma)/TABWeCrit;

        scalar y1 = this->y() - Wetmp;
        scalar y2 = this->yDot()/omega;

        // update distortion parameters
        scalar c = cos(omega*dt);
        scalar s = sin(omega*dt);
        scalar e = exp(-rtd*dt);
        y2 = (this->yDot() + y1*rtd)/omega;

        this->y() = Wetmp + e*(y1*c + y2*s);
        if (this->y() < 0)
        {
            this->y() = 0.0;
            this->yDot() = 0.0;
        }
        else
        {
            this->yDot() = (Wetmp - this->y())*rtd + e*omega*(y2*c - y1*s);
        }
    }
    else
    {
        // reset distortion parameters
        this->y() = 0;
        this->yDot() = 0;
    }
}


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

template<class ParcelType>
Foam::SprayParcel<ParcelType>::SprayParcel(const SprayParcel<ParcelType>& p)
:
    ParcelType(p),
    d0_(p.d0_),
    position0_(p.position0_),
    liquidCore_(p.liquidCore_),
    KHindex_(p.KHindex_),
    y_(p.y_),
    yDot_(p.yDot_),
    tc_(p.tc_),
    ms_(p.ms_),
    injector_(p.injector_),
    tMom_(p.tMom_),
    user_(p.user_)
{}


template<class ParcelType>
Foam::SprayParcel<ParcelType>::SprayParcel
(
    const SprayParcel<ParcelType>& p,
    const polyMesh& mesh
)
:
    ParcelType(p, mesh),
    d0_(p.d0_),
    position0_(p.position0_),
    liquidCore_(p.liquidCore_),
    KHindex_(p.KHindex_),
    y_(p.y_),
    yDot_(p.yDot_),
    tc_(p.tc_),
    ms_(p.ms_),
    injector_(p.injector_),
    tMom_(p.tMom_),
    user_(p.user_)
{}


// * * * * * * * * * * * * * * IOStream operators  * * * * * * * * * * * * * //

#include "SprayParcelIO.C"


// ************************************************************************* //
SprayParcel.C (13,830 bytes)   

dkxls

2013-08-07 12:40

reporter   ~0002377

If this is implemented, bug #950 should then fix the last critical issue for spray simulation:
http://www.openfoam.org/mantisbt/view.php?id=950

user2

2013-08-09 13:12

  ~0002384

Thanks for the report - fixed by commit b592ac

Issue History

Date Modified Username Field Change
2013-08-05 09:47 dkxls New Issue
2013-08-06 11:36 user2 Priority high => normal
2013-08-06 11:36 user2 Severity crash => major
2013-08-06 11:37 user2 Note Added: 0002370
2013-08-06 11:37 user2 Status new => resolved
2013-08-06 11:37 user2 Fixed in Version => 2.2.x
2013-08-06 11:37 user2 Resolution open => fixed
2013-08-06 11:37 user2 Assigned To => user2
2013-08-06 17:36 dkxls Note Added: 0002374
2013-08-06 17:36 dkxls Status resolved => feedback
2013-08-06 17:36 dkxls Resolution fixed => reopened
2013-08-06 17:36 dkxls File Added: liquidMixtureProperties.C
2013-08-06 17:36 dkxls File Added: liquidMixtureProperties.H
2013-08-06 17:37 dkxls File Added: SprayParcel.C
2013-08-07 12:35 dkxls Tag Attached: Evaporation
2013-08-07 12:35 dkxls Tag Attached: Lagrangian
2013-08-07 12:35 dkxls Tag Attached: spray
2013-08-07 12:35 dkxls Tag Attached: sprayEngineFoam
2013-08-07 12:35 dkxls Tag Attached: sprayFoam
2013-08-07 12:40 dkxls Note Added: 0002377
2013-08-07 12:40 dkxls Status feedback => assigned
2013-08-09 13:12 user2 Note Added: 0002384
2013-08-09 13:12 user2 Status assigned => resolved
2013-08-09 13:12 user2 Resolution reopened => fixed