Numerical Analysis of Critical Reynolds Number of Wavy Taylor Vortex Flow with Changing Aspect Ratio ()
1. Introduction
Taylor vortex flow is one of the important vortex flows that have been studied since its classic study made by G. I. Taylor in 1923 [1] . State of the flow between inner and outer cylinders of a rotating dual cylinder transits from Couette flow to Taylor vortex flow and to wavy Taylor vortex flow as the Reynolds number changes and as an aspect ratio increases [2] [3] [4] . The aspect ratio is defined to be the ratio of the gap between inner and outer cylinders to cylinder height. The Taylor vortex flow is the one where torus flows, called cells, are stacked. The wavy Taylor vortex flow is the one where every vortex of the Taylor vortex flow is subjected to the axial cyclic oscillation resulting in wavy and unstable flow over time [5] [6] .
Since Taylor’s study, Taylor vortex flow has been studied by many researchers, and the complexity of the flow has been clarified. Unsteady flow (e.g. Taylor vortex flow) causes an unstable change in the physical quantities that characterize the flow. Pacheco et al. [7] showed experimentally that in small aspect-ratio, Taylor-Couette flows have a band in the parameter space where rotating waves become steady non-axisymmetric solutions via infinite-period bifurcations. Martin et al. [8] showed that imposing axial flow in the annulus and radial flow through the cylindrical walls in a Taylor Couette system alters the stability of the flow. To analyze these unsteady flows, authors focused on quantitative values such as mean energy [9] . The kinetic energy and enstrophy for flows with different final modes are compared.
The critical Reynolds number when the flow state transits from the Taylor vortex flow to wavy Taylor vortex flow has been found by a visualization experiment of the flow when the aspect ratio is small. In the numerical analysis of the transition to the wavy Taylor vortex flow, the critical value, at which the wavy Taylor vortexes are generated by the changes of the aspect ratio, has been identified by obtaining changes in the kinetic energy due to axial velocity in every cell. However, the analysis has not yet extended until the changes in the kinetic energy converge. It should be noted that comparison between result of the experiment and of the numerical analysis is insufficient. In this study, the Taylor vortex flow in the normal mode, when the cylinder ends are fixed, is numerically analyzed based on chaos theory methodology. The purpose of the numerical analysis is to identify the critical Reynolds number when the flow state transits from the Taylor vortex flow to wavy Taylor vortex flow at various aspect ratios and at various number of cells.
2. Numerical Method
The Reynolds number Re is a dimensionless number, defined as the ratio of the inertial force in the equation of motion to the viscous force, and is given by
$\mathrm{Re}=VD/\nu $ (1)
Here, V is a representative velocity, given as the rotation speed of the inner cylinder (
${r}_{\text{in}}\omega $ ), D is the width of the gap between the inner and outer cylinders, given by the difference between their representative radii as r_{out} −r_{in}, and
$\nu $ is the kinematic viscosity of the fluid. r_{out} = 30 cm and r_{in}, = 20 cm, and the radius ratio is 0.666. The physical parameters are made dimensionless by using the gap as the representative length and the velocity of the inner cylinder as the representative velocity. The aspect ratio Γ is defined as the ratio of the cylinder length L to the gap width D and is given by
$\Gamma =L/D$ (2)
The governing equations are the unsteady three-dimensional incompressible Navier-Stokes equation with cylindrical coordinates (r, θ, z) and the continuity equation.
We use both SOR and ILUCGS methods to solve Poisson’s equation for pressure. The stress-free boundary condition was used for the upper-end wall and the stationary (non-slip) condition is used for the lower end wall. We applied Neumann conditions based on the momentum equation for pressure. As the initial condition, all velocity components are zero. Mixed solution of water and glycerin is assumed to be the working fluid, and its dynamic viscosity is 6.0 × 10^{–6} m^{2}/s. For the discretization method, we apply the QUICK method for convection terms, the second-order central difference method for the other space integration, and Euler’s method for time integration. Grids are staggered and equidistant in each direction. The number of grid points is 41 in the radial direction, and the number of grid points in the axial direction is proportionally adjusted so that it becomes 41 for the aspect ratio of 1.0. The number of grid points in the circumferential direction is 74. In order to examine the validity of the number of grid points, we analyzed Taylor vortex flow using several types of grids under various numerical conditions, and concluded that there are no differences among the modes that are finally formed, the formation of modes up to the final mode, and the manner of decay of the vortexes.
3. Mode Analysis Method
In the preceding studies on the numerical analysis, the changes in kinetic energy of axial velocity in every cell have been used to identify whether the flow state is the Taylor vortex flow or wavy Taylor vortex flow. This identification is possible because the Taylor vortex flow is defined to be the wavy Taylor vortex flow when the Taylor vortex flow becomes non-linear and unstable due to the cyclic axial oscillation. The criterion for this identification is the amplitude of the kinetic energy changes, and the flow state is identified to be the wavy Taylor vortex flow when the amplitude in the neighborhood of the estimated critical Reynolds number reaches or exceeds this criterion. When the amplitude is less than (does not reach) this criterion, the flow state is identified to be the Taylor vortex flow. The mode is decided based on the amplitude when TS is 2,000,000 (at TS = 2,000,000). The critical Reynolds number when the flow transits from the Taylor vortex flow to wavy Taylor vortex flow are decided. The computation for a relatively short period of time when TS is 2,000,000 (at TS = 2,000,000) reveals that the flow between inner and outer cylinders at around the critical Reynolds number is unsteady, and the computation is not extended to the point where the kinetic energy changes converge. Therefore, the flow state may have been decided to be the wavy Taylor vortex flow based on the kinetic energy changes even when the Taylor vortex flow is unsteady. The criterion of the energy amplitude is decided by the trend of energy amplitude changes in the neighborhood of the estimated critical Reynolds number, and thus the wavy Taylor vortex flow may exist even when the energy amplitude is less than the criterion. It will be difficult to decide the criterion that can clearly identify the Taylor vortex flow from wavy Taylor vortex flow or vice versa.
In this study, axial velocity changes in the Taylor vortex flow is expressed with an attractor in the chaos theory by taking the above into considerations. The mode of finally stabled flow and convergence of energy changes are decided respectively based on the changes and convergence of attractor orbit.
3.1. Mode Identification
The flow state of the Taylor vortex flow bifurcates from one state to another state as the rotation speed of the inner cylinder increases. The flow state, reached when the rotation speed of the inner cylinder is gradually increased while keeping the flow stable, is called the primary mode. The flow state, reached when the rotation speed is rapidly increased until the speed reaches the specific speed, is called the secondary mode. In the case where both ends of the cylinder are fixed, the primary mode corresponds to the normal mode and the secondary mode has the normal and anomalous modes. When the flow has even number of vortexes rotating from the outer cylinder towards the inner cylinder on both ends, the mode of the flow is called the normal mode. When the number of vortexes rotating from the inner cylinder towards the outer cylinder on both ends or on the one end is odd or even, the mode of the flow is called the anomalous mode. In the numerical analysis, the Reynolds number is increased from 0 (zero) to the target value at constant increase rate and at the specified number of time steps. Thus, the increase rate of Reynolds number is decided by the number of time steps, and the mode of the flow and the number of cells may bifurcate depending on the increase rate. The condition(s) of the increase rate to make the mode bifurcate is not clarified.
Therefore, the number of time steps is arbitrarily set in this study, and the vector diagram of the flow velocity across the axial cross-section between inner and outer cylinders is drawn at this number of time steps in order to see if the number of cells expected in this analysis is obtained and to see if the normal mode of the Taylor vortex flow is generated.
After checking that the number of cells expected in this analysis and the condition(s) to generate the normal mode vortexes, the numerical analysis is run and is kept running until the energy fluctuation in the cells converges. The attractor is created based on the axial velocity changes in the Taylor vortex flow which is located in the lowest layer of the gap between inner and outer cylinders. The axial velocity is analyzed at points ① through ④ in Figure 1. In this figure, MZ is the number of nodes in axial direction and C_{N} is the number of cells.
3.2. Creation of Attractor
The attractor means an “orbit on which the physical quantity in the phase space converges”, and when the attractor shows signs of remaining the constant orbit,
Figure 1. Points at which axial velocity is analyzed.
it can be deemed that changes in physical quantity have been converged. The attractor changes its shape with varying parameter. The parameter of the flow in the gap between inner and outer cylinders of the rotating dual cylinder is the Reynolds number (Re). The shape of the attractor varies with increase of Reynolds number and it follows the process described with Ruelle-Takens-Newhouse scenario shown in Figure 2. When the flow has Reynolds number less than Re0, it is the linear flow that does not fluctuate with time and the attractor converges at one point called the equilibrium point or fixed point. When the Reynolds number of the flow reaches Re0, the flow becomes non-linear which oscillates at frequency (ω_{1}), and the attractor converges into the circular orbit called the limit cycle. When the Reynolds number reaches Re1, the flow becomes the one that has two frequencies (ω_{1} and ω_{2}), and the attractor takes the form of T2 torus. When the Reynolds number exceeds Re2, the flow becomes to have one additional frequency and the attractor takes the form of T3 torus. When the Reynolds number exceeds Rec, the flow becomes turbulent, and the attractor becomes the chaotic attractor. In this state, the attractor has a wide spectrum.
In this study, the functions, W (w (t), w (t + τ), and w (t + 2τ)) that can draw the attractor in the three-dimensional coordinate system are embed as shown in Figure 3 by incorporating the delay time (τ) in the axial velocity changes (w (t)). The attractor based on the axial velocity changes is drawn as shown in Figure 4 by plotting the function (W) on the phase space of the three-dimensional coordinate system. In drawing the attractor, the delay time (τ) is set to be “τ = f/4 ” which is equal to one cycle (f) of the autocorrelation function of the axial velocity changes (w (t)) as shown in Figure 5. The autocorrelation function is derived by obtaining the first-order autocorrelation coefficient (r) calculated by the next equation. When the axial velocity does not have constant period and the delay time is unobtainable, t = 1 is set to be the delay time in creating the attractor.
Figure 2. Ruelle-Takens-Newhouse scenario.
Figure 3. Embedding of functions W (w (t), w (t + τ), and w (t + 2τ)) into three-dimensional coordinate system.
Figure 4. Attractor drawn in three-dimensional coordinate system.
$r=\frac{{\displaystyle \underset{i=1}{\overset{N-1}{\sum}}\left({x}_{i}-\stackrel{\xaf}{x}\right)\left({x}_{i\text{+}1}-\stackrel{\xaf}{x}\right)}}{{\displaystyle \underset{i=1}{\overset{N-1}{\sum}}{\left({x}_{i}-\stackrel{\xaf}{x}\right)}^{2}}}$ (3)
3.3. Mode Determination with Attractor
The mode that can finally and stably exist under the condition(s) is decided based on the created attractor. In the case where the vortex mode is the Taylor
Figure 5. Autocorrelation function (r (t)) of axial velocity (w (t)).
vortex flow, the final axial velocity is constant because the Taylor vortex flow does not fluctuate with time. Thus, the attractor converges at the fixed point. In the case where the vortex mode is the wavy Taylor vortex flow, the attractor converges in the limit cycle or in the quasi-torus orbit because the final axial velocity oscillates at one frequency or multiple frequencies. Taking these behaviors of the attractor into account, the mode that can finally exist in the flow meeting the computation conditions is decided.
In this study, the attractor is created based on the axial velocity change at every 1,000,000 of the calculation Time Step (TS) (TS = 1,000,000), and the convergence of the attractor is monitored to decide if the changes in vortex energy converge and to decide the mode. Since a huge amount of time is required in the numerical computation, the number of time steps is limited to 10,000,000 (TS = 10,000,000) and when the attractor does not converge at this limit, the mode is decided based on the changes of the attractor orbit. The attractors are color-coded with different color at every 1,000,000 (TS = 1,000,000) and they are superimposed to make it easy to check if the attractor orbit converges or not, and to monitor changes of the attractor. Figure 6 and Figure 7 indicate the attractor’s behavior respectively in the case of the Taylor vortex flow and wavy Taylor vortex flow.
When the attractor converges in the circular orbit or keeps expanding while forming the circular orbit until TS = 10,000,000, it can be considered that the attractor does not converge at the fixed point and the flow is decided to be the wavy Taylor vortex flow. When the attractor orbit has diminished, it can be considered that the attractor finally converges at the fixed point and the flow is decided to be the Taylor vortex flow. In practice, the axial velocity change of the Taylor vortex flow cannot converge finally and completely at the fixed point, but it keeps exhibiting the microscopic fluctuations. The attractor does not converge at the fixed point. That is, the flow is decided to be the Taylor vortex flow when the attractor does not draw the similar circular orbit.
3.4. Determination of Critical Reynolds Number
It is assumed that the accuracy of the Reynolds number is +/−10 at a certain aspect ratio (Γ), and the critical Reynolds number is defined to be the number when the state of vortexes becomes completely wavy (i.e., the wavy Taylor vortex flow). It is assumed that, in the case where the flow state is the Taylor vortex flow under the Re1 condition, the state transits from the Taylor vortex flow to the wavy Taylor vortex flow under the Re2 condition, in which the Reynolds number is greater than the number in Re1 condition by 20. At this point, the critical Reynolds number is Re2 (i.e., Rec = Re2).
4. Result of Numerical Analysis
This section shows the attractor in the neighborhood of the critical Reynolds number at the aspect ratio at which the critical Reynolds number is calculated in this study. The attractor is numerically analyzed at the point ① through ④ shown in Figure 1. The form of the attractor is different from one point to other points but the behavior of attractor’s change and the number of time steps to reach convergence is similar. The axial velocity component is largest at the point ② among these four points, and the attractor created at this point ② is shown in this section.
4.1. Normal Mode (Six Cells), Γ = 7.3
Figure 8 indicates the changes of the attractor at Re = 1080 for the number of time steps from 2 million to 5 million. The attractor orbit changes linearly showing very small orbit changes with the increase in the number of time steps. That is, the flow state is decided to be the Taylor vortex flow because the flow seldom changes with time. Figure 9 indicates the changes of the attractor at Re = 1100 for the number of time steps from 7 million to 10 million. It is shown that the attractor draws the same orbit and that it converges into the limit cycle. Therefore, the flow state is decided to be wavy (i.e., the wavy Taylor vortex flow). Based on the above observation, the critical Reynolds number (Rec) at the aspect ratio of 7.3 is determined to be 1100 (i.e., Rec = 1100).
4.2. Normal Mode (Six Cells), Γ = 7.6
Figure 10 indicates the attractor at Re = 1080 for the number of time steps from 2 million to 5 million. The attractor orbit changes linearly showing very small orbit changes with the increase in the number of time steps. That is, the flow state is decided to be the Taylor vortex flow because the flow seldom changes with time. Figure 11 indicates the changes of the attractor at Re = 1100 for the number of time steps from 6 million to 8 million. It is shown in this figure that the attractor does not draw the constant circler orbit but the orbit varies three-dimensionally with the change of the time. However, changes of the orbit with the time steps are insignificant and the orbit is virtually identical. Therefore, the flow state is decided to be wavy (i.e., the wavy Taylor vortex flow) because of the flow changes at the constant period. Based on the above observa
Figure 8. Γ = 7.3 and Re = 1080 (Number of time steps = 2 million to 5 million).
Figure 9. Γ = 7.3 and Re = 1100 (Number of time steps = 7 million to 10 million).
Figure 10. Γ = 7.6 and Re = 1080 (Number of time steps = 2 million to 5 million).
Figure 11. Γ = 7.6 and Re = 1100 (Number of time steps = 7 million to 10 million).
tion, the critical Reynolds number (Rec) at the aspect ratio of 7.6 is determined to be 1100 (i.e., Rec = 1100).
4.3. Normal Mode (Six Cells), Γ = 7.7
Figure 12 indicates the attractor at Re = 1040 for the number of time steps from 7 million to 10 million. The attractor orbit keeps diminishing and its change keeps decreasing as the time steps increase. That is, the flow state is decided to be the Taylor vortex flow because the flow seldom changes with time. Figure 13 indicates the attractor at Re = 1060 for the number of time steps from 7 million to 10 million. The flow state is decided to be wavy (i.e., the wavy Taylor vortex flow) because the attractor does not converge at the 10 million steps but the attractor keeps forming and expanding the circle. Based on the above observation, the critical Reynolds number (Rec) at the aspect ratio of 7.7 is determined to be 1060 (i.e., Rec = 1060).
5. Discussion
Figure 14 shows the distribution of the critical Reynolds numbers obtained by this numerical analysis and by the previous experiments. White points indicate the Reynolds numbers obtained by the numerical analysis, and black points do the Reynolds numbers obtained by the experiments. Circle, rectangle, triangle, and X show 2 cells, 4 cells, 6 cells, and 8 cells, respectively. This numerical analysis for 6 cells is carried out respectively at the aspect ratio between 7.3 and 7.7 (i.e., Γ = 7.3-7.7), and the analysis for 8 cells is carried out respectively at the aspect ratio between 6.0 and 6.2 (i.e., Γ = 6.0-6.2). The critical Reynolds number is calculated for 6 cells at the aspect ratio of 7.3, 7.6, and 7.7 (i.e., Γ = 7.3, 7.6, and 7.7). Comparison of distribution of the calculated and experimental values reveals that, at each cell, the critical Reynolds numbers are almost identical from the Reynolds number at small aspect ratio to the vicinity of the peak of Reynolds number. Even in the case of 8 cells, the critical Reynolds number that would be
Figure 12. Γ = 7.7 and Re = 1040 (Number of time steps = 5 million to 8 million).
Figure 13. Γ = 7.7 and Re = 1060 (Number of time steps = 6 million to 9 million).
Figure 14. Distribution of critical Reynolds number when the flow state transits from normal mode Taylor vortex flow to wavy Taylor vortex flow in rotating dual cylinder whose both ends are fixed.
obtained by the numerical analysis is expected to be identical to the Reynolds number obtained by the experiment, from the small aspect ratio to the vicinity of peak of Reynolds number. At any cells, the distribution of the Reynolds number after their peak is qualitatively similar between the calculated and experimental values. But, the experimental stability limit of the Taylor vortex flow is lower than the calculated stability limit. The difference between the experimental and calculated stability limit is significant at around the peak of the critical Reynolds number, and the aspect ratio at which the calculated value hits its peak is larger than the aspect ratio at which the experimental value hits its peak. It is also shown that the difference between the calculated and experimental values increases as the number of cells increases. This difference will be attributable to the pattern of changes in cell boundary plane of the wavy Taylor vortex flow at the aspect ratio after the peak and to the relaxation time of the flow used in the visualization experiment.
In the case of the Taylor vortex flow, the boundary plane between cells will finally become horizontal at a constant height. Since the boundary plane of the wavy Taylor vortex flow is usually wavy, it will be possible to distinguish the wavy Taylor vortex flow from the Taylor vortex flow by observing the cross-section of vortexes in axial direction. However, it is confirmed with the visualization experiment that the fluctuation of the boundary plane of the wavy Taylor vortex flow can have two patterns (horizontal and wavy patterns) after the peak. When the boundary plane is horizontal, it does not stay at the constant height but oscillates up and down while keeping the horizontal line. In other words, the boundary plane is wavy due to the repeated expansion and contraction of the entire cell. Normally, the cell boundary plane of the Taylor vortex flow is horizontal. When the Taylor vortex flow is unsteady, the boundary plane may have similar oscillation because the plane keeps oscillating periodically. It is presumable that, in the experiment, the wavy Taylor vortex flow and unsteady Taylor vortex flow may have been confused making the experimental critical Reynolds number lower. However, it is not clear if the periodic component of the unsteady Taylor vortex flow makes the cell boundary oscillate or not. The critical Reynolds numbers are not calculated in the neighborhood of the experimental values located in the region where the differences between the experimental and calculated values are observed, and thus the actual changes in the flow are not known. Therefore, it is necessary to review the experiment method and to redo the numerical analysis in the region outside of the neighborhood of the critical Reynolds number in order to identify the factor(s) contributing these differences.
6. Conclusions
This study has identified the critical Reynolds number when the state of the flow whirling between the inner and outer cylinders of the rotating dual cylinder transits from the Taylor vortex flow to wavy Taylor vortex flow. The numerical analysis making use of the attractor in the chaos theory has been used in this identification of the critical Reynolds number.
The calculated critical Reynolds numbers at various number of cells are almost identical to the values obtained by the visualization experiment in the region between the Reynolds number at small aspect ratio and the vicinity of the peak of Reynolds number. In the region where the aspect ratio is larger than the ratio at the peak critical Reynolds number, the distribution of the Reynolds number is qualitatively similar between the calculated and experimental values. But, the experimental stability limit of the Taylor vortex flow is lower than the calculated stability limit. The difference between the experimental and calculated stability limit is significant at around the peak of the critical Reynolds number, and the aspect ratio at which the calculated value hits its peak is larger than the aspect ratio at which the experimental value hits its peak. This difference is assumed to be attributable to the pattern of changes in cell boundary plane of the wavy Taylor vortex flow at the aspect ratio after the peak and to the relaxation time of the flow. It is necessary to review the experiment method and conditions of the numerical analysis in order to verify this assumption.
Acknowledgements
We would like to show our gratitude to members in our laboratory for sharing their pearls of wisdom with us during the course of this research, and we thank “anonymous” reviewers for helpful comments that greatly improved the manuscript.