Yang Numerical Simulation of Multiplicity and Stability of Mixed Convection in Rotating Curved Ducts

A numerical study is made on the fully developed bifurcation structure and stability of the mixed convection in rotating curved ducts of square cross-section with the emphasis on the effect of buoyancy force. The rotation can be positive or negative. The fluid can be heated or cooled. The study reveals the rich solution and flow structures and complicated stability features. One symmetric and two symmetric/asymmetric solution branches are found with seventy five limit points and fourteen bifurcation points. The flows on these branches can be symmetric, asymmetric, 2-cell, and up to 14-cell structures. Dynamic responses of the multiple solutions to finite random disturbances are examined by the direct transient computation. It is found that possible physically realizable fully developed flows evolve, as the variation of buoyancy force, from a stable steady multicell state at a large buoyancy force of cooling to the coexistence of three stable steady multicell states, a temporal periodic oscillation state, the coexistence of periodic oscillation and chaotic oscillation, a chaotic temporal oscillation, a subharmonic-bifurcation-driven asymmetric oscillating state, and a stable steady 2-cell state at large buoyancy force of heating.


INTRODUCTION
We study the fully developed bifurcation-driven multiplicity and dynamic responses of multiple solutions to finite random disturbances numerically by the finite-volume/Euler-Newton continuation and the direct transient computation for the mixed convection in ducts of square crosssection with the streamwise curvature, the spanwise rotation in either positive or negative direction, and the wall heating/cooling (Figure 1 with (R, Z, ϕ) as the radial, spanwise, and streamwise directions, resp.).A positive rotation gives rises to a Coriolis force in the cross-plane (RZ-plane) directed along positive R-direction and vice versa.
There is a host of areas where such flows and transport phenomena are of practical importance and where relevant issues are raised.For instance, sedimentation field-flow fractionation, one of the most powerful subtechniques of the field-flow fractionation (FFF) method of separation, uses rotating curved ducts to perform both separation and quantitation of particles and other colloidal-based substances.
This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
Aerosol centrifuges apply rotating curved ducts to the problem of separating airborne particles according to aerodynamic size.The fundamental economic gains associated with the increase of inlet fluid temperature and unit capacity lead to the use of rotating curved ducts as the cooling ducts of rotating power machinery such as gas/steam turbines and electric generators.The rotating curved ducts are also used in rotating heat exchangers, centrifugal material processing and material quality control, medical and chromatographic devices, and so forth.In order to calculate the pumping power needed for such devices, it is important to know the pressure drop in rotating curved ducts.Because secondary flows can enhance heat and mass transfer, knowledge of the magnitude of this effect in different ranges of operating parameters is important in designing and operating these devices.To avoid or reduce the flow-induced vibration and noise, we need to know when temporal oscillation appears.The present work can acquire a better understanding of these practical issues.
Early works on the rotating curved duct flows were constrained to two simplified limiting cases with strong or weak rotations.Ludwieg [1] developed a solution based on a momentum integral method for the isothermal flow in a square duct with a strong spanwise rotation.Miyazaki [2,3] examined the mixed convection in a curved circular/rectangular duct with spanwise rotation and wall heating by a finitedifference method.Because of the convergence difficulties with the iterative method used, Miyazaki's work was constrained to the case of weak curvature, rotation, and heating rate.As well, all the works employ a steady model for the fully developed laminar flow with a positive rotation of the duct.Since the solution is only for the asymptotic cases, the secondary flow revealed by these early works consists of only one pair of counter-rotating vortices in the cross-plane.The interaction of the secondary flow with the pressure-driven streamwise flow shifts the location of the maximum streamwise velocity away from the center of the duct and in the direction of the secondary velocity in the middle of the duct.
More comprehensive studies have been made in recent years by Wang and Cheng [4], and Daskopoulos and Lenhoff [9] for a circular tube; Matsson and Alfredsson [10,11], and Guo and Finlay [12] for a high-aspect-ratio rectangular duct; and Wang and Cheng [5,13,14,15], Wang [6,7,8], Selmi et al. [17], and Selmi and Nandakumar [16] for the square and rectangular ducts with a low-aspect ratio.All the works are for the steady fully developed flows.Wang and Cheng [4] developed an analytical solution for rotating curved flow with the effect of heating or cooling which allows to analyze the solution structure.Detailed flow structures and heat transfer characteristics were examined numerically by Wang and Cheng [5] and Wang [6,7,8].The rotating curved flows were visualized using smoke injection method by Wang and Cheng [13,14,15].Daskopoulos and Lenhoff [9] made the first bifurcation study numerically under the small curvature and the symmetry condition imposed along the tube's horizontal central plane.Matsson and Alfredsson [10] presented the first and comprehensive linear stability analysis.Matsson and Alfredsson [11] reported an experimental study, by hotwire measurements and smoke visualization, of the effect of rotation on both primary and secondary instabilities.Using a linear stability theory and spectral method, Guo and Finlay [12] examined the stability of streamwise-oriented vortices to 2D, spanwise-periodic disturbances (Eckhaus stability).Detailed bifurcation structure and linear stability of solutions were determined numerically by Selmi et al. [17] and Selmi and Nandakumar [16] without imposing the symmetric boundary conditions.
It is the relative motion between bodies that determines the performances such as friction and heat transfer characteristics.The duct rotation introduces both centrifugal and Coriolis forces in the momentum equation describing the relative motion of fluids with respect to the duct.For isothermal flows of a constant property fluid, the Coriolis force tends to produce vorticity while the centrifugal force is purely hydrostatic, analogous to the Earth's gravitational field [18].When a temperature-induced variation of fluid density occurs for nonisothermal flows, both Coriolis-and centrifugaltype buoyancy forces could contribute to the generation of the vorticity [18].These two effects of rotation either enhance or counteract each other in a nonlinear manner depending on the direction of duct rotation, the direction of wall heat flux and the flow domain.As well, the buoyancy force is proportional to the square of the rotation speed while the Coriolis force increases proportionally with the rotation speed itself [7].Therefore, the effect of system rotation is more subtle and complicated and yields new, richer features of flow and heat transfer in general, the bifurcation and stability in particular, for nonisothermal flows.While some of such new features are revealed by our recent analytical and numerical works [4,5,6,7,8], there is no known study on the bifurcation and stability of mixed convection in rotating curved ducts.
We note that all previous analytical/numerical studies of stability are limited to linear stability and some special disturbances.While the linear stability analysis is efficient in terms of the computation efforts required, it suffers three fundamental defects.First, it is not applicable to a finite disturbance.With a finite disturbance, a so-called stable solution based on linear stability may not be always stable.Second, it may not be so relevant for comparison with experiments.Because of the difficulty in controlling disturbances in experiments, experimental results of stability such as those by Matsson and Alfredsson [11] and Wang and Cheng [13,14,15] are essentially for finite random disturbances.Finally, the linear stability analysis provides no answer to the questions related to the dynamic behavior of the solutions, including how flows approach a stable solution after a disturbance, what happens to an unstable solution after a disturbance, whether all unstable solutions at a given set of parameters respond to disturbances in the same way, and whether the disturbances lead an unstable solution to the stable one at the same parameter value.Clearly, a fully transient computation is necessary to examine dynamic responses of the multiple solutions to the finite random disturbances.Such a computation is also capable of capturing the phenomena related to the transition to turbulence such as oscillation solution, periodic doubling, intermittency, and chaotic oscillation.
In previous numerical studies [9,16,17], branch stability is often determined by the stability of one point on the branch.This is partly due to the fact that the computation of the complete eigenvalue spectrum along the solution branches is a computationally expensive process and partly due to the assumption that the stability of solutions along a solution branch is unchanged without passing limit/bifurcation points in the literature.However, based on the bifurcation and stability theory, such a change in stability is possible [19].Therefore, a more detailed and careful stability analysis is desirable to observe the gain and loss of flow stability along solution branches without passing limit/bifurcation points.
The present work is a relatively comprehensive study on the bifurcation structure and stability of multiple solutions for the laminar mixed convection in a rotating curved duct of square cross-section (Figure 1).The governing differential equations in primitive variables are solved for detailed bifurcation structure by a finite-volume/Euler-Newton continuation method with the help of the bifurcation test function, the branch switching technique, and the parameterization of arc length (s) or local variable.Transient calculation is made to examine in detail the response of every solution family to finite random disturbances.The power spectra are constructed by the Fourier transformation of temporal oscillation solutions to confirm the chaotic flow.We restrict ourselves to the hydrodynamically and thermally fully developed region and 2D disturbances.So far, a detailed 3D numerical computation of flow bifurcation and stability is still too costly to conduct.A 2D model is still useful for a fundamental understanding of rotating curved duct flows.However, our assumption of fully developed flow limits our analysis to the one preserving the streamwise symmetry.There may be further bifurcation to flows that breaks this symmetry and that cannot be found in the present work.

