buoyancy augmented k-eps model
 

Dear all!

As I m interested in buoyancy driven flow I tried to adapt the standard k-eps model in order to account for buoyancy effects in the k and eps equation. I added the following source terms (see e.g. [1,2]):

_k-eqn: P

_eps-eqn: eps/k*C_4*max[P,0]

where P is calculated as follwos (y is the vertical direction; d partial derivative; g gravity):

P=g*(mu_eff/(sigma_h*rho))*drho/dy

My implementation in OF looks like this:

dimensionedScalar gMAG
(
"gMAG",
dimensionSet(0, 1, -2, 0, 0),
9.81
);
volVectorField gradRho = fvc::grad(rho_);
volScalarField P = (1/rho_)*gMAG*alphah_*mut_*gradRho.component(1);

additional source term in k-eqn:

fvm::SuSp(P/k_, k_)

additional source term in eps-eqn:

C4_*posFunct(P)*epsilon_/k_


posFunct() is defined as shown below:

volScalarField posFunct(volScalarField& P)
{
volScalarField P_pos = P;
for(volScalarField::iterator iter = P_pos.begin(); iter != P_pos.end(); ++iter)
{
if(*iter < 0)
*iter = 0;
}

return P_pos;
}

When I compare the results obtained with the modified k-eps model mentioned above and the standard k-eps model (used in chtMultiRegionFoam; set-up see http://www.cfd-online.com/Forums/ope...egionfaom.html), the differences in k, epsilon as well as in T and U are minor.

Has somebody already made experience with the implementation of a buoyancy augmented k-eps model?? What am I doing wrong here?? I greatly appreciate your comments!!

Thanks in advance,
Aram

[1] Novozhilov: "CFD Modeling of Compartment Fires"
[2] S. Nam and G. Bill: "Numerical Simulation of Thermal Plumes"

Hi all!!

In reference [1] it can be read that the buoyancy modification mentioned in the previous post do not improve the behavior of the standard k-eps model greatly. Hence, they propose a buoyancy modification based on the General Gradient Diffusion Hypothesis:

Add the following source terms:

_k-eqn: P

_eps-eqn: C1*epsilon_/(1-C4)*P/k_

where P is calculated as follwos (d partial derivative; g gravity, rhoRef reference density):

P=3*rhoRef*mu_t*g*(u'v'*drho/dx + v'v'*drho/dy + v'w'*drho/dz)/(k*rho^2*sigma_h)

I tried to implement it in the following way:

volVectorField gradRho = fvc::grad(rho_);
volVectorField diffGrad = R() & fvc::grad(rho_);
volScalarField P = (3.0/2.0)*rhoRef_*gMAG_*alphah_*mut_*diffGrad.component (1)/(k_*sqr(rho_));

additional source term in k-eqn:

+ fvm::SuSp(P/k_, k_)

additional source term in eps-eqn:

+ fvm::SuSp(C1_*(1-C4_)*P/k_, epsilon_)

The code compiles without a problem but when running a case I get bounding k and epsilon and after a while the calculation crashes. I assume a problem with the linearisation of the source terms. That s what I m currently investigating but could not figure out any solution. At the same time I m experimenting with the descretisation schemes.

I would greatly appreciate any comments!! Thx in advance,
Aram

[1] Van Maele and Merci: "Application of two buoyancy modified k-eps turbulence models to different types of buoyant plumes"

Hi!!

I was playing around with the discretisation/linearisation of the source terms and did some changes. The implementation showed below allows the simulation to run longer but crashes at time 6s (instead of 0.2s) due to bounding k (no bounding of epsilon). The so far obtained flow field looks plausible.

volScalarField P = (3.0/2.0)*rhoRef_*gMAG_*alphah_*mut_*diffGrad.component (1)/sqr(rho_);

No inclusion of k for calculating P; this is done in the source terms, so the source term for the dissipation equation is discretised in the following way:

+ fvm::SuSp(C1_*(1-C4_)*P/sqr(k_), epsilon_)

Befor discretising the source for k, I linearised it:

S_k = P/k = 2*P/k - (P/k^2)*k

and hence the source for the k equation is discretised as follows:

+ 2*P/k_ - fvm::Sp(P/sqr(k_), k_)

I would greatly appreciate your comments as I m not sure if this makes sens; at least it allowed my simulation not to blow up imediatly ...

Thx in advance!!
Aram

Dear all!!

I played around with different ways of formulating and implementing the source terms for my buoyancy modified k-eps model and was able to set up a non-crashing simulation. In the new implementation I first substituted mu_t in the source term P with its definition for the k-eps model. Then I set v'v' = k (suggested in [1]). So I get:

P = (3/2)*g*(rhoRef/rho)*(Cmu/sigma_h)*(k/epsilon)*(u'v'*drho/dx + k*drho/dy + v'w'*drho/dz)

The implemetation looks as follows:

volSymmTensorField RS = R();
volScalarField diffGrad = RS.component(1)*gradRho.component(0) + k_*gradRho.component(1) + RS.component(4)*gradRho.component(2);
volScalarField P = (3.0/2.0)*gMAG_*(rhoRef_/rho_)*alphah_*Cmu_*diffGrad;

I did not include the quotiont "k/epsilon", so it cancels out for the additional source term in the dissipation eqn.:

S_eps = C1*(1-C4)*(epsilon_/k_)*(k/epsilon)*P = C1*(1-C4)*P

I tried out two ways

+ fvm::SuSp(C1*(1-C4)*P/epsilon_, epsilon_)

and

+ C1*(1-C4)*P

the latter showed better behaviour. The additional source term for the k eqn. is added like:

+ fvm::SuSp(P/epsilon_, k_)

The simulation with the coarsest grid is running but I encounter some serious problems though:

_ a strong sensitivity with respect to the domain size. Simulations crash when domain size gets too small or too big.
_ sensitivity with respect to the initial value of epsilon (varied it by one to two orders of magnitude) where only 2.9e-5 (for my case) reuslts in a stable simulation.
_ when I refine the grid (from a grid spacing of 0.025 m to 0.01 m) the simulation crashes after 3s (my main problem now!).

I greatly appreciate any comment on my problem!! Thx in advance!!

Aram

[1] Van Maele and Merci: "Application of two buoyancy modified k-eps turbulence models to different types of buoyant plumes"

Hi!!

I m still struggling with the buoyancy modification of the standard k-eps model, even though I made some improvements. It turned out that the code gets a bit more stable when the density gradients are replaced with temperature gradients by the help of the ideal gas law. In doing that the production term reads:

P = - (3/2)*g*(rhoRef/T)*(Cmu/sigma_h)*(k/epsilon)*(u'v'*dT/dx + k*dT/dy + v'w'*dT/dz)

My implementation looks like:

volScalarField T_ = thermophysicalModel_.T();
volVectorField gradT = fvc::grad(T_);
volSymmTensorField RS = R();
volScalarField diffGrad = RS.component(1)*gradT.component(0) + k_*gradT.component(1) + RS.component(4)*gradT.component(2);

volScalarField P = (3.0/2.0)*gMAG_*(rhoRef_/T_)*alphah_*Cmu_*diffGrad;

(no inclusion of - (k/epsilon)). Additional source for the dissipation eqn.:

S_eps = - fvm::SuSp(C1_*(1.0-C4_)*P/epsilon_, epsilon_)

Additional source for the k eqn.:

S_k = - fvm::SuSp(P/epsilon_, k_)

Again my main problem is stability (bounding eps and k, than crash), specially when I refine the grid. The coarse grid runs (better now), wherease the fine crashs after a few seconds. My impression is that the source term in the dissipation eqn. plays a crucial role ....

Any comments are of great help for me!! Thx in advance and all the best,
Aram

1-could you finally solve the problem ?

2-could you share the code?

Please see attached files for the code which I used with OF-1.7.1.



Cheers,
Aram

Thank you for sharing. Have you update your code for later Openfoam version?
And how is the effect of adding buoyancy term?