Here is an answer I found on CFD-Online.

"

(1). If you evaluate the spatial terms at the known old time, and express the first order time derivative term by a two point difference between the future time and the old time, then it is explicit method, because the future value of the dependent variable can be computed directly (explicitly) from the known old values. In other words, the future values (and the flow field) depend solely on the old known values (and the old flow field). So, you can march in time this way, easily.

But if you take a large time step, you can easily step into a sink hole or man hole. The gradient in time in the explicit method is evaluated on the right hand side of the equation by the old time spatial variation. This is a real model of the real physics, except that in real world, the time step is continuous and infinitestimal. So, the finite time step in explicit method will deviate from the true path and produce inaccurate answers.

(2). In the implicit method, your are pretending that the spatial derivative term lives in the future time, which is not a right theory either. Because the future flow field is completely unknown, you are in a bigger trouble with this theory. The reason why you have a better chance to get the flow field is that, there exist "boundary conditions" in the future time. The solution theory becomes: how do we guess at the future flow field such that the future spatial derivative will equal to the two point time derivative, bounded by the future B.C's.

You adjust (iterate) the future flow field until such condition is satisfied (backward requirement, implicit method, or iterative). Since you are free to adjust the future flow variable (to absorb or re-distribute the errors), you have a better chance to obtain a good solution.(in explicit method, you are not allowed to adjust anything at all)

The problem is if the time step is large, and the future boundary condition has changed a lot, then you will be solving a different problem, linking the future and the old time directly by a two point difference, separated by a long time period. This will fail, because the boundary conditions in the future is no longer related to the old time.

That is why implicit methods with large time-steps are only advised for problems that approach steady state (i.e., fixed boundary conditions).

"

I think what they mean is, that for time-varying boundary conditions, even implicit solvers have a limit to the time-step.

Are you using some turbulence  model? Fine mesh  can capture flow natural instabilities that the coarse mesh had filtered. And the solver  may  not manege this new instabilities. 

The boundary conditions suggest the problem is time-dependent, the initial success then again that you have a Courant number limitation. These two may be linked. Also, which integration do you use? Implicit Euler should give enough diffusion to smear almost anything, but Crank-Nicolson is just on the stability limit. A good option would be a three level scheme.

Thanks Mr Ville. The problem in determination of courant number lies in finding velocity. As the problem is of oscillating nature, mass flow rate and velocity changes accordingly. Also mass flow flow rate keep on increasing as the cooling produces in cooler. Which velocity should I take to determine courant number? Can you give some more details as I stuck in this problem since two months. Thanks in advance.

Courant number is directly proportional to the velocity. Since you intend to find the maximum value of the courant number in the simulation, it would probably be best to use the maximum value of velocity that occurs in the entire domain. This is assuming the grid is uniform.

For a non uniform grid, the maximum value of u/(dx^2) would give the point where the courant number is maximum. Am I thinking along the right lines, Mr. Viitanen and Dr. Ricardo?

Sincerely, Saumitra, I do not know. it can be...But, as the CFL limit C>1 so the flow turns turbulent. This can be a possibility and the cause of divergence. 

Other is your solver: the P-v coupling  scheme it is important  and for  point (or line) solver method as Gauss-Seidel you need that this condition be respected:

Sum Aneighbors/Ap

How about using a local time stepping in the inner loop? That is, define the time pseudo time step from the local physics. In addition, reformulating the Courant condition such that the mass flux appears there explicitly may then provide enough relaxation. The aforementioned maximum value by Saumitra should solve the issue, but is likely uneconomical.  Also, you can try pseudolinearising the equations so that the diagonal has clearly more weight. The requirement Ricardo gave is the Scarborough criterion for a Gauss-Seidel method to converge.  Note that the instability issue may well have to do with the (physical) problem itself. Then there is no other option but to use sufficient resolution.

If you have a strong convective term, you should think about the discretization scheme. The convective term is a common source of numerical instabilities and select the right scheme, a second order upwind scheme for instance, usually improves considerably the numerical stability

Dear Paresh, what kind of non convergence do you observe? It could be possible that, due to the kind of  finite volume formulation, a hourglass-type instability occurs. It appears with a chessboard like pressure field and is associated to zero-energy spurious modes that amplifies during the time-marching solution. Let me kmow. Good luck!

Hello Giuseppe Vairo Sir.

I am using following details with respect to my code. System consists of two main components.

Regenerator and Pulse Tube. Length of Pulse tube: 230 mm and Regenerator Length: 210 mm. I am doing 1-D coding with grid size of regenerator as 500 and that of pulse tube equl to 350. Thus, dx_regenerator = 0.000420, dx_pulsetube = 0.000657, time step size: 0.0125. Method used: Gauss Siedel with under-relaxation factor of 0.20. 

You can observe this by viewing following video, how I am getting repeated values of pressure and temperature at a point (cold end of pulse tube, a point of cooling). Video shows Pn = New pressure at cold end. I am using Staggered grid for velocity/mass flow rate.  

Ok Peresh, but have you a pressure field distribution at a certain time step? If yes, you can check if the instantaneous pressure distribution, when the solution is deteriorate, experiences a chessboard-like patterns. If this is the matter, possible solutions can be found by considering different strategies for treting convective terms. Let me kmow. Bye

time step does not matter in a mathematical sense, it does matter in a physical sense because if you use too large a time step you may miss important phenomena that happen in shorter physical time , particularly as  you have  oscillating p and T. Try using relaxation factors that are proportional to the shortest residence time ,I.e max.expected velocity and min cell size.

