Time-marching scheme diverges
 

Hey guys!

I'm currently working on the simulation of the flow behind a supersonic nozzle. I tried it for 3 weeks with different versions of the Maccormack scheme, finally I thought I have the answer to my problem by using the Maccormack scheme with the equations in conservation form.

Now my problem is, that the calculation blows up :(
For the calculation of the time steps I use the CFL criterion by keeping the Courant number constant. I use the following equation:

delta_t = C*delta_x/(a + V)

Where C is the Courant number, a the speed of sound and V the magnitude of the velocity.

The values for delta_t are around 10^-5 at the beginning but decrease to values around 10^-147.
I think that the cause for the small values is, that my x component of the velocity becomes really large, but I don't know why. The location of those high values is right after the nozzle exit near the symmetry plane. (see my attached picture)

I also tried to make my simulation more stable by adding artificial viscosity, but it made it even worse...

Does anyone of you have an hint, what could be wrong?
I would be really grateful if someone could help me.

Really quickly, with what you explain, I would expect your time step to go so low because there is divergence in your equations that lead velocity or pressure (or both) to go sky high. The problem can therefore be anywhere other then simply the time-scheme.

Are you sure your mesh is ok?
Are your boundary conditions appropriate?
Are they implemented correctly?
Is your initial condition valid?

There is plenty of possibility of failure for such a CFD code, and one can not simply help you without more details.

Congrats for writing your own code, it's a great way of learning!

As Alex has pointed out there could be many places where the code could go wrong.

These points I would suggest whiled debugging

1) Give equal pressures at inlet and exit, with which you can see where the mistakes is.
This condition should ensure there is no flow as the pressrues are equal, if some quantity is increasing then there is a bug in that area.

2) There should not be any problem with supersonic boundary conditions as you copy (extrapolate) the data .


This might help you

again welcome to CFD coding domain and the world of debugging

Regards

Thank you for the tip!

I will definitely try to put the same pressure at inlet and outlet.

I will let you know if it helped.

Cheers!

Hy guys, I'm back ;)

So I checked my code by putting the input static pressure equal to the output static pressure, furthermore I put the inlet velocity to zero. With this setting my flow field doesn't change, hence the velocity at every node keeps zero and all the other values correspond to the surrounding values.

When only setting the inlet pressure to zero, but not the inlet velocity
something changes...

I also checked my code of the maccormack technique by calculating the first iteration by hand. So my implementation of this scheme should be alright.

I still have the same problem, that the velocity becomes to high and hence the value of the time steps decreases rapidly.

Now I have a question concerning the initial conditions. Is it okay to take the surrounding conditions as an initial condition at every node, except of course those at the nozzle exit?

Thanks for your help!

The answer really depends on what you are trying to calculate. If you are looking for the stationary solution (flow is either steady-state, or oscillating periodically) I strongly suggest you initialize your flow to a physical solution that respect your boundaries and that you slowly change your boundaries until achieved the condition you want to simulate. Doing so, you avoid having great jump in flow variables.

The easiest way to do that in my opinion, is to initialize the flow at zero speed, with 0 speed at your inlet. You already tested that condition, so it is working. Then each time-step, you can increase the velocity at the inlet by a small amount, until reaching the mach number you want to simulate.

If you initialize the flow at the velocity inlet, at the very first iteration, the correction for velocity components near your wall are going to be of great magnitude in order to achieve the no-slip condition.

The other way around might be to initialize with at full velocity with slip wall, and to increase gradually a slip ratio... But I never tested that, and it does not make a lot of physical sense.

Thank you Alex, your approach seems reasonable.
I set the initial Mach number at the nozzle exit to zero and increased it by a small amount each iteration, at the same time I calculate the new values of the pressure, temperature and density at the nozzle exit, since they depend on the Mach number.

But it got even worse with this technique, my time step decreases much more rapidly than before...

I dont't know what else I should try to make it work.

I enclosed a picture of the density distribution after 1000 time steps and a Courant number of 0.1. The total expired time is only 3.34e-4. My calculations always looked similar to that one.

Why don't check your code first in a more controlled situation? Does it work in a plane channel? Does it work is smooth convergent-divergent nozzle?

