[alex_rubel]

Hi
I'm new with OpenFoam and I done a bit of CFD with CFX before.
I doing a quite simple study about a submarine (without appendages) in OF 1.6.
I'm using simpleFoam but for the moment I get quite bad results, and convergence problems.
My inlet speed is 3.5 m/s and my mesh is around 6millions elements with a Y+ around 10. (this Mesh gives me very good results in CFX) The checkMesh is OK :

Mesh stats
    points:           6238133
    faces:            18602584
    internal faces:   18491528
    cells:            6182352
    boundary patches: 4
    point zones:      0
    face zones:       0
    cell zones:       0

Overall number of cells of each type:
    hexahedra:     6182352
    prisms:        0
    wedges:        0
    pyramids:      0
    tet wedges:    0
    tetrahedra:    0
    polyhedra:     0

Checking topology...
    Boundary definition OK.
    Point usage OK.
    Upper triangular ordering OK.
    Face vertices OK.
    Number of regions: 1 (OK).

Checking patch topology for multiply connected surfaces ...
    Patch               Faces    Points   Surface topology                  
    SURFACE_HULL        38368    38370    ok (closed singly connected)      
    SURFACE_SIDE        46288    46464    ok (non-closed singly connected)  
    SURFACE_OUTLET      17600    17689    ok (non-closed singly connected)  
    SURFACE_INLET       8800     8889     ok (non-closed singly connected)  

Checking geometry...
    Overall domain bounding box (-7.696 -8.712 -8.712) (21.78 8.712 8.712)
    Mesh (non-empty, non-wedge) directions (1 1 1)
    Mesh (non-empty) directions (1 1 1)
    Boundary openness (-5.04721e-17 2.96585e-16 3.75431e-17) OK.
    Max cell openness = 7.73145e-15 OK.
    Max aspect ratio = 496.544 OK.
    Minumum face area = 1.83581e-08. Maximum face area = 0.376178.  Face area magnitudes OK.
    Min volume = 6.63038e-12. Max volume = 0.102174.  Total volume = 6329.87.  Cell volumes OK.
    Mesh non-orthogonality Max: 61.2197 average: 9.26162
    Non-orthogonality check OK.
    Face pyramids OK.
    Max skewness = 2.01934 OK.

Mesh OK.

I'm using kOmegaSST model.
I use potentialFoam then lunch the simpleFoam without turbulences and finally set the turbulences on.

I set the following relaxation factors

relaxationFactors
{
    p               0.3;
    U               0.7;
    k               0.5;
    omega           0.5;
    R               0.7;
    nuTilda         0.7;
}

I'm also using kqRWallFunction, omegaWallFunction and nutWallFunction in the hull BC.


I got a lot of question with no answers :
1/ Am I using a good scheme (potential,turb off, turb on) to solve my problem ?
2/ Should I use WallFunction or fixedValue (0) for my hull BC (omega,k,nut) ?
3/ Should I use nNonOrthogonalCorrectors ?
4/ Is my relaxation factors good ?
5/ Do you have any idea to help me ?


Tanks a lot
Alex

[vkrastev]

Can you post your fvSchemes and fvSolution dictionaries?

V.

[maddalena]

Hi Alex,
a good way to enter in the OpenFOAM world of controlling any solver is this thread: http://www.cfd-online.com/Forums/ope...ime-steps.html. I found it really inspirating (also because it helped me to solve quite a lot of problems ). V collaborate to that thread as well.

mad

[alex_rubel]

thanks for theese answers. I'm not a work (It's 9pm in Australia) so i can't try right know.
I gonna post my files tomorow and try your thread.
Alex

[alex_rubel]

I posted my system file in the previous post

I red plenty of interesting threads this morning, but my understanding of english is not perfect so I really need your help for the fvScheme and fvSolution files.
For the BC I found information which say that because my Y+ is between 1<Y+<30 and I'm using High Re model I should set :
k-> zeroGradient
Omega -> WF
nut -> nutSpalartAllamarasWallFunction

What do you think about these settings ?

[vkrastev]