GOVERNING PARAMETERS AND NUMERICAL ALGORITHM
Consideration is given to a hydrodynamically and thermally fully developed laminar flow of viscous fluid in a square duct with the streamwise curvature, the spanwise rotation, and the wall heating or cooling at a constant heat flux (Figure 1).The geometry is toroidal, and hence finite pitch effect is not considered.The rotation can be positive or negative at a constant angular velocity.The duct is streamwisely and peripherally uniformly heated or cooled with a uniform peripheral temperature.The properties of the fluid, with the exception of density, are taken to be constant.The usual Boussinesq approximation is used to deal with the density variation.The gravitational force is negligible compared with the centrifugal and Coriolis forces.Consider a noninertial toroidal coordinate system (R, Z, ϕ) fixed to the duct rotating with a constant angular velocity about the O Z axis, as shown in Figure 1.We may obtain the governing differential equations, in the form of primitive variables, governing fully developed mixed convection based on conservation laws of mass, momentum, and energy.The boundary conditions are nonslip and impermeable, streamwise uniform wall heat flux and peripherally uniform wall temperature at any streamwise position.The proper scaling quantities for nondimensionalization are chosen based on our previous experience [5].The formulation of the problem is on full flow domain without imposing symmetric boundary conditions to perform a thorough numerical simulation.The readers are referred to Wang and Cheng [5] for the details of mathematical formulation of the problem.
The dimensionless governing equations contain five dimensionless governing parameters: one geometrical parameter σ (the duct curvature ratio defined by a/R c , the ratio of duct width/height a over the radius of the curvature R c , representing the degree of curvature), one thermophysical parameter Pr (the Prandtl number representing the ratio of momentum diffusion rate to that of the thermal diffusion), and three dynamical parameters Dk, L1, and L2 defined by Wang and Cheng [5].The pseudo-Dean number Dk is the ratio of the square root of the product of inertial and centrifugal forces to the viscous force and characterizes the effect of inertial and centrifugal forces.L1 represents the ratio of the Coriolis force over the centrifugal force, characterizing the relative strength of Coriolis force over the centrifugal force.L2 is the ratio of the buoyancy force over the centrifugal force and represents the relative strength of the buoyancy force.A positive (negative) value of L1 is for the positive (negative) rotation.A positive (negative) value of L2 indicates the wall heating (cooling).In the present work, we set σ = 0.02 (typically used in cooling systems of rotor drums and conductors of electrical generators) and Pr = 0.7 (a typical value for air) to study the effects of three dynamical parameters on the multiplicity and stability.While results regarding the effects of Dk and L1 are also available, we focus on the effects of L2 at Dk = 300 and L1 = 28 in the present paper due to limited space.
The governing differential equations are discretized by the finite-volume method to obtain discretization equations.The discretization equations are solved for parameter dependence of velocity, pressure, and temperature fields by the Euler-Newton continuation method with the solution branches parameterized by L2, the arc length, or the local variable.The starting points of our continuation algorithms are the three solutions at Dk = 300, L1 = 28, and L2 = 0 from our study of the effects of Dk and L1.The bifurcation points are detected by the test function developed by Seydel [19].The branch switching is made by a scheme approximating the difference between branches proposed by Seydel [19].The dynamic responses of multiple solutions to the 2D finite random disturbances are examined by the direct transient computation.The readers are referred to [20,21,22,23] for the numerical details and the check of grid dependence and accuracy.The computations are carried out on the supercomputer SP2 of the University of Hong Kong.

