[matt.mech.eng]

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

[alberto]

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).

[matt.mech.eng]

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

[alberto]

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,

[matt.mech.eng]

Originally Posted by alberto

Hi,

twoPhaseEulerFoam implements the turbulence mixture model of Gosman (see H. Rusche thesis).The implementation is done completely in the code, where you can find also the implementation of the wall functions.

To implement a new model:

- Remove the current model equations, or add a switch to decide which turbulence model you want to use. Equations are in twoPhaseEulerFoam/turbulenceModel/kEpsilon.H. In the same directory you find the headers where wall functions are implemented.

- Assuming all turbulence models you want to try rely on the hypothesis of turbulence viscosity, you do not have to deeply change the structure of the solver. Simply change the parts where the turbulent and effective viscosities are computed.
For explicit closures, this is all what you need to do.

- For models involving transport equations, you should code them, following the example of the k-eps equations already in the code.

As a side note, to compile twoPhaseEulerFoam, you need to run ./Allwmake because there are additional classes to be compiled.

Best,
Alberto

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

[alberto]

Originally Posted by matt.mech.eng

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,

[matt.mech.eng]

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

[alberto]

Originally Posted by matt.mech.eng

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?

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.

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.

What reference are you using?

Best,

[matt.mech.eng]

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

[alberto]

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,

[matt.mech.eng]

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.

[alberto]

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,

[matt.mech.eng]

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

[alberto]

Could you check the values of turbulent viscosity?

[matt.mech.eng]

turbulent viscosity was taken from Paraview..

[alberto]

This seems quite normal.

[matt.mech.eng]

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

[alberto]

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,

[Lapo]

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

[jiejie]

Originally Posted by raagh77

Hi Alberto

Am back to this thread !

I just finished the 2D simulation of bubble column, becker case (with turbulence being on)

I would like to discuss the results here

.> The velocity variations at three positions after time averaging was quite compariable
.> K-epsilon model predicted higher turbulent viscosity (nutb)
.> Timeperiod of the lateral bubble movement around 70seconds which was far more greater than what specified in the paper (paper says 16 to 20seconds)

This increase in timeperiod is either because of high prediction of nutb by k-epsilon or due to the relaxation parameters used for twoPhaseEulerFoam solver (because I remember you saying with "realaxtion parameters its difficult to obtain time accruate predictions") ??

Regards
Raghavendra

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