Well I have read through H.Rusche thesis and the k-epsilon model used in the thesis is Gosman's model. I read your earlier posts as well as the bubbleFoam wiki and I can see that Gosmans model is not used in bubbleFoam. My question is why are the source terms for k and epsilon, S_k and S_eps not included in bubbleFoam?

Also my understanding is that Gosman's model is used in twoPhaseeulerFoam?

Kind Regards

Matt

They seem to implement a simpler form of the model than the one from Gosman (in Rusche's thesis they mention that Gosman's model might cause numerical problems).

Thank you for the responses Alberto your help is always much appreciated.

I think it may be worth my while to have a look more closely at Gosman's model to hopefully determine what his source terms represent. I am currently looking into running the bubbleFoam/twoPhaseEulerFoam bubble column tutorial with turbulence and I cannot seem to get it working, just trying to broaden my understanding of the k-epsilon model. I think I have made the appropriate changes to fvSolutions and fvSchemes (added the variables to be solved and the new laplacian terms for the turbulence equations). However my solution blows up after very few timesteps.
I am at the moment trying to get the tutorial running with turbulence before I move on to my specific tank geometry (tank with parabolic profile)

Matt

Do you have parts of the system where alpha = 1? If so, the turbulence model implemented in twoPhaseEulerFoam/bubbleFoam is unstable in that limit, since beta -> 0.

You might want to consider, depending on your system, a mixture formulation, or to make some change to the implementation to manage the case beta -> 0 more robustly. If you need info on this, let me know.

Best,

Hi Alberto,

I noticed that you have already mentioned the method for making changes to the solver...

Are bubbleFoam and twoPhaseEuler the same in terms of turbulence models? Have these 2 solvers undergone much change since previous versions?

There are areas of my domain where alpha is equal to 1, so I will definately need to make some modifications here.

Kind Regards

Matt

New developments in terms of multi-fluid models came with 2.x, but not in terms of turbulence modeling.

For the limit of beta -> 0, you can probably define a minimum value of alpha in some term of the equation (see what they do in compressibleTwoPhaseEulerFoam for the drag), or drop the cells where beta < small from the solution, and fix the value of the variable there (risky approach in terms of robustness).

Best,

Hi Alberto,

I have finally attempted some modifications to the bubbleFoam solver to implement the mixture turbulence model. I have a question relating to the implementation of the first divergence term in the equations.. I am not sure what this term means?

Also I have not included the source terms yet as they are not implemented in bubblefoam as standard, but I would eventually like to implement these once I get it working without them.

Kind Regards

Matt

I am not sure about what equations you are considering. The "standard" mixture turbulent model is usually identical to the single-phase turbulence model, but considering the properties of the mixture. In other words:

rho -> rhoMix
U -> UMix
mu -> muMix (molecular viscosity)

and the resulting k and epsilon belong to the mixture. Please note that the same consideration applies to wall functions.

The divergence term div(rho_mix*U_mix*k_mix) is the convective term. In bubbleFoam U_mix is simply the face velocity of the mixture phi, and rho_mix is "rho = alpha*rhoa + beta*rhob", which should be available already.

What reference are you using?

Best,

sorry about that I didn't even mention.. I am using Rusche PhD thesis

the equations I am referring to are equations 3.66 and 3.67

I have circled the div terms which I am unsure about in the attached image.

Cheers

Matt

Those terms are what is in the code already:

fvm::div(phib, k)

and

fvm::div(phib, epsilon).

The next divergence term has to be discretized as

-fvm::Sp(fvc::div(phib), k)

Best,

I have added the extra term as you suggested and re-compiled but now I get the following error message:

--> FOAM FATAL ERROR:
incompatible dimensions for operation
[epsilon[0 2 -4 0 0 0 0] ] - [epsilon[0 0 -1 0 0 0 0] ]

From function checkMethod(const fvMatrix<Type>&, const fvMatrix<Type>&)
in file /opt/openfoam201/src/finiteVolume/lnInclude/fvMatrix.C at line 1278.

First, let me correct one typo: I wrote "Su" instead of "Sp".

If you use

fvm::div(phib, epsilon) - fvm::Sp(fvc::div(phib), epsilon)

it should work out, since the units are:

[phi] = m/s
[epsilon] = m^2/s^3

so the first term is m^2/s^4 (accounting for the units introduced by the divergence). The second term gives

[div(phib)] = 1/s
[epsilon] = m^2/s^3

so again m^2/s^4 as before.

Best,

Thanks a lot Alberto,

I have compiled the solver again, so basically equations 3.66 and 3.67 (Rusche) are implemented as my turbulence model without the source terms.

The solver now runs on the bubbleColumn tute, after halving the time step, for the full 20 sec however after about 10 sec of the solution the flow takes on the shape in the image and does not change.. Should I be using a different scheme for the k-epsilon terms?

I'm not sure what is causing this.
The settings are the same as the bubbleColumn tute with the following changes:
-added the k-epsilon terms to the fvSchemes dict
-halved time step to 0.001s
-turned correctAlpha to yes

Kind Regards


Matt

Could you check the values of turbulent viscosity?

turbulent viscosity was taken from Paraview..

This seems quite normal.

Hi Alberto,

So I have made some modifications to the createRASTurbulence.H file, namely the alphaEps coefficient (I believe this is the Schmidt number in the dissipation equation).
It is set to 0.76923 originally in the solver however Rusche uses 1.3.

My problem is that once I change it to 1.3 and recompile the solver blows up.. I have tried running the solver on a case that I had already run without changing alphaEps and it did not solve.

I also have not been able to figure out why I am having problems I had posted earlier with the flow! I am guessing that I am missing something in the code. I only changed parts of the section of kEpsilon.H as follows:


// Dissipation equation
fvScalarMatrix epsEqn
(
fvm::ddt(epsilon)
+ fvm::div(phib, epsilon)
- fvm::Sp(fvc::div(phi), epsilon)
- fvm::laplacian
(
nuEffb/alphaEps, epsilon,
"laplacian(DepsilonEff,epsilon)"
)
==
C1*G*epsilon/k
- fvm::Sp(C2*epsilon/k, epsilon)
);

#include "wallDissipation.H"

epsEqn.relax();
epsEqn.solve();

epsilon.max(dimensionedScalar("zero", epsilon.dimensions(), 1.0e-15));


// Turbulent kinetic energy equation
fvScalarMatrix kEqn
(
fvm::ddt(k)
+ fvm::div(phib, k)
- fvm::Sp(fvc::div(phi), k)
- fvm::laplacian
(
nuEffb/alphak, k,
"laplacian(DkEff,k)"
)
==
G
- fvm::Sp(epsilon/k, k)
);


Does this look ok to you?

Kind Regards


Matt

There was a typo that I corrected in my previous post: fvm::Sp(fvc::div(phi), k) should be fvm::Sp(fvc::div(phib), k). Same for epsilon: you should use phib.

Also, you don't need to re-compile the code to change alphaEps. Just specify it in the appropriate dictionary in the constant/RASProperties dictionary.

Best,

Hi, just a little question: in the latest implementation of twoPhaseEulerFoam i see the term above discussed:

fvm::Sp(fvc::div(phib),k) and the same for eps equation.

My question is: does this term have any importance for incompressible simulations?

In the previous versions of the code it isn't implemented in the incompressible version.

Thanks,

Lapo

Hi Raghavendra

I am useing twoPhaseEulerFoam to simulate the Becker's case as well. I am wondering if you fixed your increased bubble plume oscillation period issue?

Thanks
Jie