| JD_PM |
June 24, 2022 06:35 |
Implementing the temperature equation in interPhaseChangeFoam (OpenFOAM 8 version)
I want to implement a simplified version of the temperature equation (see here: https://ibb.co/RzjbBnj D_T is the effective turbulent thermal diffusivity).
I am using version 8 of OpenFOAM
This is how I proceeded:
I copy-pasted the interPhaseChangeFoam solver (path $FOAM_SOLVERS/multiphase/interPhaseChangeFoam) to my working directory and then changed its name to interPhaseChangeThermalFoam.
I want to keep the default cavitation models (just as a side note: I implemented the ZGB cavitation model, besides the three default ones) available so I can't just erase the phaseChangeTwoPhaseMixtures.
There is a solver called compressibleInterFoam ($FOAM_SOLVERS/multiphase/compressibleInterFoam) which contains a folder called twoPhaseMixtureThermo. The idea was to combine this folder with phaseChangeTwoPhaseMixtures. To do so I copy-pasted the .C and .H files of twoPhaseMixtureThermo to phaseChangeTwoPhaseMixtures.
Once I did that I renamed phaseChangeTwoPhaseMixtures as twoPhaseMixtureThermo.
Next I attach the relevant files regarding twoPhaseMixtureThermo
twoPhaseMixtureThermo/Make/files
Code:
twoPhaseMixtureThermo/phaseChangeTwoPhaseMixture.C
twoPhaseMixtureThermo/phaseChangeTwoPhaseMixtureNew.C
Kunz/Kunz.C
Merkle/Merkle.C
SchnerrSauer/SchnerrSauer.C
ZGB/ZGB.C
twoPhaseMixtureThermo.C
LIB = $(FOAM_USER_LIBBIN)/libtwoPhaseMixtureThermo
twoPhaseMixtureThermo/Make/options
Code:
EXE_INC = \
-I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
-I$(LIB_SRC)/twoPhaseModels/twoPhaseMixture/lnInclude \
-I$(LIB_SRC)/twoPhaseModels/interfaceProperties/lnInclude \
-I$(LIB_SRC)/transportModels/lnInclude \
-I$(LIB_SRC)/twoPhaseModels/incompressibleTwoPhaseMixture/lnInclude \
-I$(LIB_SRC)/twoPhaseModels/immiscibleIncompressibleTwoPhaseMixture/lnInclude \
-I$(LIB_SRC)/finiteVolume/lnInclude
LIB_LIBS = \
-lfluidThermophysicalModels \
-limmiscibleIncompressibleTwoPhaseMixture \
-lspecie \
-ltwoPhaseMixture \
-linterfaceProperties \
-lfiniteVolume
Next I attach the relevant files regarding interPhaseChangeThermalFoam
interPhaseChangeThermalFoam.C
Code:
#include "fvCFD.H"
#include "dynamicFvMesh.H"
#include "CMULES.H"
#include "subCycle.H"
#include "interfaceProperties.H"
#include "phaseChangeTwoPhaseMixture.H"
#include "kinematicMomentumTransportModel.H"
#include "pimpleControl.H"
#include "fvOptions.H"
#include "CorrectPhi.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
int main(int argc, char *argv[])
{
#include "postProcess.H"
#include "setRootCaseLists.H"
#include "createTime.H"
#include "createDynamicFvMesh.H"
#include "createDyMControls.H"
#include "initContinuityErrs.H"
#include "createFields.H"
#include "initCorrectPhi.H"
#include "createUfIfPresent.H"
turbulence->validate();
#include "CourantNo.H"
#include "setInitialDeltaT.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Info<< "\nStarting time loop\n" << endl;
while (pimple.run(runTime))
{
#include "readDyMControls.H"
// Store divU from the previous mesh so that it can be mapped
// and used in correctPhi to ensure the corrected phi has the
// same divergence
volScalarField divU("divU0", fvc::div(fvc::absolute(phi, U)));
#include "CourantNo.H"
#include "alphaCourantNo.H"
#include "setDeltaT.H"
runTime++;
Info<< "Time = " << runTime.timeName() << nl << endl;
// --- Pressure-velocity PIMPLE corrector loop
while (pimple.loop())
{
if (pimple.firstPimpleIter() || moveMeshOuterCorrectors)
{
mesh.update();
if (mesh.changing())
{
gh = (g & mesh.C()) - ghRef;
ghf = (g & mesh.Cf()) - ghRef;
if (correctPhi)
{
// Calculate absolute flux
// from the mapped surface velocity
phi = mesh.Sf() & Uf();
#include "correctPhi.H"
// Make the flux relative to the mesh motion
fvc::makeRelative(phi, U);
}
mixture.correct();
if (checkMeshCourantNo)
{
#include "meshCourantNo.H"
}
}
}
divU = fvc::div(fvc::absolute(phi, U));
surfaceScalarField rhoPhi
(
IOobject
(
"rhoPhi",
runTime.timeName(),
mesh
),
mesh,
dimensionedScalar(dimMass/dimTime, 0)
);
#include "alphaControls.H"
#include "alphaEqnSubCycle.H"
mixture.correct();
#include "UEqn.H"
#include "TEqn.H" /adding temperature equation
// --- Pressure corrector loop
while (pimple.correct())
{
#include "pEqn.H"
}
if (pimple.turbCorr())
{
turbulence->correct();
}
}
runTime.write();
Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
<< " ClockTime = " << runTime.elapsedClockTime() << " s"
<< nl << endl;
}
Info<< "End\n" << endl;
return 0;
}
TEqn.H
Code:
#include "fvCFD.H"
#include "dynamicFvMesh.H"
#include "CMULES.H"
#include "subCycle.H"
#include "interfaceProperties.H"
#include "phaseChangeTwoPhaseMixture.H"
#include "kinematicMomentumTransportModel.H"
#include "pimpleControl.H"
#include "fvOptions.H"
#include "CorrectPhi.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
int main(int argc, char *argv[])
{
#include "postProcess.H"
#include "setRootCaseLists.H"
#include "createTime.H"
#include "createDynamicFvMesh.H"
#include "createDyMControls.H"
#include "initContinuityErrs.H"
#include "createFields.H"
#include "initCorrectPhi.H"
#include "createUfIfPresent.H"
turbulence->validate();
#include "CourantNo.H"
#include "setInitialDeltaT.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Info<< "\nStarting time loop\n" << endl;
while (pimple.run(runTime))
{
#include "readDyMControls.H"
// Store divU from the previous mesh so that it can be mapped
// and used in correctPhi to ensure the corrected phi has the
// same divergence
volScalarField divU("divU0", fvc::div(fvc::absolute(phi, U)));
#include "CourantNo.H"
#include "alphaCourantNo.H"
#include "setDeltaT.H"
runTime++;
Info<< "Time = " << runTime.timeName() << nl << endl;
// --- Pressure-velocity PIMPLE corrector loop
while (pimple.loop())
{
if (pimple.firstPimpleIter() || moveMeshOuterCorrectors)
{
mesh.update();
if (mesh.changing())
{
gh = (g & mesh.C()) - ghRef;
ghf = (g & mesh.Cf()) - ghRef;
if (correctPhi)
{
// Calculate absolute flux
// from the mapped surface velocity
phi = mesh.Sf() & Uf();
#include "correctPhi.H"
// Make the flux relative to the mesh motion
fvc::makeRelative(phi, U);
}
mixture.correct();
if (checkMeshCourantNo)
{
#include "meshCourantNo.H"
}
}
}
divU = fvc::div(fvc::absolute(phi, U));
surfaceScalarField rhoPhi
(
IOobject
(
"rhoPhi",
runTime.timeName(),
mesh
),
mesh,
dimensionedScalar(dimMass/dimTime, 0)
);
#include "alphaControls.H"
#include "alphaEqnSubCycle.H"
mixture.correct();
#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"
<< " ClockTime = " << runTime.elapsedClockTime() << " s"
<< nl << endl;
}
//Adding temperature equation
/*solve
(
fvm::ddt(T) + fvm::div(phi,T) == fvm::laplacian(DT,T)
);*/
Info<< "End\n" << endl;
return 0;
}
createFields.H
Code:
//read the thermal diffusivity constant from transportProperties
Info<< "Reading field p_rgh\n" << endl;
volScalarField p_rgh
(
IOobject
(
"p_rgh",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
Info<< "Reading field U\n" << endl;
volVectorField U
(
IOobject
(
"U",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
//scalar field for the temperature
Info<< "Reading field T\n" << endl;
volVectorField T
(
IOobject
(
"T",
runTime.timeName(),
mesh,
IOobject::MUST_READ, // Must read T from the latest time directory
IOobject::AUTO_WRITE // Data will be written according to controlDict
),
mesh
);
#include "createPhi.H"
Info<< "Creating phaseChangeTwoPhaseMixture\n" << endl;
autoPtr<phaseChangeTwoPhaseMixture> mixturePtr
(
phaseChangeTwoPhaseMixture::New(U, phi)
);
phaseChangeTwoPhaseMixture& mixture = mixturePtr();
volScalarField& alpha1(mixture.alpha1());
volScalarField& alpha2(mixture.alpha2());
const dimensionedScalar& rho1 = mixture.rho1();
const dimensionedScalar& rho2 = mixture.rho2();
// Need to store rho for ddt(rho, U)
volScalarField rho
(
IOobject
(
"rho",
runTime.timeName(),
mesh,
IOobject::READ_IF_PRESENT
),
alpha1*rho1 + alpha2*rho2
);
rho.oldTime();
// Construct incompressible turbulence model
autoPtr<incompressible::momentumTransportModel> turbulence
(
incompressible::momentumTransportModel::New(U, phi, mixture)
);
#include "readGravitationalAcceleration.H"
#include "readhRef.H"
#include "gh.H"
volScalarField p
(
IOobject
(
"p",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
p_rgh + rho*gh
);
label pRefCell = 0;
scalar pRefValue = 0.0;
setRefCell
(
p,
p_rgh,
pimple.dict(),
pRefCell,
pRefValue
);
if (p_rgh.needReference())
{
p += dimensionedScalar
(
"p",
p.dimensions(),
pRefValue - getRefCellValue(p, pRefCell)
);
p_rgh = p - rho*gh;
}
mesh.setFluxRequired(p_rgh.name());
mesh.setFluxRequired(alpha1.name());
#include "createFvOptions.H"
However, when I compile twoPhaseMixtureThermo I get the following error
g++: error: I/software/alternate/opensource/openfoam.user/fc30/8/OpenFOAM-8/src/twoPhaseModels/interfaceProperties/lnInclude: No such file or directory
I do not get why, as lnInclude exists within my twoPhaseMixtureThermo directory.
There is an old thread asking the same question ( https://www.cfd-online.com/Forums/op...hangefoam.html) but it did not solve my issue |
