[Tibo99]

Hi foamers!

I try to implement a Cmu constant which varies in term of 'z' in the domain for the k-epsilon turbulence model. So, in other words, I try to implement a local Cmu in the domain.

I looked for some code example on the web and here on CFD-Online, but had no luck.

Is there someone that could point me out toward something in order to be able to achieve that?

Thank for your help!

Best Regards,

[geth03]

i think that one is possible.

you would need to create a new turbulence model which inherits from the
k-epsilon turbulence model.

then you can overwrite/re-implement the correct()-function, bc here k and epsilon equations are solved. you can create a new Cmu, which is calculated based on the 'z' in your domain. adjust the k-epsilon equation afterwards.

keep in mind that this modification demands advanced programming skills.

i don't know any tutorial or whatsoever that shows how to do that.

[Tibo99]

Thank you very much for replying!

I think your're right about how complex this could be.

But, I end up with the following attempt.

I obviously copy the whole structure for the k-e model and name it kEpsilonABL.

And then, I saw in the literature that Cmu(z) = u*^4/k(z)^2 (see picture in attachment).

So, What I've done is, I switched the variable 'Cmu_' in the code for 'uStar_' and then assign the unit with 'dimensionSet' in order to be consistent with the unit, like so:

 uStar_
    (
        dimensioned<scalar>::lookupOrAddToDict
        (
            "uStar",
            this->coeffDict_,
            dimensionSet(0, 1, -1, 0, 0, 0, 0),
            0.511
        )
    ),

And then, I added to following line

const volScalarField Cmu_ = pow(uStar_,4)/pow(k_,2);

here,

void kEpsilonABL<BasicTurbulenceModel>::correctNut()
{
    const volScalarField Cmu_ = pow(uStar_,4)/pow(k_,2);
    this->nut_ = Cmu_*(sqr(k_)/epsilon_);
    this->nut_.correctBoundaryConditions();
    fv::options::New(this->mesh_).correct(this->nut_);

    BasicTurbulenceModel::correctNut();
}

Finnaly, I compile everything.

Everything compute and converge. The shape of the profile for the tubulent kinetic energy 'k' is good, but don't have the right results. The values are about 3x greater of what the analytic results give. Obviously, since the value for 'k' are not good, the values on the velocity profile 'Ux' and 'epsilon' are not there, but the shape of the profiles are good though.

Any suggestions? Do you think this could or should work?

Best Regards!

[geth03]

i think your approach is ok.

i assume, from your 1-D dependancy of Cmu_, that your flow problem is also 1-D. is that correct?

[Tibo99]

Hi!

Yes, you're right.

I started with this approach to make the BC work with the model and matching the analytics solution.

By confirming this, I'll step in with a 2D simulation and I guess it should work.

But, like I said, the result are no good in 1D. It doesn't match with the analytics results by using the code's line I just showed in my last post.

I can't find where is the issue by doing this that way.....

Thank for replying though!

Regards,

[mAlletto]

Did you check that the cmu used in the boundary condition for nut is the same as the one computed in your modification

[Tibo99]

Good question!

Since I didn't get any error when I launched the simulation, I was assuming that the 'Cmu' use for 'nut' in the boundary conditions was the same (global variable).

I modify it in the espilon boundary condition, but not in nut. Don't know why!!!!!!! No rational answer!!! lol

Since 'k' is use in 'nut', I'll add a constant 'ustar' and then compute 'Cmu' as I did in the turbulence model but in the loop.

Let see if this will work........

Thank again!!!

[Tibo99]

Here is the results of the simulation. Obviously, it didn't work. It look like something happen at the bottom of the 'k' profile but, other than that, the shape of every profile seems to follow the right law.

See file in attachment for the BC applied. A source term has been apply on epsilon and k (see Parente: https://www.researchgate.net/publica...n_of_ABL_flows).

I added the following code line in the nut file:

const scalar Cmu25 = pow025(pow(ustar_,4)/pow(k[celli],2));

in the loop here:

forAll(nutw, facei)
    {
        label celli = patch().faceCells()[facei];
		
	const scalar Cmu25 = pow025(pow(ustar_,4)/pow(k[celli],2));
		
        scalar uStar = Cmu25*sqrt(k[celli]);
        scalar yPlus = uStar*y[facei]/nuw[facei];

        scalar Edash = (y[facei] + z0_[facei])/z0_[facei];

        nutw[facei] =
            nuw[facei]*(yPlus*kappa_/log(max(Edash, 1+1e-4)) - 1);

        if (debug)
        {
            Info<< "yPlus = " << yPlus
                << ", Edash = " << Edash
                << ", nutw = " << nutw[facei]
                << endl;
        }
    }

Any suggestions? Any thought?

Thank you very much for your tips!

Regards,

[mAlletto]

From a first look you use two ustar. One you priscibe and one computed from k. This is in contradiction. Than k ist too high. Did ou check the model constants in the epsilon equation

[Tibo99]

Make sense about the 'ustar' so, since I didn't had any error message before, I came back to the original state for the 'nut' file with 'Cmu'.

Really weird that I end up with the same results though.

Talking about the constant in epsilon equation, I'll try to explain quickly why I use these constant. First, my goal is to replicate the results of what Yang(2009) got in his paper (https://www.sciencedirect.com/scienc...67610508001815) by using his constant with the method developed by Parente (Cmu varies) (https://www.researchgate.net/publica...n_of_ABL_flows)

The constant are (Yang,2009):
z0 = 0.000225m
u*= 0.511 m/s
C1e=1.5
C2e=2.51
Cmu=0.028 (will varies in our case)
sigmaK=1.67
sigmaEpsilon=2.51
C1=-0.17 (k profile constant)
C2=1.62 (k profile constant)

Since Parente use a different version of 'k' profile as Yang use, the 'k' profile constant become : C1=-0.414, C2=3.945. Parente doesn't show in his paper the constant of the transport equation. That why I used Yang data.

Parente use a source term on epsilon with the constant C1 and C2(see attachment). The notation is the same as used in the 'k' profile, which is confusing because, when I developed the source term on epsilon myself, these constant are Ce1 and Ce2 instead (see picture). I make Cmu varies in this source term as well.

However, when I use the epsilon constant in the source term, it unstable. It doesn't converge. But if I use the 'k' profile constant in the source term, it give the result I showed in my post #8.

About the BC for epsilon, I canceled these line:

//const scalar Cmu25 = pow025(Cmu_);
//const scalar Cmu75 = pow(Cmu_, 0.75);

And added these line in the loop (see file in attachment):

const scalar Cmu25 = pow025(pow(ustar_,4)/pow(k[celli],2));
const scalar Cmu75 = pow((pow(ustar_,4)/pow(k[celli],2)), 0.75);

Again, thank you very much for taking the time to read my post and interact with me.

Best regards and stay safe,

[cunconbkhp]

Hi Tibo99!

I'm trying to do exacly the same thing with Openfoam version 2212. I'm following the work of Salazar. Hi did some modification from Parent's work. I'm looking for a way to add the source term Se(z) like in your picture. Howerver i haven't found the good way to define z in the code and include it in tmp<fvScalarMatrix> epsEqn. That's great if you could give me some hints. Thanks in advance.