Diverging result for Temperature field in interFoam
 

Hi all..

I added an energy equation to interFoam, applied to a falling film flow on a cylindrical surface. The result shows a temperature field growing unbounded. I have no idea what could be wrong...maybe the implementation itself or my choice of scheme. Has anyone tried to add an energy equatioin to interFoam? Are there any peculiar problems I should be aware of?

Please help...
Thanks

How do you calculate the face fluxes in the temperature equation? They should be consistent with the fluxes calculated for the momentum equation.

Hi thanks for the response.

The energy equation I added is:

fvm::ddt(rho*cp,T)
+ fvm::div(rhoPhi*cpf,T)
- fvm::laplacian(kappaf,T)

rhoPhi comes from the momentum equation, cpf and kappaf are the cell face values of cp (specific heat) and kappa (thermal conductivity) respectively.

First I calculated face values of gamma (vol fraction of phase 1) as gammaf :
surfaceScalarField gammaf = min(max(fvc::interpolate(gamma), scalar(0)), scalar(1));

then I used this value to compute face values for cp and kappa which I called cpf and kappaf as follows:

surfaceScalarFiled cpf = gammaf*cp1 + (scalar(1) - gammaf)*cp2
surfaceScalarFiled kappaf = gammaf*kappa1 + (scalar(1) - gammaf)*kappa2.

Something tells me that my implementation is suspect though. I am not just quite sure...

Thanks.

It depends on how you calculate your rhoPhi. It is updated after the solution of the gammaEqn. You should update your energy flux (rhof*Cpf*phi) to be consistent with this volume flux, i.e. not doing a simple linear interpolation of Cp.
Another hint: when you are dealing with different materials (having different diffusion coefficients for energy), the you must apply harmonic face interpolation)

Thanks for the reply Edin.

Regarding rhoPhi, I use the values directly from the updates given by the gammaEqn with no modifications whatsoever. Now two questions,

1.) How do I update the values for Cp without doing a simple linear interpolation? Should I employ a hihger order scheme for this?

2.) Please could you explain how harmonic face interpolation is implemented? Is there a dedicated scheme for this in OpenFOAM?

Thanks.

Hi Edin,

I have reviewed the results closely and its apparent that the instability results from the sharp change in thermal properties at the interface. I noticed that the temperature changes sharply at the region where gamma changes from its liquid value to its value for the gas phase. And as the computation marches in time, the instability propagates within the computation domain. So it may come down to a handling the change in thermal properties at the interphase. Does anyone have a clue?

Thanks.

Hi, Ovie, I also doing the same job just like you. The follows are my revision of the interFoam, named interFoamI. I am using the solver to do calculation. I am not sure whether it is OK. Correct me if there are something wrong. Good luck. Chiven

-------------------------------------------------------------------------------------------------------------


copy $WM_PROJECT_DIR/src/transportModels to $WM_PROJECT_USER_DIR/src/transportModels

Revise the file of twoPhaseMixture.C
rho1_(nuModel1_->viscosityProperties().lookup("rho")),
rho2_(nuModel2_->viscosityProperties().lookup("rho")),
at this place, add:
Pr1_(nuModel1_->viscosityProperties().lookup("Pr")),
Pr2_(nuModel2_->viscosityProperties().lookup("Pr")),
beta1_(nuModel1_->viscosityProperties().lookup("beta")),
TRef1_(nuModel1_->viscosityProperties().lookup("TRef")),

nuModel1_->viscosityProperties().lookup("rho") >> rho1_;
nuModel2_->viscosityProperties().lookup("rho") >> rho2_;
at this place, add:
nuModel1_->viscosityProperties().lookup("Pr") >> Pr1_;
nuModel2_->viscosityProperties().lookup("Pr") >> Pr2_;
nuModel1_->viscosityProperties().lookup("beta") >> beta1_;
nuModel1_->viscosityProperties().lookup("TRef") >> TRef1_;

Revise the file of twoPhaseMixture.H
dimensionedScalar rho1_;
dimensionedScalar rho2_;
at this place, add:

dimensionedScalar Pr1_;
dimensionedScalar Pr2_;
dimensionedScalar beta1_;
dimensionedScalar TRef1_;

//- Return const-access to phase1 density
const dimensionedScalar& rho1() const
{
return rho1_;
}
//- Return const-access to phase2 density
const dimensionedScalar& rho2() const
{
return rho2_;
};
at this place, add:
//- Return const-access to phase1 Prandtl number
const dimensionedScalar& Pr1() const
{
return Pr1_;
}
//- Return const-access to phase2 Prandtl number
const dimensionedScalar& Pr2() const
{
return Pr2_;
};
//- thermal expansion coefficient
const dimensionedScalar& beta1() const
{
return beta1_;
}
//- reference temperature
const dimensionedScalar& TRef1() const
{
return TRef1_;
}

