Dear Sirs,

I know OpenFOAM is tunned to achieve better performance than comercial codes in good meshes, and I can gladly see it do so.

However I am having some problems with convergence with kepsilon model (also tried LRR, and a v2f I implemented myself - still testing). I dont think my mesh is that bad, but I would like to tunne my parameters so that I get a more robust solver.

Does anyone out there has any guidelines on schemes and parameters towards a more robust solver?

I have observed some issues regarding kepsilon, for instance:

1- if I set the tolerance for k and epsilon equation to 10^-5, my solution blows up. Setting it to 0.001 corrects it, but it limits my accuracy to 0.001, where I would like it to be much better.

2- velocity usually converge faster, while pressure residuals stagnate at some level. Is there a way to modify those parameters so that pressure converges faster?

3- I tried to use other solvers, since I thought I might be having trouble because of coupling between the equations (mostly in my v2f implementation). Is it right to say a GaussSeidel solver is more robust when there is a strong coupling than a BICCG? I set number of sweeps to 1 and relative tolerance to 0.999 to guarantee that I was "mixing well" the equations (one sweep for each variable). Is there a misconception? (I know it is slower). Do you think BICCG could do as a good job as GaussSeidel in strongly coupled systems, while still being much faster? (Would you limit the number of iterations the same way I did for GS?)

I apreciate any comments and help.
Luiz

1) I have not seen this behaviour before and without knowing more about your case and setup I don't think I can help much except to offer these general guidelines:

First it's a very good idea to start by running potentialFoam, this is a good way of checking the BCs as well as generating a sensible starting U field.

What you really need next is a decent estimate for the inital turbulence property fields but we haven't yet written such a utility so what we tend to do is start the run with very low under-relaxation on these fields, say 0.05 or even lower. After a few iterations this can be raised to a normal level. If you are running transient you might find it useful to start by running simpleFoam first to evolve the turbulence fields to a sensible distribution.

After doing these things you should find you can have the tolerance on k as tight as you like.

2) We don't normalise the residuals the same way as other codes and one consequence is that p may appear not to converge adequately. However, we find that so long as U has converged well it doesn't matter that the p residual is higher. The best thing is to look at the way the fields are evolving.

3) I don't think there is any reason to use GaussSeidel over BICCG, the later is a much more efficient solver and we have never experienced problems with it.

P.S. I have just finished an initial implementation of the realizable k-epsilon model which I will include in the next release. Would you like to contribute your v2f implementation?

Sure! I was already ready to contribute, as soon as I clean out some coments and some dirty implementations. Also it did not passed my tests (mostly because of problems in convergence), but if you want I can send you right now.

Do you want now?

Sincerely yours,
Luiz

It's OK, it can wait until you are happy with it. If you could make it consistent in style with the models already in FOAM it will make my life easier because I am very particular about programming style and will go through it before I commit it to the release. I know I am a bit anal about programming style, it's just the way I am and as a consequence the way FOAM is :-)

I am doing my best. I used your kepsilon as a model for the classes, styles, etc...

I am still trying to get OpenFOAM fully functional in my computer. It used to work with FOAM, but when I migrated to OpenFOAM, I get segmentation fault. (The compilation was ok, but when I run the case, I got this message).

As soon as I solve this problem I will test your suggestions...

Thanks again,
Luiz.

PS- I had a (memory) problem with my password and did not have the opportunity to wish you guys a Marry Xstmas and HNY, as well as thank you for this wonderful code.

Hi

I am a new user of OpenFOAM and would like to learn how to hook my own module to it. Earlier I implemented v2f model in Fluent using UDF and it works fine. I would appreciate if you could share your v2f-OpenFOAM implementation with me and help me out in compilation procedure. Currently I am validating various turbulence models implemented in OpenFOAM with that in Fluent for backward facing step, wall and impinging jet cases. I would also like to know how to invoke energy equation, so that I can validate for heat transfer cases.