SOLUTION STRUCTURE
The bifurcation structure is shown in Figure 2 for L2 values from −20 up to 70 at σ = 0.02, Pr = 0.7, Dk = 300, and L1 = 28.In Figure 2, the radial velocity component u at r = 0.9 and z = 0.14 (where the flow is sensitively dependent on L2) is used as the state variable, enabling the most clear visualization of all solution branches.Here r = R/a, z = Z/a.Three solution branches, labeled by AS1, AS2, and S3, respectively, are found.Here, S stands for symmetric solutions with respect to the horizontal central plane z = 0, and AS indicates that the branch has both symmetric and asymmetric solutions.Branch AS1 has sixty nine limit points labeled by AS1 1 to AS1 69 , eleven bifurcation points connecting its subbranches denoted by AS1 AS1−1 to AS1 AS1−11 , two bifurcation points connecting itself to AS2 labeled by AS1 AS2−1 and AS1 AS2−1 , and one bifurcation point connecting itself to S3 denoted by AS1 S3 .Branch AS2 has four limit points labeled by AS2 1 to AS2 4 .Branch S3 is a symmetric solution branch and has two limit points S3 1 and S3 2 .The location of fourteen bifurcation points and seventy five limit points is available in [23].To visualize the details of branch connectivity and some limit/bifurcation points, the locally enlarged state diagrams are also shown in Figure 2. As Figure 2 is only 1D projection of 12 400 dimensional solution branches, all intersecting points except fourteen bifurcation points should not be interpreted as connection points of branches.
For a large L2 value (L2 ≤ −14.5 or L2 ≥ 63.1), the buoyancy force dominates the mixed convection.There is unique flow and temperature field for a specified value of L2 in these two ranges.Figure 3 illustrates the secondaryflow patterns, the streamwise velocity profiles and temperature profiles at L2 = −17 and L2 = 65, respectively.In the figure, the stream function ψ, streamwise velocity w, and temperature t are normalized by their corresponding maximum absolute values ψ max , w max , and t max .A vortex with a positive (negative) value of the secondary-flow stream function indicates a counter-clockwise (clockwise) circulation.The readers are referred to [5] for a detailed discussion of the flow structures shown in Figure 3  effects on the flow resistance and heat transfer in particular.
For any L2 value in −14.5 < L2 < 63.1, however, we can have multiple solutions.Figure 4 shows typical secondaryflow patterns of six solutions (thirty-nine solutions in total) at L2 = −11.7.It is observed that the nonlinear competition of driven forces leads to not only a rich solution structure but also complicated flow structures.Therefore, the mixed convection in rotating curved ducts is much more complicated than that available in the literature.