1) There's no need to use potentialFoam or start the simulation with turbulence switched off: if the numerical settings are appropriate (and mesh is not a disaster) the run should start normally with turbulence switched on.

2) The BC's at the wall are substantially correct (for k you can also use kqRWallFunction, but is the same as zeroGradient). You can try also nutWallFunction for nut (it's a simplified wall function compared with the nutSpalartAllmaras one)

3) In the fvSchemes file there's an error (divergence and laplacian schemes are for epsilon instead of omega): I dont' think this is the file you employed for your runs, as the simulation should not have started if the fvSchemes entries were not consistent with the turbulent quantities involved in your case. However, if you change epsilon to omega this file is also quite ok as a first attempt (it is not the best choice if you are searching for accuracy, but it should work also on not so good meshes)

4) In the fvSolution you have to lower a lot the tolerance (NOT the relTol) parameter for all the quantities: something like 10^-12 for p and 10^-10 for the other quantities should be ok. Also (but this is not so fundamental) you can add 1 or 2 nonOrthogonalCorrectors. Finally, the PCG solver for the pressure is ok, but as your mesh is quite big probably you could have a faster solution using GAMG instead (have a look in some of the tutorials or here in the forum to see how to set properly this kind of solver: of course, tolerance and relTol parameters should be the same as for PCG)

Good luck

V.

[alex_rubel]

Hi, i ran a simulation this week end with the following settings :
k->zeroGradient
Omega->WF
nut->SpalartAllmarasWF
I ran the problem with the turbulence and the high speed and I get a quite good convergence using the fvSchemes and fvSolution attached.
But I get 70% of error for my Cd,
You said that my files are good for the convergence but might be improved to get a better accuracy. Can you help me for that
Thanks
Alex

[alex_rubel]

i try the GAMG solver today and It works very well, after 1000 TStep I got the same Residuals "quality" and it reduce the ExecutionTime from 53674s (PCG) from 19778s (GAMG)
I red very interesting post about the nCellsInCoarsestLevel parameter and if I understand evrything I should set it up to sqrt(MeshNbrofElement) ? I curently using 10.
If you can give me advice of the relTol parameter it would help me a lot.
I may ask you too much but I hope you can help me
Alex

[vkrastev]

Originally Posted by alex_rubel

Hi, i ran a simulation this week end with the following settings :
k->zeroGradient
Omega->WF
nut->SpalartAllmarasWF
I ran the problem with the turbulence and the high speed and I get a quite good convergence using the fvSchemes and fvSolution attached.
But I get 70% of error for my Cd,
You said that my files are good for the convergence but might be improved to get a better accuracy. Can you help me for that
Thanks
Alex

OK, once you have stability let's try to search for accuracy!

1) In fvSchemes, try first to set div(phi,U) to

Gauss linearUpwindV cellMDLimited Gauss linear 1;

2) If this change causes stability issues, put the limiter also on grad(U) (NOT on grad(p)), by setting:

grad(U) cellMDLimited Gauss linear 1;

3) If you still have stability troubles (but I don't think this will be the case), put a limiter on the laplacianSchemes too, by setting:

default Gauss linear limited 0.5;

4) If you will be successful in using the linearUpwind scheme for div(phi,U) (either with or without the additional limiting in points 2 and 3), you can try to extend it also on the convection terms for k and omega (but to my knowledge and experience the biggest improvement in accuracy tend to come from the momentum convection term dicretization choice), by setting:

div(phi,k) Gauss linearUpwind cellMDLimited Gauss linear 1;

and the same for omega.

Hope this helps!

V.

[vkrastev]

Originally Posted by alex_rubel

i try the GAMG solver today and It works very well, after 1000 TStep I got the same Residuals "quality" and it reduce the ExecutionTime from 53674s (PCG) from 19778s (GAMG)
I red very interesting post about the nCellsInCoarsestLevel parameter and if I understand evrything I should set it up to sqrt(MeshNbrofElement) ? I curently using 10.
If you can give me advice of the relTol parameter it would help me a lot.
I may ask you too much but I hope you can help me
Alex

