[Dorian1504]

Hello,

I would like to modify the nacaAirfoil tutorial case for sonicFoam in a steady state case. Is it possible to directly calculate the steady state without computing the transient regime with sonicFoam?
I naively tried to change the ddtSchemes from Euler to steadyState, as I can see on rhoSimpleFoam. However, sonicFoam return this error :

Courant Number mean: 0.000174184 max: 0.948172
#0 Foam::error:rintStack(Foam::Ostream&) at ??:?
#1 Foam::sigFpe::sigHandler(int) at ??:?
#2 ? in "/lib/x86_64-linux-gnu/libc.so.6"
#3 Foam::divide(Foam::Field<double>&, Foam::UList<double> const&, Foam::UList<double> const&) at ??:?
#4 Foam:perator/(Foam::UList<double> const&, Foam::UList<double> const&) at ??:?
#5 Foam::diagonalSolver::solve(Foam::Field<double>&, Foam::Field<double> const&, unsigned char) const at ??:?
#6 Foam::fvMatrix<double>::solveSegregated(Foam::dict ionary const&) at ??:?
#7 Foam::fvMatrix<double>::solve(Foam::dictionary const&) at ??:?
#8 Foam::fvMatrix<double>::solve() at ??:?
#9 ? at ??:?
#10 __libc_start_main in "/lib/x86_64-linux-gnu/libc.so.6"
#11 ? at ??:?
Floating point exception (core dumped)

What should I do to transform a transient case into steady state case ?
Moreover, I use sonicFoam because I heard that rhoSimpleFoam doesn't provide good result for high velocity. Is it true, why ? I would like to simulate a mach 3 simulation.


Thank you for your help.

[Dorian1504]

No one ?

I've also tried to run keeping euler parameter, however after few iterations it diverges.

Time = 8e-08

Courant Number mean: 0.000149386 max: 0.609122
diagonal:  Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
PIMPLE: iteration 1
DILUPBiCG:  Solving for Ux, Initial residual = 6.13615e-09, Final residual = 1.85653e-12, No Iterations 1
DILUPBiCG:  Solving for Uy, Initial residual = 5.16652e-06, Final residual = 2.86941e-10, No Iterations 1
DILUPBiCG:  Solving for Uz, Initial residual = 9.55783e-08, Final residual = 1.32978e-11, No Iterations 1
DILUPBiCG:  Solving for e, Initial residual = 0.0202665, Final residual = 2.08286e-10, No Iterations 3
DILUPBiCG:  Solving for p, Initial residual = 0.00535035, Final residual = 2.12454e-13, No Iterations 4
diagonal:  Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
time step continuity errors : sum local = 1.82758e-16, global = -2.99055e-18, cumulative = -3.59721e-16
DILUPBiCG:  Solving for p, Initial residual = 7.05603e-09, Final residual = 6.17777e-13, No Iterations 2
diagonal:  Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
time step continuity errors : sum local = 2.1755e-16, global = -8.59089e-18, cumulative = -3.68312e-16
DILUPBiCG:  Solving for epsilon, Initial residual = 2.03893e-09, Final residual = 3.88331e-13, No Iterations 1
DILUPBiCG:  Solving for k, Initial residual = 5.73229e-06, Final residual = 1.00964e-12, No Iterations 3
bounding k, min: -1516.73 max: 186866 average: 1271.06
ExecutionTime = 92.55 s  ClockTime = 96 s

forceCoeffs forces write:
    Cm    = -0.0115058
    Cd    = 0.021749
    Cl    = 0.00377122
    Cl(f) = -0.00962021
    Cl(r) = 0.0133914

Time = 1.2e-07

Courant Number mean: 0.000149386 max: 179.954
diagonal:  Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
PIMPLE: iteration 1
DILUPBiCG:  Solving for Ux, Initial residual = 5.81952e-09, Final residual = 1.14377e-12, No Iterations 1
DILUPBiCG:  Solving for Uy, Initial residual = 5.2802e-06, Final residual = 3.18083e-12, No Iterations 2
DILUPBiCG:  Solving for Uz, Initial residual = 9.16697e-08, Final residual = 9.40797e-10, No Iterations 1
DILUPBiCG:  Solving for e, Initial residual = 0.000202044, Final residual = 7.08963e-10, No Iterations 5