Checking the code in a more controlled situation is not so easy, then I would have to use a transformed grid, which would complicate it even more...
I thought about maybe simulating the supersonic flow over a ramp and check if the compression shock corresponds to the analytical solution.

Now I have a question about my nozzle exit boundary condition. As Alex suggested it, I increase the values at the nozzle exit until they reach the values that I want. This also means that for the first iterations my inflow will be subsonic. That implies, that I need another boundary condition. For the supersonic case I can fix all values at the nozzle exit, since two characteristics plus the streamlin go inside the domain, which means I can fix three values. For the subsonic case this means I should only fix two values.

Am I right with my assumption?
If yes, how should I achieve that? Seems to be really chaotic...

Your solution just look does not steady yet. It's normal for 1000 iteration with CFL=0.1

Yes that's true, my solution hasn't reached a steady state yet. But the problem is, that when I do more iterations it starts to diverge...

I got some good advice. It is likely that my problem is numerical instability, since my values of the velocity for example become NAN at high iterations. So maybe some rearrangement of my equations will help.

I also tried to apply my code to a more easy geometry, a pipe, like FMDenaro suggested. But also in this simple case it diverges.

Your geometry does not seem to me particularly difficult, rather it did not accurately represents your original problem with the nozzle. In addition usually in 90% the problems are associated with the incorrect boundary conditions.

PS: Prof. Denaro is absolutely right when he says that you need to test the new code to a simple problem with a known analytical solution.

I tried to make a quick analysis of your problem in the proposed geometry. I really do not know the parameters that you specified on the BC 1) so there I set the speed of M = 2.5 and the static pressure same as ambient pressure. In this formulation, for McCormack really have a hard work :) ( see video longitudal velocity and density) For correlation with your task CFL=0.1, snapshots make each 1000 iterations

Hey SergeAS,

can you tell me which boundary condition you applied to BC1. Only fixed Mach number? And what about the other values, do you simply extrapolate them?

I fixed the values of the Mach number, the density and the temperature at BC1. I thought because the flow at BC1 is already supersonic I can do it like that.

BC1 is Dirichlet BC (all parameters fixed)

PS: According at the indicated place of occurrence of numerical instability and the character transient of solutions you are likely to have faced with the carbuncle phenomenon which is typical for the chosen numerical scheme. In your place, I would have used more recent shock capturing schemes. You can certainly try to reduce the Courant number (of course under the condition that your numerical scheme is implemented without errors).

For validation your code at first I would recommend two classical problems with shock capturing:

- Sod shock tube problem (1D)
- Oblique shock problem (2D)

Hy,

I would like to check my code for example with the 2D oblique shock problem, but I think this would just produce even more workload, without helping me with my problem...

The reason why I used the Maccormack technique was, that I have no knowledge of other techniques and it seems to be very easy. The implementation of the technique and the boundary conditions was no problem.
I'm just struggling with the instability.

The values of internal energy become zero after some iterations, always near my corner (BC1 to BC2). I tried to interpolate the value at the corner from the neighboring nodes, I also tried to change the distribution of the nozzles exit velocity to a tanh distribution, so that the value of velocity at the corner is zero.

Maybe you can help me with that corner.

What type of BC you set on the BC2? No-slip wall ?

In the corner points (where there are crossing different BC), I would recommend to set that type of BC (of 2), which imposes more restrictions, such as Dirichlet BC in this particular case. This usually increases the stability of the solutions near these points

Tends enthalpy to zero usually indicates a problem of correct resolution the rarefaction wave. In this case I would recommend to test your code on "Sod shock tube" problem.

Hy,

concerning the corner, the boundary condition that I apply to the node at the corner is a free-slip BC.

I managed to solve the problem with the instability near the corner. Now I add artificial viscosity, which made my code much more stable.

But it still diverges at very high numbers of iterations.

The most obvious problem is, that after a coupe of thousand iterations the velocity component in x direction becomes negative all over the mesh. With the time the velocity becomes even more negative. This is specially the case for the overexpanded nozzle.
The underexpanded case still doesn't work...