the rms  seems 'freeze' . Well, it is very common with this method.  If is not concerning a coding errors, It is regarding the error frequencies. ( see Non Newmann stability analyses:  http://en.wikipedia.org/wiki/Von_Neumann_stability_analysis). To circumvent this issue,  the multigrid method are applied. Coarse meshes have a fast decrement of residues. The multigrid approach acts with a 'virtual coarsening' , or in  other words, only in mathematical sense. Theres is no geometric refinement or coarsening.  

This issues is characteristic of Jacobi, Gauss Seidel and SOR methods. Maybe its all right with your code,  but the solver approach is not robust to hadle with advective dominance, or  complex source terms. In my first codes, I've experienced the same thing. 

http://en.wikipedia.org/wiki/Von_Neumann_stability_analysis

I also think like 'Markatos'.

Why don't you change the solver. Try to use the TDMA solver (1D,2D, 3D) or SIP solver(2D or 3D).

given that your problem is oscilating in nature, the time step is limited to at most one tenth of a period of time or even smaller. It's nothing to do with implicit or explicit.

Some implicit scheme can be conditionally stable. The behaviour  of your scheme shows instability. Either choose small time step or change the scheme. Certainly advection in the equation is the source of instability.

Actually I tried to reduce the time step size (dt) for a given value of grid size (dx) so that my CFL should be lower than 1. In my case mass flow rate at steady state is far more (due to cooling in cooler), than initial value. Thus CFL increases as the simulation proceeds. So, I need to keep CFL as low as possible so that CFL at steady state cannot exceed more than 1. But dt value cannot be lowered to certain value for a given grid size and When I try to lower it (say, 5e-5 sec for 100x100 grid), it gives me SEGMENTATION FAULT (CORE DUMPED) error. I am defining all the variables/data types as 'double'

Please share information if you know more about this. 

Forget CFL , as you claim to solve implicitly. I am certain you have a mistake in your boundary conditions , especially at the inlet mass flow rate.What do u mean  "due to cooler"?

In pulse tube 'cooler', the mass flow rate keep on increasing due to cooling process occurs at cold end of pulse tube. This is why velocity keep on increasing. 

Of course. So at every time step you update the mass flow rate.Don't you?

I would recommend the Hestenes-Stiefel or conjugate gradient algorithm over the Gauss-Seidel. I have successfully used both, but with the CG algorithm you can do some clever updating of the other parameters. There are quite a few scholarly articles on applications of the CG method.  It's almost always faster too.

Yes. Mass flow rate gets increased and so CFL also increases. Today I observed that mass flow when increased beyond  9 gram per seconds, code gets instability and CFL also oscillates between two values. 

Forget CFL and increase relaxation (or reduce time step) and let me know.

R u using relaxation of the false-time -step type  or the conventional one?

Properties are changing during the iteration. If they are not properly controlled the simulation will blow up quickly. Use under relaxation to control the fluid properties. Better to use direct solver for solving the matrix too.

Recently I found that when I neglect pressure drop by not considering dp/dx i.e. momentum equation, I am getting convergence till final steady state. That means, Instability lies in momentum equation.

Please look at the image (attached herewith). In the equation, I have to provide pressure boundary condition at the left/warm end of regenerator. The sinusoidal pressure waveform is applied at the left/warm end (by rotary valve).

Now, I want to apply this information as boundary in the governing equation (6) shown in the image. What I did in present code is, 1st node pressure is calculated from boundary pressure (sinusoidal) by use of governing equation (6)

i.e. Pnew[1]=Pboundary - 2*(0.5*dx)*fF*mf1*mf1/(Dh*rho*A2) where, mf1=mass flow at boundary. And pressure at subsequent nodes are calculated from Pnew[1].

When I implement this as boundary condition, It works for some grid size but not for all. When code gets instable, pressure value at given node, oscillates between two fixed values

What could be the error in boundary implementation?

It is dangerous to prescribe pressure (it is a relative variable, except if you are compressible, and can adjust at any level as it is the p difference that matters not its actual value). Use as BC the known mass flow rate ( for the continuity eqn) and the velocities (for the momentum eqns).

If you insist to use P , fix it at a reference point so all calculated p"s are relative to that one and cannot oscillate. Or use a penalty function and prescribe P OUTSIDE your domain , but this is not easily implemented

Are your physical properties calculated as function of Pressure? According to Eq  7-8 the pressure and mass flow can change the signal as function of dx. this can arise non physical properties values as negative densities, viscosities etc. Look for strange values in the physical properties. 

Yes, All thermo-physical properties are calculated as function of pressure and temperature both (With polynomial regression equation of 7th).

You may fix the Courant number to some value less than 1. This will calculate your time-step automatically according to your velocity. Then, the simulation can be stable.

Look for  unphysical values . ANSYS CFX  uses  a EOS but generate a table with these properties  limited for some values in pressure and temperature and interpolate with a surface response. Sometimes, during the solver iterations, pressure and/or temperature extrapolate the initial data and the interpolation  is made, but the program  warning  to avoid that and recommends  change the limits  of pressure and temperature if the situation persists. If  this happens generally   a divergence or overflow error happens and the solver fail. 

If this can happen with a code as ANSYS CFX, it can happen with your also. Try to monitor these properties and see if they are Ok.

Grid and time step are linked through the Courant condition. You cannot change one without changing the other.