Introducing a scalar field for buoyancy production in k-epsilon
 

Hi dear FOAMers!
I want to add the gradient of a scalar filed (like S) for calculation the production term due to buoyancy in k and epsilon equations. So I can make the new turbulence model by defining the coefficients and gravitational acceleration which are needed for calculating the buoyancy production. But I have no idea how should I do the procedure to make the modified k-epsilon model read the scalar values from the case?!?
In brief:
production due to buoyancy=Gb=beta*C3e*g*(dS/dy)

Problem is dS/dy , S is the scalar filed which is advected in the domain by velocity.

Best REGARDS:)

nothing?:confused::confused:

There already is a buoyant kEpsilon model : https://github.com/OpenFOAM/OpenFOAM...yantKEpsilon.H. If, after checking out its implementation, it is not what you want it should be a good starting point for implementing your own variation.

Caelan

Hi Clealan!
Thank you for your help. I am already dealing with that. My another question is about the thermal expansion coefficient which is appearing in buoyancy tern (G) called beta. Is there any equivalent for this beta for scalar flux?? Because most of the literature is devoted to thermal plume and beta is adopted for buoyancy due to non-zero temperature gradient field.
Regards Sina.

G=(beta)*gi*(nut/Prt)*(DT/Dxi).
T is temperature or scalar filed. But beta is only defined for temperature not any scalar

const volScalarField& S_ = U_.mesh().lookupObject<volScalarField>("S");

if you have defined the scalarField somewhere you can retrieve it with the above function

Hi Michael Alletto!
thank you for your answer yes exactly i figured this out. but do you have any idea about what i asked in my response to another friend Clealan:


My another question is about the thermal expansion coefficient which is appearing in buoyancy tern (G) called beta. Is there any equivalent for this beta for scalar flux?? Because most of the literature is devoted to thermal plume and beta is adopted for buoyancy due to non-zero temperature gradient field.
Regards Sina.

G=(beta)*gi*(nut/Prt)*(DT/Dxi).
T is temperature or scalar filed. But beta is only defined for temperature not any scalar

Which scalar are you interested in? The term accounting for the buoyancy production in the k equation comes into the equation when the buoyancy is present also in the momentum equation. If one derives the equation for k from the momentum equation with the term accounting for the buoyancy one gets the terrm rho* beta * gi * ui'T'. So if you try to derive the k equation from the momentum equation with your scalar you should see which constant you need instead of beta.

Dear Michael!
Thank you for your response again, I understand the general formulation of buoyancy and its relation with momentum equations. But the question is still remaining. My scalar filed is Salinity which is entering the domain due to momentum as a jet flow. So the density is affected by salinity with the Millero-Poission eqn of state (P=f(S,T)). So what should I do for beta, when I am dealing with salinity concentration as a scalar of my domain?
Regards