--> FOAM FATAL ERROR: 
Maximum number of iterations exceeded

    From function Foam::scalar Foam::species::thermo<Thermo, Type>::T(Foam::scalar, Foam::scalar, Foam::scalar, Foam::scalar (Foam::species::thermo<Thermo, Type>::*)(Foam::scalar, Foam::scalar) const, Foam::scalar (Foam::species::thermo<Thermo, Type>::*)(Foam::scalar, Foam::scalar) const, Foam::scalar (Foam::species::thermo<Thermo, Type>::*)(Foam::scalar) const) const [with Thermo = Foam::hConstThermo<Foam::perfectGas<Foam::specie> >; Type = Foam::sensibleInternalEnergy; Foam::scalar = double; Foam::species::thermo<Thermo, Type> = Foam::species::thermo<Foam::hConstThermo<Foam::perfectGas<Foam::specie> >, Foam::sensibleInternalEnergy>]
    in file /opt/OpenFOAM/OpenFOAM-4.x/src/thermophysicalModels/specie/lnInclude/thermoI.H at line 66.

FOAM aborting

#0  Foam::error::printStack(Foam::Ostream&) at ??:?
#1  Foam::error::abort() at ??:?
#2  Foam::hePsiThermo<Foam::psiThermo, Foam::pureMixture<Foam::constTransport<Foam::species::thermo<Foam::hConstThermo<Foam::perfectGas<Foam::specie> >, Foam::sensibleInternalEnergy> > > >::calculate() at ??:?
#3  Foam::hePsiThermo<Foam::psiThermo, Foam::pureMixture<Foam::constTransport<Foam::species::thermo<Foam::hConstThermo<Foam::perfectGas<Foam::specie> >, Foam::sensibleInternalEnergy> > > >::correct() at ??:?
#4  ? at ??:?
#5  __libc_start_main in "/lib/x86_64-linux-gnu/libc.so.6"
#6  ? at ??:?
Aborted (core dumped)

I do not know where the problem could come from, the 3D boundary conditions, the schemes? I've tried to keep as much as possible the structure of the tutorial case.

The checkMesh report:

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Create time

Create polyMesh for time = 0

Time = 0

Mesh stats
    points:           2227663
    faces:            5752919
    internal faces:   5381269
    cells:            1789298
    faces per cell:   6.22266
    boundary patches: 4
    point zones:      0
    face zones:       0
    cell zones:       0

Overall number of cells of each type:
    hexahedra:     1463941
    prisms:        68406
    wedges:        0
    pyramids:      0
    tet wedges:    658
    tetrahedra:    0
    polyhedra:     256293
    Breakdown of polyhedra by number of faces:
        faces   number of cells
            4   32968
            5   31963
            6   63479
            7   1582
            8   754
            9   74606
           10   532
           11   1114
           12   37023
           13   119
           14   1122
           15   10968
           18   63

Checking topology...
    Boundary definition OK.
    Cell to face addressing 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
                     box    60000    60300  ok (non-closed singly connected)
                   inlet     5000     5151  ok (non-closed singly connected)
                  outlet     5000     5151  ok (non-closed singly connected)
                 airfoil   301650   354595      ok (closed singly connected)

Checking geometry...
    Overall domain bounding box (-20 -10 -5) (20 10 5)
    Mesh has 3 geometric (non-empty/wedge) directions (1 1 1)
    Mesh has 3 solution (non-empty) directions (1 1 1)
    Boundary openness (2.06255e-16 -1.53285e-16 2.26207e-16) OK.
    Max cell openness = 3.56498e-16 OK.
    Max aspect ratio = 4.89722 OK.
    Minimum face area = 3.1772e-08. Maximum face area = 0.0400983.  Face area magnitudes OK.
    Min volume = 8.75771e-12. Max volume = 0.00801009.  Total volume = 7999.91.  Cell volumes OK.
    Mesh non-orthogonality Max: 55.825 average: 10.378
    Non-orthogonality check OK.
    Face pyramids OK.
    Max skewness = 3.28533 OK.
    Coupled point location match (average 0) OK.

Mesh OK.

End

[sheaker]

Hello. Im not sure but it looks like the only compressible steady-state solver is rhoSimpleFoam (openFoam 1.6-ext):
Description
Steady-state SIMPLE solver for laminar or turbulent RANS flow of
compressible fluids.


I can see in Your log:
bounding k, min: -1516.73 max: 186866 average: 1271.06

It usually means that Your inlet k and epsilon are far away from real values. // Not usually but it may cause bounding.

Also You should keep an eye on Your mesh. High non-orthogonality is a common problem of divergence.

Have a nice day.
Sheaker

[Dorian1504]

Thank you for your reply !

