[tyskies]

Hello everyone,


I just implemented the EB model of Manceau, R., & Hanjalić, K. (2002) following the lines of Javadi, A. (2016). However, I'm having trouble to validate my implementation with the channel flow at Retau=395... I attached my results with the precursor simulation (left) and the simulation with EB model (right).

kOmegaSST_results.jpg EBRSM_results.jpg

Note that the u+-profile is good wtih the precursor simulation (not with the EB model) and that the flow accelerate in the centerline with the new turbulence model. This is confirmed when we look at the pressure gradient in the log of both simulations:
0.985 for the precursor simulation and 0.689 with the other.

Does anyone faced those kind of issue ?
Does anybody knows how to solve this ?

Thanks.


NB: I'm using OpenFOAM v2012.

[HPE]

Good to double check your implementation: discrepancy in fvc::grad(U)

Did the EBRSM converge? Did you run the simulation long enough?

[tyskies]

Originally Posted by HPE

Good to double check your implementation: discrepancy in fvc::grad(U)

Did the EBRSM converge? Did you run the simulation long enough?

Thanks for your answer HPE !

Since yesterday:
1) I've double checked my implementation as you suggested.
2) I've compared my implementation with SSG model because the part where gradU is used is the same as in EB model.
To me, all seems good! Still attached my model and the test case:

validation_ebrsm.zip

And yes, the results that I show are converged...

Do you have any other ideas ?
Thanks a lot for the help!

[HPE]

- uncomment "correctNut()" func.
- run the case with "simpleFoam" with "consistent=true"

May be the following is equivalent to what you have, but the new is less confusing to me
- volTensorField Omega(skew(gradU)); => volTensorField Omega(gradU.T() - gradU); // remember OpenFOAM's gradU is actually the tranpose of the usual Jacobian notation
- twoSymm(b & Omega) => twoSymm(b & Omega.T())

hope these help

[tyskies]

Originally Posted by HPE

- uncomment "correctNut()" func.
- run the case with "simpleFoam" with "consistent=true"

May be the following is equivalent to what you have, but the new is less confusing to me
- volTensorField Omega(skew(gradU)); => volTensorField Omega(gradU.T() - gradU); // remember OpenFOAM's gradU is actually the tranpose of the usual Jacobian notation
- twoSymm(b & Omega) => twoSymm(b & Omega.T())

hope these help

Thanks for your answer !


I've modified the code the way you suggested and re-run the case with pisoFoam and simpleFoam. Unfortunately (or not!), I get the same results with both solvers:

pisoFoam_EBRSM_results.jpg simpleFoam_EBRSM_results.jpg

EDIT: The horizontal axis of the top right corner figure is y+, not (x2/delta).


I think the problem is due to the way the solver adjusts the pressure gradient because it drops quickly at the beginning of the simulation...
Have you ever read a thread on those kind of issues with the meanVelocityForce option ?


Another time thanks for your help!

[HPE]

I am confident that "meanVelocityForce" is fine.
How many iterations did you run in simpleFoam case?

[tyskies]

Originally Posted by HPE

I am confident that "meanVelocityForce" is fine.
How many iterations did you run in simpleFoam case?

Yes, actually I agree with you on meanVelocityForce (thinking of it during the night!).

I ran 100 000 iterations.
It looks like the Reynolds stresses and the U+-profile converged at the 25 000th iteration.
Moreover, the pressure gradient/wall shear stress is 0.7185 at the end of the simulation whereas it should be around 0.98.

[tyskies]

Originally Posted by HPE

I am confident that "meanVelocityForce" is fine.
How many iterations did you run in simpleFoam case?

I've done some tests and when I remove the boundary condition on epsilon I get the same results. Everything is working as if the boundary condition on epsilon isn't applied.

As my programming skills in OpenFoam are limited, have you any idea how I can improve this part ?
Moreover, what does the following command epsEqn.ref().boundaryManipulate ?


Thanks for your help!

[tyskies]

Originally Posted by HPE

I am confident that "meanVelocityForce" is fine.
How many iterations did you run in simpleFoam case?

I've been able to obtain correct results with the following lines for the epsilon boundary condition:

    const fvPatchList& patches = this->mesh_.boundary();
    forAll(patches, patchi)
    {
        const fvPatch& curPatch = patches[patchi];
        if (isA<wallFvPatch>(curPatch))
        {
          //const scalarField& nuw = nu().boundaryField()[patchi];
          forAll(curPatch, facei)
          {
              label celli = curPatch.faceCells()[facei];
              epsilon_[celli] = 2.0 * 0.002532 * k_[celli] / 0.0003 / 0.0003 ;  // / (nearWallDist(this->mesh_).y()[facei] * nearWallDist(this->mesh_).y()[facei]); // this->y_[patchi][facei] / this->y_[patchi][facei];
          }
       }
    }

Results.jpg


Now, I need to properly implement the following boundary condition: epsilon = 2 * nu * k_ / y / y

However, when I use the following lines:

    const fvPatchList& patches = this->mesh_.boundary();
    forAll(patches, patchi)
    {
        const fvPatch& curPatch = patches[patchi];
        if (isA<wallFvPatch>(curPatch))
        {
          forAll(curPatch, facei)
          {
              label celli = curPatch.faceCells()[facei];
              // epsilon_[celli] = (1/0.0003/0.0003) * nuw[facei] * k_[celli] ;  // / (nearWallDist(this->mesh_).y()[facei] * nearWallDist(this->mesh_).y()[facei]); // this->y_[patchi][facei] / this->y_[patchi][facei];
              epsilon_[celli] = 2.0 * nu_[celli] * k_[celli] / curPatch[facei].y() / curPatch[facei].y();
           }
       }
    }

I get an error :

error: no match for ‘operator[]’ (operand types are ‘const Foam::fvPatch’ and ‘Foam::label’ {aka ‘int’})
  412 |               epsilon_[celli] = 2.0 * nu_[celli] * k_[celli] / curPatch[facei].y() / curPatch[facei].y();

What am I doing wrong ?
Thanks for your help!