- Email: [email protected]
Convection in a low Prandtl number fluid
,. . . . . . . .
CRYSTAL G R O W T H
ELSEVIER
Journal of Crystal Growth 162 (1996) 95-106
Convection in a low Prandtl number fluid A.C. Skeldon a,* ,l, D.S. Riley a K.A. Cliffe b a Department of Theoretical Mechanics, University Park, Nottingham NG7 2RD, UK b AEA Technology, B424.4 Harwell Laboratory, Didcot, Oxon OXll ORA, UK
Received 21 August 1995; accepted 13 October 1995
Abstract
As a simple two-dimensional model of convection in the liquid phase during crystal growth using the Bridgman technique, we consider the fluid flow in a shallow rectangular cavity heated from one side. For an aspect ratio of 4 (aspect ratio = length/heigh0, previous numerical studies have shown the existence of two types of oscillatory solution. The first of these has been well-studied for a range of Prandtl numbers from 0 to 0.015. However, the form of the second has been observed only in a time-dependent study for a Prandtl of zero. We locate the new Hopf bifurcation which gives rise to this latter solution and study its dependence on Prandtl number. Before the onset of oscillations, various steady corotating multi-cell solutions are found depending on the aspect ratio of the cavity. We compute how the transition between three corotating cells and two corotating cells takes place for changing aspect ratio. Understanding can be facilitated by comparison with analogous problems. We consider tilting the cavity to the B6nard configuration, where the fluid is heated from below rather than from the side. In this case a no-flow solution exists where heat is transferred by conduction alone. This solution becomes unstable to counter-rotating cells at a symmetry-breaking bifurcation. We study how this symmetrybreaking bifurcation disconnects as the cell is tilted and the solutions evolve into the side-wall cavity solutions. In addition we trace the saddle-node and Hopf bifurcations found in the side-wall heated problem for changing tilt and reveal that they also exist in the B6nard convection limit; disconnected solutions in the B6nard problem have not been studied previously.
1. Introduction
The growth of crystals from a melt is an important process for the semiconductor industry. The crystal quality is determined by many factors, from interfacial dynamics at the molecular level to the effect of melt convection on the bulk distribution of dopants. The difficulty o f the problem leads to the consideration of simplified models which generally concentrate on one aspect: in this paper we focus on
* Corresponding author. Fax: +44 171 477 8000. Present address: Department of Mathematics, City University, Northampton Square, London ECIV 0HB, UK. Elsevier Science B.V. SSDI 0022-0248(95)00923-X
the convective instabilities that occur in the melt. F o l l o w i n g a similar rationale, Hurle [1] studied convection in a shallow boat of rectangular cross-section heated from the side as a simple model of the fluid dynamics associated with the Bridgman process o f crystal growth. He found that, if a small temperature difference was applied across the melt, then the flow was steady. However, on increasing the temperature difference, the flow became oscillatory. A number of experimental studies, for example Miiller and W i e h e l m [2], showed that these oscillations were linked to " s t r i a t i o n s " which appear in the final crystal and which are detrimental to the performance of semi-conductor devices made from the crystal.
96
A.C. Skeldon et aL /Journal of Crystal Growth 162 (1996) 95-106
Much of the theoretical work on this configuration to date focuses on an infinite shallow cavity in which it is assumed that the flow is parallel and that there is an imposed, constant, thermal flux along the cell. The governing equations can then be solved exactly, resulting in a cubic polynomial for the vertical velocity dependence and a quintic polynomial for the corresponding temperature profile. Stability analyses of this "Hadley cell" solution have been performed, first by Hart [3,4] and Gill [5], and later refined by Kuo and Korpela [6] and by Laure and Roux [7]. These have shown that the parallel-flow solution is stable at low Grashof numbers, but is destabilised as the Grashof number is increased. In the case of rigid, insulating upper and lower boundaries and for a Prandtl number less than 0.034, the resulting flow consists of stationary corotating cells whose axes are aligned transverse to the imposed temperature gradient. The critical Grashof number for this instability asymptotes to about 8000 as the Prandtl number is reduced to zero. This transition is a shear instability, analogous to the shear instability of a plane parallel flow found by Birikh [8]. For a Prandtl number greater than 0.034, the first instability of the steady parallel flow is to an oscillatory roll solution, where now the axes of the rolls are parallel to the imposed temperature gradient (longitudinal rolls). A qualitatively similar picture holds if the upper and lower boundaries are both conducting, or if the upper boundary is stress free. Cormack et al. [9] considered a shallow cavity, but with a finite length. They showed that in the central region of this cavity it was reasonable to assume a parallel flow solution providing Gr2Pr2/A3 ~ 3 × 105, where A >> 1 is the aspect ratio of the cavity. Note that for low Prandtl numbers (Pr = 0.02) and for more moderate aspect ratios (A 4) this implies that the parallel solution is a reasonable solution providing Gr ,~ 220 000. Much of the related numerical work has considered rectangular cavities filled with a low Prandtl number fluid and heated from the side. In the main these have concentrated on Prandtl numbers of either 0 or 0.015 and an aspect ratio of 4, which were the values used in a benchmarking exercise at a GAMM workshop [10]. The contributions to this workshop considered rigid lower boundary and end walls, and either a rigid or stress-free surface. Qualitatively
different flows were observed with the different fluid boundary conditions. For rigid walls, the flow is steady and parallel at low Grashof number, but as the Grashof number is increased, this is replaced by a steady solution with three corotating cells. On increasing the Grashof number further, this three-cell solution is in turn replaced by an oscillatory solution at a Hopf bifurcation. During one cycle, the oscillatory solution interchanges between a one-cell state and a three-cell state. Winters [11] traced the path of this Hopf bifurcation and showed that decreasing the aspect ratio or increasing the Prandtl number inhibits the onset of oscillation, causing the critical Grashof number to increase. For the case of zero Prandtl number, the dependence on the aspect ratio was confirmed by Okada and Ozoe [12]. The GAMM workshop also highlighted the presence of a second steady solution consisting of two corotating cells; Winters showed that this arose at a saddle-node bifurcation. Pulicani et al. [13] computed the time-dependent solutions for increasing Gr (see also Ref. [14]). As a result of the GAMM workshop, they conclude that zero Prandtl number calculations are representative of small Prandtl number calculations (Pr--~ 10 -2) and thus restrict their study to this case. This has the advantage that the temperature field is then unaffected by the flow and can be solved analytically, leaving only the fluid equations to be solved numerically. As in the GAMM workshop, they find that, on increasing the Grashof number, there is first a transition from one steady cell to three steady corotating cells, and then a further transition to an oscillatory solution. In addition, they find that on further increase of the Grashof number, this periodic solution is first replaced by a quasi-periodic motion, then by a second type of periodic solution; the quasi-periodic motion exhibits some hysteresis. In this paper we extend the above results. We show that if the aspect ratio is increased to 5.5 then it is the two corotating cell solution which is preferred, that is, it is the two cell solution which is reached by quasi-static increase of the Grashof number from zero. The three-cell solution becomes disconnected and we compute the exchange mechanism to demonstrate how this transition occurs. Winters [1 l] found one path of Hopf bifurcations and mentions in passing the presence of a second: we have computed the
A.C. Skeldon et al./Journal of Crystal Growth 162 (1996)95-106
second Hopf bifurcation in detail. For aspect ratio 4, this new Hopf bifurcation occurs at a critical Grashof number greater than that computed by Winters. However, we trace both types of Hopf bifurcation and show that an exchange of stability between the two Hopf bifurcations can occur. The presence of this second Hopf bifurcation helps to explain the results obtained by Pulicani et al. One question that naturally arises in the present study is: how does the side-wall heating problem relate to the B~nard problem where a fluid layer is heated from below, and, in particular, how does the disconnected branch arise? A transition from one configuration to the other can be made by tilting the cavity. In contrast to the side-wall heating problem, the Btnard problem has a no-flow solution where the fluid is stationary and the heat is transferred by conduction alone. For small Grashof numbers, the conduction solution is stable, however as the Grashof number is increased, this solution becomes unstable through a pitchfork bifurcation producing a number of counter-rotating cells. This transition is buoyancy driven and the number of cells produced is determined by the aspect ratio. Cliffe and Winters [15] showed that as a square cell is tilted, the leading pitchfork bifurcation to a unicellular flow unfolds: there is a smooth transition to a stable unicellular flow representing the natural circulation. A disconnected solution arises at a saddle-node bifurcation with a stable branch representing anomalous un-natural circulation. Tracing the saddle-node bifurcation for increasing tilt showed that for the case Pr = 1 it asymptoted to a fixed angle of tilt, 22 ° and consequently, no analogous solutions exist in the side-wall heating problem. Riley and Winters [16] showed that in a similar problem where the cavity is filled with a porous medium (the Lapwood problem), the disconnected solution branch again does not exist in the side-wall limit, although they found a more complicated sequence of transitions involving isola formations. These studies motivated us to look further at the relationship between the B~nard problem and the side-wall heating problem. We compute the bifurcation structure for a tall thin cavity filled with a low Prandtl number fluid and heated from below and show how this structure is altered as the cavity is tilted. In addition we trace the two types of Hopf bifurcation and the saddle-node bifurcation which
97
gives rise to stable corotating cells in the side-wall heated problem for changing tilt. In the next section the equations of motion for the problem are stated. The numerical methods used are discussed in Section 3 which then leads to the results, presented in Section 4. The results are discussed and our conclusions drawn in Section 5.
2. Problem formulation A rectangular cavity of height H and length L filled with a low Prandtl number fluid is considered. One side of the cavity is maintained at a constant temperature T h, whilst the opposite side is kept at a lower temperature, Tc. The remaining two boundaries are assumed insulating. All the walls are taken to be impermeable and non-slip. In order to relate the side-wall heating problem to the B~nard problem, the angle between the side of the cell and the gravity vector, 0 is defined such that when 0 = 0 the gravity vector and the temperature gradient are parallel,the Btnard configuration; when 0 = 7r/2, the side-wall heating problem is recovered. We assume that the convection is described by the Navier-Stokes and energy equations subject to the Boussinesq approximation and non-dimensionalise using H, G r l / 2 o / H , H Z / v , T h - T~, Grpov2/H 2 as the length,velocity, time, temperature and pressure scales respectively where Gr = g ~ H 3 (T h - Tc)/v 2. This results in the four non-dimensional equations:
Ou Ov --+--=0, Ox Oy Ou Ot + Grl/2
( OU v OU I U~x + -~y]
02u 02u am = -Gr~/2~x +-+-+ Grf/2Tcos 0 OX 2
Oy 2
- - + Gr ~/2 u - - + Ot Ox -~y ] Op 02 v 02v = -- Grl/20---y + --Ox2 + --Oy2 + Grl/ZT sin O, --+Gr
Ot
~/2 u - - + v - -
Ox
Oy
=
--+
Ox 2
cgy2 J"
A.C. Skeldon et aL /Journal of Crystal Growth 162 (1996) 95-106
98
Along the non-slip and insulating boundaries:
u=v=O,
3T --=0,
Oy
+fp ~
whilst on the isothermal and non-slip boundaries:
u=v=O,
T= - ±2 '
on the cold side, and u=v=0,
T=½,
on the opposite warm side. When the angle of tilt is non-zero the equations of motion together with the boundary conditions result in a centrosymmetric problem in that, if the origin of our coordinate system is taken as the centre of the cell, the problem is invariant under the transformation: x~-x, y~-y, u~-u, v--*-v, T-~ - r .
(2)
When the angle of tilt is zero, that is, when the fluid is heated from below, the problem is in addition invariant under the reflection symmetry:
x--~ - x ,
y--~y,
u-~ - u ,
v--~v,
T--~ - T .
(3)
3. Numerical methods
We consider the pressure-integrated weak form for Eq. (1), that is:
L
+Grl/2
+ 7y)
- Gr~/2Tcos O]
+
dxdy=0.
Where, f , , f~, f r and fp are test functions. These test functions are chosen so that there is no contribution from the boundary integral. A standard Galerkin formulation of the finite element method was adopted in which the test functions were chosen from the same set of functions as the basis functions. The temperature and velocities were represented as a sum of quadratic basis functions whilst linear discontinuous basis functions were used for the pressure. The rectangular domain is divided into rectangular nine-node elements. In order to capture any rapid changes which occur at the boundaries of the cavity, the spacing of the nodes is graduated so that there are more elements near the walls. Grids of 30 X 12, 40 X 16 and 60 X 24 elements were used. The coarsest grid gave qualitatively correct results, however, there was a large variation in the accuracy to which bifurcation points were found, depending on the values of the parameters and the particular type of bifurcation point studied. At best, the coarsest grid gave critical parameter values which were to within 1% of those found by using the finest grid, but at worst it computed the value of the critical parameters only to within 40%. In particular, accurate location of the second type of Hopf bifurcation at Prandtl numbers of above 0.01 was only achieved with the finer grids. Most of the results presented below were computed using the 60 X 24 element grid. The resulting discretised equations are of the form:
Ox
M - - + f ( x , , ~ ) --0. -
Grl/2POf"~x + Of. Ou Ox Ox
--Gr~/2Dox
+
Of. Ou Oy Oy
aL a v + af~ av Ox ax Oy ay
+--
dt
+ GrV~
Ox
(5)
Here, M is the mass matrix, x is the solution vector and t~ is a vector of the control parameters, that is, the Prandtl number, the Grashof number, the aspect ratio and the tilt angle. We seek steady-state solutions to this problem by solving:
f ( x , a ) =0,
+A ~
(4)
(6)
using Newton's method, which is guaranteed to converge providing a good enough initial guess to the
A.C. Skeldon et al./ Journal of Crystal Growth 162 (1996) 95-106
solution is made. For low Grashof numbers, we find x - - 0 is a sufficiently good guess, but to reach higher Grashof numbers, a path-following approach is employed. That is, the Grashof number is incremented in small steps; at each step Newton's method is used to find a solution, using the solution at the previous step to give an initial estimate for the new solution. Both stable and unstable solutions can be followed, but this method fails at saddle-nc~le points. This, however, is readily overcome by parameterising the solution branch with arclength, as originally suggested by Keller [17]. Only stable solutions are physically realisable and thus monitoring changes of stability of the solutions is also important. Generically, changes of stability occur either when a single real eigenvalue passes through zero or when a complex conjugate pair of eigenvalues crosses the imaginary axis. In the former case, a steady-state bifurcation occurs and in the latter, a Hopf bifurcation. In each case an "extended" system of necessary conditions can be solved to locate the bifurcation point. For example, in order to find a saddle-node bifurcation we follow Jepson and Spence [18], and solve the extended system: =0,
99
along the branch of solutions the bifurcation parameter increases and then decreases (or vice versa), and secondly there is a change in the number of eigenvalues with positive real part. As the determinant of the Jacobian equals the product of the eigenvalues, this second feature is signalled by a corresponding change in the sign of determinant. Furthermore, since we compute the Jacobian in order to use Newton's method and calculating the determinant is computationally inexpensive, this second feature can be checked with little additional computational time. When a Hopf bifurcation occurs a complex conjugate pair of eigenvalues crosses the imaginary axis. Hence, although a change of stability may occur, there is no corresponding change in the sign of the determinant of the Jacobian. Computing all the eigenvalues of the Jacobian is computationally intensive, prohibitively so in this case. Instead we use a new method developed by Cliffe et al. [21]. Rather than compute the entire spectrum of eigenvalues, this method seeks to compute only the "most dangerous" i.e. those of smallest real part: it is these eigenvalues which determine when bifurcations occur. By tracking these eigenvalues a good estimate of where Hopf bifurcations occur can be made.
= o,
l(~b) - 1 = 0.
(7)
Here, ~b is the eigenvector of the Jacobian matrix, /x(x0,~t), and the last equation in Eq, (7) is a normalisation condition for ~. There is a corresponding extended system which gives necessary conditions for a Hopf bifurcation to occur (cf. Jepson [19] and Griewank and Reddien [20]). In both cases the extended systems are solved using Newton's method and paths of bifurcation points are followed in parameter space in the same way that paths of regular solution points are followed. As for the location of regular solution points, in order for Newton's method to converge to a solution of the extended system, a good enough initial guess for the bifurcation is required. Our experience suggests that in practice, the basin of convergence can be very small and hence a very good guess is needed. This does not pose a problem for the location of saddle-node bifurcations since these have two readily monitored characteristics. Firstly, stepping
4. Results 4.1. Side-wall heating problem
Consistent with previous studies, at small Grashof number and A = 4, we find that the solution is parallel over a large percentage of the cavity. As the Grashof number is increased, the region over which the flow is parallel gradually shrinks; further increasing the Grashof number results in the growth of two smaller cells at either end of the cavity and gives rise to the existence of three corotating cells. These transitions are the remnant of the steady transverse shear instability found by Hart. The parallel-flow solution in an infinite layer has a continuous translational symmetry. At Hart's instability this symmetry is broken and the resulting stable corotating cells have only a discrete translational symmetry. The finite cavity does not have the translational symmetry of the infinite layer, and hence this bifurcation is
100
A.C. Skeldon et aL /Journal of Crystal Growth 162 (1996) 95-106
unfolded into a smooth evolution to corotating cells, where the cells rotate such that the flow passes up the hot wall and down the cold wall. If the aspect ratio is increased to 5.5, instead of a transition from a parallel flow to one cell and then to three cells, there is a smooth evolution to two cells, as shown in Fig. 1 - this is qualitatively different from the behaviour at lower aspect ratio. We have computed the mechanism by which this change occurs and the results are summarised in Fig. 2. This figure displays a bifurcation set showing the computed path of saddle-node bifurcation as a function of aspect ratio and Grashof number. The point labelled H is a hysteresis point and occurs at an aspect ratio A n = 5.417. A transcritical bifurcation point, labelled T occurs at A T = 5.422. Cubic splines have been used to join the computed points (indicated by diamonds). In order to aid the interpretation of this figure, schematic bifurcation diagrams for four different values of the aspect ratio are shown in Fig. 3. Each diagram shows a measure of the existing steady solutions against the Grashof number using the convention that stable solution branches are marked with a solid line and unstable solutions with a dashed line.
Fig. 1. The smooth evolution to a stable two cell solution for A = 5.5. (a) Gr = 10000, (b) Gr = 20000, (c) Gr = 40000, (d) Gr = 60000, (e) Gr = 100000, and (f) Gr = lg0000.
54,,,~0
i
I
I
I
I
I 5.42
t 5.425
I 5.43
t 5.435
54O0O 53500 '~
53000
~ r~
52600 52O00 51500 5.41
, 5.415
5.44
Aspect ratio Fig. 2. Bifurcation set showing the paths o f saddle-node bifurcations.
Fig. 3a shows the solution structure for increasing Grashof number and A < An. At this aspect ratio and for low Grashof number only one solution branch exists, that is, the one/three cell solution branch. At high Grashof number, there are in addition two further solutions (one of which is stable and one of which is unstable), arising at a saddle-node bifurcation point. The stable solution consists of two corotating cells and is the disconnected solution found by the participants of the GAMM workshop. On increasing the aspect ratio to A = A H, a kink forms in the primary branch, so that for A n < A < A T (Fig. 3b) there are two further saddle-node bifurcations. At this aspect ratio the disconnected branch is qualitatively unchanged. Increasing the aspect ratio to A T = 5.422, the disconnected branch joins the primary branch at a transcritical bifurcation point as shown in Fig. 3c. On increasing the aspect ratio still further the transcritical bifurcation point unfolds so that the two cell solution becomes the primary branch and the one/three cell solution becomes the disconnected branch. Note that the entire exchange mechanism occurs in a very small region of parameter space. This type of exchange mechanism was originally suggested by Benjamin [22], and has been explored experimentally and numerically by Cliffe and Mullin [23] for the Taylor-Couette problem. In general, it is to be expected that increasing the aspect ratio further will change the primary branch to ever
A.C. Skeldon et al./ Journal of Crystal Growth 162 (1996) 95-106
greater numbers of cells through a sequence of similar exchange processes. We now turn to the consideration of Hopf bifurcations that cause the steady bifurcation structure to become unstable. At the GAMM workshop, the participants demonstrated the existence of an oscillatory state characterized by an oscillation between one and three cells. Winters computed the Hopf bifurcation at which this originated. The dangerous eigenvalue method enables us to recognise the existence and, in turn, locate a second Hopf bifurcation along the primary branch for A = 4. The form of the resulting
101
(a)
(a)
4,
't
Gr
(b) Fig. 4. (a) Streamlines and (b) isotherms, for eight equally spaced time intervals in one cycle of the oscillation created in the anti-centrosymmetric Hopf bifurcation at A = 4, Gr = 300370, Pr = 0.015. Gr
(c)
~_
_
time-dependent solution can be constructed, since locally: x--- x0 + e[ gg COS(tOt) +~:I sin(tO/)],
Gr
Gr
tbxee ce}}s
Fig. 3. Qualitative bifurcation diagrams showing the exchange process between two and threc ceUs. Stable solutions are indicated with a solid line and unstable solutions with a dotted line. (a) A < A n , ( b ) A n < A < A - r , ( c ) A = A . r , ( d ) A > A.r.
(8)
where x o is the steady-state solution at the bifurcation point and ~R and ~I are respectively the real and imaginary parts of the eigenvector associated with the Hopf bifurcation point. In Fig. 4a and 4b the streamlines and isotherms are plotted for eight equally spaced points around one cycle: for this and all similar figures presented in this paper time increases anti-clockwise. For the streamlines the perturbation caused by the periodic oscillation has been added to the base flow and shows that the left and the right cells are alternately weak and then strong. For small Prandtl numbers, and the Grashof numbers under consideration here, conduction dominates advection. As a consequence, the perturbation to the temperature field as a result of the periodic motion is small,
102
A.C. Skeldon et aL /Journal of Crystal Growth 162 (1996) 95-106
and it is more illuminating t o view the disturbance to the temperature field alone (Fig. 4b). The two types of oscillation created at the two Hopf bifurcation points have different frequencies and different symmetries. That computed by Winters [ l 1] and others is such that at each point in the cycle the streamlines and the isotherms retain the centrosymmetry of the problem and we refer to this Hopf bifurcation as Hopf c. For the second type of periodic oscillation computed here and shown in Fig. 4 the streamlines and isotherms are not centrosymmetric at any point in time but points in the cycle at time t are related centrosymmetrically with points at t + T / 2 , where T is the period of the oscillation. We refer to this as the anti-centrosymmetric Hopf bifur-
/ I i /
/
./
e/
/
.~ ,"
,t,"
_. . . . .
Gr Fig. 6. One possible Hopf/Hopf unfolding. 600000 -
-
500000
4OOO0O =300000
2ooooo
100000 0
0
I
I
0.005
0.01
I
0.015
I
0.02
I
0.025
0.03
Prandtl number 4OO 350
(b)
I
I
I
I
30O
~. 25o =
200
15o 100 5o o
I
• 0.005
I
I
0.01
0.015
I
0.02
I
0.025
0.03
Prmidtl number Fig. 5. A = 3.5: (a) Bifurcation set showing the paths of the centrosymmelric Hopf bifurcation (Hopf¢) and the anti-centrosymmetric Hopf bifurcation (Hopfa). (b) Frequency Dependence of the two oscillations produced at the Hopf bifurcations.
cation, Hopf~. In both types of oscillations, hotter and colder patches of fluid are advected by the flow. The anti-centrosymmetric Hopf is similar in nature to the second type of periodic oscillation computed by Pulicani et al. [13] for Pr = 0. For A = 4 and Pr = 0.015, the second Hopf bifurcation occurs at a critical Grashof number which is greater than for the first Hopf bifurcation, and the resulting time-dependent solutions are necessarily unstable. However, Fig. 5a shows paths of both types of Hopf bifurcation for varying Prandfl number and for an aspect ratio of A = 3.5. At P r - 0.01 there is a Hopf-Hopf bifurcation point and the order in which the two Hopf bifurcations occur interchanges. In the three-dimensional parameter space spanned by Pr, Gr and A there exists a line of Hopf-Hopf interaction points. Local to such points the dynamical phenomena can be very rich (cf. Van Gils, Krupa and Langford [24] and Golubitsky and Langford [25]): quasi-periodicity and chaos are common features. One typical local scenario,that would explain the results obtained by Pulicani et al., is shown in Fig. 6. Note that although the two Hopf bifurcations are close for Prandtl numbers close to zero, this is no longer true for a Prandtl number of 0.02 relevant to liquid gallium or mercury and the dynamical behaviour observed for Pr = 0 may well be qualitatively different from that for higher Prandtl numbers.
A.C. Skeldon et al./Journal of Crystal Growth 162 (1996) 95-106
Changing the thermal boundary conditions from insulating to conducting does not qualitatively change any of the above results; indeed in the limit of zero Prandtl number conducting and insulating boundaries are equivalent. For small but non-zero Prandtl number conducting boundaries have a destabilising effect and result in lower critical Grashof numbers for the transitions to oscillatory behaviour.
600000
5OOOOO
4 4ooooo r~
200000 1OOOOO
When 0 = 0 in Eq. (1), the gravity vector and the temperature gradient are parallel and, at low Grashof number, a stable solution exists where the fluid is stationary and heat is transferred by conduction alone.
/ /""
,
~,
o~
Gr
/
'~ 3ooooo ,J=
4.2. Tilt p r o b l e m
(b )
103
S
../
////"
,,,¢ .,
j
j'
///
t
\". t.h:: ce'l, - - . . !~'o ~l,ja
Fig. 7. Schematic bifurcation diagrams for A = 4, Pr = 0.015 (a) no tilt (h) small tilt. Computed streamlines and isotherms are shown for typical points, stable branches are indicated with a solid line and unstable branches with a dotted line. In (a) the bifurcation to one cell occurs at Gr = 147000, to two cells at Gr = 193000 and to three cells at Gr = 293 000.
0
0
I 5
I 10
I 15
I 2O
25
Angle of tilt
Fig.
8. Path of saddle-node bifurcations for A = 4, Pr = 0.015 and for increasing tilt.
When the buoyancy force is sufficient to overcome the viscous forces in the fluid, the conduction solution becomes unstable. For the case of tall thin cavities of aspect ratio 4 and a Prandtl number of 0.015, the first transition is to a single convection c e l l ' a t a Grashof number of 147 300. At higher Grashof numbers, the conduction solution is further destabilised to different numbers of stacked counterrotating cells. The first three such transitions are shown in the form of a schematic bifurcation diagram in Fig. 7a. Again, stable solutions are indicated with a solid curve while unstable solutions are shown with a dashed curve. The bifurcations from the conducting solution are pitchfork bifurcations since they break one of the symmetries of the equations. If the bifurcation is to an odd number of cells it breaks the u p / d o w n symmetry (cf. Eq. (3)) of the equations whereas if it is to an even number of cells the bifurcation breaks the centrosymmetry of the equations (cf. Eq. (2)). If the cell is tilted, the equations still have centrosymmetry and hence the transitions to even numbers of cells remain as pitchfork bifurcations. However, the equations no longer have u p / d o w n symmetry and thus the pitchfork bifurcations to odd numbers of cells become disconnected. For small angles of tilt, the pitchfork bifurcations are replaced by a smooth evolution to an oddnumber of cells rotating in one direction and a disconnected solution of cells rotating in the opposite direction, as shown schematically in Fig. 7b. The disconnected
104
A.C. Skeldon et aL /Journal of Crystal Growth 162 (1996) 95-106
solution arises at a saddle-node bifurcation. For a square cavity containing a fluid of unit Prandtl number, Cliffe and Winters [15] traced the saddle-node bifurcation as a function of tilt. They found that the saddle-node bifurcation asymptoted to a fixed angle of 22 ° indicating that the disconnected solution does not exist in the side-wall heated limit. Results of a similar calculation for A = 4 and P r = 0 . 0 1 5 are presented in Fig. 8. This shows that the critical Grashof number for the saddle-node bifurcation again increases sharply with increase in tilt. We have not followed the saddle-node bifurcation above a tilt angle of 8°, but by this stage the critical Grashof number is already 600000. In the side-wall heated
limit the disconnected solutions would consist of convection cells where the flow is up the cold wall and down the hot wall, so it is not surprising that their existence is greatly inhibited. In contrast, many of the other features which we have computed for a shallow cavity heated from the side are manifested in the tall cavity heated from below. Specifically, on increasing the Grashof number in the side-wall heated problem for an aspect ratio of 4 and a Prandtl number of 0.015 the flow evolved smoothly from a one-cell to a three-cell state before being destabilised at the centrosymmetric Hopf bifurcation. A second Hopf bifurcation was found further along the, now unstable, one/three cell solu-
Pf~ . . . . . .
~
Hopf~
Gr Fig. 9. Schematic bifurcation diagram showing the Prandtl number 0.015 and heated from below.
transitions along the one-cell solution branch in a tall thin cavity filled with a fluid of
A.C. Skeldon et al./Journal of Crystal Growth 162 (1996) 95-106
tion branch, this time at the anti-centrosymmetric Hopf bifurcation. For the same Prandtl number we have computed a very similar bifurcation sequence for the one-cell solution branch in the Bdnard problem, as shown schematically in Fig. 9. Streamlines and isotherms are shown for representative points. Fig. 9 illustrates the fact that the one-cell solution smoothly evolves to a three corotating cell solution, this is analogous to the shear instability found in the side-wall heated problem and is to be contrasted with the three counter-rotating cells which bifurcate from the conduction solution in a buoyancy-driven instability (cf. Fig. 7). Along the one/three cell solution branch there are again two Hopf bifurcations with similar structures to the Hopf bifurcations which occur in the side-wall heating problem. As for Fig. 4, the streamlines and isotherms for the resulting oscil-
800000
|
1
i
i
i
t (a~ I~ ~ r
i
i
i
!
Hopfc bifurcations -0-Hopf bifurcations -+--
Lk
105
latory motion are plotted for eight equally spaced points around one cycle. Again, for clarity, only the perturbation to the temperature field is shown. Note that it is now the anti-centrosymmetric oscillatory solution that bifurcates first. The fact that these two Hopf bifurcations are indeed analogous to those found in the side-wall limit is illustrated in Fig. 10a for a Prandtl number of 0.015, where the two curves of Hopf bifurcations have been traced as a function of tilt. The exchange of stability between the two types of Hopf bifurcation occurs at a tilt angle of 0 -- 57 °. As discussed above, this can lead to complicated local dynamical behaviour. For an aspect ratio of 4 in the side-wall heated cavity there exists a disconnected two corotating cell solution arising at a saddle-node bifurcation. The third line drawn on Fig. 10a shows the path of this saddle-node bifurcations for increasing tilt, From this we see that the stable disconnected two cell solution also exists in the Bdnard problem, the streamlines and isotherms in this limit are shown in Fig. 10b. Again the corotating cell solution should be contrasted with the counter-rotating two-cell solution produced in the buoyancy driven instability and shown in Fig. 7.
~" 200000
5. Conclusion
0
I
I
I
1
0
10
20
30
I
1
=
I
I
I
40
50
60
70
80
90
Angle of tilt
(b) f r f J
J P J
r
J
Fig. 10. (a) Bifurcation curves as a function of tilffor Pr = 0.015. The saddle-node bifurcation curve refers to the saddle-node bifurcation which gives rise to two corotating cells. (b) Disconnected two corotating cell solution in the Bdnard limit.
We have shown that in a cavity heated from the side containing a low Prandtl number fluid there is a change in the qualitative form of the primary branch as the aspect ratio is increased from 4 to 5.5. We have computed the exchange mechanism which results in this qualitative transformation of the solution. In general, increasing the aspect ratio further is expected to result in additional changes to the structure of the primary solution. Such changes are of fundamental importance since if a fixed mass of fluid is gradually solidified from one end, even if the Grashof number is low enough so that oscillations do not occur, there could be a sequence of catastrophic changes from one flow state to another. In the case of P r = 0 . 0 1 5 , A = 4 the GAMM workshop and related calculations found that the steady three-cell solution became unstable at a Hopf bifurcation point. We have located a second Hopf bifurcation and calculated its dependence on the
106
A.C. Skeldon et al. /Journal of Crystal Growth 162 (1996) 95-106
Prandtl number. In both cases, increasing the Prandtl number inhibits the oscillation. For low Prandtl number, the second Hopf bifurcation can exchange stability with the first Hopf at a Hopf-Hopf point. This transition helps to explain the origin of the quasiperiodic behaviour computed by Pulicani et al. [13] for Pr = 0. A comparison of the bifurcation structure of the shallow, side-wall heated cavity with convection in a tall thin cavity heated from below has been made. Whilst the buoyancy driven instabilities of the latter case do not occur in the side-wall limit, both problems exhibit similar sequences of shear instabilities. Acknowledgements The authors would like to acknowledge the support of the Engineering and Physical Sciences Research Council (GR/K41311) (A.C.S.) and the Corporate Research Programme of AEA Technology (K.A.C.). The authors are grateful to Dr T. Mullin and Professor D.T.J. Hurle for originally suggesting the problem to them and to Professor K.H. Winters for useful discussions of the problem. We thank the referee for constructive comments on an earlier draft of this paper. References [1] D.T.J. Hurle, Phil. Mag. 13 (1966) 305. [2] A. M[fller and M. Wiehelm, Z. Naturf, A 19 (1964) 254.
[3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13]
[14] [15] [16] [17]
[18] [19] [20] [21] [22] [23] [24] [25]
J. Hart, J. At. Sci. 29 (1972) 687. J. Hart, J. Fluid Mech. 132 (1983) 271. A.E. Gill, J. Fluid Mech. 64 (1974) 577. H.P. Kuo and S.A. Korpela, Phys. Fluids 31 (1988) 33. P. Laure and B. Roux, CRAS Ser. II 305 (1987) 1137. R.V. Birikh, PMM (1966) 30 432. D.E. Cormack, L.G. Leal and J. Imberger, J. Fluid Mech. 65 (1974) 209. B. Roux, Ed., Notes on Numerical Fluid Mechanics, Vol. 27 (1990). K.H. Winters, Int. J. Numer. Methods Eng. 25 (1988) 401. K. Okada and J. Ozoe, J. Crystal Growth 126 (1993) 330. J.P. Pulicani, E.C. Del Arco, A. Randriamampianina, P. Bontoux and R. Peyret, Int. J. Numer. Methods Fluids 10 (1990) 481. E.C. Del Arco, J.P. Pulicani and A. Randriamampianina, CRAS Ser. I1 309 (1989) 1869. K.A. Cliffe and K.H, Winters, J. Comp. Phys. 54 (1984) 531. D.S. Riley and K.H. Winters, J, Fluid Mech. 215 (1990) 309. H.B. Keller, Numerical Solution of Bifurcation and Nonlinear Eigenvalue Problems, in: Applications of Bifurcation Theory, Ed. P.H. Rabinowitz (Academic Press, New York, 1977) p. 359. A. Jepson and A. Spence, SIAM J. Numer. Anal. 22 (1985) 347. A.D. Jepson, PhD Thesis, Part 2, California Institute of Technology ( 1981). A. Griewank and G. Reddien, IMA J. Numer. Anal. 3 (1983) 295. K.A. Cliffe, T.J. GarraU and A. Spence, Ad. Cornp. Math. 1 (1993) 337. T.B. Benjamin, Proc. Roy. Soc. 359 (1978) 1. K.A. Cliffe and T. Mullin, J. Fluid Mech. 153 (1985) 243. S.A. van Gils, M. Krupa, and W.F. Langford, Nonlinearity 3 (1990) 825. M. Golubitsky and W.F. Langford, J. Diff. Eq. 41 (1981) 375.
CRYSTAL G R O W T H
ELSEVIER
Journal of Crystal Growth 162 (1996) 95-106
Convection in a low Prandtl number fluid A.C. Skeldon a,* ,l, D.S. Riley a K.A. Cliffe b a Department of Theoretical Mechanics, University Park, Nottingham NG7 2RD, UK b AEA Technology, B424.4 Harwell Laboratory, Didcot, Oxon OXll ORA, UK
Received 21 August 1995; accepted 13 October 1995
Abstract
As a simple two-dimensional model of convection in the liquid phase during crystal growth using the Bridgman technique, we consider the fluid flow in a shallow rectangular cavity heated from one side. For an aspect ratio of 4 (aspect ratio = length/heigh0, previous numerical studies have shown the existence of two types of oscillatory solution. The first of these has been well-studied for a range of Prandtl numbers from 0 to 0.015. However, the form of the second has been observed only in a time-dependent study for a Prandtl of zero. We locate the new Hopf bifurcation which gives rise to this latter solution and study its dependence on Prandtl number. Before the onset of oscillations, various steady corotating multi-cell solutions are found depending on the aspect ratio of the cavity. We compute how the transition between three corotating cells and two corotating cells takes place for changing aspect ratio. Understanding can be facilitated by comparison with analogous problems. We consider tilting the cavity to the B6nard configuration, where the fluid is heated from below rather than from the side. In this case a no-flow solution exists where heat is transferred by conduction alone. This solution becomes unstable to counter-rotating cells at a symmetry-breaking bifurcation. We study how this symmetrybreaking bifurcation disconnects as the cell is tilted and the solutions evolve into the side-wall cavity solutions. In addition we trace the saddle-node and Hopf bifurcations found in the side-wall heated problem for changing tilt and reveal that they also exist in the B6nard convection limit; disconnected solutions in the B6nard problem have not been studied previously.
1. Introduction
The growth of crystals from a melt is an important process for the semiconductor industry. The crystal quality is determined by many factors, from interfacial dynamics at the molecular level to the effect of melt convection on the bulk distribution of dopants. The difficulty o f the problem leads to the consideration of simplified models which generally concentrate on one aspect: in this paper we focus on
* Corresponding author. Fax: +44 171 477 8000. Present address: Department of Mathematics, City University, Northampton Square, London ECIV 0HB, UK. Elsevier Science B.V. SSDI 0022-0248(95)00923-X
the convective instabilities that occur in the melt. F o l l o w i n g a similar rationale, Hurle [1] studied convection in a shallow boat of rectangular cross-section heated from the side as a simple model of the fluid dynamics associated with the Bridgman process o f crystal growth. He found that, if a small temperature difference was applied across the melt, then the flow was steady. However, on increasing the temperature difference, the flow became oscillatory. A number of experimental studies, for example Miiller and W i e h e l m [2], showed that these oscillations were linked to " s t r i a t i o n s " which appear in the final crystal and which are detrimental to the performance of semi-conductor devices made from the crystal.
96
A.C. Skeldon et aL /Journal of Crystal Growth 162 (1996) 95-106
Much of the theoretical work on this configuration to date focuses on an infinite shallow cavity in which it is assumed that the flow is parallel and that there is an imposed, constant, thermal flux along the cell. The governing equations can then be solved exactly, resulting in a cubic polynomial for the vertical velocity dependence and a quintic polynomial for the corresponding temperature profile. Stability analyses of this "Hadley cell" solution have been performed, first by Hart [3,4] and Gill [5], and later refined by Kuo and Korpela [6] and by Laure and Roux [7]. These have shown that the parallel-flow solution is stable at low Grashof numbers, but is destabilised as the Grashof number is increased. In the case of rigid, insulating upper and lower boundaries and for a Prandtl number less than 0.034, the resulting flow consists of stationary corotating cells whose axes are aligned transverse to the imposed temperature gradient. The critical Grashof number for this instability asymptotes to about 8000 as the Prandtl number is reduced to zero. This transition is a shear instability, analogous to the shear instability of a plane parallel flow found by Birikh [8]. For a Prandtl number greater than 0.034, the first instability of the steady parallel flow is to an oscillatory roll solution, where now the axes of the rolls are parallel to the imposed temperature gradient (longitudinal rolls). A qualitatively similar picture holds if the upper and lower boundaries are both conducting, or if the upper boundary is stress free. Cormack et al. [9] considered a shallow cavity, but with a finite length. They showed that in the central region of this cavity it was reasonable to assume a parallel flow solution providing Gr2Pr2/A3 ~ 3 × 105, where A >> 1 is the aspect ratio of the cavity. Note that for low Prandtl numbers (Pr = 0.02) and for more moderate aspect ratios (A 4) this implies that the parallel solution is a reasonable solution providing Gr ,~ 220 000. Much of the related numerical work has considered rectangular cavities filled with a low Prandtl number fluid and heated from the side. In the main these have concentrated on Prandtl numbers of either 0 or 0.015 and an aspect ratio of 4, which were the values used in a benchmarking exercise at a GAMM workshop [10]. The contributions to this workshop considered rigid lower boundary and end walls, and either a rigid or stress-free surface. Qualitatively
different flows were observed with the different fluid boundary conditions. For rigid walls, the flow is steady and parallel at low Grashof number, but as the Grashof number is increased, this is replaced by a steady solution with three corotating cells. On increasing the Grashof number further, this three-cell solution is in turn replaced by an oscillatory solution at a Hopf bifurcation. During one cycle, the oscillatory solution interchanges between a one-cell state and a three-cell state. Winters [11] traced the path of this Hopf bifurcation and showed that decreasing the aspect ratio or increasing the Prandtl number inhibits the onset of oscillation, causing the critical Grashof number to increase. For the case of zero Prandtl number, the dependence on the aspect ratio was confirmed by Okada and Ozoe [12]. The GAMM workshop also highlighted the presence of a second steady solution consisting of two corotating cells; Winters showed that this arose at a saddle-node bifurcation. Pulicani et al. [13] computed the time-dependent solutions for increasing Gr (see also Ref. [14]). As a result of the GAMM workshop, they conclude that zero Prandtl number calculations are representative of small Prandtl number calculations (Pr--~ 10 -2) and thus restrict their study to this case. This has the advantage that the temperature field is then unaffected by the flow and can be solved analytically, leaving only the fluid equations to be solved numerically. As in the GAMM workshop, they find that, on increasing the Grashof number, there is first a transition from one steady cell to three steady corotating cells, and then a further transition to an oscillatory solution. In addition, they find that on further increase of the Grashof number, this periodic solution is first replaced by a quasi-periodic motion, then by a second type of periodic solution; the quasi-periodic motion exhibits some hysteresis. In this paper we extend the above results. We show that if the aspect ratio is increased to 5.5 then it is the two corotating cell solution which is preferred, that is, it is the two cell solution which is reached by quasi-static increase of the Grashof number from zero. The three-cell solution becomes disconnected and we compute the exchange mechanism to demonstrate how this transition occurs. Winters [1 l] found one path of Hopf bifurcations and mentions in passing the presence of a second: we have computed the
A.C. Skeldon et al./Journal of Crystal Growth 162 (1996)95-106
second Hopf bifurcation in detail. For aspect ratio 4, this new Hopf bifurcation occurs at a critical Grashof number greater than that computed by Winters. However, we trace both types of Hopf bifurcation and show that an exchange of stability between the two Hopf bifurcations can occur. The presence of this second Hopf bifurcation helps to explain the results obtained by Pulicani et al. One question that naturally arises in the present study is: how does the side-wall heating problem relate to the B~nard problem where a fluid layer is heated from below, and, in particular, how does the disconnected branch arise? A transition from one configuration to the other can be made by tilting the cavity. In contrast to the side-wall heating problem, the Btnard problem has a no-flow solution where the fluid is stationary and the heat is transferred by conduction alone. For small Grashof numbers, the conduction solution is stable, however as the Grashof number is increased, this solution becomes unstable through a pitchfork bifurcation producing a number of counter-rotating cells. This transition is buoyancy driven and the number of cells produced is determined by the aspect ratio. Cliffe and Winters [15] showed that as a square cell is tilted, the leading pitchfork bifurcation to a unicellular flow unfolds: there is a smooth transition to a stable unicellular flow representing the natural circulation. A disconnected solution arises at a saddle-node bifurcation with a stable branch representing anomalous un-natural circulation. Tracing the saddle-node bifurcation for increasing tilt showed that for the case Pr = 1 it asymptoted to a fixed angle of tilt, 22 ° and consequently, no analogous solutions exist in the side-wall heating problem. Riley and Winters [16] showed that in a similar problem where the cavity is filled with a porous medium (the Lapwood problem), the disconnected solution branch again does not exist in the side-wall limit, although they found a more complicated sequence of transitions involving isola formations. These studies motivated us to look further at the relationship between the B~nard problem and the side-wall heating problem. We compute the bifurcation structure for a tall thin cavity filled with a low Prandtl number fluid and heated from below and show how this structure is altered as the cavity is tilted. In addition we trace the two types of Hopf bifurcation and the saddle-node bifurcation which
97
gives rise to stable corotating cells in the side-wall heated problem for changing tilt. In the next section the equations of motion for the problem are stated. The numerical methods used are discussed in Section 3 which then leads to the results, presented in Section 4. The results are discussed and our conclusions drawn in Section 5.
2. Problem formulation A rectangular cavity of height H and length L filled with a low Prandtl number fluid is considered. One side of the cavity is maintained at a constant temperature T h, whilst the opposite side is kept at a lower temperature, Tc. The remaining two boundaries are assumed insulating. All the walls are taken to be impermeable and non-slip. In order to relate the side-wall heating problem to the B~nard problem, the angle between the side of the cell and the gravity vector, 0 is defined such that when 0 = 0 the gravity vector and the temperature gradient are parallel,the Btnard configuration; when 0 = 7r/2, the side-wall heating problem is recovered. We assume that the convection is described by the Navier-Stokes and energy equations subject to the Boussinesq approximation and non-dimensionalise using H, G r l / 2 o / H , H Z / v , T h - T~, Grpov2/H 2 as the length,velocity, time, temperature and pressure scales respectively where Gr = g ~ H 3 (T h - Tc)/v 2. This results in the four non-dimensional equations:
Ou Ov --+--=0, Ox Oy Ou Ot + Grl/2
( OU v OU I U~x + -~y]
02u 02u am = -Gr~/2~x +-+-+ Grf/2Tcos 0 OX 2
Oy 2
- - + Gr ~/2 u - - + Ot Ox -~y ] Op 02 v 02v = -- Grl/20---y + --Ox2 + --Oy2 + Grl/ZT sin O, --+Gr
Ot
~/2 u - - + v - -
Ox
Oy
=
--+
Ox 2
cgy2 J"
A.C. Skeldon et aL /Journal of Crystal Growth 162 (1996) 95-106
98
Along the non-slip and insulating boundaries:
u=v=O,
3T --=0,
Oy
+fp ~
whilst on the isothermal and non-slip boundaries:
u=v=O,
T= - ±2 '
on the cold side, and u=v=0,
T=½,
on the opposite warm side. When the angle of tilt is non-zero the equations of motion together with the boundary conditions result in a centrosymmetric problem in that, if the origin of our coordinate system is taken as the centre of the cell, the problem is invariant under the transformation: x~-x, y~-y, u~-u, v--*-v, T-~ - r .
(2)
When the angle of tilt is zero, that is, when the fluid is heated from below, the problem is in addition invariant under the reflection symmetry:
x--~ - x ,
y--~y,
u-~ - u ,
v--~v,
T--~ - T .
(3)
3. Numerical methods
We consider the pressure-integrated weak form for Eq. (1), that is:
L
+Grl/2
+ 7y)
- Gr~/2Tcos O]
+
dxdy=0.
Where, f , , f~, f r and fp are test functions. These test functions are chosen so that there is no contribution from the boundary integral. A standard Galerkin formulation of the finite element method was adopted in which the test functions were chosen from the same set of functions as the basis functions. The temperature and velocities were represented as a sum of quadratic basis functions whilst linear discontinuous basis functions were used for the pressure. The rectangular domain is divided into rectangular nine-node elements. In order to capture any rapid changes which occur at the boundaries of the cavity, the spacing of the nodes is graduated so that there are more elements near the walls. Grids of 30 X 12, 40 X 16 and 60 X 24 elements were used. The coarsest grid gave qualitatively correct results, however, there was a large variation in the accuracy to which bifurcation points were found, depending on the values of the parameters and the particular type of bifurcation point studied. At best, the coarsest grid gave critical parameter values which were to within 1% of those found by using the finest grid, but at worst it computed the value of the critical parameters only to within 40%. In particular, accurate location of the second type of Hopf bifurcation at Prandtl numbers of above 0.01 was only achieved with the finer grids. Most of the results presented below were computed using the 60 X 24 element grid. The resulting discretised equations are of the form:
Ox
M - - + f ( x , , ~ ) --0. -
Grl/2POf"~x + Of. Ou Ox Ox
--Gr~/2Dox
+
Of. Ou Oy Oy
aL a v + af~ av Ox ax Oy ay
+--
dt
+ GrV~
Ox
(5)
Here, M is the mass matrix, x is the solution vector and t~ is a vector of the control parameters, that is, the Prandtl number, the Grashof number, the aspect ratio and the tilt angle. We seek steady-state solutions to this problem by solving:
f ( x , a ) =0,
+A ~
(4)
(6)
using Newton's method, which is guaranteed to converge providing a good enough initial guess to the
A.C. Skeldon et al./ Journal of Crystal Growth 162 (1996) 95-106
solution is made. For low Grashof numbers, we find x - - 0 is a sufficiently good guess, but to reach higher Grashof numbers, a path-following approach is employed. That is, the Grashof number is incremented in small steps; at each step Newton's method is used to find a solution, using the solution at the previous step to give an initial estimate for the new solution. Both stable and unstable solutions can be followed, but this method fails at saddle-nc~le points. This, however, is readily overcome by parameterising the solution branch with arclength, as originally suggested by Keller [17]. Only stable solutions are physically realisable and thus monitoring changes of stability of the solutions is also important. Generically, changes of stability occur either when a single real eigenvalue passes through zero or when a complex conjugate pair of eigenvalues crosses the imaginary axis. In the former case, a steady-state bifurcation occurs and in the latter, a Hopf bifurcation. In each case an "extended" system of necessary conditions can be solved to locate the bifurcation point. For example, in order to find a saddle-node bifurcation we follow Jepson and Spence [18], and solve the extended system: =0,
99
along the branch of solutions the bifurcation parameter increases and then decreases (or vice versa), and secondly there is a change in the number of eigenvalues with positive real part. As the determinant of the Jacobian equals the product of the eigenvalues, this second feature is signalled by a corresponding change in the sign of determinant. Furthermore, since we compute the Jacobian in order to use Newton's method and calculating the determinant is computationally inexpensive, this second feature can be checked with little additional computational time. When a Hopf bifurcation occurs a complex conjugate pair of eigenvalues crosses the imaginary axis. Hence, although a change of stability may occur, there is no corresponding change in the sign of the determinant of the Jacobian. Computing all the eigenvalues of the Jacobian is computationally intensive, prohibitively so in this case. Instead we use a new method developed by Cliffe et al. [21]. Rather than compute the entire spectrum of eigenvalues, this method seeks to compute only the "most dangerous" i.e. those of smallest real part: it is these eigenvalues which determine when bifurcations occur. By tracking these eigenvalues a good estimate of where Hopf bifurcations occur can be made.
= o,
l(~b) - 1 = 0.
(7)
Here, ~b is the eigenvector of the Jacobian matrix, /x(x0,~t), and the last equation in Eq, (7) is a normalisation condition for ~. There is a corresponding extended system which gives necessary conditions for a Hopf bifurcation to occur (cf. Jepson [19] and Griewank and Reddien [20]). In both cases the extended systems are solved using Newton's method and paths of bifurcation points are followed in parameter space in the same way that paths of regular solution points are followed. As for the location of regular solution points, in order for Newton's method to converge to a solution of the extended system, a good enough initial guess for the bifurcation is required. Our experience suggests that in practice, the basin of convergence can be very small and hence a very good guess is needed. This does not pose a problem for the location of saddle-node bifurcations since these have two readily monitored characteristics. Firstly, stepping
4. Results 4.1. Side-wall heating problem
Consistent with previous studies, at small Grashof number and A = 4, we find that the solution is parallel over a large percentage of the cavity. As the Grashof number is increased, the region over which the flow is parallel gradually shrinks; further increasing the Grashof number results in the growth of two smaller cells at either end of the cavity and gives rise to the existence of three corotating cells. These transitions are the remnant of the steady transverse shear instability found by Hart. The parallel-flow solution in an infinite layer has a continuous translational symmetry. At Hart's instability this symmetry is broken and the resulting stable corotating cells have only a discrete translational symmetry. The finite cavity does not have the translational symmetry of the infinite layer, and hence this bifurcation is
100
A.C. Skeldon et aL /Journal of Crystal Growth 162 (1996) 95-106
unfolded into a smooth evolution to corotating cells, where the cells rotate such that the flow passes up the hot wall and down the cold wall. If the aspect ratio is increased to 5.5, instead of a transition from a parallel flow to one cell and then to three cells, there is a smooth evolution to two cells, as shown in Fig. 1 - this is qualitatively different from the behaviour at lower aspect ratio. We have computed the mechanism by which this change occurs and the results are summarised in Fig. 2. This figure displays a bifurcation set showing the computed path of saddle-node bifurcation as a function of aspect ratio and Grashof number. The point labelled H is a hysteresis point and occurs at an aspect ratio A n = 5.417. A transcritical bifurcation point, labelled T occurs at A T = 5.422. Cubic splines have been used to join the computed points (indicated by diamonds). In order to aid the interpretation of this figure, schematic bifurcation diagrams for four different values of the aspect ratio are shown in Fig. 3. Each diagram shows a measure of the existing steady solutions against the Grashof number using the convention that stable solution branches are marked with a solid line and unstable solutions with a dashed line.
Fig. 1. The smooth evolution to a stable two cell solution for A = 5.5. (a) Gr = 10000, (b) Gr = 20000, (c) Gr = 40000, (d) Gr = 60000, (e) Gr = 100000, and (f) Gr = lg0000.
54,,,~0
i
I
I
I
I
I 5.42
t 5.425
I 5.43
t 5.435
54O0O 53500 '~
53000
~ r~
52600 52O00 51500 5.41
, 5.415
5.44
Aspect ratio Fig. 2. Bifurcation set showing the paths o f saddle-node bifurcations.
Fig. 3a shows the solution structure for increasing Grashof number and A < An. At this aspect ratio and for low Grashof number only one solution branch exists, that is, the one/three cell solution branch. At high Grashof number, there are in addition two further solutions (one of which is stable and one of which is unstable), arising at a saddle-node bifurcation point. The stable solution consists of two corotating cells and is the disconnected solution found by the participants of the GAMM workshop. On increasing the aspect ratio to A = A H, a kink forms in the primary branch, so that for A n < A < A T (Fig. 3b) there are two further saddle-node bifurcations. At this aspect ratio the disconnected branch is qualitatively unchanged. Increasing the aspect ratio to A T = 5.422, the disconnected branch joins the primary branch at a transcritical bifurcation point as shown in Fig. 3c. On increasing the aspect ratio still further the transcritical bifurcation point unfolds so that the two cell solution becomes the primary branch and the one/three cell solution becomes the disconnected branch. Note that the entire exchange mechanism occurs in a very small region of parameter space. This type of exchange mechanism was originally suggested by Benjamin [22], and has been explored experimentally and numerically by Cliffe and Mullin [23] for the Taylor-Couette problem. In general, it is to be expected that increasing the aspect ratio further will change the primary branch to ever
A.C. Skeldon et al./ Journal of Crystal Growth 162 (1996) 95-106
greater numbers of cells through a sequence of similar exchange processes. We now turn to the consideration of Hopf bifurcations that cause the steady bifurcation structure to become unstable. At the GAMM workshop, the participants demonstrated the existence of an oscillatory state characterized by an oscillation between one and three cells. Winters computed the Hopf bifurcation at which this originated. The dangerous eigenvalue method enables us to recognise the existence and, in turn, locate a second Hopf bifurcation along the primary branch for A = 4. The form of the resulting
101
(a)
(a)
4,
't
Gr
(b) Fig. 4. (a) Streamlines and (b) isotherms, for eight equally spaced time intervals in one cycle of the oscillation created in the anti-centrosymmetric Hopf bifurcation at A = 4, Gr = 300370, Pr = 0.015. Gr
(c)
~_
_
time-dependent solution can be constructed, since locally: x--- x0 + e[ gg COS(tOt) +~:I sin(tO/)],
Gr
Gr
tbxee ce}}s
Fig. 3. Qualitative bifurcation diagrams showing the exchange process between two and threc ceUs. Stable solutions are indicated with a solid line and unstable solutions with a dotted line. (a) A < A n , ( b ) A n < A < A - r , ( c ) A = A . r , ( d ) A > A.r.
(8)
where x o is the steady-state solution at the bifurcation point and ~R and ~I are respectively the real and imaginary parts of the eigenvector associated with the Hopf bifurcation point. In Fig. 4a and 4b the streamlines and isotherms are plotted for eight equally spaced points around one cycle: for this and all similar figures presented in this paper time increases anti-clockwise. For the streamlines the perturbation caused by the periodic oscillation has been added to the base flow and shows that the left and the right cells are alternately weak and then strong. For small Prandtl numbers, and the Grashof numbers under consideration here, conduction dominates advection. As a consequence, the perturbation to the temperature field as a result of the periodic motion is small,
102
A.C. Skeldon et aL /Journal of Crystal Growth 162 (1996) 95-106
and it is more illuminating t o view the disturbance to the temperature field alone (Fig. 4b). The two types of oscillation created at the two Hopf bifurcation points have different frequencies and different symmetries. That computed by Winters [ l 1] and others is such that at each point in the cycle the streamlines and the isotherms retain the centrosymmetry of the problem and we refer to this Hopf bifurcation as Hopf c. For the second type of periodic oscillation computed here and shown in Fig. 4 the streamlines and isotherms are not centrosymmetric at any point in time but points in the cycle at time t are related centrosymmetrically with points at t + T / 2 , where T is the period of the oscillation. We refer to this as the anti-centrosymmetric Hopf bifur-
/ I i /
/
./
e/
/
.~ ,"
,t,"
_. . . . .
Gr Fig. 6. One possible Hopf/Hopf unfolding. 600000 -
-
500000
4OOO0O =300000
2ooooo
100000 0
0
I
I
0.005
0.01
I
0.015
I
0.02
I
0.025
0.03
Prandtl number 4OO 350
(b)
I
I
I
I
30O
~. 25o =
200
15o 100 5o o
I
• 0.005
I
I
0.01
0.015
I
0.02
I
0.025
0.03
Prmidtl number Fig. 5. A = 3.5: (a) Bifurcation set showing the paths of the centrosymmelric Hopf bifurcation (Hopf¢) and the anti-centrosymmetric Hopf bifurcation (Hopfa). (b) Frequency Dependence of the two oscillations produced at the Hopf bifurcations.
cation, Hopf~. In both types of oscillations, hotter and colder patches of fluid are advected by the flow. The anti-centrosymmetric Hopf is similar in nature to the second type of periodic oscillation computed by Pulicani et al. [13] for Pr = 0. For A = 4 and Pr = 0.015, the second Hopf bifurcation occurs at a critical Grashof number which is greater than for the first Hopf bifurcation, and the resulting time-dependent solutions are necessarily unstable. However, Fig. 5a shows paths of both types of Hopf bifurcation for varying Prandfl number and for an aspect ratio of A = 3.5. At P r - 0.01 there is a Hopf-Hopf bifurcation point and the order in which the two Hopf bifurcations occur interchanges. In the three-dimensional parameter space spanned by Pr, Gr and A there exists a line of Hopf-Hopf interaction points. Local to such points the dynamical phenomena can be very rich (cf. Van Gils, Krupa and Langford [24] and Golubitsky and Langford [25]): quasi-periodicity and chaos are common features. One typical local scenario,that would explain the results obtained by Pulicani et al., is shown in Fig. 6. Note that although the two Hopf bifurcations are close for Prandtl numbers close to zero, this is no longer true for a Prandtl number of 0.02 relevant to liquid gallium or mercury and the dynamical behaviour observed for Pr = 0 may well be qualitatively different from that for higher Prandtl numbers.
A.C. Skeldon et al./Journal of Crystal Growth 162 (1996) 95-106
Changing the thermal boundary conditions from insulating to conducting does not qualitatively change any of the above results; indeed in the limit of zero Prandtl number conducting and insulating boundaries are equivalent. For small but non-zero Prandtl number conducting boundaries have a destabilising effect and result in lower critical Grashof numbers for the transitions to oscillatory behaviour.
600000
5OOOOO
4 4ooooo r~
200000 1OOOOO
When 0 = 0 in Eq. (1), the gravity vector and the temperature gradient are parallel and, at low Grashof number, a stable solution exists where the fluid is stationary and heat is transferred by conduction alone.
/ /""
,
~,
o~
Gr
/
'~ 3ooooo ,J=
4.2. Tilt p r o b l e m
(b )
103
S
../
////"
,,,¢ .,
j
j'
///
t
\". t.h:: ce'l, - - . . !~'o ~l,ja
Fig. 7. Schematic bifurcation diagrams for A = 4, Pr = 0.015 (a) no tilt (h) small tilt. Computed streamlines and isotherms are shown for typical points, stable branches are indicated with a solid line and unstable branches with a dotted line. In (a) the bifurcation to one cell occurs at Gr = 147000, to two cells at Gr = 193000 and to three cells at Gr = 293 000.
0
0
I 5
I 10
I 15
I 2O
25
Angle of tilt
Fig.
8. Path of saddle-node bifurcations for A = 4, Pr = 0.015 and for increasing tilt.
When the buoyancy force is sufficient to overcome the viscous forces in the fluid, the conduction solution becomes unstable. For the case of tall thin cavities of aspect ratio 4 and a Prandtl number of 0.015, the first transition is to a single convection c e l l ' a t a Grashof number of 147 300. At higher Grashof numbers, the conduction solution is further destabilised to different numbers of stacked counterrotating cells. The first three such transitions are shown in the form of a schematic bifurcation diagram in Fig. 7a. Again, stable solutions are indicated with a solid curve while unstable solutions are shown with a dashed curve. The bifurcations from the conducting solution are pitchfork bifurcations since they break one of the symmetries of the equations. If the bifurcation is to an odd number of cells it breaks the u p / d o w n symmetry (cf. Eq. (3)) of the equations whereas if it is to an even number of cells the bifurcation breaks the centrosymmetry of the equations (cf. Eq. (2)). If the cell is tilted, the equations still have centrosymmetry and hence the transitions to even numbers of cells remain as pitchfork bifurcations. However, the equations no longer have u p / d o w n symmetry and thus the pitchfork bifurcations to odd numbers of cells become disconnected. For small angles of tilt, the pitchfork bifurcations are replaced by a smooth evolution to an oddnumber of cells rotating in one direction and a disconnected solution of cells rotating in the opposite direction, as shown schematically in Fig. 7b. The disconnected
104
A.C. Skeldon et aL /Journal of Crystal Growth 162 (1996) 95-106
solution arises at a saddle-node bifurcation. For a square cavity containing a fluid of unit Prandtl number, Cliffe and Winters [15] traced the saddle-node bifurcation as a function of tilt. They found that the saddle-node bifurcation asymptoted to a fixed angle of 22 ° indicating that the disconnected solution does not exist in the side-wall heated limit. Results of a similar calculation for A = 4 and P r = 0 . 0 1 5 are presented in Fig. 8. This shows that the critical Grashof number for the saddle-node bifurcation again increases sharply with increase in tilt. We have not followed the saddle-node bifurcation above a tilt angle of 8°, but by this stage the critical Grashof number is already 600000. In the side-wall heated
limit the disconnected solutions would consist of convection cells where the flow is up the cold wall and down the hot wall, so it is not surprising that their existence is greatly inhibited. In contrast, many of the other features which we have computed for a shallow cavity heated from the side are manifested in the tall cavity heated from below. Specifically, on increasing the Grashof number in the side-wall heated problem for an aspect ratio of 4 and a Prandtl number of 0.015 the flow evolved smoothly from a one-cell to a three-cell state before being destabilised at the centrosymmetric Hopf bifurcation. A second Hopf bifurcation was found further along the, now unstable, one/three cell solu-
Pf~ . . . . . .
~
Hopf~
Gr Fig. 9. Schematic bifurcation diagram showing the Prandtl number 0.015 and heated from below.
transitions along the one-cell solution branch in a tall thin cavity filled with a fluid of
A.C. Skeldon et al./Journal of Crystal Growth 162 (1996) 95-106
tion branch, this time at the anti-centrosymmetric Hopf bifurcation. For the same Prandtl number we have computed a very similar bifurcation sequence for the one-cell solution branch in the Bdnard problem, as shown schematically in Fig. 9. Streamlines and isotherms are shown for representative points. Fig. 9 illustrates the fact that the one-cell solution smoothly evolves to a three corotating cell solution, this is analogous to the shear instability found in the side-wall heated problem and is to be contrasted with the three counter-rotating cells which bifurcate from the conduction solution in a buoyancy-driven instability (cf. Fig. 7). Along the one/three cell solution branch there are again two Hopf bifurcations with similar structures to the Hopf bifurcations which occur in the side-wall heating problem. As for Fig. 4, the streamlines and isotherms for the resulting oscil-
800000
|
1
i
i
i
t (a~ I~ ~ r
i
i
i
!
Hopfc bifurcations -0-Hopf bifurcations -+--
Lk
105
latory motion are plotted for eight equally spaced points around one cycle. Again, for clarity, only the perturbation to the temperature field is shown. Note that it is now the anti-centrosymmetric oscillatory solution that bifurcates first. The fact that these two Hopf bifurcations are indeed analogous to those found in the side-wall limit is illustrated in Fig. 10a for a Prandtl number of 0.015, where the two curves of Hopf bifurcations have been traced as a function of tilt. The exchange of stability between the two types of Hopf bifurcation occurs at a tilt angle of 0 -- 57 °. As discussed above, this can lead to complicated local dynamical behaviour. For an aspect ratio of 4 in the side-wall heated cavity there exists a disconnected two corotating cell solution arising at a saddle-node bifurcation. The third line drawn on Fig. 10a shows the path of this saddle-node bifurcations for increasing tilt, From this we see that the stable disconnected two cell solution also exists in the Bdnard problem, the streamlines and isotherms in this limit are shown in Fig. 10b. Again the corotating cell solution should be contrasted with the counter-rotating two-cell solution produced in the buoyancy driven instability and shown in Fig. 7.
~" 200000
5. Conclusion
0
I
I
I
1
0
10
20
30
I
1
=
I
I
I
40
50
60
70
80
90
Angle of tilt
(b) f r f J
J P J
r
J
Fig. 10. (a) Bifurcation curves as a function of tilffor Pr = 0.015. The saddle-node bifurcation curve refers to the saddle-node bifurcation which gives rise to two corotating cells. (b) Disconnected two corotating cell solution in the Bdnard limit.
We have shown that in a cavity heated from the side containing a low Prandtl number fluid there is a change in the qualitative form of the primary branch as the aspect ratio is increased from 4 to 5.5. We have computed the exchange mechanism which results in this qualitative transformation of the solution. In general, increasing the aspect ratio further is expected to result in additional changes to the structure of the primary solution. Such changes are of fundamental importance since if a fixed mass of fluid is gradually solidified from one end, even if the Grashof number is low enough so that oscillations do not occur, there could be a sequence of catastrophic changes from one flow state to another. In the case of P r = 0 . 0 1 5 , A = 4 the GAMM workshop and related calculations found that the steady three-cell solution became unstable at a Hopf bifurcation point. We have located a second Hopf bifurcation and calculated its dependence on the
106
A.C. Skeldon et al. /Journal of Crystal Growth 162 (1996) 95-106
Prandtl number. In both cases, increasing the Prandtl number inhibits the oscillation. For low Prandtl number, the second Hopf bifurcation can exchange stability with the first Hopf at a Hopf-Hopf point. This transition helps to explain the origin of the quasiperiodic behaviour computed by Pulicani et al. [13] for Pr = 0. A comparison of the bifurcation structure of the shallow, side-wall heated cavity with convection in a tall thin cavity heated from below has been made. Whilst the buoyancy driven instabilities of the latter case do not occur in the side-wall limit, both problems exhibit similar sequences of shear instabilities. Acknowledgements The authors would like to acknowledge the support of the Engineering and Physical Sciences Research Council (GR/K41311) (A.C.S.) and the Corporate Research Programme of AEA Technology (K.A.C.). The authors are grateful to Dr T. Mullin and Professor D.T.J. Hurle for originally suggesting the problem to them and to Professor K.H. Winters for useful discussions of the problem. We thank the referee for constructive comments on an earlier draft of this paper. References [1] D.T.J. Hurle, Phil. Mag. 13 (1966) 305. [2] A. M[fller and M. Wiehelm, Z. Naturf, A 19 (1964) 254.
[3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13]
[14] [15] [16] [17]
[18] [19] [20] [21] [22] [23] [24] [25]
J. Hart, J. At. Sci. 29 (1972) 687. J. Hart, J. Fluid Mech. 132 (1983) 271. A.E. Gill, J. Fluid Mech. 64 (1974) 577. H.P. Kuo and S.A. Korpela, Phys. Fluids 31 (1988) 33. P. Laure and B. Roux, CRAS Ser. II 305 (1987) 1137. R.V. Birikh, PMM (1966) 30 432. D.E. Cormack, L.G. Leal and J. Imberger, J. Fluid Mech. 65 (1974) 209. B. Roux, Ed., Notes on Numerical Fluid Mechanics, Vol. 27 (1990). K.H. Winters, Int. J. Numer. Methods Eng. 25 (1988) 401. K. Okada and J. Ozoe, J. Crystal Growth 126 (1993) 330. J.P. Pulicani, E.C. Del Arco, A. Randriamampianina, P. Bontoux and R. Peyret, Int. J. Numer. Methods Fluids 10 (1990) 481. E.C. Del Arco, J.P. Pulicani and A. Randriamampianina, CRAS Ser. I1 309 (1989) 1869. K.A. Cliffe and K.H, Winters, J. Comp. Phys. 54 (1984) 531. D.S. Riley and K.H. Winters, J, Fluid Mech. 215 (1990) 309. H.B. Keller, Numerical Solution of Bifurcation and Nonlinear Eigenvalue Problems, in: Applications of Bifurcation Theory, Ed. P.H. Rabinowitz (Academic Press, New York, 1977) p. 359. A. Jepson and A. Spence, SIAM J. Numer. Anal. 22 (1985) 347. A.D. Jepson, PhD Thesis, Part 2, California Institute of Technology ( 1981). A. Griewank and G. Reddien, IMA J. Numer. Anal. 3 (1983) 295. K.A. Cliffe, T.J. GarraU and A. Spence, Ad. Cornp. Math. 1 (1993) 337. T.B. Benjamin, Proc. Roy. Soc. 359 (1978) 1. K.A. Cliffe and T. Mullin, J. Fluid Mech. 153 (1985) 243. S.A. van Gils, M. Krupa, and W.F. Langford, Nonlinearity 3 (1990) 825. M. Golubitsky and W.F. Langford, J. Diff. Eq. 41 (1981) 375.