[makaveli_lcf]

Hi everyone!

I discovered following while simulating buoyant flow with 1.7.0 and 1.7.1 versions (Ubuntu 10.04 LTS packages from apt):

1. Comparable to 1.7.0 in version 1.7.1 kappat field is introduced in TEqn for incompressible SIMPLE and PIMPLE solvers.

v1.7.0:

{
    volScalarField kappaEff
    (
        "kappaEff",
        turbulence->nu()/Pr + turbulence->nut()/Prt
    );

    fvScalarMatrix TEqn
    (
        fvm::ddt(T)
      + fvm::div(phi, T)
      - fvm::laplacian(kappaEff, T)
    );

    TEqn.relax();
    TEqn.solve(mesh.solver(T.select(finalIter)));

    rhok = 1.0 - beta*(T - TRef);
}

v1.7.1:

{
    kappat = turbulence->nut()/Prt;
    kappat.correctBoundaryConditions();

      volScalarField kappaEff("kappaEff", turbulence->nu()/Pr + kappat);

    fvScalarMatrix TEqn
    (
        fvm::ddt(T)
      + fvm::div(phi, T)
      - fvm::laplacian(kappaEff, T)
    );

    TEqn.relax();
    TEqn.solve(mesh.solver(T.select(finalIter)));

    rhok = 1.0 - beta*(T - TRef);
}

2. For kappat field Jayatille wall function is introduced, but it has a bug in kappat calculation which is fixed in v2.0.1:

diff for kappatJayatillekeWallFunctionFvPatchScalarField.C between 1.7.1 and 2.0.1 results in

210c210
<     const scalarField& nuw = rasModel.nu().boundaryField()[patchI];
---
>     const scalarField& nuw = rasModel.nu()().boundaryField()[patchI];
238c238
<             scalar kt = nu*(yPlus/(Prt_/kappa_*log(E_*yPlusTherm) + P) - 1/Pr);
---
>             scalar kt = nu*(yPlus/(Prt_*(log(E_*yPlus)/kappa_ + P)) - 1/Pr);
264c264,268
< makePatchTypeField(fvPatchScalarField, kappatJayatillekeWallFunctionFvPatchScalarField);
---
> makePatchTypeField
> (
>     fvPatchScalarField,
>     kappatJayatillekeWallFunctionFvPatchScalarField
> );

I noticed this issue because got a floating point operation error actually coming from
log(E_*yPlusTherm) where yPlusTherm can be equal zero.

You cal also notice that parenthesis around (log(E_*yPlus)/kappa_ + P)) result in different expression
than Prt_/kappa_*log(E_*yPlusTherm) + P.

Hope it would be helpful for someone))))

Cheers,
A.