Revise the file fo files
LIB = $(FOAM_USER_LIBBIN)/libincompressibleTransportModelsI

Then to compile the incompressibleTransportModelsI
wmake libso

Copy the files in the directory of "OpenFOAM-1.6/applictions/solvers/multiphase/" to the directory of
"$WM_PROJECT_USER_DIR//applictions/solvers/multiphase/"
In the directory of "$WM_PROJECT_USER_DIR//applictions/solvers/multiphase/", perform the following revisions.
cp -r interFoam interFoamI
cd interFoamI
mv interFoam.C interFoamI.C
multiphase/interFoamI> vi interFoamI.C
#include "UEqn.H"
at this place, add:
#include "TEqn.H"

interFoamI/Make> vi files
interFoamI.C
EXE = $(FOAM_USER_APPBIN)/interFoamI

interFoamI/Make> vi options
EXE_INC = \
-I$(LIB_SRC)/transportModels \
-I$(WM_PROJECT_USER_DIR)/src/transportModels/incompressible/lnInclude \
-I$(LIB_SRC)/transportModels/interfaceProperties/lnInclude \
-I$(LIB_SRC)/turbulenceModels/incompressible/turbulenceModel \
-I$(LIB_SRC)/finiteVolume/lnInclude
EXE_LIBS = \
-L$(FOAM_USER_LIBBIN) \
-linterfaceProperties \
-lincompressibleTransportModelsI \
-lincompressibleRASModels \
-lincompressibleLESModels \
-lfiniteVolume

make a new file named TEqn.H:
{
volScalarField kappaEff
(
"kappaEff",
turbulence->nut()/Pr
);
fvScalarMatrix TEqn
(
fvm::ddt(T)
+ fvm::div(phi, T)
- fvm::laplacian(kappaEff, T)
);
TEqn.relax();
TEqn.solve();
rho == alpha1*rho1*(1.0-beta1*(T-TRef1)) + (scalar(1) - alpha1)*rho2;

}


multiphase/interFoamI> vi createFields.H


Info<< "Reading thermophysical properties\n" << endl;
Info<< "Reading field T\n" << endl;
volScalarField T
(
IOobject
(
"T",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);

const dimensionedScalar& rho1 = twoPhaseProperties.rho1();
const dimensionedScalar& rho2 = twoPhaseProperties.rho2();
at the place, add:
const dimensionedScalar& Pr1 = twoPhaseProperties.Pr1();
const dimensionedScalar& Pr2 = twoPhaseProperties.Pr2();
const dimensionedScalar& beta1 = twoPhaseProperties.beta1();
const dimensionedScalar& TRef1 = twoPhaseProperties.TRef1();

// Need to store rho for ddt(rho, U)
volScalarField rho
(
IOobject
(
"rho",
runTime.timeName(),
mesh,
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
),
alpha1*rho1*(1.0-beta1*(T-TRef1)) + (scalar(1) - alpha1)*rho2,
alpha1.boundaryField().types()
);
rho.oldTime();
Herein, add:
// Need to store Pr
volScalarField Pr
(
IOobject
(
"Pr",
runTime.timeName(),
mesh,
IOobject::READ_IF_PRESENT
),
alpha1*Pr1 + (scalar(1) - alpha1)*Pr2
);
Pr.oldTime();

multiphase/interFoamI> vi alphaEqnSubCycle.H
add this line in the end:
rho == alpha1*rho1*(1.0-beta1*(T-TRef1)) + (scalar(1) - alpha1)*rho2;
Pr == alpha1*Pr1 + (scalar(1) - alpha1)*Pr2;

wmake

Ovie, I think you need to work in Cp in rhoPhi, it needs to be consistent with the volume flux calculate in gammaEqn.H. By that I mean, Cp needs to be taken into account at that point in the calculation.

Chiven, why do you incorporate "turbulence->nut()/Pr"? I wouldn't expect to see this in a solver like interFoam. Why don't you define a kappa for both fluids?

Thanks Feijooos.
However, I am not quite sure I understand what you meant by working in Cp in rhoPhi. Is there some other way to compute Cp for the "mixture" other than interpolating using the values of gamma and the constant cp values for the individual phases? Or are you suggesting that cp somehow has to be incoporated in the gammaEqn when solving for gamma? Please be kind to clarify.

Chiven, thanks for the insight. I am trying to review the suggested code. It looks reasonable though.

Thanks all..

Dear feijooos, thank you for the comments. Maybe you are right, and i shall accept your suggestions, and do some calculations.
best regards,
Chiven

Dear Ovie, my revision of the codes have NOT been validated. I am doing some calculations, but it is calculating very slowly. I am not sure whether it is OK. If I have a new progress, I shall post it. Good luck.
Chiven

@ovie
1) do not update Cp alone, you don' need it. simply update the whole energy flux as well as rho*Cp instead, using the same way as mass flux and rho are updated.
2) you may simply use Gauss harmonic for the diffusion term, or implement it on your own like 1.0/fvc::interpolate(1/lambda).

