[Alan]

Hi all,

I am trying to use sonicFoam for some inviscid simulations, and my test case (supersonic flow over a wedge) does not seem to be following the analytical solution.

I am simulating a 14.03 deg wedge at Mach 7, sea level atmospheric conditions. If you look at the image below, it shows that the shock sits lower than the analytical solution (the line source).



I can adjust the thermophysical properties get the shock to sit in the right position, but this does not give me much faith in other solutions! My thermophysical properties are:

thermoType      ePsiThermo<pureMixture<constTransport<specieThermo<hConstThermo<perfectGas>>>>>;

mixture         air 1 28.96 1005 0 0 0.715;

If I adjust my cp I can get the shock to sit in the right position. I have also tried a finer mesh, but there is not change in shock position, so I do not believe it is a meshing issue.

Does anyone have any ideas as to what could be causing this? I have simulated the same wedge in other solvers (rhoCentralFoam and rhopSonicFoam) and the results match the analytical theory.

Thanks in advance for any help,
Alan

[alberto]

Hello, I'm sorry to answer with a question: what does rhoSonicFoam does?

The major difference I can see between the solvers you used is that rhopSonicFoam and rhoCentralFoam solve for conserved variables, while sonicFoam does not.

I hope this additional test can help to clarify.

Best,

[Alan]

No problems at all Alberto, your help is much appreciated. Here is an image of the wedge simulated in rhoSonicFoam, as you can see shock lies in the correct position.



Thanks again,
Alan

[alberto]

Interesting.

To exclude a possibility: in your run in rhoCentralFoam, did you set thermoPhysical properties to the same data you report for air (some tutorials use the normalised gas)?

Best,

[Alan]

Hi Alberto,

If I run the wedge in rhoCentralFoam with the same thermophysical properties as my sinicFoam case (i.e. not normalised) then the shock matches the theory. It appears to be a problem with the sonicFoam solver, as it gives the same (incorrect) shock position in both normalised and regular thermophysical properties.

Thanks again for your help, it is much appreciated.

[alberto]

If there is no bug, it might be related to how the equations are discretized, as I said before, but it's just a guess.

I'd suggest to use codes with the conservative form of the equations for shock capturing.

Best,

[KrisT]

I tried sonicFoam on OF-1.3, shock tube case in tutorial. It gives the wrong result (as said before, probably due to the fact that scheme is not conservative).

I checked against exact solution, obtained from C-code in the link below

http://num.math.uni-goettingen.de/kn...ro_ss2006.html

Haven't checked later versions (but that would be easy to do).

[Alan]

So is there a way to fix this? sonicFoam seems to be the obly stable solver for my problem, the others seem to get anomalies in the solution that diverge

[KrisT]

Originally Posted by Alan

So is there a way to fix this? sonicFoam seems to be the obly stable solver for my problem, the others seem to get anomalies in the solution that diverge

What kind of problems do you get with rhoCentralFoam? Have you tried varying CFL on that solver?

[Alan]

Well, in rhoCentralFoam I am unable to progress the solution without it crashing within a few iterations (getting a maximum number of iterations exceeded error). This is the error I get:

Create time

Create mesh for time = 0.00264945

Reading thermophysical properties

Selecting thermodynamics package ePsiThermo<pureMixture<constTransport<specieThermo<hConstThermo<perfectGas>>>>>


