compressibleInterFilmFoam: Experimental VoF solver supporting VoF<->film transfer

This commit is contained in:
Henry Weller
2017-06-27 15:55:43 +01:00
parent 9a4851faf9
commit d74f354f0c
30 changed files with 1744 additions and 146 deletions

View File

@ -0,0 +1,4 @@
VoFPatchTransfer/VoFPatchTransfer.C
compressibleInterFilmFoam.C
EXE = $(FOAM_APPBIN)/compressibleInterFilmFoam

View File

@ -0,0 +1,43 @@
EXE_INC = -ggdb3 \
-I. \
-I.. \
-I../../VoF \
-I../twoPhaseMixtureThermo \
-I$(LIB_SRC)/transportModels/compressible/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
-I$(LIB_SRC)/transportModels/twoPhaseMixture/lnInclude \
-I$(LIB_SRC)/transportModels/interfaceProperties/lnInclude \
-I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
-I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/thermophysicalProperties/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/thermophysicalFunctions/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/SLGThermo/lnInclude \
-I$(LIB_SRC)/regionModels/regionModel/lnInclude \
-I$(LIB_SRC)/regionModels/surfaceFilmModels/lnInclude \
-I$(LIB_SRC)/dynamicMesh/lnInclude \
-I$(LIB_SRC)/meshTools/lnInclude \
-I$(LIB_SRC)/dynamicFvMesh/lnInclude \
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
EXE_LIBS = \
-ltwoPhaseMixtureThermo \
-ltwoPhaseSurfaceTension \
-lcompressibleTransportModels \
-lfluidThermophysicalModels \
-lspecie \
-ltwoPhaseMixture \
-ltwoPhaseProperties \
-linterfaceProperties \
-lturbulenceModels \
-lcompressibleTurbulenceModels \
-lSLGThermo \
-lsurfaceFilmModels \
-lsurfaceFilmDerivedFvPatchFields \
-ldynamicMesh \
-lmeshTools \
-ldynamicFvMesh \
-lfiniteVolume \
-lfvOptions

View File

@ -0,0 +1,36 @@
{
fvScalarMatrix TEqn
(
fvm::ddt(rho, T) + fvm::div(rhoPhi, T)
- fvm::Sp(contErr, T)
- fvm::laplacian(mixture.alphaEff(turbulence->mut()), T)
+ (
(
fvc::div(fvc::absolute(phi, U), p)
+ fvc::ddt(rho, K) + fvc::div(rhoPhi, K)
)()() - contErr*K
)
*(
alpha1/mixture.thermo1().Cv()
+ alpha2/mixture.thermo2().Cv()
)()()
==
fvOptions(rho, T)
+ pos(Srho)
*(
surfaceFilm.Sh()()/mixture.thermo1().Cp()()
+ surfaceFilm.Tref*Srho
)
);
TEqn.relax();
fvOptions.constrain(TEqn);
TEqn.solve();
fvOptions.correct(T);
mixture.correctThermo();
mixture.correct();
}

View File