@ Edin,

Thanks for the insight. I would review the suggested approach and see how it impacts the convergence of my result.

@ Edin,
The temperature field converges now. Thanks.

For those interested, the modifications I made are as follows:

1.) twoPhaseMixture.H file:

add the following for specific heat and Prandtl number:
dimensionedScalar cp1_;
dimensionedScalar cp2_;
dimensionedScalar Pr1_;
dimensionedScalar Pr2_;

Add the following functions to return the new data members added to the class.
//- Return const-access to phase1 Prandtl #
const dimensionedScalar& Pr1() const
{
return Pr1_;
}

//- Return const-access to phase2 Prandtl #
const dimensionedScalar& Pr2() const
{
return Pr2_;
};

//- Return const-access to phase1 specific heat
const dimensionedScalar& cp1() const
{
return cp1_;
}

//- Return const-access to phase2 specific heat
const dimensionedScalar& cp2() const
{
return cp2_;
};

Add the function kappaf():
//- Return the face-interpolated thermal conductivity
tmp<surfaceScalarField> kappaf() const;



2.) twoPhaseMixture.C file:

modify the constructor to include the following lines for the new data members:
cp1_(nuModel1_->viscosityProperties().lookup("cp")),
cp2_(nuModel2_->viscosityProperties().lookup("cp")),
Pr1_(nuModel1_->viscosityProperties().lookup("Pr")),
Pr2_(nuModel2_->viscosityProperties().lookup("Pr")),

add definition of kappaf() function as follows:
tmp<surfaceScalarField> twoPhaseMixture::kappaf() const
{
surfaceScalarField alpha1f =
min(max(fvc::interpolate(alpha1_), scalar(0)), scalar(1));

return tmp<surfaceScalarField>
(
new surfaceScalarField
(
"kappaf",
alpha1f*rho1_*cp1_*(1/Pr1_)*fvc::interpolate(nuModel1_->nu())
+ (scalar(1) - alpha1f)*rho2_*cp2_*(1/Pr2_)*fvc::interpolate(nuModel2_->nu())
)
);
}


Modify the read() function as follows:
bool twoPhaseMixture::read()
{
if (transportModel::read())
{
if
(
nuModel1_().read(subDict(phase1Name_))
&& nuModel2_().read(subDict(phase2Name_))
)
{
nuModel1_->viscosityProperties().lookup("rho") >> rho1_;
nuModel2_->viscosityProperties().lookup("rho") >> rho2_;
nuModel1_->viscosityProperties().lookup("cp") >> cp1_;
nuModel2_->viscosityProperties().lookup("cp") >> cp2_;
nuModel1_->viscosityProperties().lookup("Pr") >> Pr1_;
nuModel2_->viscosityProperties().lookup("Pr") >> Pr2_;

return true;
}
else
{
return false;
}
}
else
{
return false;
}
}

Compile using Allwmake. You can execute Allwmake from:
~/OpenFOAM/OpenFOAM-1.5/src/transportModels$

My modified solver I called interTempFoam.
I made the following modifications.

3.) CreateFields.H:
After the expression:
const dimensionedScalar& rho2 = twoPhaseProperties.rho2();
:
Add:
const dimensionedScalar& cp1 = twoPhaseProperties.cp1();
const dimensionedScalar& cp2 = twoPhaseProperties.cp2();

At the end of intialization of volScalarField P:
Add the lines:

// this is used in the unsteady term in the heat equation
// initialization does not matter since this would be recomputed at the end // of gammaEqnSubcycle.H file

Info<< "Reading / calculating rho*cp\n" << endl;
volScalarField rhoCp
(
IOobject
(
"rho*Cp",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE
),
gamma*rho1*cp1 + (scalar(1) - gamma)*rho2*cp2,
gamma.boundaryField().types()
);
rhoCp.oldTime();