Thanks and regards

GS

Hello Henry,

I am running simpleFoam with realizable k-epsilon(RKE) model to validate backward facing step case. RKE code is the same what you have written and the modifications required to make it visible in the GUI have been made. Now I can select it from the FoamX GUI. While running RKE, I am finding it difficult to get properly developed flow and results are not good even after running it for 12000 iterations. As compared to RKE standard k-epsilon(SKE) model takes just 1000 iteratios. Also, I have looked at both the codes i.e. RKE as well as SKE, following differences have been observed:

RKE epsilon and k equations:
------------------------------------------
tmp<fvscalarmatrix> epsEqn
(
fvm::ddt(epsilon_)
+ fvm::div(phi_, epsilon_)
- fvm::laplacian(DepsilonEff(), epsilon_)
==
C1*magS*epsilon_ + boundarySource
- fvm::Sp
(
C2*epsilon_/(k_ + sqrt(nu()*epsilon_)) + boundaryCentral,
epsilon_
)
);

epsEqn().relax();
solve(epsEqn);
bound(epsilon_, epsilon0_);


// Turbulent kinetic energy equation

tmp<fvscalarmatrix> kEqn
(
fvm::ddt(k_)
+ fvm::div(phi_, k_)
- fvm::laplacian(DkEff(), k_)
==
G - fvm::Sp(epsilon_/k_, k_)
);
-------------------------------------------

SKE epsilon and k equations:

tmp<fvscalarmatrix> epsEqn
(
fvm::ddt(epsilon_)
+ fvm::div(phi_, epsilon_)
- fvm::Sp(fvc::div(phi_), epsilon_)
- fvm::laplacian(DepsilonEff(), epsilon_)
==
C1*G*epsilon_/k_ + boundarySource
- fvm::Sp(C2*epsilon_/k_ + boundaryCentral, epsilon_)
);

epsEqn().relax();
solve(epsEqn);
bound(epsilon_, epsilon0_);


// Turbulent kinetic energy equation

tmp<fvscalarmatrix> kEqn
(
fvm::ddt(k_)
+ fvm::div(phi_, k_)
- fvm::Sp(fvc::div(phi_), k_)
- fvm::laplacian(DkEff(), k_)
==
G
- fvm::Sp(epsilon_/k_, k_)
);
----------------------------------------

Why do we have terms "- fvm::Sp(fvc::div(phi_), epsilon_)" and "- fvm::Sp(fvc::div(phi_)", k_) in SKE but not in RKE. I have looked at RNG model also, and these terms are absent there as well.

In the file "realizableKE.H" we have :

---------------------------------------------
tmp<volscalarfield> rCmu
(
const volTensorField& gradU,
const volScalarField& S2,
const volScalarField& magS
);

tmp<volscalarfield> rCmu
(
const volTensorField& gradU
);

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

Why do we need to define "gradU" two times. Kindly excuse me if I am asking stupid questions. I am new to C++ and OOPS.

"Cmu" in RKE is computed as "rCmu", therefore do we really need to read "Cmu"?

The other constants like "A0", "C2", "alphak" and "alphaEps" are passed correctly.

Regards
GS

Hello Henry,

I am running simpleFoam with realizable k-epsilon(RKE) model to validate backward facing step case. RKE code is the same what you have written and the modifications required to make it visible in the GUI have been made. Now I can select it from the FoamX GUI. While running RKE, I am finding it difficult to get properly developed flow and results are not good even after running it for 12000 iterations. As compared to RKE standard k-epsilon(SKE) model takes just 1000 iteratios. Also, I have looked at both the codes i.e. RKE as well as SKE, following differences have been observed:

RKE epsilon and k equations:
------------------------------------------
tmp<fvscalarmatrix> epsEqn
(
fvm::ddt(epsilon_)
+ fvm::div(phi_, epsilon_)
- fvm::laplacian(DepsilonEff(), epsilon_)
==
C1*magS*epsilon_ + boundarySource
- fvm::Sp
(
C2*epsilon_/(k_ + sqrt(nu()*epsilon_)) + boundaryCentral,
epsilon_
)
);