A ok I got it. Ok the beta in the momentum equtions comes from the Boussinesq Approximation, i.e. you linearize the density around a a given value and expand it as a function of temperature. The same could be done for the salinity: you can expend it around a given salinity value (see http://www.o3d.org/eas-ocean-modelin...3-EqMotion.pdf)

so i assume you can substitude the coefficient of termal expansion beta by the coeffcient of salinity concentration, i.e



this is what I assume after a short search about change of density as a function of salinity.

Hope this is usefull

Hi!
Thank you for your research and quick response. I got your hints. Actually I am testing the results of this approach now. So for every specific scalar in a domain one should go after finding the relation between rho and that specific scalar to find a similar manner of beta? Is that true?
Regards, Sina.

In some manner. You should really write down the derivation of the transport equation for the turbulent kinetic energy equation considering buoyancy to understand if everything is correct (see e.g. https://www.cambridge.org/core/books...A39631674EC3C9).

In principle what you do if you have a source term in the momentum equation: you decompose it into mean and fluctuation and multiply it with the velocity velocity fluctuation and you end up something like this:

thnaks a lot. I will go further with your helps and hints and I will notice the results here again.
Regards Sina

Just one question out of curiosity:

What are the applications you are doing the modeling?

of coures. I am working on neagitve buoyant jets which are buoaynat due to salinity of the jet.
can I ask another question ? I want to make a if statement with < and > :
if ( a> ...)

do...
else
do...
but the error is concerned about "<" and ">" and value which I entroduced for comparison :a> 0 (the zero)! how come?
I am doing this in the .C of my solver.

the last part of the error is :
' Foam:: tmp<Foam::GeomtericField<double, Foam::fvpatchField, Foam::volMesh> >' is not derived from 'const Foam::Pair <Type>'
Regards Sina

Form the short code you sent i presume you want to compare a volumetric viele with a float. But hard to say without and code.

More specific:



if (U.component(1) > 0)

var1=1.0;
else

var1=0.5;




and I want to consider U.component(1) as Uy (U.component(1)=Uy) of velocity field.
Regards Sina

What you are doing ist comparing a volScalarField with a float.

You should write

But with a bit oft efford you will find a lot oft threads whitch deal with the topic oft Looping over fields in OF.

Michael

I tried this now and I also checked the synax you wrote in other posts and discussions about forAll loops, the synax is correct but I get two error:
in brief:
..... has no member named 'size'
for (Foam::label i=; i<(list).size(); i++)


....in expansion of macro 'forAll'
forAll (U.component(1) , Celli)
^

...no match for operator []' ...
if (U.component(1)[Celli] > 0)
^

(***this ^ is under the [ ***)
Regards Sina

Maybe
forAll (U , Celli)

U[Celli].component(0)

Thanks a lot. That was correct. solver compiled with no error:):)

Greeting!
During running a test case with my modified solver (pimpleFoam) and modified LRR turbulence model, I get an error :symbol lookup error:Undefined symbol .....
and it is referring to dimensioned Vector
I just defined g (9.81) as a dimensionedVedtor in LRR.C and it compiled with no error. But it seems the solver is unfamiliar with that!
I should say that the solver is working without problem with other cases.
Exactly before calculation the R (Reynolds stress tensor) during the run, it jumps out. Inside the LRR.C I exactly before Reynolds stress eqns added that g (9.81) definition to use it as a vector in calculation of buoyancy term (Gb) to add that in to the Reynolds strees eqns.
Regards
Sina

The problem changed to another thing:
no longer previous error! But I can't switch the wallReflection off now! The new modified LRR model is working properly, but it can't distinguish that the wallReflection is switched off! Is it coming from compilation procedure?

I am now using Gcoef in buoyantKE. I am wondering the meaning of this->alpha_ in Gcoef and g should include the direction or not?
In epsilonsource/ksource:

return -fvm::Susp(...)

If it return -.... in code, G should be -G*... But in equation, it should be +G*...

I am confused about it. Thank you.

alpha is the phase fraction. So this is used for two phase flows.


g is a vector and it has a direction: https://github.com/OpenFOAM/OpenFOAM...Epsilon.C#L104
And it is used in scalar products.


SuSp does add a source implizitly to an solver equation if it is positive, because this increases the stability of the solution. This means it is added to the diagonal of the matrix to be solved. If it is negative, it will be added explicitly, so to the "right" side of the equation.


So it makes a difference, if you put something like -fvm::SuSp(-G, k) or +fvm::SuSp(G, k) to an equation, because it checks the sign of the first argument.

In this case isn't alpha thermal diffusivity? See eg : https://github.com/OpenFOAM/OpenFOAM...lDiffusivity.H.

Caelan

If you look for tutorial, where buoyantKEpsilon is used, you find them in:
tutorials/multiphase/driftFluxFoam/RAS/mixerVessel2D (dahl, tank3D)
These are multiphase cases, where alpha is used.


But to make sure, this is correct, I patched the buoyantKEpsilon/buoyantKEpsilon.C to output alpha dimension:

compare this to the dimensions of the thermal diffusivity:

https://github.com/OpenFOAM/OpenFOAM...usivity.H#L116

Good to know -- thanks for checking.

Caelan

Hi,

I don't think that the alpha in compressibleTurbulenceModels is the phase fraction.

alpha will be the phase fraction for the phaseCompressibleTurbulenceModels.
https://github.com/OpenFOAM/OpenFOAM...ulenceModels.C

However, for incompressible or compressible TurbulenceModels, alpha will be just a geometricOneField.
https://github.com/OpenFOAM/OpenFOAM...moModels.C#L30
https://github.com/OpenFOAM/OpenFOAM...tricOneField.H

In incompressibleTurbulenceModels, rhoField is also just a geometricOneField (uniform value of 1, without dimension).
When we use these models in multiphase-solver such as interFoam, density of fluid is not considered in the turbulenceModel. rho is 1 at everywhere.

Sofisticated and elegant source-code of turbulence models makes the meanings of alpha and rho in turbulence models a little unclear, I think.
We can get the both compressible and incompressible models from the same code in return for this.