[Oke'e]

Please can you help me with an explanation of what rhok is?

[makaveli_lcf]

http://en.wikipedia.org/wiki/Boussin...%28buoyancy%29

[Oke'e]

Thanks for the prompt response. I know about the boussinesque approximation. My question is because of a case I'm setting up in OpenFOAM. I am simulating stratification in a cylindrical tank. I have modified the 'hotroom' example under the bouyantBoussinesqPimpleFoam solver and changed the entries in the 'transportProperties' file to those of water. When the run the case with the original entries (which I suspect to be for air), I get realistic results - the stratification is obvious. But when I change the properties - Pr, beta, nu, to those of water the results when viewed show no signs of stratification. Could you be of help please?

[makaveli_lcf]

Hello Howard!

It is always better when you ask detailed questions. Because for the obvious ones you can find answers yourself and it is more useful for you to learn finding answers.

Would you please post you pictures which do not satisfy you and also the parameters which you used?

Cheers,
A.

[makaveli_lcf]

Do you mean simulating compressible flow? Of couse for water (almost incompressible liquid) there is no such stratification as for e.g. air.

[Oke'e]

Hello Dr. Alexander,

Please pardon the delay in my response. I travelled out of town without my laptop and was without access to computer facilities.

The changes I made in the transportProperties file were

nu : changed from 1e-05 to 4.5e-07 (m2sec-1)
beta: changed from 3e-03 to 5e-04 (K-1)
Pr: changed from 0.9 to 3.0
Prt: changed from 0.7 to 2.3.
Tref was left unchanged.

The RASModel keyword in the RASProperities file was changed to laminar.

In the 0/p file the wall pressures were set with the $internalField entry, while the p_rgh value was set to constant 0 in the 0/p_rgh file.

You'll see from the attached picture that there is no evidence of stratification. Stratification is obtainable in incompressible flows (the 'buoyantBoussinesqPimpleFoam' solver is for "Transient solver for buoyant, turbulent flow of incompressible fluids" according to the user guide.) This is the basis for stratified hot water storages which is what I'm trying to study.

Howard

[chiragomshanti]

Dear Dr.Alexander,
In the file "kappatJayatillekeWallFunctionFvPatchScalarField.C" kappatw shouldn't be kappat? I searched much, but didn't find any use of kappatw.
With Regards,
Chirag Patel

[mabinty]

dear all,

i wonder why the JayatillekeWallFunction is not available for the the calculation of alphat in the compressible RANS version. for now, in the compressible version alphat is evaluated as mut/Prt at the wall. i m thinking of implementing the procedure of the JayatillekeWallFunction for compressible RANS in alphatWallFunctions.

does it make sens or is alphatw = mutw/Prt resonable anyway??

@Chirag Patel

kappatw is taken from the expression

scalarField& kappatw = *this;

thanks in advance for your input!

cheers,
aram

[mabinty]

hi,

I though about something like this:

// Populate boundary values
scalarField& alphatw = *this;

forAll(alphatw, faceI)
{
       [ ... ]

        if (yPlus > yPlusTherm)
        {
            scalar mu = muw[faceI];
            scalar at = mu*(yPlus/(Prt_*(log(E_*yPlus))/kappa_ + P) - 1/Pr);
            alphatw[faceI] = max(0.0, at);
        }
        else
        {
            alphatw[faceI] = 0.0;
        }

in turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/alphatWallFunctions/alphatWallFunction/alphatWallFunctionFvPatchScalarField.C (OF 1.7.1)

cheers,
Aram

[makaveli_lcf]

hi Aram!

you are right) But I found an implementation in version OF 2.10:

see http://foam.sourceforge.net/docs/cpp...ce.html#l00195

00195 void alphatJayatillekeWallFunctionFvPatchScalarField::updateCoeffs()
00196 {
00197     if (updated())
00198     {
00199         return;
00200     }
00201 
00202     const label patchI = patch().index();
00203 
00204     const RASModel& rasModel = db().lookupObject<RASModel>("RASProperties");
00205 
00206     const scalar Cmu25 = pow025(Cmu_);
00207 
00208     const scalarField& y = rasModel.y()[patchI];
00209 
00210     const scalarField& muw = rasModel.mu().boundaryField()[patchI];
00211 
00212     const scalarField& alphaw = rasModel.alpha().boundaryField()[patchI];
00213     scalarField& alphatw = *this;
00214 
00215     const tmp<volScalarField> tk = rasModel.k();
00216     const volScalarField& k = tk();
00217 
00218     const fvPatchVectorField& Uw = rasModel.U().boundaryField()[patchI];
00219     const scalarField magUp(mag(Uw.patchInternalField() - Uw));
00220     const scalarField magGradUw(mag(Uw.snGrad()));
00221 
00222     const scalarField& rhow = rasModel.rho().boundaryField()[patchI];
00223     const fvPatchScalarField& hw =
00224         rasModel.thermo().h().boundaryField()[patchI];
00225 
00226     // Heat flux [W/m2] - lagging alphatw
00227     const scalarField qDot((alphaw + alphatw)*hw.snGrad());
00228 
00229     // Populate boundary values
00230     forAll(alphatw, faceI)
00231     {
00232         label faceCellI = patch().faceCells()[faceI];
00233 
00234         scalar uTau = Cmu25*sqrt(k[faceCellI]);
00235 
00236         scalar yPlus = uTau*y[faceI]/(muw[faceI]/rhow[faceI]);
00237 
00238         // Molecular Prandtl number
00239         scalar Pr = muw[faceI]/alphaw[faceI];
00240 
00241         // Molecular-to-turbulent Prandtl number ratio
00242         scalar Prat = Pr/Prt_;
00243 
00244         // Thermal sublayer thickness
00245         scalar P = Psmooth(Prat);
00246         scalar yPlusTherm = this->yPlusTherm(P, Prat);
00247 
00248         // Evaluate new effective thermal diffusivity
00249         scalar alphaEff = 0.0;
00250         if (yPlus < yPlusTherm)
00251         {
00252             scalar A = qDot[faceI]*rhow[faceI]*uTau*y[faceI];
00253             scalar B = qDot[faceI]*Pr*yPlus;
00254             scalar C = Pr*0.5*rhow[faceI]*uTau*sqr(magUp[faceI]);
00255             alphaEff = A/(B + C + VSMALL);
00256         }
00257         else
00258         {
00259             scalar A = qDot[faceI]*rhow[faceI]*uTau*y[faceI];
00260             scalar B = qDot[faceI]*Prt_*(1.0/kappa_*log(E_*yPlus) + P);
00261             scalar magUc = uTau/kappa_*log(E_*yPlusTherm) - mag(Uw[faceI]);
00262             scalar C =
00263                 0.5*rhow[faceI]*uTau
00264                *(Prt_*sqr(magUp[faceI]) + (Pr - Prt_)*sqr(magUc));
00265             alphaEff = A/(B + C + VSMALL);
00266         }
00267 
00268         // Update turbulent thermal diffusivity
00269         alphatw[faceI] = max(0.0, alphaEff - alphaw[faceI]);
00270 
00271         if (debug)
00272         {
00273             Info<< "    uTau           = " << uTau << nl
00274                 << "    Pr             = " << Pr << nl
00275                 << "    Prt            = " << Prt_ << nl
00276                 << "    qDot           = " << qDot[faceI] << nl
00277                 << "    yPlus          = " << yPlus << nl
00278                 << "    yPlusTherm     = " << yPlusTherm << nl
00279                 << "    alphaEff       = " << alphaEff << nl
00280                 << "    alphaw         = " << alphaw[faceI] << nl
00281                 << "    alphatw        = " << alphatw[faceI] << nl
00282                 << endl;
00283         }
00284     }
00285 
00286     fixedValueFvPatchField<scalar>::updateCoeffs();
00287 }
00288

I didn't go into details, but I think you can use it. It would be nice if you comment on what is the differece between what you suggest (I also would write similar code).

Cheers,
Alex

[mabinty]

hi Alex,

thanks for the hint (still using OF 1.7). I had a look on the implementation of OF 2.1. Three main differences in the why how i proposed it:

1) alphat isnt set to zero for yPlus < yPlusLam, but I couldn t find a reference so far were the equation comes from
2) The BC calculates alphaEff and updates alphat by subtracting alpha (laminar)
3) For yPlus > yPlusLam the calculation of alphaEff differs from the version I proposed by an additional term in the denominator which contains velocities (which I couldn t find out yet where it comes from):