Actually It is not true that there is an "ideal choice" for the nCellsInCoarsestLevel parameter, as the best choice (in terms of solver's efficiency trade off between number of iterations necessary to reach a given convergence level and mean time required for each iteration to be completed) strongly depend on the dimension of your problem and on the type of Multigrid solver in question (you can read some very exhaustive explanations about it here http://www.cfd-online.com/Forums/ope...e-steps-2.html , see in particular pages 2 and 3 of the tread an search for the posts coming from arjun). Apart from this, nCellsInCoarsestLevel from 10 to 50 should be ok for your case.
About the relTol parameters, It's my opinion that using too low values in a steady-state and underrelaxed solution practice (as is the SIMPLE algorithm in the simpleFoam solver) is not useful at all and produces only a big waste of time for each iteration...However, a minimum convergence criterion should be provided within the single iteration, thus I think that your values are quite ok (0.05 for p and 0.1 for the other quantities). If you want to make some tests, try to lower them to 0.001 for p and 0.01 for the other quantities, but I don't believe this will produce significant improvements in the convergence pattern of your residuals.

Hope this helps, too!

V.

[alex_rubel]

Thanks, I gonna try with that.
Do you have information about the nCellsInCoarsestLevel parameter ? I red plenty of different explanation about this. Apparently it's set how many cellsare going to be solve by the same equation but I'm not sure at all
Alex

Edit
Tanks again for your explanation it's a lot more clear in my head

[alex_rubel]

Hi,
I get some good improvement for my simulation with a Y+ of 30 (KomegaSST with wallFunctions) using

divSchemes
{
    default         none;
    div(phi,U)      Gauss linearUpwindV Gauss linear;
    div(phi,k)      Gauss linearUpwind Gauss linear;
    div(phi,omega)  Gauss linearUpwind Gauss linear;
    div((nuEff*dev(grad(U).T()))) Gauss linear;

I get Cd = 0.00142 instead of 0.00175 before. The Cdexp = 0.001
I tried your settings but it leads to bounding problem for k and omega and an equivalent result.

1) In fvSchemes, try first to set div(phi,U) to

Gauss linearUpwindV cellMDLimited Gauss linear 1;

2) If this change causes stability issues, put the limiter also on grad(U) (NOT on grad(p)), by setting:

grad(U) cellMDLimited Gauss linear 1;

3) If you still have stability troubles (but I don't think this will be the case), put a limiter on the laplacianSchemes too, by setting:

default Gauss linear limited 0.5;

4) If you will be successful in using the linearUpwind scheme for div(phi,U) (either with or without the additional limiting in points 2 and 3), you can try to extend it also on the convection terms for k and omega (but to my knowledge and experience the biggest improvement in accuracy tend to come from the momentum convection term dicretization choice), by setting:

div(phi,k) Gauss linearUpwind cellMDLimited Gauss linear 1;

Even if my Cderror is still quite high I really like to find imporvement for a small Y+ (1 to 10) simulation and without wallFunctions and this schemes


g

radSchemes
{
    default         Gauss linear;
    grad(p)         Gauss linear;
    grad(U)         Gauss linear;
// upwind 1st order , linear = second
}

divSchemes
{
    default         none;
    div(phi,U)      Gauss upwind;
    div(phi,k)      Gauss upwind;
    div(phi,omega) Gauss upwind;
    div(phi,epsilon) Gauss upwind;    
    div((nuEff*dev(grad(U).T()))) Gauss linear;
}

laplacianSchemes
{
    default         Gauss linear corrected;

}

I get convergence issues
For the WallFunctions I set zeroGradient for k, omegaWallFunction and nutSpalartAllmarasWallFunctions do I have to to set a very small fixedValue instead of zeroGradient, and what about nut ?
I tried with zerGradient instead of omegaWallFunction but I didn't see any difference

Thanks Alex

ps : I used Gauss linearUpwind Gauss linear; instead of Gauss linearUpwind cellMDLimited Gauss linear 1; I know that cellMDLimited is a modifier for the gradient scheme Gauss linear but that's all is you have more information it will be great.
pps : I hop I 'll be able to help you back one day

