Hi,

I am trying to write a Fortran 90 program to solve the problem of axisymmetric Taylor-Couette flow between concentric cylinders, in the infinite cylinder approximation. I am using a 2nd order, finite-difference, vorticity-stream function formulation in cylindrical coordinates computing the azimuthal velocity on top of circular Couette flow. I use Crank-Nicolson time-stepping on the linear terms and Adams-Bashforth on the nonlinear terms with operator factorisation.

However, I am having a problem that is frustratingly difficult to track down. As a test I compute the velocity growth rates against known results, but an increase in the radial and/or axial truncation parameters does not cause the growth rate to converge to some value as you would expect. Changing the truncation always causes order 1 differences in growth rate. In fact, it appears that the growth rate is tending to zero as I increase the truncation. However, the code always predicts the correct critical Reynolds number for the onset of Taylor vortices, no matter how good or bad the truncation, and the general shapes of all the variables are correct.

You would expect that a decrease in the truncation would cause a worsening growth rate.

The code is nonlinear but I have linearised it to help debugging, as well as trying other parameters (radius ratio, Reynolds number etc.). From the tests I have done it seems that the basic underlying diffusion equations (without nonlinear terms) as well as the Poisson equation on their own are correct. I am pretty sure there are no typos having checked the code plenty of times and the time-step is not too large.

Does anybody have any ideas as to what the problem could be here? Has anyone come across anything like this in any other problems? Any suggestions at all are welcome.

I hope I haven't been to vague in describing the above but I can post further details or parts of the code if that would help,

Thanks in advance,

Anthony.

A possible area to look at may be the boundary conditions you apply to your intermediate predicted variables. Are they second order estimates for the new time step? It is a fairly common problem area for predictor-corrector schemes which leads to a degraded or inconsistent accuracy at the boundary. Depending on the influence of the boundary, this may not become apparent except at fine resolutions.

Thanks for the reply,

All boundary conditions are 2nd order, but as far as I can tell this is not where the problem lies. With operator factorisation, there is no intermediate time step, and this was an area I have looked very closely at, especially the boundary conditions for vorticity, which depend on the stream function.

Actually, I think I may have partly tracked down the problem. It seems to be something to do with solving the Poisson equation for the stream function. The way I solve this is with the ADI method, so I only have to solve a tridiagonal system first in the radial direction, and then in the axial direction using the Thomas algorithm. Unfortunately, I have to iterate this procedure to find the solution within some desired tolerance.

Depending on the initial seed I give to the problem, (sometimes as small as 10^(-12) if I want to check growth rates) the growth rate of the solution follows a pattern: very quickly it reaches a constant value, which is what I have been taking to be the actual growth rate, but then after a certain amount of time (not within the initial transient) it becomes very erratic, looking like a numerical instability, before eventually settling down to the true growth rate. Just at this point there is a very slight change in the exponential behaviour of the solution which causes a log plot to have a steeper gradient.

Looking at the profile of where the program spends most of its time, 70% is spent in the Thomas algorithm for the Poisson equation, and this is far too much. I am now thinking that this is where the problem lies, so I will probably have to reformulate the Poisson equation as a linear system matrix*psi=vorticity and try multigrid methods or cyclic reduction which unfortunately means much larger matrices and vectors (I don't think I can use the Fourier transform because of the variable coefficients due to the cylindrical radius).

If anyone has any other ideas let me know.

Thanks again for your help,

Anthony.