STABILITY OF MULTIPLE SOLUTIONS
Recognizing that there is no study on dynamic responses of multiple solutions to finite random disturbances in the literature, a relatively comprehensive transient computation is made to examine the dynamic behavior and stability of typical steady solutions with respect to four sets of finite random disturbances with d = 4%, 10%, 15%, and 40%, respectively.The results presented in this paper are those obtained from the disturbance with d =10% unless otherwise stated.
Seven subranges are identified with each having distinct dynamic responses to the finite random disturbances.The first ranges from L2 = −20 to L2 = −14.5,where the finite random disturbances lead all steady solutions at any fixed L2 to a steady symmetric multicell state on AS1 a with the same L2.The second covers the range −14.5 < Dk ≤ −13.6,where there is coexistence of three stable steady symmetric multicell states.In the third subrange −13.6 < L2 ≤ 12.1, all steady solutions evolve to a temporal periodic solution.The fourth subrange is from L2 = −12.1 to L2 = −11.5,where the solutions response to the finite random disturbances in the form of either periodic oscillation or chaotic oscillation.There is the coexistence of periodic and chaotic oscillations.In the fifth subrange −11.5 < L2 ≤ −10.5, all steady solutions evolve to a temporal chaotic solution.The next subrange −10.5 < L2 ≤ −10.2 serves as a transition between the chaotic oscillation and the stable steady 2-cell flow.The solutions response to the finite random disturbances in the form of subharmonic-bifurcation-driven asymmetric oscillation.In the last subrange L2 > −10.2, the finite random disturbances lead all steady solutions at any fixed L2 to a stable steady symmetric 2-cell state on AS1 s with the same L2.
Stable steady symmetric multicell state: −20 ≤ L2 ≤ −14.5 Figure 5 typifies the responses of solutions on AS1 to the finite random disturbances.In the figure, the deviation of velocity components from their initial steady values is plotted against the nondimensional time τ at (0.9, 0.14), (0.94, 0.1), and (0.96, 0.06) for L2 = −15.We plot both radial (u-) and spanwise (v-) velocity components for the first point (0.9, 0.14) while only u-velocity component is shown for the last two points.To facilitate the comparison, we use these four velocity components (either velocity itself or derivation velocity from its initial steady value) in all figures illustrating dynamic responses of the multiple solutions to finite random disturbances.It is observed that all deviation velocities vanish after a short period of time.The flows and temperature profiles return to their initial steady symmetric multicell state similar to that shown in Figure 3a.
Coexistence of three stable steady symmetric multicell states: −14.5 < L2 ≤ −13.6 Figure 6 illustrates the typical response of solutions a, d, and e at L2 = −14 labeled in Figure 2b to the finite random disturbances.It shows that the finite random disturbances lead eventually these solutions to the steady solutions a, c, and e.There is coexistence of three stable steady symmetric multicell states a, c, and e.
Temporal periodic solution: −13.6 The dynamic response of solution b at L2 = −12.1 (Figure 2b) is shown in Figure 7a.The finite random disturbances here lead the solution to a temporal periodic state with a period of 0.065.Some typical secondary-flow patters are detailed in Figure 7b within one period of τ.We clearly observe the temporal oscillations between symmetric 12-cell flows and symmetric 14-cell flows.A detailed study on dynamic responses of the other solutions at L2 = −12.1 and the comparison of flow and temperature fields within one period show that the finite random disturbances lead all the solutions at the same L2 to the same periodic oscillation.
A similar dynamic evolution pattern exists for all cases with different values of L2.This signals the similarity of flow and temperature fields within one period among the periodic states for different values of L2 in the range −13.6 < Dk ≤ −12.1.Our detailed examination of flow and temperature fields has confirmed this and shown that the flow structures in Figure 7b are typical for all L2 in this range.