Maximum number of iterations exceeded#0  Foam::error::printStack(Foam::Ostream&) in "/home/aharrland/OpenFOAM/OpenFOAM-1.6/lib/linuxGccDPOpt/libOpenFOAM.so"
#1  Foam::error::abort() in "/home/aharrland/OpenFOAM/OpenFOAM-1.6/lib/linuxGccDPOpt/libOpenFOAM.so"
#2  Foam::ePsiThermo<Foam::pureMixture<Foam::constTransport<Foam::specieThermo<Foam::hConstThermo<Foam::perfectGas> > > > >::calculate() in "/home/aharrland/OpenFOAM/OpenFOAM-1.6/lib/linuxGccDPOpt/libbasicThermophysicalModels.so"
#3  Foam::ePsiThermo<Foam::pureMixture<Foam::constTransport<Foam::specieThermo<Foam::hConstThermo<Foam::perfectGas> > > > >::ePsiThermo(Foam::fvMesh const&) in "/home/aharrland/OpenFOAM/OpenFOAM-1.6/lib/linuxGccDPOpt/libbasicThermophysicalModels.so"
#4  Foam::basicPsiThermo::addfvMeshConstructorToTable<Foam::ePsiThermo<Foam::pureMixture<Foam::constTransport<Foam::specieThermo<Foam::hConstThermo<Foam::perfectGas> > > > > >::New(Foam::fvMesh const&) in "/home/aharrland/OpenFOAM/OpenFOAM-1.6/lib/linuxGccDPOpt/libbasicThermophysicalModels.so"
#5  Foam::basicPsiThermo::New(Foam::fvMesh const&) in "/home/aharrland/OpenFOAM/OpenFOAM-1.6/lib/linuxGccDPOpt/libbasicThermophysicalModels.so"
#6  main in "/home/aharrland/OpenFOAM/OpenFOAM-1.6/applications/bin/linuxGccDPOpt/rhoCentralFoam"
#7  __libc_start_main in "/lib/tls/i686/cmov/libc.so.6"
#8  _start at /usr/src/packages/BUILD/glibc-2.9/csu/../sysdeps/i386/elf/start.S:122


    From function specieThermo<thermo>::T(scalar f, scalar T0, scalar (specieThermo<thermo>::*F)(const scalar) const, scalar (specieThermo<thermo>::*dFdT)(const scalar) const) const
    in file /home/dm2/henry/OpenFOAM/OpenFOAM-1.6/src/thermophysicalModels/specie/lnInclude/specieThermoI.H at line 68.

FOAM aborting

Here is an image I get of the solution at the point of it crashing:



This is what the simulation looks like in rhopSonicFoam:



Note the 'stagnation' points at the top of the compression wedge. These slowly over the time of the solution go to zero, causing the time step (with adjust time step on) to go to a very large number.

Here is the simulation in sonicFoam:



This gives a nice stable solution, although with the incorrect shock angle.

I have a well refined, structured mesh, so I do not believe this to be the problem. This is also only a simple case which I need to expand to much more complex geometries, so the apparent stability of sonicFoam is quite appealing to me.

Thanks for your help KrisT.

[Alan]

Originally Posted by KrisT

What kind of problems do you get with rhoCentralFoam? Have you tried varying CFL on that solver?

No matter how small I make the CFL it still crashes.

[alberto]

Originally Posted by Alan

No matter how small I make the CFL it still crashes.

In rhoCentralFoam, what is giving troubles to you seems to be the energy equation, given the error message you obtain. Is it converging before the crash?

Are you performing inviscid runs in rhoCentralFoam? It seems to me the behaviour at the walls is different.

Additionally, could you show a zoom of the mesh with the solution under it (Surface with edges in paraview) around the tip of the triangle at the exit (where rhoCentralFoam has a sort of blue circle in the solution).

Best,
Alberto

[Alan]

Hi Alberto,

Im not sure what you mean by 'is it converging before the crash?'. Do you mean is the energy equation converging, or is the entire solution converging? I am solving this case as inviscid, as far as I am aware rhoCentralFoam is only an inviscid solver? My wall BC's are slip for the velocity file, and zeroGradient for p and T.

Here is the image that you asked for (note the flow is from left to right, so this is at the inlet, not the outlet):



I have also attached a similar picture around the entry to the isolator, which is the area of flow I believe is giving the solvers all the trouble.



Thanks again for your help Alberto, it is very much appreciated.

Alan

[alberto]

Originally Posted by Alan

Hi Alberto,

Im not sure what you mean by 'is it converging before the crash?'. Do you mean is the energy equation converging, or is the entire solution converging?

Yes, I'm asking if the energy equation is converging.

I am solving this case as inviscid, as far as I am aware rhoCentralFoam is only an inviscid solver? My wall BC's are slip for the velocity file, and zeroGradient for p and T.

No, rhoCentralFoam can solve for both the inviscid and the viscous flow. Check the code (see flag "inviscid").

Here is the image that you asked for (note the flow is from left to right, so this is at the inlet, not the outlet):

The grid seems good. Maybe it won't make any difference, but could you avoid having that tip touching the inlet conditions?

Best,
Alberto

[KrisT]

Maybe I'm stupid here, but what happens if you ramp up the Mach number?

Would you get a stable solution for a lower Mach number? Is it possible to specify varying (increasing) inlet Mach number?

/j

[immortality]

Hi Alan
although a somewhat old thread.but I saw the same problem when a moving shock goes near the wall.residuals of U.x and h increases,although crashing doesn't occur.I use maxCo too.
has had anyone this problem too to share?
well,then sonicFoam solves only non-conservative form of equations.that was a very good information you gave to me!