[vkrastev]

Originally Posted by alex_rubel

Hi,
I get some good improvement for my simulation with a Y+ of 30 (KomegaSST with wallFunctions) using

divSchemes
{
    default         none;
    div(phi,U)      Gauss linearUpwindV Gauss linear;
    div(phi,k)      Gauss linearUpwind Gauss linear;
    div(phi,omega)  Gauss linearUpwind Gauss linear;
    div((nuEff*dev(grad(U).T()))) Gauss linear;

I get Cd = 0.00142 instead of 0.00175 before. The Cdexp = 0.001
I tried your settings but it leads to bounding problem for k and omega and an equivalent result.
Even if my Cderror is still quite high I really like to find imporvement for a small Y+ (1 to 10) simulation and without wallFunctions and this schemes


g

radSchemes
{
    default         Gauss linear;
    grad(p)         Gauss linear;
    grad(U)         Gauss linear;
// upwind 1st order , linear = second
}

divSchemes
{
    default         none;
    div(phi,U)      Gauss upwind;
    div(phi,k)      Gauss upwind;
    div(phi,omega) Gauss upwind;
    div(phi,epsilon) Gauss upwind;    
    div((nuEff*dev(grad(U).T()))) Gauss linear;
}

laplacianSchemes
{
    default         Gauss linear corrected;

}

I get convergence issues
For the WallFunctions I set zeroGradient for k, omegaWallFunction and nutSpalartAllmarasWallFunctions do I have to to set a very small fixedValue instead of zeroGradient, and what about nut ?
I tried with zerGradient instead of omegaWallFunction but I didn't see any difference

Thanks Alex

ps : I used Gauss linearUpwind Gauss linear; instead of Gauss linearUpwind cellMDLimited Gauss linear 1; I know that cellMDLimited is a modifier for the gradient scheme Gauss linear but that's all is you have more information it will be great.
pps : I hop I 'll be able to help you back one day

I'm happy you have reached some improvements in your results, but let me point out some additional questions:

1) From your last post it is not totally clear what you want to do now to further improve your computed Cd value...However, if you want to pass from wall functions modeling to a direct solution of the boundary layer, you cannot definitely leave the first order upwind scheme on the convection terms, as accuracy is a milestone for the b.l. direct solving approach.

2) If you had a better behavior without the gradient limiters turned on, this is probably due to your hexa-mesh (hexa-meshes exhibits usually a more stable behavior then the tetra-dominant one in conjunction with higher-order interpolation schemes)

3) To directly solve the b. l. with the SST model you'll need to:

-Have the first nodes away from the wall inside the y+<1 limit
-Set k at the wall to fixedValue 10^-10 (or, in general, to a very low value)
-Set omega at the wall to fixedValue <a very high value> (it should be something like 10^6, but there is a formula usually used in the literature which depends on 1/y^2, where y is the the absolute distance of the first nodes: I think you can find this formula somewhere here in the forum or at least in some paper related to the SST model).
-Set nut at the wall to calculated (an initial "value uniform 0" shoul be ok)

Of course, lowering the y+ will probably lead to near wall cells with very high aspect ratios, but you should try to remain below the value of about 500 which you had in our y+=10 case (thus you have to reduce the cell's dimension also in the streamwise direction). In any case, as a first attempt you can try the same numerical setup which gave to you the best results in the WF configuration.

About the cellMDLimited option, it is a muldidimensional (MD) gradient limiter (it means that the limiter is imposed separately in the three dimensions), but honestly still I don't know what exactly are the limiting criteria adopted (as in OF there are three other kind of limiters: cellLimited, faceMDLimited, and faceLimited).

Good luck with your new runs!

V.

[alex_rubel]

Hi V,

1) From your last post it is not totally clear what you want to do now to further improve your computed Cd value...However, if you want to pass from wall functions modeling to a direct solution of the boundary layer, you cannot definitely leave the first order upwind scheme on the convection terms, as accuracy is a milestone for the b.l. direct solving approach.

