[jp3g12]

Greetings Foamers!

I am simulating a symmetric airfoil with a Reynolds number of 6,000,000 million and a zero angle of attack in OpenFOAM. The flow velocity is 87 m/s with standard density and kinematic viscosity (also the chord of the airfoil is 1 m). I am using the realizable k-epsilon turbulence model but I am struggling a bit with it. I have simulated this same airfoil at the same conditions using Spalart-Allmaras and SST and the results I got were close to experimental/computational results. The mesh worked with SA and SST although it has very high aspect ratio cells on the wake of the airfoil.

Just as an example of how 'crappy' the results are, I have attached the residuals for 1st order upwind. k and epsilon seems to oscillate a lot. I am not using wall functions because my mesh is very fine near the wall (approximately achieves about a wall y+ of around 0.2). Instead I am imposing a zeroGradient BC on the wall for k and epsilon. The boundary conditions for k, epsilon and nut were computed using the following relations from Spalart and Rumsey:

nut = 2e-7 *U
k = 1e-6 *U^2
epsilon = 4.5e-7*U^3

I don't know exactly what is going wrong with this case. It may be the fvSchemes I am using. In general, just to be safe, I run calculations for a few iterations with

    div(phi,U)          Gauss upwind;
    div(phi,epsilon)    Gauss upwind;
    div(phi,k)          Gauss upwind;

to then move to second order schemes

 div(phi,U)      Gauss linearUpwind grad(U);  
 div(phi,epsilon) Gauss linearUpwind grad(epsilon);
 div(phi,k)      Gauss linearUpwind grad(k);

But as it can be seen, I don't even get nice residuals with 1st order, and when I try 2nd order it just blows up!!!

Could anyone look at my case and point out possible mistakes or things a more experience user would do?

https://www.dropbox.com/sh/krabswgrufc2cm4/bP7FCRefq4

Thank you very much for the help!

James

[linnemann]

Hi

First off if these the correct dimensions of the mesh

Overall domain bounding box (-484.457 -1 -507.806) (501 0 507.806)

It seems really large, 1km long/high.

[jp3g12]

Originally Posted by linnemann

Hi

First off if these the correct dimensions of the mesh

Overall domain bounding box (-484.457 -1 -507.806) (501 0 507.806)

It seems really large, 1km long/high.

Yes that is the domain. It is a NASA validation case, and the domain lenght is about 1 km long 9they did that to avoid any interference with the domain boundaries). That should not be a problem however because I have been simulating this case succesfully for the SA and SST models. Furthermore, I used the same settings I posted for the airFoil2D tutorial and I got the same 'crappy' results.

As I said, there has to be something wrong with the numerical schemes or with the boundary conditions, but as far as I am concerned, the BCs I am using are pretty standard. I am aware that k-epsilon turbulence model is very tricky when it comes to the epsilon equation...

Can you see anything fundamentally wrong with my set up?

Thank you for your help!

James

[FelixL]

Hello James,


The realizable k-e-model in OpenFOAM is not a low-Re turbulence model. Therefore, you can not use it without wall functions. This is the reason why your results are bad.

Use the low-Re k-e-models implemented in OpenFOAM if you want to try out this class of turbulence models.


Greetings Felix

[jp3g12]

Hi Felix,

I will give it a shot. I will use the Launder Sharma turbulence model. My only question now, is why can't I use the realizable KE turbulence model without using wall functions even if I am resolving the entire BL?

Thanks,

James

[FelixL]

Hello James,


epsilon based turbulence models become unstable and yield unphysical results when integrated down to the wall through the viscous sublayer. The only way to overcome this is either to use wall functions (-> High-Re models) or to add damping functions into the turbulence model (-> Low-Re models).

For more details, please read Turbulence Modeling for CFD by D.C. Wilcox.


Greetings
Felix

[jp3g12]

Hello Felix,

So, I have tried to simulate the case using the LaunderSharma turbulence model. The residuals for the first order case look pretty good as it can be seen. I would like to point out that I had to change the fvSolution from what I posted earlier to:

 
   p
    {
        solver          PCG;
        preconditioner  DIC;
        tolerance       1e-16;
        relTol          0.001;
    }
    U
    {
        solver          PBiCG;
        preconditioner  DILU;
        tolerance       1e-15;
        relTol          0.001;
    }
    k
    {
        solver          PBiCG;
        preconditioner  DILU;
        tolerance       1e-8;
        relTol          0.001;
    }
    epsilon
    {
        solver          PBiCG;
        preconditioner  DILU;
        tolerance       1e-8;
        relTol          0.001;
    }