@ -0,0 +1,325 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2017 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 "VoFPatchTransfer.H"
#include "twoPhaseMixtureThermo.H"
#include "addToRunTimeSelectionTable.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
namespace regionModels
{
namespace surfaceFilmModels
{
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
defineTypeNameAndDebug(VoFPatchTransfer, 0);
addToRunTimeSelectionTable(transferModel, VoFPatchTransfer, dictionary);
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
VoFPatchTransfer::VoFPatchTransfer
(
surfaceFilmModel& film,
const dictionary& dict
)
:
transferModel(type(), film, dict),
deltaFactorToVoF_
(
coeffDict_.lookupOrDefault<scalar>("deltaFactorToVoF", 1.0)
),
deltaFactorToFilm_
(
coeffDict_.lookupOrDefault<scalar>("deltaFactorToFilm", 0.5)
),
alphaToVoF_
(
coeffDict_.lookupOrDefault<scalar>("alphaToVoF", 0.5)
),
alphaToFilm_
(
coeffDict_.lookupOrDefault<scalar>("alphaToFilm", 0.1)
),
transferRateCoeff_
(
coeffDict_.lookupOrDefault<scalar>("transferRateCoeff", 0.1)
)
{
const polyBoundaryMesh& pbm = film.regionMesh().boundaryMesh();
patchIDs_.setSize(pbm.size());
if (coeffDict_.found("patches"))
{
const wordReList patchNames(coeffDict_.lookup("patches"));
const labelHashSet patchSet = pbm.patchSet(patchNames);
Info<< " applying to patches:" << nl;
label pidi = 0;
forAllConstIter(labelHashSet, patchSet, iter)
{
const label patchi = iter.key();
patchIDs_[pidi++] = patchi;
Info<< " " << pbm[patchi].name() << endl;
}
patchIDs_.setSize(pidi);
patchTransferredMasses_.setSize(pidi, 0);
}
else
{
Info<< " applying to all patches" << endl;
forAll(patchIDs_, patchi)
{
patchIDs_[patchi] = patchi;
}
patchTransferredMasses_.setSize(pbm.size(), 0);
}
if (!patchIDs_.size())
{
FatalErrorInFunction
<< "No patches selected"
<< exit(FatalError);
}
}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
VoFPatchTransfer::~VoFPatchTransfer()
{}
// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
void VoFPatchTransfer::correct
(
scalarField& availableMass,
scalarField& massToTransfer
)
{
NotImplemented;
}
void VoFPatchTransfer::correct
(
scalarField& availableMass,
scalarField& massToTransfer,
scalarField& energyToTransfer
)
{
// Do not correct if no patches selected
if (!patchIDs_.size()) return;
const scalarField& delta = film().delta();
const scalarField& rho = film().rho();
const scalarField& magSf = film().magSf();
const polyBoundaryMesh& pbm = film().regionMesh().boundaryMesh();
const twoPhaseMixtureThermo& thermo
(
film().primaryMesh().lookupObject<twoPhaseMixtureThermo>
(
twoPhaseMixtureThermo::dictName
)
);
const volScalarField& heVoF = thermo.thermo1().he();
const volScalarField& TVoF = thermo.thermo1().T();
volScalarField CpVoF = thermo.thermo1().Cp();
const volScalarField& rhoVoF = thermo.thermo1().rho()();
const volScalarField& alphaVoF = thermo.alpha1();
forAll(patchIDs_, pidi)
{
const label patchi = patchIDs_[pidi];
label primaryPatchi = -1;
forAll(film().intCoupledPatchIDs(), i)
{
const label filmPatchi = film().intCoupledPatchIDs()[i];
if (filmPatchi == patchi)
{
primaryPatchi = film().primaryPatchIDs()[i];
}
}
if (primaryPatchi != -1)
{
scalarField deltaCoeffs
(
film().primaryMesh().boundary()[primaryPatchi].deltaCoeffs()
);
film().toRegion(patchi, deltaCoeffs);
scalarField hp(heVoF.boundaryField()[primaryPatchi]);
film().toRegion(patchi, hp);
scalarField Tp(TVoF.boundaryField()[primaryPatchi]);
film().toRegion(patchi, Tp);
scalarField Cpp(CpVoF.boundaryField()[primaryPatchi]);
film().toRegion(patchi, Cpp);
scalarField rhop(rhoVoF.boundaryField()[primaryPatchi]);
film().toRegion(patchi, rhop);
scalarField alphap(alphaVoF.boundaryField()[primaryPatchi]);
film().toRegion(patchi, alphap);
scalarField Vp
(
film().primaryMesh().boundary()[primaryPatchi]
.patchInternalField(film().primaryMesh().V())
);
film().toRegion(patchi, Vp);
const polyPatch& pp = pbm[patchi];
const labelList& faceCells = pp.faceCells();
// Accumulate the total mass removed from patch
scalar dMassPatch = 0;
forAll(faceCells, facei)
{
const label celli = faceCells[facei];
scalar dMass = 0;
if
(
delta[celli] > 2*deltaFactorToVoF_/deltaCoeffs[facei]
|| alphap[facei] > alphaToVoF_
)
{
dMass =
transferRateCoeff_*delta[celli]*rho[celli]*magSf[celli];
massToTransfer[celli] += dMass;
energyToTransfer[celli] += dMass*film().hs()[celli];
}
if
(
alphap[facei] > 0
&& delta[celli] > 0
&& delta[celli] < 2*deltaFactorToFilm_/deltaCoeffs[facei]
&& alphap[facei] < alphaToFilm_
)
{
dMass =
-transferRateCoeff_*alphap[facei]*rhop[facei]*Vp[facei];
massToTransfer[celli] += dMass;
energyToTransfer[celli] += dMass*hp[facei];
}
availableMass[celli] -= dMass;
dMassPatch += dMass;
}
patchTransferredMasses_[pidi] += dMassPatch;
addToTransferredMass(dMassPatch);
}
}
transferModel::correct();
if (writeTime())
{
scalarField patchTransferredMasses0
(
getModelProperty<scalarField>
(
"patchTransferredMasses",
scalarField(patchTransferredMasses_.size(), 0)
)
);
scalarField patchTransferredMassTotals(patchTransferredMasses_);
Pstream::listCombineGather
(
patchTransferredMassTotals,
plusEqOp<scalar>()
);
patchTransferredMasses0 += patchTransferredMassTotals;
setModelProperty<scalarField>
(
"patchTransferredMasses",
patchTransferredMasses0
);
patchTransferredMasses_ = 0;
}
}
void VoFPatchTransfer::patchTransferredMassTotals
(
scalarField& patchMasses
) const
{
// Do not correct if no patches selected
if (!patchIDs_.size()) return;
scalarField patchTransferredMasses
(
getModelProperty<scalarField>
(
"patchTransferredMasses",
scalarField(patchTransferredMasses_.size(), 0)
)
);
scalarField patchTransferredMassTotals(patchTransferredMasses_);
Pstream::listCombineGather(patchTransferredMassTotals, plusEqOp<scalar>());
forAll(patchIDs_, pidi)
{
const label patchi = patchIDs_[pidi];
patchMasses[patchi] +=
patchTransferredMasses[pidi] + patchTransferredMassTotals[pidi];
}
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace surfaceFilmModels
} // End namespace regionModels
} // End namespace Foam
// ************************************************************************* //

View File

@ -0,0 +1,144 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2017 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::regionModels::surfaceFilmModels::VoFPatchTransfer
Description
Transfer mass between the film and the VoF in the continuous phase.
SourceFiles
VoFPatchTransfer.C
\*---------------------------------------------------------------------------*/
#ifndef VoFPatchTransfer_H
#define VoFPatchTransfer_H
#include "transferModel.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
namespace regionModels
{
namespace surfaceFilmModels
{
/*---------------------------------------------------------------------------*\
Class VoFPatchTransfer Declaration
\*---------------------------------------------------------------------------*/
class VoFPatchTransfer
:
public transferModel
{
// Private member functions
//- Disallow default bitwise copy construct
VoFPatchTransfer(const VoFPatchTransfer&);
//- Disallow default bitwise assignment
void operator=(const VoFPatchTransfer&);
protected:
//- Factor of the cell height above which the film is transferred
// to the VoF
scalar deltaFactorToVoF_;
//- Factor of the cell height below which the VoF may be transferred
// to the film
scalar deltaFactorToFilm_;
//- VoF limit above which all of the film is transferred to the VoF
scalar alphaToVoF_;
//- VoF limit below which the VoF may be transferred to the film
scalar alphaToFilm_;
//- Transfer rate coefficient
scalar transferRateCoeff_;
//- List of patch IDs at which the film is removed
labelList patchIDs_;
//- Transferred mass for each patch at which the film is removed
scalarField patchTransferredMasses_;
public:
//- Runtime type information
TypeName("VoFPatchTransfer");
// Constructors
//- Construct from surface film model
VoFPatchTransfer(surfaceFilmModel& film, const dictionary& dict);
//- Destructor
virtual ~VoFPatchTransfer();
// Member Functions
//- Correct
virtual void correct
(
scalarField& availableMass,
scalarField& massToTransfer
);
//- Correct kinematic and thermodynamic transfers
virtual void correct
(
scalarField& availableMass,
scalarField& massToTransfer,
scalarField& energyToTransfer
);
//- Accumulate the total mass injected for the patches into the
// scalarField provided
virtual void patchTransferredMassTotals
(
scalarField& patchMasses
) const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace surfaceFilmModels
} // End namespace regionModels
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -0,0 +1,148 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2017 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/>.
Application
compressibleInterFoam
Description
Solver for 2 compressible, non-isothermal immiscible fluids using a VOF
(volume of fluid) phase-fraction based interface capturing approach.
The momentum and other fluid properties are of the "mixture" and a single
momentum equation is solved.
Turbulence modelling is generic, i.e. laminar, RAS or LES may be selected.
\*---------------------------------------------------------------------------*/
#include "fvCFD.H"
#include "CMULES.H"
#include "EulerDdtScheme.H"
#include "localEulerDdtScheme.H"
#include "CrankNicolsonDdtScheme.H"
#include "subCycle.H"
#include "rhoThermo.H"
#include "twoPhaseMixture.H"
#include "twoPhaseMixtureThermo.H"
#include "turbulentFluidThermoModel.H"
#include "SLGThermo.H"
#include "surfaceFilmModel.H"
#include "pimpleControl.H"
#include "fvOptions.H"
#include "fvcSmooth.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
int main(int argc, char *argv[])
{
#include "postProcess.H"
#include "setRootCase.H"
#include "createTime.H"
#include "createMesh.H"
#include "createControl.H"
#include "createTimeControls.H"
#include "createFields.H"
#include "createAlphaFluxes.H"
#include "createSurfaceFilmModel.H"
#include "createFvOptions.H"
volScalarField& p = mixture.p();
volScalarField& T = mixture.T();
const volScalarField& psi1 = mixture.thermo1().psi();
const volScalarField& psi2 = mixture.thermo2().psi();
filmModelType& surfaceFilm = tsurfaceFilm();
turbulence->validate();
if (!LTS)
{
#include "readTimeControls.H"
#include "CourantNo.H"
#include "setInitialDeltaT.H"
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Info<< "\nStarting time loop\n" << endl;
while (runTime.run())
{
#include "readTimeControls.H"
if (LTS)
{
#include "setRDeltaT.H"
}
else
{
#include "CourantNo.H"
#include "alphaCourantNo.H"
#include "setDeltaT.H"
}
runTime++;
Info<< "Time = " << runTime.timeName() << nl << endl;
surfaceFilm.evolve();
// --- Pressure-velocity PIMPLE corrector loop
while (pimple.loop())
{
#include "alphaControls.H"
#include "compressibleAlphaEqnSubCycle.H"
volScalarField::Internal Srho(surfaceFilm.Srho());
contErr -= posPart(Srho);
#include "UEqn.H"
#include "TEqn.H"
// --- Pressure corrector loop
while (pimple.correct())
{
#include "pEqn.H"
}
if (pimple.turbCorr())
{
turbulence->correct();
}
}
runTime.write();
Info<< "ExecutionTime = "
<< runTime.elapsedCpuTime()
<< " s\n\n" << endl;
}
Info<< "End\n" << endl;
return 0;
}
// ************************************************************************* //

View File

@ -0,0 +1,6 @@
Info<< "\nConstructing surface film model" << endl;
typedef regionModels::surfaceFilmModels::surfaceFilmModel filmModelType;
SLGThermo slgThermo(mesh, mixture.thermo1());
autoPtr<filmModelType> tsurfaceFilm(filmModelType::New(mesh, g));

View File

@ -0,0 +1,129 @@
{
volScalarField rAU("rAU", 1.0/UEqn.A());
surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU));
volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p_rgh));
surfaceScalarField phiHbyA
(
"phiHbyA",
fvc::flux(HbyA)
+ fvc::interpolate(rho*rAU)*fvc::ddtCorr(U, phi)
);
MRF.makeRelative(phiHbyA);
surfaceScalarField phig
(
(
mixture.surfaceTensionForce()
- ghf*fvc::snGrad(rho)
)*rAUf*mesh.magSf()
);
phiHbyA += phig;
// Update the pressure BCs to ensure flux consistency
constrainPressure(p_rgh, U, phiHbyA, rAUf, MRF);
tmp<fvScalarMatrix> p_rghEqnComp1;
tmp<fvScalarMatrix> p_rghEqnComp2;
if (pimple.transonic())
{
surfaceScalarField phid1("phid1", fvc::interpolate(psi1)*phi);
surfaceScalarField phid2("phid2", fvc::interpolate(psi2)*phi);
p_rghEqnComp1 =
fvc::ddt(rho1) + fvc::div(phi, rho1) - fvc::Sp(fvc::div(phi), rho1)
+ correction
(
psi1*fvm::ddt(p_rgh)
+ fvm::div(phid1, p_rgh) - fvm::Sp(fvc::div(phid1), p_rgh)
);
deleteDemandDrivenData(p_rghEqnComp1.ref().faceFluxCorrectionPtr());
p_rghEqnComp1.ref().relax();
p_rghEqnComp2 =
fvc::ddt(rho2) + fvc::div(phi, rho2) - fvc::Sp(fvc::div(phi), rho2)
+ correction
(
psi2*fvm::ddt(p_rgh)
+ fvm::div(phid2, p_rgh) - fvm::Sp(fvc::div(phid2), p_rgh)
);
deleteDemandDrivenData(p_rghEqnComp2.ref().faceFluxCorrectionPtr());
p_rghEqnComp2.ref().relax();
}
else
{
#include "rhofs.H"
p_rghEqnComp1 =
pos(alpha1)
*(
(
fvc::ddt(alpha1, rho1) + fvc::div(alphaPhi1*rho1f)
- (fvOptions(alpha1, mixture.thermo1().rho())&rho1)
)/rho1
- fvc::ddt(alpha1) - fvc::div(alphaPhi1)
+ (alpha1*psi1/rho1)*correction(fvm::ddt(p_rgh))
) - surfaceFilm.Srho()/rho1;
p_rghEqnComp2 =
pos(alpha2)
*(
(
fvc::ddt(alpha2, rho2) + fvc::div(alphaPhi2*rho2f)
- (fvOptions(alpha2, mixture.thermo2().rho())&rho2)
)/rho2
- fvc::ddt(alpha2) - fvc::div(alphaPhi2)
+ (alpha2*psi2/rho2)*correction(fvm::ddt(p_rgh))
);
}
// Cache p_rgh prior to solve for density update
volScalarField p_rgh_0(p_rgh);
while (pimple.correctNonOrthogonal())
{
fvScalarMatrix p_rghEqnIncomp
(
fvc::div(phiHbyA)
- fvm::laplacian(rAUf, p_rgh)
);
solve
(
p_rghEqnComp1() + p_rghEqnComp2() + p_rghEqnIncomp,
mesh.solver(p_rgh.select(pimple.finalInnerIter()))
);
if (pimple.finalNonOrthogonalIter())
{
p = max(p_rgh + (alpha1*rho1 + alpha2*rho2)*gh, pMin);
p_rgh = p - (alpha1*rho1 + alpha2*rho2)*gh;
dgdt =
(
alpha1*(p_rghEqnComp2 & p_rgh)
- alpha2*(p_rghEqnComp1 & p_rgh)
);
phi = phiHbyA + p_rghEqnIncomp.flux();
U = HbyA
+ rAU*fvc::reconstruct((phig + p_rghEqnIncomp.flux())/rAUf);
U.correctBoundaryConditions();
fvOptions.correct(U);
}
}
// Update densities from change in p_rgh
mixture.thermo1().correctRho(psi1*(p_rgh - p_rgh_0));
mixture.thermo2().correctRho(psi2*(p_rgh - p_rgh_0));
rho = alpha1*rho1 + alpha2*rho2;
// Correct p_rgh for consistency with p and the updated densities
p_rgh = p - rho*gh;
p_rgh.correctBoundaryConditions();
K = 0.5*magSqr(U);
}

View File

@ -87,10 +87,7 @@ surfaceScalarField rhoPhi
fvc::interpolate(rho)*phi
);
volScalarField dgdt
(
pos0(alpha2)*fvc::div(phi)/max(alpha2, scalar(0.0001))
);
volScalarField dgdt(alpha1*fvc::div(phi));
// Construct compressible turbulence model
autoPtr<compressible::turbulenceModel> turbulence

View File

@ -63,8 +63,8 @@ Foam::twoPhaseMixtureThermo::twoPhaseMixtureThermo
thermo1_ = rhoThermo::New(U.mesh(), phase1Name());
thermo2_ = rhoThermo::New(U.mesh(), phase2Name());
thermo1_->validate(phase1Name(), "e");
thermo2_->validate(phase2Name(), "e");
// thermo1_->validate(phase1Name(), "e");
// thermo2_->validate(phase2Name(), "e");
correct();
}