Timemarching scheme diverges
1 Attachment(s)
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 timescheme.
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! 
Quote:
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 timestep, 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 noslip 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. 
1 Attachment(s)
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.34e4. 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 convergentdivergent 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. 
Quote:
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. 
Quote:
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. 
Quote:
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 freeslip 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? 
Quote:
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. Quote:
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. 
All times are GMT 4. The time now is 18:08. 