alphaEff = mu*(yPlus/(Pr_t*log(E *yPlus)/kappa_ + 0.5*rhow*uTau*(Pr_t*sqr(mag(Up)) + (Pr_t - Pr)*sqr(mag(Uc)))))

For yPlus < yPlusLam a similar term can be found. I was searching for a reference on that but couldn t find one so far. Any comments on that are appreciated!! I think I will give it a try in OF 1.7 anyway.

cheers,
Aram

[makaveli_lcf]

Hi,

those additional terms are for compressible flow, thus there are none of them in incompressible BC version. Check for the full version:

http://my.fit.edu/itresources/manual...pdf/flu_th.pdf, (page 124)

Sorry for the wrong link to the Theory Guide)

Cheers!

[mabinty]

That s great! Got it now!!

Thanks again!!

Cheers,
Aram

[vahid.najafi]

Hi every foamers.

I want to add k(kinetic turbulence energy) in one solver.

for this work, added next line in the code of solver:

const volScalarField &k=U_.db().lookupObject<volScalarField>("k")
and wmake was Successfully .

but How can I Understand that this k is the same with k(kinetic turbulence energy)???Is it true?????or not????

Because I replase another word in this :

const volScalarField &M=U_.db().lookupObject<volScalarField>("M")
but nothing error was not occured!!!!!!??????

[makaveli_lcf]

Could the lookup method just give you a null pointer, then no error is detected

[vahid.najafi]

Thanks Dr. Alexander A. Vakhrushev for your answer.

its means this work not true???this k is not kinetic turbulence energy???
if your answer yes!!!Please suggest me one way to Correction this code with k(kinetic turbulence energy)!!!

[makaveli_lcf]

Could you show more information, e.g.post your solver (if it is not secret) to understand your problem. Actually it should be ok what yo use to get k field.

[vahid.najafi]

Thanks for your answer Dr.
it is not secret!!!
my solver is in Appendix.I want to added kinetic turbulence energy in Directory:
OpenFOAM/vahid 2.0.1/run/myinterPhaseChangeFoam/phaseChangeTwoPhaseMixtures/SchnerrSauer/SchnerrSauer.c
in file:SchnerrSauer.c

