No shock in airfoil 0012 case despite of Mach number exceeds 1
 

Hi Foamers,
I'm pretty new to OpenFoam and CFD in general. I managed to get good results for low subsonic speeds on the NACA 0012. Now I'm trying to go transsonic for my thesis.

I'm struggling to get a shock visible on an NACA 0012 at freestream mach number of 0.82. First I decided to use rhoSimpleFoam. Correct me if I'm wrong but there are no non-reflecting boundary conditions for this solver. So I went on and tried the following:

• rhoPimpleFoam
• kOmegaSST turbulence model with Wallfunction
• waveTransmissive boundary conditions for in- and outlets
• adjustTimeStep with maxCo 0.95

Mesh:

• the has 125426 cells
• Mesh non-orthogonality Max: 79.17, average: 16.17
• Picture see below

Solver:

• relaxation: p. 0.4; rho 0.6; "(U|h|k|epsilon|omega).*" 0.7;
• for p : solver PCG;
  preconditioner DIC;
• for "(rho|U|h|k|epsilon|omega)" : solver PBiCG;
  preconditioner DILU;

Results:

• Yplus is a little too low with around 14 but increasing it doesn't solve the problem
• k, and Omega stay constant throughout the domain. The logfile says nothing about solving for k or omega which is very strange. I'm using waveTransmissve for them as well
• Mach number exceeds 1 but no shock occurs. Below Mach number and Isobars are plotted
• I have a flow acceleration in some cells behind the airfoil see picture below

After I wasn't lucky with rhoPimpleFoam, I also tried sonicFoam but didn't succeed either to get a shock. Neither increasing the velocity nor changing schemes or increasing domain size helped.

I would very much appreciate some critics or advise how I could go on, as I'm running out of Ideas.

Thanks a lot

Use rhoCentralFoam. That is good for shock capturing and compressible flows.

I agree with praveen. With Ma > 1 you should not use a pressure-based solver.

Thanks a lot for the hint rhoCentralFoam really works much better than the others for my case. Even if I don't understand why, as I also tried
- rhoPimpleFoam (which should be density based right ?)
- sonicFoam
- rhoSimpleFoam

I didn't used rhoCentralFoam in the first place as the manual suggests that it does not consider turbulence. But it does.

I'll put my results in here if they're finished. Thanks lot !!

Hi Alberto, could you please justify this statement? To my knowledge, pressure-based solvers CAN handle transonic flows with shocks, and this has been established quite a long time ago (see for instance "Turbulent Transonic Flow Simulation Using a Pressure-Based Method", Y. G. Lai, R. M. C. SO and A. J. Przekwas, International Journal of Engineering Science, vol 33, issue 4, pp 469-483, 1995). I'm not a density-based solver expert, so of course I might be wrong, but I think that the "supremacy" of density-based approaches for this kind of problems is mainly due to the higher computational efficiency in the solution of Riemann-type problems, rather than a matter of accuracy. What's your opinion about this?

Regards

V.

Hi,

yes, they can manage "high-Mach" number flows. The literature on "all-speeds" pressure based solvers is very long, and they all basically rely on similar algorithms.

In practice, compressible methods allow better numerical schemes to be implemented, which reflect the physics of the problem, and naturally deal with some of the numerical difficulties presented by this kind of flow taking advantage of the mathematical nature of the equations.

P.S. I am not an expert either, but recent experience showed that using a pressure based code to solver high-speed compressible flows is a waste of time. Way too many convergence and stability problems compared to a density-based one.

Best,

rhoPimpleFoam is pressure-based compressible and unsteady.

Please report it on http://www.openfoam.com/mantisbt/main_page.php so they can fix that.

Best,

