Abstract: Unsteady solutions of the thermal flows through a curved rectangular duct of aspect ratio 2 is investigated numerically by using the spectral method over a wide range of the Dean number, 0≤Dn≤1000, with a temperature difference between the vertical sidewalls for the Grashof number, 1000≤Gr≤1500. The outer wall of the duct is heated while the inner wall is cooled. In the present study, unsteady solutions of the thermal flows for a single case of the Grashof number Gr = 1500 is investigated in detail and complete flow behavior, covering the wide range of Gr, is shown by a phase diagram. It is found that the steady flow turns into chaotic flow through periodic or multi-periodic flows for large or moderate Gr; for small Gr, however, the steady flow turns into chaotic state through various flow instabilities, if Dn is increased.
INTRODUCTION
Flows in a curved duct have attracted considerable attention not only because of their ample applications in chemical, mechanical, civil, biomechanical or biological engineering but also because of the physically interesting features under the action of the centrifugal force caused by the curvature of the duct. Curved diffusing passages are extremely used in many engineering applications, such as in air conditioning systems, refrigeration, heat exchangers, ventilators, centrifugal pumps and blade-to-blade passages in modern gas turbines. Blood flow in veins and arteries is another example of curved duct flows. One of the interesting phenomena of the flow through a curved duct is the bifurcation of the flow because generally there exist many steady solutions due to channel curvature. Studies of the flow through a curved duct have been made, experimentally or numerically, for various shapes of the cross section. For example, a circle (Dennis and Ng, 1982; Yanase et al., 1989), a semi-circle (Nandakumar and Masliyah, 1982), an oval (Kao, 1992) and a rectangle (Ligrani and Niver, 1988; Finlay and Nandakumar, 1990; Yanase et al., 2002). Recently, Yanase et al. (2005b) performed numerical simulations of non-isothermal flows (0≤Gr≤1000) through a curved rectangular duct of aspect ratio 2, where they obtained many branches of steady solutions and addressed the time-dependent behavior of the unsteady solutions.
Time dependent analysis of fully developed curved duct flows was initiated by Yanase and Nishiyama (1988) for a rectangular cross section and by Yanase et al. (1989) for a circular cross section in connection with the bifurcation diagram of steady solutions. In both the studies they investigated unsteady solutions for the case where dual solutions exist. The time-dependent behavior of the flow in a curved rectangular duct of large aspect ratio was investigated by Yanase et al. (2002) numerically. They performed time-evolution calculations of the unsteady solutions with and without symmetry condition and found that periodic oscillations appear with symmetry condition while aperiodic time variation without symmetry condition. Recently, Wang and Yang (2005) performed numerical as well as experimental investigations of periodic oscillations for the fully developed flow in a curved square duct. Flow visualization in the range of Dean numbers from 50 to 500 was conducted in their experiment. They showed, both experimentally and numerically, that a temporal oscillation takes place between symmetric/asymmetric 2-cell and 4-cell flows when there are no stable steady solutions. However, the time-dependent behavior of the flows through a curved square duct was investigated by Mondal et al. (2006) where they performed time-evolution calculations of the solutions and showed the transitional behavior of the unsteady solutions.
A remarkable characteristic of the flow through a curved duct is to enhance thermal exchange between two differentially heated vertical sidewalls, because it is possible that the secondary flow may convey heat and then increases heat flux between two sidewalls (Chandratilleke and Nursubyakto, 2003). Recently, Yanase et al. (2005) studied the bifurcation structure as well as the effects of secondary flows on convective heat transfer for moderate Grashof numbers (Gr≤1000). Very recently, Mondal et al. (2007) extended the study of Yanase et al. (2005) for larger Grashof numbers (1000≤Gr≤1500) and studied the flow characteristics. However, they did not perform unsteady solutions of the thermal flows. This study is, therefore, an extension of Mondal et al. (2007) with a view to study the nonlinear behavior of the unsteady solutions for large Grashof numbers.
The present research shows numerical results of the unsteady solutions through a curved rectangular duct with differentially heated vertical sidewalls for the large Grashof number. Complete flow behavior, covering a wide range of the Grashof number, is also presented by a phase diagram.
MATHEMATICAL FORMULATIONS
Consider a viscous incompressible fluid streaming through a curved duct with a constant curvature. The cross section of the duct is a rectangle with width 2d and height 2h. It is assumed that the outer wall of the duct is heated while the inner one is cooled. The temperature of the outer wall is T0+<ΔT and that of the inner wall is T0 - <ΔT, where <ΔT>0. The x, y and z axes are taken to be in the horizontal, vertical and axial directions, respectively. It is assumed that the flow is uniform in the z direction and that it is driven by a constant pressure gradient G along the center-line of the duct as shown in Fig. 1.
The variables are non-dimensionalized by using the representative length d and representative velocity
where v is the kinematic viscosity of the fluid. We introduce the non-dimensional variables defined as:
where, u, v and w are the velocity components in the x, y and z directions, respectively; t is the non-dimensional time, P the non-dimensional pressure, <δ the non-dimensional curvature defined as <δ = d/L, L being the radius of the duct curvature and temperature is nondimensionalized by <ΔT. Henceforth, all the variables are nondimensionalized if not specified. In the above method of nondimensionlization, the variables with prime denote the dimensional quantities.
| Fig. 1: | Coordinate system of the curved rectangular duct |
Since the flow field is uniform in the z-direction, the sectional stream function <ψ is introduced as:
| (1) |
A new coordinate variable ý is introduced in the y direction as y = lý, where l = h/d is the aspect ratio of the cross section. From now on, y denotes ý for the sake of simplicity. The basic equations for w, <ψ and T are then derived from the Navier-Stokes equations and the energy equation with the Boussinesq approximation as:
| (2) |
| (3) |
| (4) |
where,
| (5) |
The Dean number Dn, the Grashof number Gr and the Prandtl number Pr which appear in Eq. 2-4 are defined as:
| (6) |
where <μ, <γ, <κ and g are the viscosity, the coefficient of
thermal expansion, the coefficient of thermal diffusivity and the gravitational
acceleration, respectively.
The rigid boundary conditions for w and <ψ are used as:
| (7) |
and the temperature T is assumed to be constant on the walls as:
| (8) |
NUMERICAL METHOD
In order to solve the Eq. 2-4 numerically, the spectral method is used. This method is thought to be the best numerical method for solving the Navier-Stokes equations as well as the energy equation (Gottlieb and Orszag, 1977). By this method the variables are expanded in a series of functions consisting of the Chebyshev polynomials. (Mondal et al., 2007).
In this study, the resistance coefficient <λ is used as the representative quantity of the flow state which is also defined in our earlier study (Mondal et al., 2007). In the present study, numerical calculations are carried out for the curvature <δ = 0.1 over a wide range of the Dean number 0≤Dn≤1000 for the Grashof number 1000≤Gr≤1500 for l = 2. The resistance coefficient <λ is used to discriminate the steady solution branches and to pursue the time evolution of the unsteady solutions.
| Fig. 2: | Steady solution branches for Gr = 1500 and 100≤Dn≤1000 |
RESULTS
Steady solutions: Using the path continuation technique with various initial guesses as discussed by Mondal (2006), we obtained five branches of asymmetric steady solutions in our previous research (Mondal et al., 2007) for the curvature <δ = 0.1 over a wide range of the Dean number 0≤Dn≤1000 and the Grashof number 1000≤Gr≤1500. A bifurcation diagram of steady solutions, for example and for convenience, is shown in Fig. 2 for 100≤Dn≤1000 and Gr = 1500 using <λ, the representative quantity of the solutions. The steady solution branches are named the 1st steady solution branch (1st branch, long dashed line), the 2nd steady solution branch (2nd branch, thin solid line), the 3rd steady solution branch (3rd branch, thick solid line), the 4th steady solution branch (4th branch, dash dotted line) and the 5th steady solution branch (5th branch, dotted line), respectively. The steady solution branches are distinguished by the nature and number of secondary flow vortices appearing in the cross section of the duct as discussed in Mondal et al. (2007). Stability characteristics of the flows were also conducted in that research. However, time-dependent behavior of the unsteady solutions was left to investigate which is conducted in the present study.
Unsteady solutions for Gr =1500: In order to study the nonlinear behavior of the unsteady solutions and furthermore to determine the transitional behavior of the unsteady solutions, time evolution calculations are performed. Though the present study covers unsteady solutions over a wide range of Gr (1000≤Gr≤1500), in the present study, however, we discuss the results of time-evolution calculations for Gr = 1500 only and complete unsteady flow behavior, covering the wide range of Gr, is discussed in the next section and is shown by a phase diagram in Fig. 8.
| Fig. 3: | Contours of secondary flow (top) and temperature profile (bottom) for Dn = 100 and 150 at t = 10 |
Time evolutions of <λ for Dn≤105 and 144≤Dn≤165 show that the value of <λ quickly approaches steady state. The reason is that, in our earlier research (Mondal et al., 2007), we found that among five branches of steady solutions, only the first branch, which existed throughout the whole range of the Dean number, was linearly stable in two different intervals of the Dean number, 0≤Dn≤105 and 144≤Dn≤165. Therefore, time-evolution results are consistence with the stability results. Figure 3 shows secondary flow pattern and temperature profiles for Dn = 100 and Dn = 150 at Gr = 1500 when t = 10. In Fig. 3, to draw the contours of <ψ and T we use the increments <Δψ = 0.6 and <ΔT = 0.2, respectively. The same increments of <ψ and T are used for all the figures in this study, unless specified. The right-hand side of each duct box of <ψ and T is in the outside direction of the duct curvature. In the figures of the secondary flow, solid lines (<ψ≥0) show that the secondary flow is in the counter clockwise direction while the dotted lines (<ψ<0) in the clockwise direction. Similarly, in Fig. 3 of the temperature field, solid lines are those for T≥0 and dotted ones for T<0. As shown in Fig. 3, the unsteady solution for Dn≤105 is a single-vortex solution while that for 144≤Dn≤165 a two-vortex solution. The secondary flows are asymmetric with respect to the horizontal center plane y = 0. The reason is that, heating the outer wall causes deformation of the secondary flow and yields asymmetry of the flow.
| Fig. 4: | The results for Gr = 1500 and Dn =125. (a) Time evolution
of <λ and the value of <λ for the first steady solution for 0≤t≤40
and (b) secondary flows (top) and temperature profiles (bottom) for one
period of oscillation at 32.0≤t≤35.45 |
Then, in order to see the unsteady flow behavior, obtained for 144≤Dn≤165 where the solution is linearly unstable on the first branch, time evolution calculations are performed at Dn = 125. Figure 4a shows time-evolution result for Dn = 125, where it is seen that the flow oscillates multi-periodically. In Fig. 4a, to observe the relationship between the periodic solution and the steady state, the value of <λ for the steady solution branch at Dn = 125 is also shown by a straight line using the same kind of line, which was used in the bifurcation diagram in Fig. 2. As shown in Fig. 4a, the periodic solution at Dn = 125 oscillates in the region above the 1st steady solution branch. It shows that the 1st branch plays a role of an envelope of this periodic oscillation. To observe the periodic change of the flow patterns and temperature distributions, contours of secondary flow and temperature profiles are shown in Fig. 4b, for one period of oscillation at 32.0≤t≤35.45, where it is seen that the multi-periodic oscillation at Dn = 125 is a two-vortex solution with one large vortex dominating the other one.
| Fig. 5: | The results for Gr = 1500 and Dn = 500. (a) Time evolution
of <λ and the values of <λ for the steady solutions for 0≤t≤10
and (b) secondary flows (top) and temperature profiles (bottom) for 8.50≤t≤8.66 |
Time evolution of <λ is then performed for Dn = 500 as shown in Fig. 5a where, the values of <λ for the steady solution branches are also plotted by straight lines. As shown in Fig. 5a, the flow oscillates periodically which takes place in the region between the first and third steady solution branches. Initial condition independence has also been examined using the initial condition from another steady solution branch and it is found that the periodic oscillation drifts in the same place. Secondary flow patterns and temperature distributions, for one period of oscillation at 850≤t≤8.66, are shown in Fig. 5b. It is found that the secondary flow becomes complete two-vortex solution as the Dean number is increased. Time evolution of <λ for Dn = 600 is shown in Fig. 6a, where the values of <λ for the steady solution branches are also plotted by drawing straight lines in order to see the relationship between the steady and unsteady solutions calculated at the same Dean number. As seen in Fig. 6a, the flow oscillates irregularly with the large windows of quasi-periodic oscillations, which suggests that the flow is chaotic.
| Fig. 6: | The results for Gr = 1500 and Dn = 600. (a) Time evolution
of <λ and the values of <λ for the steady solutions for 0≤t≤20
and (b) secondary flows (top) and temperature profiles (bottom) for 14.0≤t≤16.0 |
It is found that the chaotic solution at Dn = 600 fluctuates around <λ = 0.2122 on the upper branch of the 4th steady solution at Dn = 600. Contours of secondary flow and temperature profile for Dn = 600 at 14.0≤t≤16.0 are shown in Fig. 6b, where it is seen that the chaotic oscillation at Dn = 600 is a two- and four-vortex solution.
Next, time evolution of <λ together with the values of <λ for the steady solution branches, indicated by straight lines, are shown in Fig. 7a for Dn = 1000, where it is found that the flow is chaotic. As seen in Fig. 7a, the unsteady solution at Dn = 1000 oscillates above all the steady solution branches and the upper part of the 3rd branch, which has the maximum <λ (<λ≈ 0.15927)) at Dn = 1000, looks like an attractor of this chaotic solution.
| Fig. 7: | The results for Gr = 1500 and Dn = 1000. (a) Time evolution
of <λ and the values of <λ for the steady solutions for 0≤t≤10
and (b) secondary flows (top) and temperature profiles (bottom) for 6.0≤t≤8.0 |
The chaotic solution for Dn = 1000 is called a strong chaos but that for Dn = 600 a weak chaos (Mondal et al., 2006), because the chaotic solution at Dn = 600 is still trapped by the steady solution branches but that for Dn = 1000 tends to get away from them. In order to observe the change of the flow patterns and temperature distributions, contours of secondary flow and temperature profiles for Dn = 1000 are shown in Fig. 7b, where the increments <Δψ = 1.2 and <ΔT = 0.4 are used to draw the contours of secondary flow and temperature profile, respectively. As seen in the secondary flow patterns, the chaotic solution at Dn = 1000 is composed of four-vortex solutions only.
Unsteady solutions in the Dn-Gr plane: Here, the distribution of the steady, periodic and chaotic solutions, obtained by the time evolution computations, is presented by a phase diagram as shown in Fig. 8 in the Dean number versus Grashof number (Dn-Gr) plane for 0≤Dn≤1000 and 1000≤Gr≤1500. In Fig. 8, the circles indicate stable steady solutions, crosses periodic solutions and triangles chaotic solutions. As shown in Fig. 8, the steady flow turns into chaotic flow through various flow instabilities. For Gr<1250, the periodic solution appears in three different intervals of the Dean number in the scenario 'steady → periodic → steady → periodic → steady → periodic → chaotic', if Dn is increased.
| Fig. 8: | Distribution of the time-dependent solutions in the Dean number
vs. Grashof number (Dn-Gr) plane for 0≤Dn≤1000 and 1000≤Gr≤1500
(•) steady-state solution, (χ) periodic solution and (<Δ)
chaotic solution |
For larger Gr (Gr>1250), however, the periodic solution occurs in two distinct intervals of Dn and the flow undergoes 'steady → periodic → steady → periodic → chaotic', if Dn is increased. It is found that the flow characteristics drastically change at Gr ≈ 1250 and the chaotic solution occurs for Dn ≥ 600 whether Gr is small or large.
DISCUSSION
In the present study, a numerical simulation of fully developed two-dimensional (2-D) flow of viscous incompressible fluid through a curved rectangular duct is presented and a brief discussion on the plausibility of applying 2-D calculations to study curved duct flows will be given in this section. It has been shown by many experimental and numerical (Mondal et al., 2006) investigations that curved duct flows easily attain asymptotic fully developed 2-D states (uniform in the main flow direction) at most 270° from the inlet. Wang and Yang (2005) showed that even periodic flows can be analyzed by 2-D calculation. They showed that for an oscillating flow, there exists a close similarity between the flow observation at 270° and 2-D calculation. In fact the periodic oscillation, observed in the cross section, was a traveling wave advancing in the downstream direction. The same phenomenon is also obtained in the present numerical calculation. Therefore, it is found that 2-D calculations can predict the existence of three dimensional traveling wave solutions as an appearance of 2-D periodic oscillation.
In Fig. 8 we presented the bifurcation diagram where it is shown that there exists a region of stable steady solutions in lowest Dn region and oscillating solutions appear if Dn is increased. If Dn is increased further, a region of stable steady solutions again appears. The flow then turns into chaotic region through periodic solutions if Dn is increased further and the transition process from periodic to chaotic oscillation can be predicted by 2-D analysis as shown by Yamamoto et al. (1995). The transition to chaos of the periodic oscillation, obtained by the 2-D calculation in the present study, may correspond to destabilization of traveling waves in the curved duct flows like that of Tollmien-Schlichting waves in a boundary layer. Our 2-D analysis, therefore, may contribute to the study of curved duct flows by giving a complete outline for not only fully developed but also developing curved duct flows.