Then I have a question regarding the Sod shock tube. Should I hold the boundary values at the left and right end of the tube fixed over time? What will this test tell me? Only, how accurate my code is?

How changes the pressure in your computation domain?
At the limit, velocity of the gas must be directed away from areas of high pressure to low pressure. Otherwise you possible mess up signs in your code.

This test shows how your numerical scheme is able to resolve discontinous in the form of waves of compression and rarefaction.

In addition to this test, you can also find some useful tests (Matlab code) in the appendix to the book of Pieter Wesseling "Principles of Computional Fluid Dynamics".

Hy,

the direction of my velocity is from higher to lower pressure.

I did the simulation for two different outlet (not the nozzle exit!) boundary conditions.
One where I simply extrapolated all values from the interior. In this case my pressure at the outlet is the biggest (one order higher than the specified pressure for the surroundings!) and decreases in the direction to the nozzle exit, which means, that I have everywhere negative velocities.

The other condition was with a pressure outlet BC, where I have the pressure of the surroundings at the outlet. In this case the pressure increases from the outlet in direction to the nozzle.

Another weird thing is, that sometimes the pressure between two cells suddenly becomes negative(!) or has a much lower value so that the sign of the velocity at that node is different from the others.

Another issue is, that the density becomes negative(!) after a couple of thousand iterations.

1) In your case outlet BC should be d[rho,U,V,h]/dx=0
2) Negative values ​​of pressure, temperature and density suggest that your solution is unstable, In this case for chosen numerical method I would suggest to reduce the time step or CFL
...or use more recent numerical methods.

PS: extrapolation is bad way, because in this case you are replacing the physical laws of conservation of computational tricks

Hy,

alright, then I will try it with zero gradient instead of extrapolation.

Does the reduction of the CFL really make the code much more stable. I understand the reason for CFL and of course it should not be to big and in my case I see that I definitely shouldn't choose Courant numbers above 0.1. But I mean, is there a certain point where the code becomes stable?

Another question is, how does the spacing of the grid influence the stability?
My grid has a uniform spacing in x and y direction of about 5mm. Is that small enough? Or does the overall size of the domain influence the stability.

As I already stated, I use the Maccormack technique because it is pretty simple. In the 1970s they calculated stuff like the space shuttle with that, so why shouldn't it work in 2013. But maybe you could tell me some other easy to program more modern techniques. Maybe I'll give it a try :)

I really appreciate your help SergeAS!

The grid spacing will influence the stability in terms of how much will interact different BC. For example in the problems of supersonic aerodynamics bow shock wave can interact with the inlet boundary, moreover mesh size affects the resolution accuracy shocks. As to the Courant number should be noted that it is calculated on the basis of the conditions of linear stability, while the problem with the shock waves are very nonlinear. If we talk about the McCormack method, then it has both pluses and minuses. The most significant downside I just would like to mention a strong dispersion error in the high-gradient area and hence stability problems in these areas.
Personally, I do not use the McCormack method for about 15 years.

Ok SergeAS you convinced me.
I will try it with a different technique. Although my current code might work properly, I can't really check it because the calculation would take days...

Something is confusing me, I read a little bit about the modern shock capturing methods and they are always applied to equations that are discretized with the finite volume technique.
Is it necessary in my case to use the finite volume discretization? My geometry is simple and I of course use a structured grid for it.
Is there also a modern shock capturing technique that can be used with a finite difference discretization?

Could someone suggest me a good technique for unsteady hyperbolic PDEs?

Cheers!

FV discretizations starts from the weak form of the hyperbolic equations, for this reason they are used in presence of non-regular solutions (shock).
Try a reading to the book of LeVeque

Which book do you mean the book Numerical Methods for Conservation Laws?

What does weak form of the equations mean?

For the FDM I used the Euler equations in strong conservation form (this is how Anderson calls them).
Can't I just discretize this equation with the FVM for example by using a second order upwind scheme?
I enclosed how I would discretize the Euler equation in strong conservation form. Would this discretization be correct and could I solve it simply explicitly for every node?

have a reading of http://books.google.it/books/about/F...UC&redir_esc=y