So, if I well understood I need to change the inlet boundary conditions for k and epsilon in order to be closer of real values? I've just kept the value writtent on the tutorial case.

Ok, I will follow your advice and try to increase the mesh quality. Is there an efficient way to do it ? Like changing the number of non-orthogonal correctors I guess.

Dorian

[sheaker]

Setting additional correctors loops could possibly solve the problem.

I am also trying to simulate an airfoil but in my case it is an incompressible transient flow (later I will try a compressible transient solver).
I have prepared mesh for NACA 65(2)-415 with wall function (y+35). I cannot check it in simulation because my laptop is down but I think the mesh is pretty good and should work well. (See attachments.)
Here is a case I will use with my mesh and own settings just when my laptop back from service: https://github.com/traviscarrigan/OpenFOAM-NACA2412

If You are interested in calculating drag and lift coefficient then Your mesh is a key to right results.

If You want to, I can give You a simple c++ code that calculates k, epsilon, y+ and some parameters for Salome Viscous Layer 2D.


Have a nice day.
Sheaker

[Dorian1504]

I tried to perform a simulation with a much more simple geometry than an airfoil and it seemed to run much better. Thus I'm almost sure the problem is my mesh as you said.

I'm still a beginner in meshing, I've used snappyHexMesh to create mine. You said you have used salome mesher for your case, is there any advantages by using salome instead of snappyHexMesh for a 3D-case ? I take any advice which can help me to create a suitable mesh.

Thank you for your case, I will study it too. Yes I'm well interested by the calculation of aerodynamics coefficients.

I would be glad to see your code for the calculation of the parameters, everything is good to increase my understanding of CFD. I hope to be able to modify it to match my own case.

Thank you very much for your help !

[sheaker]

I haven't use snappyHexMesh so I can tell You nothing about advantages and disadvantages of this tool. I think that snappyHexMesh cannot create proper wall function. Salome allow You to use a lot of hypothesis for meshing 3D,2D and 1D elements while snappyHexMesh isn't so customizable.
You said You are interested in 3D case so it isn't an airfoil but a wing, right? I think that meshing a 3D wing with a wall function would be extremely hard. It took me a few days to found out a way to create a 2D airfoil mesh. If I were You I would start from 2D airfoil mesh trying to match values of lift and drag coefficient with NACA data sheet. (That is exactly what I'm trying to do.) Then, when You successful simulate an airfoil You could start a 3D case but...
My mesh is about 180 000 elements. I believe a 3D case could need 200 times more elements and it is 36 000 000. (It is just my coarse calculation.) Are You sure You could afford so big simulation?

About mesh itself. It is a long way to create mesh in Salome and I am unable to give You a setp-by-setp solution. But take those hints:
- Chose Your parameters important for creating a geometry and mesh. Use attached code.txt. You will need velocity, density, kinematic viscosity (first the base, then the exponent), reference length (chord).
You will get Re number, Turbulence intensity, k and epsilon. Then wall function calculator for salome will start.
- Chose more parameters. You will need to set y+ (for example 35)
You will get thickness of the first layer for y+35,
- Estimate total thickness of viscous layers (for example 100 mm), set stretch factor ratio (for example 1.1, stretch factor ratio is a ratio between two neighboring cells).
You will get proper total thickness of viscous layers and important Scale Factor parameter (scale factor is a ratio between first and last cell).