I couldn't make the other one work even for 1st order. Now when I switch from Gauss upwind in divSchemes for U, k and epsilon to Gauss linearUpwind for these quantities, the solution does not seem to converge anymore. The drag values and lift are consequently very far from what could be expected for a NACA 0012 at a zero angle of attack. I'm still using the same BCs that I was using before (so zeroGradient for both epsilon and k at the wall). Hence my case set-up in the zero folder is still the same that I posted earlier and the only thing that I have changed is the fvSolution file to what I have posted over here and the turbulence model to a low Reynolds number one.

Any suggestions?

Thanks for the help!!

James

[FelixL]

Hello James,


your boundary conditions are wrong. At the wall k and epsilon both equal zero, when using Low-Re models. See, e.g., http://www.cfd-online.com/Wiki/Low-Re_k-epsilon_models

Usually its better to set those quantities to low nonzero values, for example 1e-16. This avoids div by zero errors.


Greetings
Felix

P.s.: this information is also found in the book I mentioned. So I strongly advise you to read it, if you want to work with the different turbulence models.

[jp3g12]

Thanks Felix,

I somewhat got confused with other post that recommended using zero gradient. Also I tried setting there k and epsilon to zero but it would give me a floating point exception (due to the singularity in epsilon equation)... I will try to take a look at Wilcox but it is not a easy book to find !

Will try these settings!

Thanks for the help,

James

[jp3g12]

Hi Felix,

yet another question about this !!!

I implemented all the changes you suggested and now k and epsilon are zero in the wall (k, epsilon = 1e-16). I presume the rest of the boundary conditions are correct since they have remained untouched from the SST simulations I was able to complete succesfully a while back.

The changes however do not improve the convergence at all, and even for Gauss Upwind schemes I see no signs of convergence (the pressure residuals do not drop at all). I simulated the airFoil2D geometry with the same settings and it appeared to work very well.

Now I have no idea of what could be wrong with the set up!

Thanks,

James

[FelixL]

Hello James,


it's probably the tolerance settings for k and epsilon. As you can see in the residual plot, these quantities "converge" very quickly and OpenFOAM doesn't solve those equations anymore after a few iterations. However, residuals of order 1e-8 don't necessarily mean that the k and epsilon fields are correct, since they depent on the flow field which still is being solved.

Reduce the tolerance for k and epsilon to something like 1e-16 and see if it helps.


Greetings
Felix

[jp3g12]

Hi Felix,

As far as I am aware, I have the tolerance of every quantity set to 1e-16. Furthermore for k and epsilon I have a minIter of 1 (I don't know if it works though but I assume it does because I don't get any errors). Could it be the reltol?

Thank you,

James

[FelixL]

Hello James,


post #7 tells me something different about your tolerance settings? I really mean "tolerance", not "relTol". Also the residual plot confirms the 1e-8-setting for k and Epsilon.

minIter doesn't work in the official OpenFOAM distributions. It has no effect.


Greetings
Felix

[jp3g12]

Hi Felix,

Well I hate to say this but I seem to be (still) struggling with this. So I implemented all of your suggestions and the solution does not converge at all even using 1st order upwind.

I have revised the BCs over and over and, as long as I am aware, they should be correct now. I have also tried different numerical schemes and relaxation factors and they did not improve the situation. Furthermore in some cases I even get negative drag values, which bothers me a lot. I have re-checked the mesh and it is the same mesh I used for SA and SST (which worked just fine)...

I am not able to see anything really wrong with my set up. I attached it, in case someone wants to take a look at it. Recommendations are welcome!

Thank you very much for your help!

James

[linnemann]

Hi

Been trying to get your case to run with my usual tricks.

Nothing so far works and I always trace the problem to the same Cells in the mesh.

I have no solution to this problem.

[jp3g12]

Hi Linnemann,

yeah I have been working on it a lot and I don't seem to be able to find a solution either. It is very weird because that same mesh worked on SA and SST. Obviously the problem may be due to the fact that K-Epsilon is more sensitive to BCs near the wall than the other two turbulence models.

Thanks for the help

James

[linnemann]

Hi again

decided to try one more time and being really radical with the relaxation factors so I started out with.

        U               0.3;
        k               0.01;
        epsilon         0.01;

Then after 100 iterations

        U               0.3;
        k               0.1;
        epsilon         0.1;

After 200 iterations

        U               0.4;
        k               0.1;
        epsilon         0.1;

After 300 iterations

        U               0.4;
        k               0.2;
        epsilon         0.2;

After 400 iterations

        U               0.4;
        k               0.3;
        epsilon         0.3;

After 500 iterations

        U               0.5;
        k               0.3;
        epsilon         0.3;

After 600 iterations

        U               0.5;
        k               0.4;
        epsilon         0.4;

So far it has been running up til 2000itt.

The residuals are really not stable and the forces are fluctuating quite a bit.

But it runs and it does not crash and the velocity looks reasonable