Coexistence of periodic and chaotic oscillations:
The coexistence of periodic and chaotic oscillations is shown by the dynamic response of four solutions at L2 = −11.7 to the finite random disturbances in Figure 8.While solutions b and c in Figure 4 response to the finite random disturbances in the form of periodic oscillation, solutions d and e evolve to chaotic oscillations.The typical secondaryflow patterns for the periodic oscillating solutions are similar to those in Figure 7b and are characterized by symmetric  2b and solution d : u(0.9, 0.14) = 2.493; and (c) e labeled in Figure 2b and solution e : u(0.9, 0.14) = − 16.52.patterns of twelve cells at least.Those for the chaotic oscillation are shown in Figure 9 and are featured by symmetric and asymmetric flows of six cells at most.Furthermore, the selection of periodic or chaotic oscillation is also dependent on disturbances.10%.The power spectra of four velocity temporal series in Figure 10 are constructed by the Fourier transformation.They contain the broadband noise, indicating the flow being chaotic [19].Their sensitivity to the initial conditions further confirms this.
Subharmonic-bifurcation-driven asymmetric oscillation: The transition from the stable steady flow in L2 > −10.2 to the chaotic oscillation in −11.5 < L2 ≤ −10.4 is characterized by subharmonic-bifurcation-driven asymmetric oscillation and shown by the variation of dynamic responses of solutions as L2 decreases from −10.2 to −10.35 (Figure 11).The power spectra of the velocity temporal series in Figure 11 are constructed by the Fourier transformation and shown in Figure 12.At L2 = −10.2, the oscillation is periodic and is from the Hopf bifurcation.The frequency f 1 in Figure 12a is universal for all velocity components, pressure, and temperature at all points in the flow domain.bifurcation doubles the period of the temporal oscillation at L2 = −10.3,generating another frequency f 3 that is a quarter of f 1 (Figures 11c and 12c).At L2 = −10.35, the oscillation is aperiodic-like (Figure 11d); the spectrum has many closely spaced peaks (Figure 12d), forecasting the onset of chaos.
Stable steady symmetric 2-cell state: −10.2 < L2 ≤ 70 Figure 13 typifies the responses of solution a at L2 = 50 (Figure 2a) to the finite random disturbances.The disturbances lead the solution a to the solution c whose secondaryflow pattern is shown in Figure 13.An extensive transient computation concludes that the solution on AS1 s is the unique stable state in −10.2 < L2 ≤70.

CONCLUDING REMARKS
The governing differential equations from the conservation laws for the mixed convection in rotating curved ducts are discretized by the finite-volume method to obtain discretization equations, a set of nonlinear algebraic equations.The discretization equations are solved for parameter dependence of flow and temperature fields by the Euler-Newton continuation with the solution branches parameterized by L2, the arc length, or the local variable.The bifurcation points are detected by the test function.The branch switching is made by a scheme approximating the difference between branches.One symmetric and two symmetric/asymmetric solution branches are found with fourteen bifurcation and seventy five limit points.Both solution and flow structures are much more richer than those available in the literature.
The dynamic responses of multiple solutions to the 2D finite random disturbances are examined by the direct transient computation.The finite random disturbances are found to lead the steady solutions to a stable steady multicell state in −20 < Dk ≤ −14.5, the coexistence of three stable steady multicell states in −14.5 < L2 ≤ −13.6, a temporal periodic oscillation in −13.6 < Dk ≤ −12.1, the coexistence of periodic and chaotic oscillating states in −12.1 < L2 ≤ −11.5, a chaotic oscillation in −11.5 < L2 ≤ −10.5, a subharmonicbifurcation-driven asymmetric oscillation in −10.5 < L2 ≤ −10.2, and a stable steady 2-cell state in −10.2 < L2 ≤ 70.

Figure 1 :
Figure 1: Physical problem and coordinating system.