For Your geometry You need total thickness of viscous layers.
- You could create a sketch like in attachment nr 1. Close-up on tailing edge in attachment nr 2. As You can notice there is an airfoil which is offset of distance of total thickness of viscous layers and then some links are created. You could create it in freeCAD (it doesn't allow to offset, I think. I used AutoCAD.) and export as *.iges file. Then import to geometry module in salome. Don't forget to create background.
- You will need to create groups of edges (see attachment nr 3) for every part of Your geometry.

In mesh module use NetGen 1D-2D for background. For viscous layers parts use button "assign a set of hypothesis: 2D automatic quadriangulation" and then add sub-meshes. Sub-meshes based on geometry groups You have created before. Use wire discretisation with local length hypothesis on group 3 and wire discretisation with nb. segments hypothesis on group 1 and group 2. Also on group 1 and group 2 click button on the right from the hypothesis menu to set the settings. Use number of segments from wall function calculator. Use type of distribution to scale distribution and set scale factor from wall function calculator. (Sometimes You need to revers edges.) Always create new hypothesis instead of using existing one.

- Mesh every part separately and then create a compound mesh.

- Somehow it is important to create geometry from one rectangular face. I mean I have to create special cut_tools for every single part of my geometry and cut the final geometries parts from this one rectangular face.

- Create flat 2D mesh, create groups of edges with names of patches for openFoam. Then extrude mesh in Z direction (1 cell depth) to create 2D mesh for openFoam. Remove groups of edges. Change names in groups of faces.

- Export as UNV. In openFoam use ideasUnvToFoam mesh.unv and transformPoints -scale '(0.001 0.001 0.001)'.

That was my procedure for 2D airfoil mesh but I'm not an expert and I'm not sure if this is a good way to create mesh.

I'm sorry for chaotic reply - my language is very complicated.

I hope this will help You.

Have a nice day.
Sheaker

[Dorian1504]

Thank you four all these details, I realize I was really naive as to my approach and I have under estimated the difficulties to run a simulation (including meshing).

In fact, I'm interested in a 3d-case, I didn't know the difference between airfoil and wing until now, I still need to improve my english. I have to succeed to perform accurate simulations in order to run an optimization algorithm using the aerodynamic coefficients calculated by the openFoam. So, I will have to create suitable mesh for every airfoils created by my generator and check their coefficients. Run a 3D-case should be ok with the equipment of my laboratory.

Starting directly by a 3d wing seems to be pretentious, then I will begin with the 2d-airfoil.

If you agreed, is there a way to communicate with you and sharing our results or problems?

Thanks again

Dorian

[sheaker]

I wonder why non expert answers You. I think there should be many "foamers" with actual experience in calculating aerodynamical coefficients.

I found something that maybe will be useful for us. Haven't read it yet but it looks like a simple airfoil calculations:
https://www.hpc.ntnu.no/display/hpc/...l+Calculations

I have three more things in my mind:
1. Indifferent whether you like to use wall function or not you should double-check proper y+ size. Long time ago I read somewhere that y+ of wall function for the finite volume method should be different, than for the finite element method. I'm not sure about this, but it seems reasonable. I found this but i haven't read it yet:
How to calculate y-plus?
2. Your velocity on the surface of Your wing changes with chord:
https://www.mh-aerotools.de/airfoils...s/veldist1.gif
but You should keep your y+ value constant on the entire surface of Your wing, I guess. This led me to another question: Do You have experimental data or You will be trying to fit simulation results to theoretical aerodynamics coefficients?
If You have experimental data of velocities it will be easier;
else, You should estimate those velocities and then prepare suitable mesh.
3. If You are going to simulate high angle of attack then k-epsilon turbulence model isn't good. As I remember k-epsilon loses separation point and doesn't simulate its effect. Try this, i haven't read it yet:
http://www.uobabylon.edu.iq/uobColeg...r.Emad%20).pdf

I will send You a PM with my email. It will be easier to share files.

Have a nice day.
Sheaker

[JPBLourenco]

Originally Posted by sheaker

I haven't use snappyHexMesh so I can tell You nothing about advantages and disadvantages of this tool. I think that snappyHexMesh cannot create proper wall function. Salome allow You to use a lot of hypothesis for meshing 3D,2D and 1D elements while snappyHexMesh isn't so customizable.
You said You are interested in 3D case so it isn't an airfoil but a wing, right? I think that meshing a 3D wing with a wall function would be extremely hard. It took me a few days to found out a way to create a 2D airfoil mesh. If I were You I would start from 2D airfoil mesh trying to match values of lift and drag coefficient with NACA data sheet. (That is exactly what I'm trying to do.) Then, when You successful simulate an airfoil You could start a 3D case but...
My mesh is about 180 000 elements. I believe a 3D case could need 200 times more elements and it is 36 000 000. (It is just my coarse calculation.) Are You sure You could afford so big simulation?

About mesh itself. It is a long way to create mesh in Salome and I am unable to give You a setp-by-setp solution. But take those hints:
- Chose Your parameters important for creating a geometry and mesh. Use attached code.txt. You will need velocity, density, kinematic viscosity (first the base, then the exponent), reference length (chord).
You will get Re number, Turbulence intensity, k and epsilon. Then wall function calculator for salome will start.
- Chose more parameters. You will need to set y+ (for example 35)
You will get thickness of the first layer for y+35,
- Estimate total thickness of viscous layers (for example 100 mm), set stretch factor ratio (for example 1.1, stretch factor ratio is a ratio between two neighboring cells).
You will get proper total thickness of viscous layers and important Scale Factor parameter (scale factor is a ratio between first and last cell).