Well, this is quite interesting, since now I'm working on transonic internal flows and indeed I'm facing stability issues with the pressure based sonicFoam solver (although I'm quite sure that some of them are mesh-quality related)...I'll definitely give a try to the rhoCentralFoam (with turbulence added) solver to see if things are better in an overall sense (accuracy+stability). Thank you for the interesting discussion!

Best

V.

Hi schwermetall,
It would be great if you could share your case setup as well. It would serve as a great starting point for my airfoil investigation. Thx!

Hi Foamers,
first of all thanks for the support.
Unfortunately I don't get a steady state solution. At first a I got two shocks on the upper and lower side in the region where they belong. But with increasing time my solution gets unsteady.
There is a region of low pressure that emerges from the trailing edge disturbing the complete flow field around the airfoil. This behaviour could be seen in sonicFoam as well, see first post.

What would you consider a reasonable physical time after the flow around the airfoil should reach a steady state if the conditions are as follows
- freestream velocity: 277 m/s
- airfoil length 1 m
- domain length 40 m
- Co<0.5
- 125 000 cells

I added some picture with different time steps. I don't get a steady state after after 0.14 seconds.
Very strange is, that there are two shocks on the upper side at time 0.11. After some time they occur on the lower side.

thanks a lot

Maybe it needs more iterations to reach steady state. Note that rhoCentralFoam uses global time stepping, so reaching steady state can be very slow. And moreover it uses forward euler time stepping. Atleast, some 2/3/4 stage RK scheme with local time stepping would be more robust and better for steady state problems. But this needs to be coded and is not available as a scheme. But your problem could be something else also, I cannot say.

It seems I found the problem concerning the unsteady behavior. I changed the mesh at the trailing edge, so that cells around the last point of the airfoil get less skewed. Below you can see pictures from physical time equals 0.152

Nevertheless the fluctuations in density and velocity at the leading edge remain. Does anyone have an idea what they could come from? I already changed the default divScheme from linear to upwind, but that doesn't help.
Below my fvSchemes for the rhoCentralFoam solver:

fluxScheme Kurganov;
ddtSchemes
{
default Euler;
}
gradSchemes
{
default Gauss linear;
}
divSchemes
{
default Gauss upwind;
div(tauMC) Gauss linear;
}
laplacianSchemes
{
default Gauss linear corrected;
}
interpolationSchemes
{
default linear;
reconstruct(rho) vanLeer;
reconstruct(U) vanLeerV;
reconstruct(T) vanLeer;
}
snGradSchemes
{
default corrected;
}

Grateful for any hints

From your fvSchemes file it is not clear if you are using a turbulence model or not. In addition, your mesh seems of good quality but maybe a little coarse near the wall, considering the very high speed of your case. So, what are your turbulence model, approximate max y+ value at the airfoil surface and, consequently, wall treatment?

Best

V.

hi vkrastev
I'm using RAS model with kOmega turbulence model plus wallfunction, not with resolved boundary layer. My y+ ranges between 26 and 167 (see below) . So consequently the mesh should be even a little coarser near the surface right ?

Boundary consitions are
fixedValue for Velocity at the inlet
zeroGrad for Velocity at the outlet

and
fixedValue for pressure at the outlet
zeroGrad for pressure at the inlet

No, your y+ in my opinion is ok for a wall function treatment. Instead, your boundary conditions seem not so adequate for a compressible case. I don't know if this will solve your problem, but i would try something different both for pressure and velocity. If the incoming flow is already supersonic you can try (see the nacaAirfoil tutorial inside the sonicFoam/ras/ tutorial folder):

supersonicFreeStream for velocity at the inlet
zeroGradient for pressure at the inlet

inletOutlet for velocity at the outlet
waveTransmissive for pressure at the outlet

Otherwise, if the incoming flow is subsonic,you can change the inlet conditions in:

fixedValue for velocity
waveTransmissive for pressure

Regards

V.

PS - Now I recall that your inlet condition was subsonic, so you can try directly the second option

Hi
thanks for the advice. I thought about changing the boundary condition, but the thing that kept from doing it was, that I couldn't find any pressure waves being reflected at the boundaries. Shouldn't I see at least anything coming back into the domain ?

Nevertheless I'm going to try it and I'll report the results.
By the way what do you recommend for this lInf value when using waveTransmissive? I already played with it but I'm not sure.
From what I understand it is the distance behind the boundary where the given boundary value will be reached ?

Regards

Well, actually I haven't played that much with the waveTransmissive BC's, but the fact is that the lesser the lInf the more reflective will become the BC (in the limit of lInf=0 you'll have actually a totally reflective fixed pressure condition), so maybe you can start with a value of the order of the chord lenght of the airfoil and see what happens. Regarding the pressure waves, maybe they can occur in some time steps out of your saving points, bun even if they don't there could be some backflow at the boundaries and the inletOutlet condition on velocity can (in principle) account for these backflows while the zeroGradient condition cannot. Anyway, as I said in my previous post, maybe the problem are not the BC's, but is always better to apply the most appropriate combination of them, in order to focus on other (eventual) issues.

Best

V.

Hi all
I would like start with your mesh first. make the y+ to 2.
Then you should use boundary conditions as suggestd by Vkrastev.
Infact simply Inlet BC works good so that you can state all pressure, temperture and velocity etc... at inlet and at oulet you can make use of zerogradient B.C for all.
Use sonicFoam solver this is a turbulent compressible flow. If you are using 1.6 version there are some issues with sonicFoam solver (you cannot resolve thermal boundary layer).

Density based solvers are generally used for not only capturing shocks but also their interactions and these are very sensitive. where pressure based solvers are not better for these cases.

use small time step for analysis.

Thanks
Kiran Ambilpur

Not if he is using wall functions (wall functions work bad with too low y+ values: this is a quite basic CFD rule for wall treatment of high-Re flows).

Well, putting to zeroGradient all the non-fixed quantities works usually fine for incompressible flows, but not so in compressible cases where a bit much care should be employed in the BC choice.

sonicFoam is an option, but (to my experience) it also has stability issues in transonic flows, while the density based rhoCentralFoam (with turbulence added) seem more stable (though requiring a bit smaller time step size).


This is of course desirable whatever the solver (for rhoCentralFoam maxCo < 0.5 should be fine).

Regards

V.

I already tried sonicFoam (2.0) using different boundary conditions schemes etc etc ....
But I wasn't able to get a shock visible, with that solver.

I'm using adjustTimeStep yes; with maxCo of 0.5;

So thanks for the ideas, but I already tried that.

Hi guys
I now changed the boundary conditions for velocity at the outlet to
type inletOutlet
inletValue (0 0 0)

and pressure at the outlet to
waveTransmissive with lInf 2

Doesn't seem to change anything.

I'm now trying somthing completly different. I'm using this MUSCL scheme for divSchemes. Does any of you has experience with it ?

From what I understand this schemes tries to keep fluctuations as low as possible. Which is exactly my problem. Any experience with that ??

If you have already tried upwind (and you did, right?), this was the most stable choice you can make about divSchemes. Why don't you try instead with a different turbulence model (for instance, the Spalart-Allmaras one)? Another possibility would be to resolve the boundary layer with a low-Re model, but this would increase a lot the computational time. Finally, you should also consider that the flow could be instable in nature or that it would require a long time to reach a steady-state condition.

Best

V.

Hi vkrastev,
first of all many thanks for staying with me and my problem.

Yes I tried upwind, but not for div(tauMC). I'm now trying MUSCL 0.9 for div(phi,omega) and div(phi,k).

As it takes Ages to make a run, I haven't done major changes like changing the turbulence model today, but I'll let it run over night and report tomorrow.

This is also the reason why resolving boundary layer is no option for me. To get a physical time of 0.2 seconds, my computer is running for nearly 1.5 days. My target is to get an airfoil polar in the end (like Cl over alpha, Cd over Cl ...). So if the computation time increases further for each angle of attack I have to stop there.

No problem, I'm interested in the performances of rhoCentralFoam too (though my case is quite different from yours).

div(tauMC) is not a convective term, but an explicitly calculated divergence term, so I don't think it could make some difference in the stability of the calculations (actually I'm not even sure it could be calculated in other ways than linear interpolation, as in this case there's no convection velocity involved and thus the concept of "upwind direction" becomes meaningless). Instead, div(phi,<something>) are convective terms, which are always source of stability issues (and in my experience this is especially true for convective terms involving turbulent quantities): so, if you have already tried the upwind scheme on them, then you cannot do anything more about their stability behavior.

I understand. In this case you can change the turbulence model retaining the wf treatment near the walls (running for a bit more the simulation and see what happens remains also an option).

Best

V.

Hi vkrastev

concerning div(tauMC) I now think, that I remember that Openfoam returned an error when I tried the upwind scheme. But your explanation is pretty convinient.

short Update:
For the laplacianSchemes it is possible to choose cellLimitedSchemes. If I understand the intention correct, theses Schemes reduce fluctuations by limiting the gradient between neighbouring cell centers. (found a very good explanation http://www.cfd-online.com/Forums/ope...lllimited.html )

so I changed my fvSchemes as follows:

gradSchemes
{
default none;
grad(p) cellLimited leastSquares 0.9;
grad(U) Gauss linear;
grad(rho) cellLimited leastSquares 0.9;
grad(rhoU) cellLimited leastSquares 0.9;
grad((1|psi)) Gauss linear;
grad(e) Gauss linear;
grad(sqrt(((Cp|Cv)*(1|psi)))) Gauss linear;
grad(T) Gauss linear;
grad(omega) Gauss linear;
grad(k) Gauss linear;
}

The reason was, that pressure/density are probably the ones causing trouble.

It now seems (! no long time confirmation) as if the pressure oscillations are diminishing. See pictures below

div(tauMC) must be discretized with a central scheme, since it does not introduce directional bias. Direction-biased schemes are appropriate only for convection.

Best,

I would systematically apply the same limiter (I'd suggest cellLimited leastSquares 1 to all the gradients. Notice that the number here has a different effect in comparison to the effect it has in Laplacian and snGrad! 1 does NOT mean you do not limit.)

Best,

That is my point too.

Hey Foamers,
I did as alberto suggested:
gradSchemes
{
default cellLimited leastSquares 1;
{

The result got much better after that as you can see below. But the fluctuations haven't vanished completely.
So I'll now try the Spalart-Almaras model and see what happens then.

Just a short add-on from my side: it seems (at least in my internal flow case) that Tadmor's flux scheme is more stable than Kurganov's one (probably because it's slightly more dissipative, at least as I was able to understand from Kurganov and Tadmor's original paper).

Regards

V.

Hi Foamers,
as I'm running the Spalart-Almaras model I get the following warning:

--> FOAM Warning :
From function tmp<volScalarField> SpalartAllmaras::k() const
in file SpalartAllmaras/SpalartAllmaras.C at line 262
Turbulence kinetic energy not defined for Spalart-Allmaras model. Returning zero field

I had a look at the source code, but as my c++ is worse than my Chinese I can really tell what is happening there.
The Spalart-allmaras model doesn't need the kinetic turbulent energy at all. So I don't understand why the model is creating a field k mit zero entries.

I would like to understand that. Especially because when I ran the incompressible case, I didn't get that warning. In the incompressible Spalart-Allmaras code there is a similar line than the one shown above.

That piece of code is missing in OF 1.7.1 (which is my reference version), and I really don't understand why is there in OF 2.0.0/1/x ... I think you can ignore the warning message without any consequence for your calculations. Anyway, have you tried the Tadmor scheme instead of the Kurganov one?

Regards

V.

Hi vkrastev
just changed to Tadmor, but it does not seem to change the fluctuations.

As the fluctuations are spacial, I'm thinking the problem might be connected to fvSolution. Maybe I gonna try
instead of diagonal. What do you think?

No, I don't think you should change the diagonal solver options for rho, rhoU and rhoE, because as far as I have understood a little the rhoCentralFoam algorithm, these quantities are solved explicitly via the Kurganov/Tadmor reconstrucion procedure. What about the Spalart-Allmaras model compared to the kOmega?

V.

Hi vkrastev
my actual case is running with spalart-Almaras. As you can see below, the Amplitude of the fluctuations diminished a little bit.

The position of the shock probably becomes better after some time, thats at least what normally happens.

I'm open for new ideas ... thanks a lot

I think that you have to let the simulations go ahead until at least the shock position becomes sufficiently stable for both models and then make some comparison. In the end your results appear not bad, and honestly I don't have at the moment other suggestions for smearing out those little oscillations...

Regards

V.

So after waiting for ages I reached a time of 0.186 and you're the oscillations are slowly decaying, as you can see below. So I'm going to wait a little longer.

I got some input from my supervisor and she thinks the problem I have here is a odd even decoupling. I'll have a look at the literature, see what I find out about it and if there's a way to get better results

Do you mean a pressure-velocity decoupling? This could be interesting, but honestly I can't say if it is a reasonable hypothesis...Anyway, if I understand you properly, the pressure-velocity coupling in the rhoCentralFoam solver occurs in the inviscid part of the solution, which is driven by the Tadmor/Kurganov algorithms: so the answer to the question would be inside the nature of such algorithms, which (at least to my actual knowledge of compressible flow solution procedures) are quite a mystery to me, so you'll have to face it all by yourself..Good luck!

V.

Yes it's about a pressure velocity decoupling.

:D thanks a lot for all your help until now. I'll let you if I find anything out, but maybe I can unmystify it a little bit ;-)

Have nice weekend.

regards from Munich

Hi schwermetall,

I observed something very similar to this when simulating a supersonic channel flow using rhoCentralFoam and kOmegaSST. The wall pressure at top and bottom walls of the channel is fluctuating a lot and it looks a lot like in your profiles. However, temperature and velocity profiles along the channel are fine, so it only seems to be a pressure problem.

I did as you suggested and changed to cellLimited for grad schemes: It helped, but the fluctuations are not completely gone. As I am usually using grids with y+ of around 30 and wall functions I tried to refine the grid. For y+ = 15 and wall functions the instabilities looked weaker. And when I deactivated the wall functions for a fine grid of y+ = 1 also the pressure fluctuations were gone completely. However, the results now seem to suffer quite a lot from numerical dissipation and are considerably less accurate in temperature and velocity profiles.

Did you try to refine your grid at the walls? It would be interesting to see if you can observe the same behaviour of the solver for your NACA case.