I have add (turbulence kinetic Energy or k ) in a part of this solver.
first:
added:
-I$(LIB_SRC)/turbulenceModels \
-I$(LIB_SRC)/turbulenceModels/incompressible/RAS/RASModel \

in option.


and then:
added
//.................................................. .....changed
// Construct incompressible turbulence model
autoPtr<incompressible::RASModel> turbulence
(
incompressible::RASModel::New(U, phi, twoPhaseProperties())
);
//.................................................. .....changed


in creatField.

and then:
added

#include "RASModel.H"
and:

//.................................................. .............
volScalarField turbKinEnergy = turbulence().k(); //turbKinEnergy = turbulence().k();
volScalarField turbDisEnergy = turbulence().epsilon();
in this and turbkinEnergy=k
//.................................................. .............


in myinterPhaseChangeFoam.c

and wmake is ok!!!



but
I want in this solver in directory:
/run/myinterPhaseChangeFoam/phaseChangeTwoPhaseMixtures/SchnerrSauer/SchnerrSauer.c




I added
#include "RASModel.H"

in it and:


Foam::tmp<Foam::volScalarField>
Foam:haseChangeTwoPhaseMixtures::SchnerrSauer: Coeff
(
const volScalarField& p
) const
{
volScalarField limitedAlpha1(min(max(alpha1_, scalar(0)), scalar(1)));
volScalarField rho
(
limitedAlpha1*rho1() + (scalar(1) - limitedAlpha1)*rho2()
);

return
//......I want to change it( <<k>> turbulence multiple in it):
(3*rho1()*rho2())*sqrt(2/(3*rho1()))*k()
*rRb(limitedAlpha1)/(rho*sqrt(mag(p - pSat()) + 0.01*pSat()));
//.................................................. ......
}


dont successful wmake, and seen(was not declared ):

I/opt/openfoam201/src/OSspecific/POSIX/lnInclude -fPIC -c $SOURCE -o Make/linux64GccDPOpt/SchnerrSauer.o
phaseChangeTwoPhaseMixtures/SchnerrSauer/SchnerrSauer.C: In member function â€کFoam::tmp<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> > Foam:haseChangeTwoPhaseMixtures::SchnerrSauer: Coeff(const Foam::volScalarField&) const’:
phaseChangeTwoPhaseMixtures/SchnerrSauer/SchnerrSauer.C:113: error: k() was not declared in this scope
make: *** [Make/linux64GccDPOpt/SchnerrSauer.o] Error 1


please help me.

another way that I worked:


Foam::tmp<Foam::volScalarField>
Foam:haseChangeTwoPhaseMixtures::SchnerrSauer: Coeff
(
const volScalarField& p
) const
{
volScalarField limitedAlpha1(min(max(alpha1_, scalar(0)), scalar(1)));
volScalarField rho
(
limitedAlpha1*rho1() + (scalar(1) - limitedAlpha1)*rho2()
);

return
const volScalarField &k=U_.db().lookupObject<volScalarField>("k");

(3*rho1()*rho2())*sqrt(2/(3*rho1()))*k()
*rRb(limitedAlpha1)/(rho*sqrt(mag(p - pSat()) + 0.01*pSat()));


and wmake was Successfully .

but How can I Understand that this k is the same with k(kinetic turbulence energy)???Is it true???
Because I replase another world in this :
const volScalarField &M=U_.db().lookupObject<volScalarField>("M")
but not error occured!!!!!!??????


please help me more about this problem????Thanks.

[camoesas]

Hello Everybody,

I am exhuming this thread because it is the only one were somebody talks about yPlustherm. So maybe the Pro's see their old thread poping up and somebody can give me a reply:

I am working with the above mentioned wall function. What I dont understand is the calculation of yPlustherm. The docu says it is the value of yPlus at the edge of the thermal laminar sublayer. But where does the formula come from:

00084         scalar f = ypt - (log(E_*ypt)/kappa_ + P)/Prat; 
00085         scalar df = 1.0 - 1.0/(ypt*kappa_*Prat); 
00086         scalar yptNew = ypt - f/df;

Fluent uses a similar alphajayatillekewallfunction but says that yPlusTherm is the point where the linear and the log law intersect.
So my question is: What is yPlustherm and whats the background for this formula?
So I would be gratefull if somebody could lift the fog.


Thanks!

Camoesas