openfoam/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright (C) 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 "twoPhaseSystem.H"
#include "surfaceInterpolate.H"
#include "fixedValueFvsPatchFields.H"
#include "fvcCurl.H"

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

Foam::twoPhaseSystem::twoPhaseSystem
(
    const fvMesh& mesh
)
:
    IOdictionary
    (
        IOobject
        (
            "phaseProperties",
            mesh.time().constant(),
            mesh,
            IOobject::MUST_READ_IF_MODIFIED,
            IOobject::NO_WRITE
        )
    ),

    mesh_(mesh),

    phase1_
    (
        *this,
        *this,
        wordList(lookup("phases"))[0]
    ),

    phase2_
    (
        *this,
        *this,
        wordList(lookup("phases"))[1]
    ),

    Cvm_
    (
        "Cvm",
        dimless,
        lookup("Cvm")
    ),

    Cl_
    (
        "Cl",
        dimless,
        lookup("Cl")
    ),

    drag1_
    (
        dragModel::New
        (
            subDict("drag"),
            phase1_,
            phase1_,
            phase2_
        )
    ),

    drag2_
    (
        dragModel::New
        (
            subDict("drag"),
            phase2_,
            phase2_,
            phase1_
        )
    ),

    heatTransfer1_
    (
        heatTransferModel::New
        (
            subDict("heatTransfer"),
            phase1_,
            phase1_,
            phase2_
        )
    ),

    heatTransfer2_
    (
        heatTransferModel::New
        (
            subDict("heatTransfer"),
            phase2_,
            phase2_,
            phase1_
        )
    ),

    dispersedPhase_(lookup("dispersedPhase")),

    residualPhaseFraction_
    (
        readScalar(lookup("residualPhaseFraction"))
    ),

    residualSlip_
    (
        "residualSlip",
        dimVelocity,
        lookup("residualSlip")
    )
{
    if
    (
        !(
            dispersedPhase_ == phase1_.name()
         || dispersedPhase_ == phase2_.name()
         || dispersedPhase_ == "both"
        )
    )
    {
        FatalErrorIn("twoPhaseSystem::twoPhaseSystem(const fvMesh& mesh)")
            << "invalid dispersedPhase " << dispersedPhase_
            << exit(FatalError);
    }

    Info << "dispersedPhase is " << dispersedPhase_ << endl;

    // Ensure the phase-fractions sum to 1
    phase2_.volScalarField::operator=(scalar(1) - phase1_);
}


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

Foam::tmp<Foam::volScalarField> Foam::twoPhaseSystem::rho() const
{
    return phase1_*phase1_.thermo().rho() + phase2_*phase2_.thermo().rho();
}


Foam::tmp<Foam::volVectorField> Foam::twoPhaseSystem::U() const
{
    return phase1_*phase1_.U() + phase2_*phase2_.U();
}


Foam::tmp<Foam::surfaceScalarField> Foam::twoPhaseSystem::phi() const
{
    return
        fvc::interpolate(phase1_)*phase1_.phi()
      + fvc::interpolate(phase2_)*phase2_.phi();
}


Foam::tmp<Foam::volScalarField> Foam::twoPhaseSystem::dragCoeff() const
{
    tmp<volScalarField> tdragCoeff
    (
        new volScalarField
        (
            IOobject
            (
                "dragCoeff",
                mesh_.time().timeName(),
                mesh_
            ),
            mesh_,
            dimensionedScalar("dragCoeff", dimensionSet(1, -3, -1, 0, 0), 0)
        )
    );
    volScalarField& dragCoeff = tdragCoeff();

    volVectorField Ur(phase1_.U() - phase2_.U());
    volScalarField magUr(mag(Ur) + residualSlip_);

    if (dispersedPhase_ == phase1_.name())
    {
        dragCoeff = drag1().K(magUr);
    }
    else if (dispersedPhase_ == phase2_.name())
    {
        dragCoeff = drag2().K(magUr);
    }
    else if (dispersedPhase_ == "both")
    {
        dragCoeff =
        (
            phase2_*drag1().K(magUr)
          + phase1_*drag2().K(magUr)
        );
    }
    else
    {
        FatalErrorIn("twoPhaseSystem::dragCoeff()")
            << "dispersedPhase: " << dispersedPhase_ << " is incorrect"
            << exit(FatalError);
    }

    volScalarField alphaCoeff(max(phase1_*phase2_, residualPhaseFraction_));
    dragCoeff *= alphaCoeff;

    // Remove drag at fixed-flux boundaries
    forAll(phase1_.phi().boundaryField(), patchi)
    {
        if
        (
            isA<fixedValueFvsPatchScalarField>
            (
                phase1_.phi().boundaryField()[patchi]
            )
        )
        {
            dragCoeff.boundaryField()[patchi] = 0.0;
        }
    }

    return tdragCoeff;
}