// this is used in the convective heat flux term
// Initialisation does not matter because rhoPhiCpf is reset after the
// gamma solution before it is used in the T equation.
Info<< "Reading / calculating rho*phi*cp\n" << endl;
surfaceScalarField rhoPhiCpf
(
IOobject
(
"rho*phi*cpf",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE
),
rhoPhi*cp1
);

4.) gammaEqnSubcycle.H
After the line:
rho == gamma*rho1 + (scalar(1) - gamma)*rho2;

Add:
rhoCp == gamma*rho1*cp1 + (scalar(1) - gamma)*rho2*cp2;

5.) gammaEqn.H
After the expression:
MULES::explicitSolve(gamma, phi, phiGamma, 1, 0);
rhoPhi = phiGamma*(rho1 - rho2) + phi*rho2;

Add:
rhoPhiCpf = phiGamma*(rho1*cp1 - rho2*cp2) + phi*rho2*cp2;

6.) TEqn.H
create the TEqn.H file using the following expressions:

surfaceScalarField kappaf = twoPhaseProperties.kappaf();

fvScalarMatrix TEqn
(
fvm::ddt(rhoCp, T)
+ fvm::div(rhoPhiCpf, T)
- fvm::laplacian(kappaf, T)
);

TEqn.solve();

7.) interTempFoam.C
After the lines:
#include "gammaEqnSubCycle.H"
#include "UEqn.H"

Add:
#include "TEqn.H"

I have assumed that the user knows how to create and compile/link a modified solver from an existing OpenFOAM solver. If not, there is an existing tutorial for this called: "How to add temperature to icoFoam".

One more thing, when running computations on falling film flows, you must take extra care to ensure that the heating surface is fully wetted prior to heating to ensure convergence of the temperature equation. What I normally do is run interFoam to the point where the surface is fully wetted and then map the results to my interTempFoam using a uniform field for temperature at the new start time. This helped the convergence VERY significantly.

Thanks all. If I missed something, please shout!

ovie,
this is correct, congratulations.
you are still doing linear interpolation of kappa, but I guess its influence is very small.
enjoy it.

Hi Edin,

Thanks for the review. But a quick one: which of these two expressions do you suppose would be the more appropriate expression for calculating kappaf in the function kappaf() in twoPhaseMixture.C?

alpha1f*rho1_*cp1_*(1/fvc::interpolate(1/((1/Pr1_)*(nuModel1_->nu()))))
+ (scalar(1) - alpha1f)*rho2_*cp2_*(1/fvc::interpolate(1/((1/Pr2_)*(nuModel2_->nu()))))

or

alpha1f*(1/fvc::interpolate(1/((rho1_*cp1_*(1/Pr1_)*(nuModel1_->nu())))))
+ (scalar(1) - alpha1f)*(1/fvc::interpolate(1/((rho2_*cp2_*(1/Pr2_)*(nuModel2_->nu()))))).


Or is either inappropriate?

Thanks.

hi ovie.

maybe it would be a little complicated to write it all in one expression, because face values must be obtained from cell center values of kappa which should already be some kind of averages between phases. one way would be to update kappa() at cell centers also as a weighted average, just like mu() is updated, and then evaluate the conductivity at faces simply as 1.0/fvc::interpolate(1/kappa()). the other way would be to update kappa() at cell centers as a harmonic mean, but then interpolate it linearly at faces. i think this is something you may play with for a while, you can test it on purely diffusion problems.

however, i am thinking whether the viscosity should be also evaluated using the same procedure. updating kappa() as weighted average is consistent with the diffusion coefficient for the UEqn, i.e. mu(), but the value at cell faces should definitely be evaluated using harmonic interpolation. this was shown e.g. in patankar(1980). it is clear that density is a weighted average, but i am not happy with the same definition for viscosity (and generally diffusion coefficients, like conductivity). nevertheless, as i see that all people are doing it exactly in the same way like it is in OpenFOAM, so i think the eventual errors are really small.

regards.

Ovie,

I was reading over this thread again and I saw that I missed your question on how to incorporate cp into the face fluxes. I guess you already found that out :)

I have a question for you though, why do you use face values for kappa? Why not just use cellcentered values?

Thanks,

Eelco

Hi Eelco,

Those were my thoughts precisely: to use cell-centered values for kappa in the heat equation. I dont have any logical explaination for my choice of face centered values except that in the momentum equation, I noticed that face centered values were used for mu. And I simply followed suit. I dont know how this affects the accurary and maybe convergence of the results though.

hmm interesting. I'm wondering if you can just use:

kappa == gamma*kappa1 + (scalar(1) - gamma)*kappa2;

Don't see how this would affect the solution, but you need to specify the different conductivities for both fluids.