For Your geometry You need total thickness of viscous layers.
- You could create a sketch like in attachment nr 1. Close-up on tailing edge in attachment nr 2. As You can notice there is an airfoil which is offset of distance of total thickness of viscous layers and then some links are created. You could create it in freeCAD (it doesn't allow to offset, I think. I used AutoCAD.) and export as *.iges file. Then import to geometry module in salome. Don't forget to create background.
- You will need to create groups of edges (see attachment nr 3) for every part of Your geometry.

In mesh module use NetGen 1D-2D for background. For viscous layers parts use button "assign a set of hypothesis: 2D automatic quadriangulation" and then add sub-meshes. Sub-meshes based on geometry groups You have created before. Use wire discretisation with local length hypothesis on group 3 and wire discretisation with nb. segments hypothesis on group 1 and group 2. Also on group 1 and group 2 click button on the right from the hypothesis menu to set the settings. Use number of segments from wall function calculator. Use type of distribution to scale distribution and set scale factor from wall function calculator. (Sometimes You need to revers edges.) Always create new hypothesis instead of using existing one.

- Mesh every part separately and then create a compound mesh.

- Somehow it is important to create geometry from one rectangular face. I mean I have to create special cut_tools for every single part of my geometry and cut the final geometries parts from this one rectangular face.

- Create flat 2D mesh, create groups of edges with names of patches for openFoam. Then extrude mesh in Z direction (1 cell depth) to create 2D mesh for openFoam. Remove groups of edges. Change names in groups of faces.

- Export as UNV. In openFoam use ideasUnvToFoam mesh.unv and transformPoints -scale '(0.001 0.001 0.001)'.

That was my procedure for 2D airfoil mesh but I'm not an expert and I'm not sure if this is a good way to create mesh.

I'm sorry for chaotic reply - my language is very complicated.

I hope this will help You.

Have a nice day.
Sheaker

Hi everyone.

Oskar, is it possible to have the dump file from SALOME (python script) to check how you did the boxes and mesh?
I'm trying to do the same but I'm with some problems in how I do the boxes to create the mesh.

Thanks
JoaoLourenco

[sheaker]

Hello.
I believe there is a horrible chaos inside file You ask me about.
Take this Salome project file:
https://drive.google.com/open?id=0B6...UU1V19qUXZ3d3c
It should be NACA 65(2)-415.
http://airfoiltools.com/airfoil/deta...=naca652415-il

I think it will be easier for You to preview only the right one because there are two or more bad meshes. Also many of objects are unpublished.

You were asking me for dump file:
http://wklejto.pl/txt301920
http://wklejto.pl/txt301921
I'm not proud of this files but in the end the mesh is pretty good. I was learning!

Here is checkMesh:

Mesh stats
    points:           116414
    internal points:  0
    faces:            290038
    internal faces:   134398
    cells:            77277
    boundary patches: 4
    point zones:      0
    face zones:       0
    cell zones:       0

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

Checking topology...
    Boundary definition OK.
    Cell to face addressing 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                  
    inlet               90       182      ok (non-closed singly connected)  
    outlet              90       182      ok (non-closed singly connected)  
    wall                906      1812     ok (non-closed singly connected)  
    sides               154554   116414   ok (non-closed singly connected)  

Checking geometry...
    Overall domain bounding box (-1.5 -1.5 0) (4.5 1.5 1)
    Mesh (non-empty, non-wedge) directions (1 1 0)
    Mesh (non-empty) directions (1 1 0)
    All edges aligned with or perpendicular to non-empty directions.
    Boundary openness (5.03049e-19 6.03659e-18 -4.73228e-15) OK.
    Max cell openness = 3.90892e-16 OK.
    Max aspect ratio = 15.6726 OK.
    Minumum face area = 1.8771e-09. Maximum face area = 0.141421.  Face area magnitudes OK.
    Min volume = 1.8771e-09. Max volume = 0.00727518.  Total volume = 17.9627.  Cell volumes OK.
    Mesh non-orthogonality Max: 35.748 average: 5.02129
    Non-orthogonality check OK.
    Face pyramids OK.
    Max skewness = 0.368642 OK.
    Coupled point location match (average 0) OK.

Mesh OK.

Will You be able to set rhoSimpleFoam for NACA simulation? I succeed with rhoPimpleFoam and sonicFoam solvers with kEpsilon and kOmegaSST turbulence models but not with rhoSimpleFoam.

Have a nice day.
Sheaker