Foam::tmp<Foam::volVectorField> Foam::twoPhaseSystem::liftForce
(
    const volVectorField& U
) const
{
    tmp<volVectorField> tliftForce
    (
        new volVectorField
        (
            IOobject
            (
                "liftForce",
                mesh_.time().timeName(),
                mesh_
            ),
            mesh_,
            dimensionedVector
            (
                "liftForce",
                dimensionSet(1, -2, -2, 0, 0),
                vector::zero
            )
        )
    );
    volVectorField& liftForce = tliftForce();

    volVectorField Ur(phase1_.U() - phase2_.U());

    liftForce =
        Cl_*(phase1_*phase1_.rho() + phase2_*phase2_.rho())
       *(Ur ^ fvc::curl(U));

    // Remove lift at fixed-flux boundaries
    forAll(phase1_.phi().boundaryField(), patchi)
    {
        if
        (
            isA<fixedValueFvsPatchScalarField>
            (
                phase1_.phi().boundaryField()[patchi]
            )
        )
        {
            liftForce.boundaryField()[patchi] = vector::zero;
        }
    }

    return tliftForce;
}


Foam::tmp<Foam::volScalarField> Foam::twoPhaseSystem::heatTransferCoeff() const
{
    tmp<volScalarField> theatTransferCoeff
    (
        new volScalarField
        (
            IOobject
            (
                "heatTransferCoeff",
                mesh_.time().timeName(),
                mesh_
            ),
            mesh_,
            dimensionedScalar
            (
                "heatTransferCoeff",
                dimensionSet(1, -1, -3, -1, 0),
                0
            )
        )
    );
    volScalarField& heatTransferCoeff = theatTransferCoeff();

    volVectorField Ur(phase1_.U() - phase2_.U());
    volScalarField magUr(mag(Ur) + residualSlip_);

    if (dispersedPhase_ == phase1_.name())
    {
        heatTransferCoeff = heatTransfer1().K(magUr);
    }
    else if (dispersedPhase_ == phase2_.name())
    {
        heatTransferCoeff = heatTransfer2().K(magUr);
    }
    else if (dispersedPhase_ == "both")
    {
        heatTransferCoeff =
        (
            phase2_*heatTransfer1().K(magUr)
          + phase1_*heatTransfer2().K(magUr)
        );
    }
    else
    {
        FatalErrorIn("twoPhaseSystem::heatTransferCoeff()")
            << "dispersedPhase: " << dispersedPhase_ << " is incorrect"
            << exit(FatalError);
    }

    volScalarField alphaCoeff(max(phase1_*phase2_, residualPhaseFraction_));
    heatTransferCoeff *= alphaCoeff;

    // Remove heatTransfer at fixed-flux boundaries
    forAll(phase1_.phi().boundaryField(), patchi)
    {
        if
        (
            isA<fixedValueFvsPatchScalarField>
            (
                phase1_.phi().boundaryField()[patchi]
            )
        )
        {
            heatTransferCoeff.boundaryField()[patchi] = 0.0;
        }
    }

    return theatTransferCoeff;
}


bool Foam::twoPhaseSystem::read()
{
    if (regIOobject::read())
    {
        bool readOK = true;

        readOK &= phase1_.read();
        readOK &= phase2_.read();

        lookup("Cvm") >> Cvm_;
        lookup("Cl") >> Cl_;

        return readOK;
    }
    else
    {
        return false;
    }
}


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