I recognise that my last post is not very clear, I prefer spending time in the low Y+ model than the other one even if my accuracy is not perfect. So what you say is that I have to use upwind scheme for low Y+ model ?

2) If you had a better behavior without the gradient limiters turned on, this is probably due to your hexa-mesh (hexa-meshes exhibits usually a more stable behavior then the tetra-dominant one in conjunction with higher-order interpolation schemes)

It makes sense, sometime you are happy of the long time you spent on doing a structured mesh

And thanks a lot for the advices of your third point I gonna try with that and I'll post my results asap.

I found something about a missing scalar of sqrt(2) in the KomegaSST model, I modified it (If I change the code, save it and make ./Allwmake it's ok ?), and for the moment my Cd is heading to 0.00125 instead of 0.0014 before

Alex

[vkrastev]

Originally Posted by alex_rubel

I recognise that my last post is not very clear, I prefer spending time in the low Y+ model than the other one even if my accuracy is not perfect. So what you say is that I have to use upwind scheme for low Y+ model ?

Just the opposite: in any case where you need accuracy (and the direct solving of the b. l. is such a case) you cannot use the first order upwind scheme, as it will introduce a high amount of artificial diffusivity (in simple words, the error in the solution will be significant). So you have to switch (at least for the momentum convection term div(phi,U)) to something more accurate, such as the linearUpwind scheme (but there are also other possible choices, such as the limitedLinear or the Gamma schemes).

Originally Posted by alex_rubel

I found something about a missing scalar of sqrt(2) in the KomegaSST model, I modified it (If I change the code, save it and make ./Allwmake it's ok ?), and for the moment my Cd is heading to 0.00125 instead of 0.0014 before

Yes, I've read something about this bug in the SST model implementation in OpenFOAM. As your modification to the code has been minimal, running only ./Allwmake is ok. In other cases, when the modifications are more invasive, it is safer to run first wclean all (inside the folder containing the new source files and the related Make folder) and then wmake all or the Allwmake executable (if present in the top level folder). By the way, which schemes have you used to obtain Cd=0.00125?

V.

[alex_rubel]

Hi V.

By the way, which schemes have you used to obtain Cd=0.00125?

I used :

ddtSchemes
{
    default steadyState;
}

gradSchemes
{
    default         Gauss linear;
    grad(p)         Gauss linear;
    grad(U)         Gauss linear;
//    grad(U)         cellLimited Gauss linear 1;
}

divSchemes
{
    default         none;
    div(phi,U)      Gauss linearUpwindV Gauss linear;
    div(phi,k)      Gauss linearUpwind Gauss linear;
    div(phi,omega)  Gauss linearUpwind Gauss linear;
    div((nuEff*dev(grad(U).T()))) Gauss linear;
}

laplacianSchemes
{
//    default         Gauss linear corrected;
//    default         Gauss linear limited 0.5;
    default         Gauss linear limited 0.333;
}

This settings seem to work fine for my case but if you see any improvement to do please let me know.

I ran a simulation for my low Y+ model with the following settings

nut :

SURFACE_HULL
    {
        type            calculated;
        value           uniform 0;
    }

omega :

SURFACE_HULL
    {
        type            fixedValue;
        value        uniform 10^6; //To modify
    }

k :

SURFACE_HULL
    {
        type            fixedValue;
    value         uniform 10^-10;
    }

and the following scheme :

gradSchemes
{
    default         Gauss linear;
    grad(p)         Gauss linear;
//    grad(U)         Gauss linear;
    grad(U)         cellMDLimited Gauss linear 1;
}

divSchemes
{
    default         none;
    div(phi,U)      Gauss linearUpwindV cellMDLimited Gauss linear 1;
    div(phi,k)      Gauss linearUpwind cellMDLimited Gauss linear 1;
    div(phi,omega)  Gauss linearUpwind cellMDLimited Gauss linear 1;
    div((nuEff*dev(grad(U).T()))) Gauss linear;
}

laplacianSchemes
{
//    default         Gauss linear corrected;
//    default         Gauss linear limited 0.5;
    default         Gauss linear limited 0.333;
}

I tough that the cellMDLimited modifier could help the convergence of my system but as you can see it's bad , very bad.
I used the following relaxation factors :

relaxationFactors
{
    p               0.1;
    U               0.5;
    k               0.3;
    omega           0.3;
}

I gonna try different possibility today but if you see any way to help me, it would be great. I attached the whole case if you want to have a look at it.

Alex

[vkrastev]

Originally Posted by alex_rubel

I ran a simulation for my low Y+ model with the following settings

nut :

SURFACE_HULL
    {
        type            calculated;
        value           uniform 0;
    }

omega :

SURFACE_HULL
    {
        type            fixedValue;
        value        uniform 10^6; //To modify
    }

k :

SURFACE_HULL
    {
        type            fixedValue;
    value         uniform 10^-10;
    }

and the following scheme :

gradSchemes
{
    default         Gauss linear;
    grad(p)         Gauss linear;
//    grad(U)         Gauss linear;
    grad(U)         cellMDLimited Gauss linear 1;
}

divSchemes
{
    default         none;
    div(phi,U)      Gauss linearUpwindV cellMDLimited Gauss linear 1;
    div(phi,k)      Gauss linearUpwind cellMDLimited Gauss linear 1;
    div(phi,omega)  Gauss linearUpwind cellMDLimited Gauss linear 1;
    div((nuEff*dev(grad(U).T()))) Gauss linear;
}

laplacianSchemes
{
//    default         Gauss linear corrected;
//    default         Gauss linear limited 0.5;
    default         Gauss linear limited 0.333;
}

I tough that the cellMDLimited modifier could help the convergence of my system but as you can see it's bad , very bad.
I used the following relaxation factors :

relaxationFactors
{
    p               0.1;
    U               0.5;
    k               0.3;
    omega           0.3;
}

I gonna try different possibility today but if you see any way to help me, it would be great. I attached the whole case if you want to have a look at it.

Alex

Before starting a big struggle with the fvSchemes file, can you please post the checkMesh response for your low-y+ mesh? In addition, for the best setting of omega at the wall, take a look here http://www.cfd-online.com/Forums/ope...oam-1-7-a.html

About your settings for the high-y+ mesh, I think they are pretty good (probably you will not be able to obtain better results in any case, as a consequence of the simplified wall modeling).

Regards

V.

[alex_rubel]

Hi V

About your settings for the high-y+ mesh, I think they are pretty good (probably you will not be able to obtain better results in any case, as a consequence of the simplified wall modeling).

That's what I though

There is a checkMesh in the simu7 zip file for the low Y+ model

[vkrastev]

Originally Posted by alex_rubel

There is a checkMesh in the simu7 zip file for the low Y+ model

***High aspect ratio cells found, Max aspect ratio: 2828.2, number of cells 34328
<<Writing 34328 cells with high aspect ratio to set highAspectRatioCells

Here it is the problem....I told you in one of my previous posts that you should keep a not extremely high aspect ratio for the near-wall cells, although this would constrain you to reduce the dimensions of the cells also in the streamwise and spanwise directions. At most you must keep the same max value of your high y+ case (which was around 500), but probably you'll need to reduce it further (I know that the number of cells will increase a lot, but this is a price you have to pay for the direct solution of the b. l.). So, reduce your max axpect ratio (at least to 500), set the correct b. c. for omega at the wall (see the reference in my previous post), keep all the other settings as before, and let's see what happens...

Regards

V.

[alex_rubel]

I told you in one of my previous posts that you should keep a not extremely high aspect ratio for the near-wall cells, although this would constrain you to reduce the dimensions of the cells also in the streamwise and spanwise directions.

You are totally right, I forget that, I gonna modify my mesh.

I'm doing a yaw angle study for my high Y+ model and the result are quite consistent I don't think I can improve a lot this model. I'm gonna try Kepsilon And Komega but my experiment on my case with CFX shows me that KomegaSST is the more accurate

Alex