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 bladetoblade 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 semicircle (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 nonisothermal flows (0≤Gr≤1000) through a curved rectangular duct of aspect ratio 2, where they obtained many branches of steady solutions and addressed the timedependent 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 timedependent behavior
of the flow in a curved rectangular duct of large aspect ratio was investigated
by Yanase et al. (2002) numerically. They performed timeevolution 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 2cell and 4cell flows when there are no stable steady
solutions. However, the timedependent behavior of the flows through a curved
square duct was investigated by Mondal et al. (2006) where they performed
timeevolution 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 T_{0}+<ΔT
and that of the inner wall is T_{0}  <Δ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 centerline of
the duct as shown in Fig. 1.
The variables are nondimensionalized by using the representative length d and representative velocity ,
where v is the kinematic viscosity of the fluid. We introduce the nondimensional
variables defined as:
where, u, v and w are the velocity components in the x, y and z directions,
respectively; t is the nondimensional time, P the nondimensional pressure,
<δ the nondimensional 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 zdirection, the sectional stream function <ψ is introduced as:
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 NavierStokes equations and the
energy equation with the Boussinesq approximation as:
where,
The Dean number Dn, the Grashof number Gr and the Prandtl number Pr which appear
in Eq. 24 are defined as:
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:
and the temperature T is assumed to be constant on the walls as:
NUMERICAL METHOD
In order to solve the Eq. 24 numerically,
the spectral method is used. This method is thought to be the best numerical
method for solving the NavierStokes 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, timedependent 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 timeevolution 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, timeevolution
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 righthand 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 singlevortex
solution while that for 144≤Dn≤165 a twovortex 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 timeevolution
result for Dn = 125, where it is seen that the flow oscillates multiperiodically.
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 multiperiodic
oscillation at Dn = 125 is a twovortex 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 twovortex 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 quasiperiodic 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 fourvortex 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 fourvortex solutions only.
Unsteady solutions in the DnGr 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 (DnGr) 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 timedependent solutions in the Dean number
vs. Grashof number (DnGr) plane for 0≤Dn≤1000 and 1000≤Gr≤1500
(•) steadystate 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 twodimensional (2D) flow of viscous incompressible fluid through a curved rectangular duct is presented and a brief discussion on the plausibility of applying 2D 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 2D 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 2D calculation. They showed that for an oscillating flow, there exists a close similarity between the flow observation at 270° and 2D 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 2D calculations can predict the existence of three dimensional traveling wave solutions as an appearance of 2D 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
2D analysis as shown by Yamamoto et al. (1995). The transition to chaos
of the periodic oscillation, obtained by the 2D calculation in the present
study, may correspond to destabilization of traveling waves in the curved duct
flows like that of TollmienSchlichting waves in a boundary layer. Our 2D 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.