CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > Main CFD Forum

Time-marching scheme diverges

Register Blogs Members List Search Today's Posts Mark Forums Read

Reply
 
LinkBack Thread Tools Display Modes
Old   September 18, 2013, 15:42
Default Time-marching scheme diverges
  #1
New Member
 
Obad
Join Date: Sep 2013
Posts: 29
Rep Power: 3
Obad is on a distinguished road
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.
Attached Images
File Type: jpg Skizze Randbedingungen.jpg (111.4 KB, 30 views)
Obad is offline   Reply With Quote

Old   September 18, 2013, 17:21
Default
  #2
Member
 
Join Date: Jul 2013
Posts: 37
Rep Power: 3
Alex C. is on a distinguished road
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!
Alex C. is offline   Reply With Quote

Old   September 21, 2013, 22:36
Default
  #3
Member
 
Ramesh K
Join Date: Dec 2009
Location: Bangalore
Posts: 65
Rep Power: 7
RameshK is on a distinguished road
Send a message via Yahoo to RameshK
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
RameshK is offline   Reply With Quote

Old   September 22, 2013, 18:11
Default
  #4
New Member
 
Obad
Join Date: Sep 2013
Posts: 29
Rep Power: 3
Obad is on a distinguished road
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!
Obad is offline   Reply With Quote

Old   September 27, 2013, 07:21
Default
  #5
New Member
 
Obad
Join Date: Sep 2013
Posts: 29
Rep Power: 3
Obad is on a distinguished road
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!
Obad is offline   Reply With Quote

Old   September 27, 2013, 09:35
Default
  #6
Member
 
Join Date: Jul 2013
Posts: 37
Rep Power: 3
Alex C. is on a distinguished road
Quote:
Originally Posted by Obad View Post
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.
Alex C. is offline   Reply With Quote

Old   September 28, 2013, 08:40
Default
  #7
New Member
 
Obad
Join Date: Sep 2013
Posts: 29
Rep Power: 3
Obad is on a distinguished road
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.
Attached Images
File Type: jpg density_distribution_1000_0.jpg (23.3 KB, 24 views)
Obad is offline   Reply With Quote

Old   September 30, 2013, 08:03
Default
  #8
Senior Member
 
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 1,568
Rep Power: 20
FMDenaro will become famous soon enough
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?
FMDenaro is offline   Reply With Quote

Old   October 2, 2013, 16:33
Default
  #9
New Member
 
Obad
Join Date: Sep 2013
Posts: 29
Rep Power: 3
Obad is on a distinguished road
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...
Obad is offline   Reply With Quote

Old   October 4, 2013, 11:06
Default
  #10
Member
 
SergeAS's Avatar
 
Serge A. Suchkov
Join Date: Oct 2011
Location: Moscow, Russia
Posts: 73
Blog Entries: 5
Rep Power: 5
SergeAS is on a distinguished road
Send a message via Skype™ to SergeAS
Your solution just look does not steady yet. It's normal for 1000 iteration with CFL=0.1
SergeAS is offline   Reply With Quote

Old   October 5, 2013, 13:11
Default
  #11
New Member
 
Obad
Join Date: Sep 2013
Posts: 29
Rep Power: 3
Obad is on a distinguished road
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.
Obad is offline   Reply With Quote

Old   October 5, 2013, 15:25
Default
  #12
Member
 
SergeAS's Avatar
 
Serge A. Suchkov
Join Date: Oct 2011
Location: Moscow, Russia
Posts: 73
Blog Entries: 5
Rep Power: 5
SergeAS is on a distinguished road
Send a message via Skype™ to SergeAS
Quote:
Originally Posted by Obad View Post
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.
SergeAS is offline   Reply With Quote

Old   October 5, 2013, 17:25
Default
  #13
Member
 
SergeAS's Avatar
 
Serge A. Suchkov
Join Date: Oct 2011
Location: Moscow, Russia
Posts: 73
Blog Entries: 5
Rep Power: 5
SergeAS is on a distinguished road
Send a message via Skype™ to SergeAS
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
SergeAS is offline   Reply With Quote

Old   October 9, 2013, 10:49
Default
  #14
New Member
 
Obad
Join Date: Sep 2013
Posts: 29
Rep Power: 3
Obad is on a distinguished road
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.
Obad is offline   Reply With Quote

Old   October 9, 2013, 11:24
Default
  #15
Member
 
SergeAS's Avatar
 
Serge A. Suchkov
Join Date: Oct 2011
Location: Moscow, Russia
Posts: 73
Blog Entries: 5
Rep Power: 5
SergeAS is on a distinguished road
Send a message via Skype™ to SergeAS
Quote:
Originally Posted by Obad View Post
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)
SergeAS is offline   Reply With Quote

Old   October 9, 2013, 11:58
Default
  #16
New Member
 
Obad
Join Date: Sep 2013
Posts: 29
Rep Power: 3
Obad is on a distinguished road
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.
Obad is offline   Reply With Quote

Old   October 9, 2013, 15:30
Default
  #17
Member
 
SergeAS's Avatar
 
Serge A. Suchkov
Join Date: Oct 2011
Location: Moscow, Russia
Posts: 73
Blog Entries: 5
Rep Power: 5
SergeAS is on a distinguished road
Send a message via Skype™ to SergeAS
Quote:
Originally Posted by Obad View Post
Hy,

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.
SergeAS is offline   Reply With Quote

Old   October 12, 2013, 12:03
Default
  #18
New Member
 
Obad
Join Date: Sep 2013
Posts: 29
Rep Power: 3
Obad is on a distinguished road
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?
Obad is offline   Reply With Quote

Old   October 12, 2013, 14:53
Default
  #19
Member
 
SergeAS's Avatar
 
Serge A. Suchkov
Join Date: Oct 2011
Location: Moscow, Russia
Posts: 73
Blog Entries: 5
Rep Power: 5
SergeAS is on a distinguished road
Send a message via Skype™ to SergeAS
Quote:
Originally Posted by Obad View Post
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.
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.

Quote:
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?
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".
SergeAS is offline   Reply With Quote

Old   October 12, 2013, 15:44
Default
  #20
New Member
 
Obad
Join Date: Sep 2013
Posts: 29
Rep Power: 3
Obad is on a distinguished road
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.
Obad is offline   Reply With Quote

Reply

Thread Tools
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are On
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
same geometry,structured and unstructured mesh,different behaviour. sharonyue OpenFOAM Running, Solving & CFD 13 January 2, 2013 23:40
High Courant Number @ icoFoam Artex85 OpenFOAM Running, Solving & CFD 9 January 3, 2012 09:06
Orifice Plate with a fully developed flow - Problems with convergence jonmec OpenFOAM Running, Solving & CFD 3 July 28, 2011 05:24
Convergence moving mesh lr103476 OpenFOAM Running, Solving & CFD 30 November 19, 2007 15:09
Modeling in micron scale using icoFoam m9819348 OpenFOAM Running, Solving & CFD 7 October 27, 2007 00:36


All times are GMT -4. The time now is 22:56.