CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > Software User Forums > OpenFOAM > OpenFOAM Running, Solving & CFD

Diverging result for Temperature field in interFoam

Register Blogs Members List Search Today's Posts Mark Forums Read

Like Tree45Likes

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   October 12, 2009, 10:14
Default Diverging result for Temperature field in interFoam
  #1
Member
 
Ovie Doro
Join Date: Jul 2009
Posts: 99
Rep Power: 17
ovie is on a distinguished road
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
ovie is offline   Reply With Quote

Old   October 13, 2009, 06:19
Default
  #2
Member
 
Edin Berberovic
Join Date: Mar 2009
Posts: 31
Rep Power: 17
eberberovic is on a distinguished road
How do you calculate the face fluxes in the temperature equation? They should be consistent with the fluxes calculated for the momentum equation.
eberberovic is offline   Reply With Quote

Old   October 13, 2009, 11:53
Default
  #3
Member
 
Ovie Doro
Join Date: Jul 2009
Posts: 99
Rep Power: 17
ovie is on a distinguished road
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.
Kummi likes this.
ovie is offline   Reply With Quote

Old   October 14, 2009, 03:54
Default
  #4
Member
 
Edin Berberovic
Join Date: Mar 2009
Posts: 31
Rep Power: 17
eberberovic is on a distinguished road
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)
eberberovic is offline   Reply With Quote

Old   October 14, 2009, 10:58
Default
  #5
Member
 
Ovie Doro
Join Date: Jul 2009
Posts: 99
Rep Power: 17
ovie is on a distinguished road
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.
ovie is offline   Reply With Quote

Old   October 14, 2009, 20:03
Default
  #6
Member
 
Ovie Doro
Join Date: Jul 2009
Posts: 99
Rep Power: 17
ovie is on a distinguished road
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.
ovie is offline   Reply With Quote

Old   October 14, 2009, 21:13
Default
  #7
Senior Member
 
J. Cai
Join Date: Apr 2009
Posts: 180
Rep Power: 17
chiven is on a distinguished road
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
amolrajan, shahoco and Kummi like this.
chiven is offline   Reply With Quote

Old   October 14, 2009, 21:27
Default
  #8
Member
 
Eelco Gehring
Join Date: Mar 2009
Posts: 70
Rep Power: 17
feijooos is on a distinguished road
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?
feijooos is offline   Reply With Quote

Old   October 14, 2009, 21:51
Default
  #9
Member
 
Ovie Doro
Join Date: Jul 2009
Posts: 99
Rep Power: 17
ovie is on a distinguished road
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..
ovie is offline   Reply With Quote

Old   October 14, 2009, 21:52
Default
  #10
Senior Member
 
J. Cai
Join Date: Apr 2009
Posts: 180
Rep Power: 17
chiven is on a distinguished road
Dear feijooos, thank you for the comments. Maybe you are right, and i shall accept your suggestions, and do some calculations.
best regards,
Chiven
chiven is offline   Reply With Quote

Old   October 14, 2009, 21:57
Default
  #11
Senior Member
 
J. Cai
Join Date: Apr 2009
Posts: 180
Rep Power: 17
chiven is on a distinguished road
Quote:
Originally Posted by ovie View Post

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

Thanks all..


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
chiven is offline   Reply With Quote

Old   October 15, 2009, 03:22
Default
  #12
Member
 
Edin Berberovic
Join Date: Mar 2009
Posts: 31
Rep Power: 17
eberberovic is on a distinguished road
@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).
eberberovic is offline   Reply With Quote

Old   October 15, 2009, 09:59
Default
  #13
Member
 
Ovie Doro
Join Date: Jul 2009
Posts: 99
Rep Power: 17
ovie is on a distinguished road
@ Edin,

Thanks for the insight. I would review the suggested approach and see how it impacts the convergence of my result.
ovie is offline   Reply With Quote

Old   October 15, 2009, 20:30
Default
  #14
Member
 
Ovie Doro
Join Date: Jul 2009
Posts: 99
Rep Power: 17
ovie is on a distinguished road
@ 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 is offline   Reply With Quote

Old   October 17, 2009, 08:57
Default
  #15
Member
 
Edin Berberovic
Join Date: Mar 2009
Posts: 31
Rep Power: 17
eberberovic is on a distinguished road
ovie,
this is correct, congratulations.
you are still doing linear interpolation of kappa, but I guess its influence is very small.
enjoy it.
eberberovic is offline   Reply With Quote

Old   October 17, 2009, 10:59
Default
  #16
Member
 
Ovie Doro
Join Date: Jul 2009
Posts: 99
Rep Power: 17
ovie is on a distinguished road
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.
ovie is offline   Reply With Quote

Old   October 18, 2009, 09:05
Default
  #17
Member
 
Edin Berberovic
Join Date: Mar 2009
Posts: 31
Rep Power: 17
eberberovic is on a distinguished road
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.
Elham and amolrajan like this.
eberberovic is offline   Reply With Quote

Old   October 18, 2009, 12:47
Default
  #18
Member
 
Eelco Gehring
Join Date: Mar 2009
Posts: 70
Rep Power: 17
feijooos is on a distinguished road
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
feijooos is offline   Reply With Quote

Old   October 18, 2009, 13:41
Default
  #19
Member
 
Ovie Doro
Join Date: Jul 2009
Posts: 99
Rep Power: 17
ovie is on a distinguished road
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.
ovie is offline   Reply With Quote

Old   October 18, 2009, 13:49
Default
  #20
Member
 
Eelco Gehring
Join Date: Mar 2009
Posts: 70
Rep Power: 17
feijooos is on a distinguished road
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.
feijooos is offline   Reply With Quote

Reply

Tags
interfoam energy, temperature field

Thread Tools Search this Thread
Search this Thread:

Advanced Search
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are Off
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
Getting a concentration field around a bubble in InterFoam azman OpenFOAM Running, Solving & CFD 3 June 7, 2022 05:21
Moving mesh Niklas Wikstrom (Wikstrom) OpenFOAM Running, Solving & CFD 122 June 15, 2014 07:20
Adding temperature field to InterFoam yapalparvi OpenFOAM Running, Solving & CFD 8 October 14, 2009 21:18
Problem with rhoSimpleFoam matteo_gautero OpenFOAM Running, Solving & CFD 0 February 28, 2008 07:51
Problems calculating field gh with interFoam cricke OpenFOAM Running, Solving & CFD 0 December 10, 2007 08:17


All times are GMT -4. The time now is 03:11.