epsEqn().relax();
solve(epsEqn);
bound(epsilon_, epsilon0_);


// Turbulent kinetic energy equation

tmp<fvscalarmatrix> kEqn
(
fvm::ddt(k_)
+ fvm::div(phi_, k_)
- fvm::laplacian(DkEff(), k_)
==
G - fvm::Sp(epsilon_/k_, k_)
);
-------------------------------------------

SKE epsilon and k equations:

tmp<fvscalarmatrix> epsEqn
(
fvm::ddt(epsilon_)
+ fvm::div(phi_, epsilon_)
- fvm::Sp(fvc::div(phi_), epsilon_)
- fvm::laplacian(DepsilonEff(), epsilon_)
==
C1*G*epsilon_/k_ + boundarySource
- fvm::Sp(C2*epsilon_/k_ + boundaryCentral, epsilon_)
);

epsEqn().relax();
solve(epsEqn);
bound(epsilon_, epsilon0_);


// Turbulent kinetic energy equation

tmp<fvscalarmatrix> kEqn
(
fvm::ddt(k_)
+ fvm::div(phi_, k_)
- fvm::Sp(fvc::div(phi_), k_)
- fvm::laplacian(DkEff(), k_)
==
G
- fvm::Sp(epsilon_/k_, k_)
);
----------------------------------------

Why do we have terms "- fvm::Sp(fvc::div(phi_), epsilon_)" and "- fvm::Sp(fvc::div(phi_)", k_) in SKE but not in RKE. I have looked at RNG model also, and these terms are absent there as well.

In the file "realizableKE.H" we have :

---------------------------------------------
tmp<volscalarfield> rCmu
(
const volTensorField& gradU,
const volScalarField& S2,
const volScalarField& magS
);

tmp<volscalarfield> rCmu
(
const volTensorField& gradU
);

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

Why do we need to define "gradU" two times. Kindly excuse me if I am asking stupid questions. I am new to C++ and OOPS.

"Cmu" in RKE is computed as "rCmu", therefore do we really need to read "Cmu"?

The other constants like "A0", "C2", "alphak" and "alphaEps" are passed correctly.

Regards
GS

[henry]

The RKE model is the least-well tested but worked very well for the cases I have tried including the backward-facing step.

The fvm::Sp(fvc::div(phi_), epsilon_) and fvm::Sp(fvc::div(phi_), k_) terms are optional. Note that at convergence fvc::div(phi_) = 0 so these terms only have any effect early in the solution procedure while continuity is violated due to incomplete solution of the pressure equation.

> Why do we need to define "gradU" two times

Sorry, I don't understand the question.

> "Cmu" in RKE is computed as "rCmu", therefore do we really need to read "Cmu"?

That depends on wether you want to use Cmu or rCmu in the wall-functions, this issue is not properly addressed in the papers and I opted for Cmu.

Could your convergence problems relate to the choice of outlet boundary condition? What are you currently using?

Thanks Henry for your prompt reply.

>Sorry, I don't understand the question.

In the file "realizableKE.H" we have :

---------------------------------------------
tmp<volscalarfield> rCmu
(
const volTensorField& gradU,
const volScalarField& S2,
const volScalarField& magS
);

tmp<volscalarfield> rCmu
(
const volTensorField& gradU
);

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

In the above I meant why do we need

tmp<volscalarfield> rCmu
(
const volTensorField& gradU
);

Regarding BC I am using outlet boundary condition at the exit.

Regards
GS

[henry]

Look at the calls to the two rCmu functions and what they do and you will see why there are two, one calling the other -- it saves code duplication.

Try using the inletOutlet BC at the outlet and specify (0 0 0) as the inletValue for U and whatever you are using as the inlet values for k and epsilon.