Three-dimensional bifurcations of a two-phase Rayleigh–Benard problem in a cylinder

Three-dimensional bifurcations of a two-phase

Rayleigh±Benard problem in a cylinder

C.W. Lan

_{, C.H. Wang}

Department of Chemical Engineering, National Taiwan University, Taipei 10617, Taiwan, ROC Received 3 March 2000; received in revised form 26 July 2000

Abstract

Three-dimensional (3D) bifurcations of a partially melted or solidi®ed material in a cylinder heated from below are studied numerically. Through nonlinear calculations, bifurcation diagrams are constructed for a melt of a Prandtl number of one. As the interface is ®xed, our calculated results agree reasonably well with previous calculations, but some discrepancies exist, which are further discussed through their dynamic evolutions and imperfect bifurcations of 5° tilt. As the interface is allowed to deform, the bifurcation behavior changes signi®cantly, both for the onset of con-vection and its concon-vection mode. For the initial melt aspect ratio of one, the primary bifurcation changes from su-percritical to subcritical with the increasing solid amount, and the onset mode from an axisymmetric (m0) mode to a 3D (m1) mode. Although the free interface destabilizes the conductive mode and leads to an earlier onset of convection, it may stabilize some ¯ow modes through its con®nement. Imperfect bifurcations due to a 5° tilt are further illus-trated. Ó 2001 Elsevier Science Ltd. All rights reserved.

Keywords: Three-dimensional; Bifurcation; Rayleigh±Benard; Buoyancy ¯ow; Interface

1. Introduction

In many physical problems in nature or industrial processes, thermal convection is an important mech-anism for heat and mass transfer. Particularly, as a ¯uid is heated from below, its nonlinear behavior is very rich and interesting in many scienti®c ®elds. The classic Rayleigh±Benard (RB) problem o€ers a ®rst approach to complex ¯ows as well as the transition from con-ductive to convective heat transfer modes. In many physical processes, the RB problem may further be as-sociated with freezing or melting, where the onset of thermal convection and the ¯ow modes are coupled with the deformable melt/solid interface. Such is the case for the storage of energy using the phase change of material [1,2]. In the solidi®cation processing of alloys and the growth of single crystals [3±5], it is also particularly important. For example, during crystal growth of

elec-tronic materials, the onset mode could a€ect the inter-face morphology and dopant segregation, which are crucial to crystal quality [3].

In engineering applications, the material or ¯uid is usually con®ned in a container, e.g., a box or a cylinder. A generic situation in a cylinder is shown in Fig. 1(a), where the ¯uid is heated from below. Because of the friction from the solid wall, the onset of convection and the ¯ow bifurcations depend on the boundary ditions. Indeed, as compared with the buoyancy con-vection between two parallel plates of in®nite extent, the bounded con®guration has received much less attention because of the diculty of its analytical and numerical treatment. In fact, in many engineering applications, such as crystal growth or solidi®cation in an ampoule, this con®guration is important. For this problem with only the ¯uid inside, its bifurcations (primary) have been well studied for the axisymmetric modes since the work by Liang et al. [6] and Charlson and Sani [7,8]. Detailed bifurcation diagrams were established by Yamaguchi et al. [9] through a fully nonlinear ®nite element analysis. Due to the limitation of the computation resources, its three-dimensional (3D) bifurcations were not

www.elsevier.com/locate/ijhmt

*_{Corresponding author. Tel.: 2363-3917; fax:}

+886-2-2363-3917.

E-mail address: [email protected] (C.W. Lan).

0017-9310/01/$ - see front matter Ó 2001 Elsevier Science Ltd. All rights reserved. PII: S 0 0 1 7 - 9 3 1 0 ( 0 0 ) 0 0 2 4 8 - 9

investigated until recently. Muller et al. [10] studied the ¯ow structures in a cylinder for a liquid metal (Prandtl number Pr ˆ 0:02) and water (Pr ˆ 6:7) both exper-imentally and numerically, and found that for water, the axisymmetric mode can be stable up to 10Rac, where Rac

is the primary critical Rayleigh number. However, a secondary bifurcation was found at Ra ˆ 2800 with the 3D solutions having an m2 symmetry. A latter study by Neumann [11] also obtained similar results but with more careful numerical calculations, where the mesh e€ects were discussed. Interestingly, Hardin and Sani [12] used a weakly nonlinear analysis to show, di€er-ently, that the m2 mode branched from the m0 mode supercritically at Racˆ 2800, and beyond this secondary

bifurcation the m0 mode became unstable. Such a con-troversy was commented on by Wagner et al. [13] who showed that the appearance of the m2 mode at Ra ˆ 2800 was possibly due to the numerical error of a too-coarse grid used by Neumann [11]. Nevertheless, later Wanschura et al. [14] used a mixed ®nite-di€erence/ Chebyshev collocation method to study the same prob-lem for the aspect ratio ranging from 0.9 to 1.56 and for Pr ˆ 0:02 and 1, and commented that while the m2 mode might appear but suspected that there was a small range

near Ra ˆ 2800 where the m2 mode was not stable. Despite the fact that the stability of the m0 and m2 modes are still unclear, the bifurcation structure pro-posed by Hardin and Sani [12] is perhaps the most complete to date for Pr ˆ 6:7. It is all the same dicult to make any judgment about the stability of the modes because the supercritical bifurcations are not structurally stable. A little imperfection, such as tilt or heat loss, can remove these points. Furthermore, these points are highly sensitive to numerical accuracy. So far, the re-ported results are still not mesh independent. At least, no one has proven that. Therefore, in addition to the perfect bifurcations with the ideal boundary conditions, imper-fect bifurcations due to symmetry breaking also provide useful information of which modes being more stable. However, to the best of our knowledge no such study has been conducted on this problem. The energy analysis [14] is also a useful way to investigate the stabilizing (or destabilizing) mechanisms. However, obtaining a basic state is a numerical problem. In addition to numerical accuracy, the branching structure is still dicult to ®nd. Newton's method with a continuation technique [15±17] has been proven to be an ecient way for bifurcation analyses, but it is extremely costly for 3D problems. Nomenclature

A melt aspect ratio

A surface area vector in Fig. 2(b). Cp speci®c heat

e1 unit vector in x1-direction

eg unit vector in the gravity direction

f residual vector g gravitational acceleration h melt height DH heat of fusion H domain height I ¯uxes k thermal conductivity L length

n unit normal vector Nu Nusselt number P pressure

Pr Prandtl number, vm=am

R radius

Ra thermal Rayleigh number, gbTR3…THÿ Tm†=

mmam

s arclength S surface area

S/ source term of variable /

St Stefan number, DH=Cpm…THÿ Tm†

T temperature

TC cold zone temperature

TH hot zone temperature

Tm melting point u x1_{-component of velocity} v x2_{-component of velocity} DV control volume v velocity vector w x3_{-component of velocity} x variable vector y variable vector …x; Ra†T

Greek symbols

am thermal di€usivity of melt

bT thermal expansion coecient

ni _{coordinate in computational domain}

c tilt angle

j thermal conductivity ratio …k=km†

h dimensionless temperature s dimensionless time sij singularity test function

mm kinematic viscosity

Superscripts

1 x-direction in the Cartesian coordinate 2 y-component in the Cartesian coordinate 3 z-component in the Cartesian coordinate  reduced value

Subscripts

c critical value s solid m melt

As a materials is partially melted or solidi®ed in a cylinder, as shown in Fig. 1(b), the bifurcation behavior for a two-phase RB problem could be quite di€erent. After the onset of convection, the interface shape is a coupled solution with the bifurcated mode, as show in Fig. 1(c). Such a problem is particularly interesting and important to solidi®cation and crystal growth, as well as in the applications toward energy storage problems using phase change materials. In fact, the coupled RB convection with solidi®cation between two plates with in®nitely extended domain was studied for the ®rst time by Davis et al. [18]. They found that a subcritical bi-furcation was associated with the onset of hexagonal cells for a thick solid layer. Limited by their weakly nonlinear analysis, their study was available only very near the static state. Through a fully nonlinear bifur-cation and stability analyses, Lan et al. [15] constructed the bifurcation diagrams with detailed solution con-nectivity for a two-dimensional (2D) bounded case. Additional stable solution families and the transcritical bifurcation due to the interface concavity were found. For the cylinder con®guration with an interface, Chang and Brown [16], for the ®rst time, used a fully

non-linear ®nite element analysis to construct its bifurcation diagrams. However, their calculations were based on the axisymmetric assumption, and the 3D bifurcations were not considered. They also found that the system became less stable, i.e., an earlier onset of convection, as the interface was introduced. The subcritical bifur-cation due to the interface was not observed in their study.

From the previous 2D results [15], we expect that the coupling of the 3D bifurcations and solidi®cation or melting could lead to di€erent bifurcation behavior and thus di€erent solution connections. If the solutions are structurally stable, these solutions will be of interest in practice, and they may have an impact on the design and operation. Further ®nding a subcriticality will be par-ticularly useful because it causes a catastrophic change from a stable mode to the other. Therefore, in this paper we attempt to study, for the ®rst time, the 3D bifurca-tions of the two-phase RB problem in a cylinder using a numerical bifurcation analysis. Through the fully non-linear analysis, bifurcation diagrams can be constructed. From the diagrams, the onset of the ¯ow and the possible (preferred) modes, as well as their linear

Fig. 1. Schematic of Rayleigh±Benard problems in a cylinder: (a) one-phase problem; (b) two-phase problem at a conductive mode; (c) two-phase problem at a convective mode.

stability, can be determined. The roles of the interface and the solid amount are further investigated.

In this paper, we only consider the material being pure or almost pure so that the solutal e€ects and the mushy zone can be neglected. In most energy storage or materials processing applications, e.g., the vertical Bridgman crystal growth [3], the material is con®ned in a container, i.e., with solid walls. Therefore, the use of the no-slip wall as the boundaries does not lose its gener-ality. To generalize the system, the orientation (tilt angle c) of the cylinder is arbitrary, so that the e€ect of tilting (imperfect bifurcation) can also be studied. For com-parison, we will ®rst construct the bifurcation diagram for the ®xed-interface problem. As mentioned previ-ously, the detailed stability of the convective modes is still controversial, and most of the arguments are due to numerical resolution. The 3D calculations of the bifurcation and stability analyses require tremendous computation e€ort, and it is still dicult to get mesh-independent results. Therefore, in this paper, even though we also construct the bifurcation diagrams and propose some possible stable states for the ®xed-inter-face problem, we do not attempt to settle the contro-versy. Instead, we are more interested in the solutions that appear to be more stable even under an imperfect condition. Furthermore, from there, we can investigate the role of the deformable interface on the bifurcation structures. Again, the preferred modes are proposed.

In the next section, the model description is pre-sented, followed by numerical solutions in Section 3. Section 4 is devoted to results and discussion. Conclu-sions and the comments are drawn in Section 5. 2. Model description

The schematic of the two-phase RB problem in an inclined cavity is shown in Fig. 1(c), which is repre-sented in a Cartesian coordinate (x1_{; x}2_{; x}3_{). The height}

of the cylinder is H, which can be changed in this study, and the radius is R. The e€ects of melt aspect ratio A is de®ned by Lm=R, while the initial ratio of solid and melt

heights is Ls=Lm. The temperature at the top and bottom

ends is ®xed; the lower temperature TH is higher than

the melting point (Tm) of the material inside, while the

upper temperature TC is lower. The side walls are

as-sumed adiabatic; heat loss can be considered if neces-sary. The deformable melt/solid interface height h…x1_{; x}2_†

is unknown a priori and needs to be solved simulta-neously with ®eld variables. Physical properties of the material are assumed constant, while the melt is in-compressible. Again, these assumptions can be relaxed if necessary.

Dimensionless variables are de®ned by scaling the length with R, velocity with am=R, and time with R2=am,

where am is the thermal di€usivity. The dimensionless

temperature …h† is de®ned as h  …T ÿ Tm†=…THÿ Tm†,

while hsis used for the solid phase. For the convenience

of illustration, TCis adjusted so that hsis equal to ÿ1 for

Ls=Lmˆ 1 at the upper boundary. The coordinates

(x1_{; x}2_{; x}3_{), the interface height h, and velocity}

compo-nents (u; v; w) used afterwards are all dimensionless. The Boussinesq approximation [19] is also adopted in the ¯ow calculation. This assumption is reasonable if the temperature di€erence is not too large. Since the melt/ solid density ratio is usually about one for most ma-terials, the volume change due to solidi®cation or melting is neglected. The conservation equations in a dimensionless form for the 3D incompressible laminar ¯ow of a Newtonian ¯uid and heat transfer in the melt can be described as follows:

r
v ˆ 0; …1†

ov=os ‡ v
rv ˆ ÿrP ‡ Prr2_{v ÿ Pr Ra he}

g; …2†

oh=os ‡ v
rh ˆ r2_{h; …3†}

where s; v; P, and h are the dimensionless time, velocity, pressure, and temperature, respectively. Pr is the Prandtl number and eg the unit vector of gravity orientation.

The associated dimensionless number Ra in the source term of the momentum equation is the thermal Rayleigh number, de®ned as follows:

Ra gbT…THÿ Tm†R3

ammm ;

where g is the gravitational acceleration, bTthe thermal

expansion coecients, and mmis the kinematic viscosity.

In the solid …s†, only the conductive heat transfer needs to be considered

ohs=os ˆ jsr2hs; …4†

where js ks=km is the dimensionless thermal

conduc-tivity of the solid, and it is assumed to be equal to one in this study.

The boundary conditions for temperature are set to be:

h…x1_{; x}2_{; 0† ˆ 1; h}

s…x1; x2; H=R† ˆ ÿ1; …5†

for the bottom and top surfaces. The side wall is as-sumed adiabatic,

n
rh ˆ 0; …6†

where n is the normal vector of the solid wall pointing outward. For the velocity, the no-slip boundary condi-tion is adopted

v ˆ 0: …7†

More importantly, at the melt/solid interface, the interface energy balance is imposed as

St…n
e3†oh_osÿ jsoh_ons   s‡ oh on   mˆ 0; …8†

as well as the melting-point isotherm:

h ˆ hsˆ 0: …9†

In Eq. (8), St  DH=‰Cpm…THÿ Tm†Š is the Stefan

number, where DH is the heat of fusion and Cpm the

speci®c heat of the melt. In Eq. (9), the e€ect of inter-facial energy, i.e., the Gibbs±Thompson e€ect, is ne-glected because the system dimension is usually in the order of centimeters. This e€ect can be added easily if necessary. It should be pointed out that the numerical solution presented in this report can be used for a more general case. Therefore, the assumptions made so far are simply for the ease and clarity of illustration.

The measure of convective heat transfer is through the averaged Nusselt number Nu at the melt/solid in-terface de®ned as Nu ˆ ÿ1_S Z S oh on dS; …10†

where n is the unit normal vector pointing toward the solid, and S is the area of the interface. For the case without convection (pure heat conduction), Nu ˆ 1.

3. Numerical solutions

Due to the unknown and deformed interface shape h…x1_{; x}2_{†, a ®nite volume method (FVM) [17,20] based on}

the body-®tted coordinates is adopted. Fig. 2(a) shows a portion of a sample mesh for calculation; the mesh lines de®nes the boundaries of ®nite or control volumes (CVs). A typical CV is illustrated in Fig. 2(b). The idea of FVM is simply to integrate the governing equations over every CV. For example, for each CV, after the Gauss theorem is applied, a ¯ux balance equation for the variables can be obtained with the following form: Ieÿ Iw‡ Inÿ Is‡ Itÿ Ib‡ Z DVS/dV ÿ d ds Z DV/ dV ˆ 0; …11† where Ii, i ˆ …e; w; n; s; t; b†, represents the ¯uxes of

variables across the faces of CV and / ˆ …u; v; w; P or h†; DV is the volume of CV and S/ the source term of the

governing equations. Each of the ¯uxes Iiis made of two

distinct parts, namely a convective contribution and a di€usive contribution; the mesh movement also needs to be considered in the convective part. Both terms are approximated with the central di€erence scheme; up-winding is found unnecessary here. More importantly,

CVs with zero volume are also placed on boundaries. Accordingly, the boundary conditions can be imported into the ¯ux balance equations easily. Detailed im-plementation of the scheme can be found elsewhere [17,20]. For stationary solutions, the time-dependent term is removed.

In order to locate multiple steady states for station-ary solutions and to avoid the singular Jacobian matrix at bifurcation or limiting (turning) points, the pseudo arclength continuation [21] is further adopted. The idea of continuation is to trace the solution along a branch family through an arclength s by including an additional equation for the continuation parameter p. For example, if Ra is chosen as the continuation parameter, the ad-ditional equation can be written as the following:

oxT os   s0
x…s†‰ ÿ x…s0†Š ‡ oRa_os   s0 Ra…s† ‰ ÿ Ra…s0†Š ÿ …s ÿ s0† ˆ 0; …12†

where x ˆ …u; v; w; P; h; h†T _{and s ÿ s}

0 is the step size

along the branch, where the arclength s  k…x; Ra†T_k 2.

Appending the continuation equation to the steady-state ¯ux balance equations (Eq. (11) without the tran-sient term) for all CVs in the domain, as well as the isotherm condition (Eq. (9)), leads to a set of nonlinear algebraic equations:

f…y†  0; …13†

where y ˆ ‰x; RaŠT_{. Eq. (13) is solved by Newton's}

method for all variables simultaneously. It should be pointed out that the beauty of Newton's method is that it solves all the variables globally and does not require any special interface update schemes for the unknown interface shape. During Newton's iterations, the Jaco-bian matrix ~J, de®ned as ofi=oyj, is estimated by the

forward di€erence and the Newton's linear equations are solved by the GMRES method [22]. The preconditioner for GMRES used in this study is the incomplete LU decomposition without ®ll in, i.e., the so called ILU…0† preconditioner.

Although Newton's method provides a versatile so-lution to our problem, due to the limitation of computer memory (512 MB), the ®nest mesh we can use is quite limited. The total number of nonlinear algebraic equa-tions is limited to 2  105_{. However, we have examined}

most of the standard solutions by using an iterative decoupled approach [20], and the numerical error using the mesh here is believed to be acceptable. Furthermore, although the leading eigenvalue of the Jacobian matrix from the Newton's method can be used for stability analysis [15], it is, unfortunately, still too costly to cal-culate for our 3D problem. Therefore, instead of using the eigenvalues, we have adopted a fully time-dependent calculation for the stability analysis. By providing a 3D disturbance, we can examine the time evolution of the

responses. If the original state can be restored, the mode is stable. Otherwise, the solution is not stable. Further-more, by monitoring the singularity of the Jacobian matrix through the test function sij[15,23] during

solu-tion tracking, we are able to detect the bifurcasolu-tion points.

4. Results and discussion

In the general formulation of this problem, the so-lutions depend on several parameters including Ra, c, Pr, A, Ls=Lm, etc. However, in this study we focus

mainly on the Prandtl number of one with the initial melt aspect ratio being one (Pr ˆ 1; A ˆ 1). For crystal growth or solidi®cation applications, oxide materials have a comparable Prandtl number. The unity thermal conductivity ratio (i.e., jsˆ 1) is also reasonable. The

e€ect of the solid amount (Ls=Lm) will be discussed in the

end. Although the parameters we consider here are very limited, we believe that the role of the interface for other situations may be similar.

Before calculated results are presented, several nu-merical tests on the nunu-merical scheme have been per-formed to ensure the correctness of calculations. Fig. 3 shows the comparison of calculated thermal ®elds (an m1 mode) with those obtained by Fluent [24], a com-mercial CFD package, for Ls=Lmˆ 0 (with the interface

being ®xed) at Ra ˆ 4000; Dh ˆ 0:05. As shown, good agreement is obtained for both the top-view isotherms at x3_{ˆ 0:5 and the side-view ones at x}1_{ˆ 0. Other similar}

benchmark tests can be found elsewhere [25]. The e€ect of meshes on the m1 family for the case with 5° tilt is shown in Fig. 4; the interface is free and Ls=Lmˆ 1. As

shown, for Ra < 5000, the e€ect of meshes on the solu-tion, for the mesh ®ner than M1, seems to be insigni®-cant. For stronger convection, a ®ner mesh may be necessary. However, in this study, the convection is generally very week (Ra 6 8000). Therefore, we have

Fig. 3. Comparison of the calculated thermal ®elds of an m1 mode at Ra ˆ 4000 with those obtained by Fluent.

chosen the mesh M2 for key solution families (m0 and m1 modes), while the mesh M1 for the rest of solution branches. Further re®nement of the solution is possible by using the segregated approach [20]. However, such an approach is dicult to obtain unstable branches. Fur-thermore, as mentioned previously, the stability of a solution family can be sensitive to the mesh resolution, and that is the main reason for the di€erences among some of the previous reports. In reality, due to the limitation of computing resources, it may be still dicult to obtain the so-called mesh independent bifurcation points for 3D problems. On the other hand, bifurcation structures are usually less sensitive to the numerical ac-curacy. Therefore, in this report, we do not attempt to settle the controversy among the previous studies for the ®xed-interface problem. Instead, we will illustrate its key

bifurcations, especially, when the deformable interface is introduced.

4.1. Fixed interface (A ˆ 1 and Ls=Lmˆ 0)

As the interface is ®xed and there is no-slip on the side wall, the bifurcation diagram for Nu vs Ra is shown in Fig. 5(a); the unstable branches are indicated by dashed lines. Some typical solutions at Ra ˆ 8000 are shown in Fig. 5(b). The solution structure for 2000 < Ra < 4000 is redrawn in Fig. 6, while the cor-responding test function sij along the static and m1

branches is shown below. Some typical solutions are shown in Fig. 6 as well. In Figs. 5 and 6, the solution families with a clear symmetry at x3_{ˆ 0:5 are marked by}

m0, m1, m2, or m3, respectively, for illustration. One can examine the solutions easily to get a better picture of the ¯ow structures. As shown in Fig. 5(a) or Fig. 6, the ®rst primary bifurcation is at Racˆ 2250 leading to an stable

axisymmetric family (m0). This value is in good agree-ment with previous reported values (2250 by Hardin and Sani [12] and 2200 by Neumann [11]). The dramatic change of the test function along the static line is also in good agreement with the bifurcation point shown in Fig. 5(a) or Fig. 6. This bifurcation is also the onset of convection, which is a transition from the conductive to convective modes. The second and third primary bifur-cations lead to m2 (Racˆ 2480) and m1 (Racˆ 2700)

modes, respectively. Again, they are in reasonably good agreement with previous reported values as well (Racˆ 2470; 2700, respectively [12,14]).

Beside these primary bifurcations, however, the stability of the branching families is also a major interest in practice. However, it is still quite controversial.

Fig. 5. Calculated results for the single-phase problem (A ˆ 1 and Ls=Lmˆ 0): (a) bifurcation diagram; (b) some thermal and ¯ow

structures.

Fig. 4. The e€ects of meshes on Nu for a 5° tilt with a de-formable interface.

Wanschura et al. [14] tried to settle the arguments, but no conclusions have been made to date. They also ob-tained a stable m0 solution up to about 5±10Rac for

Pr ˆ 6:7. However, our m0 solution becomes unstable at about Ra  8000, which can be identi®ed by the change of the test function. We have also used the Fluent code to test the m0 mode, and the stability limit is at about Ra ˆ 6000±7000. Therefore, our stability limit for m0 seems to be reasonable. However, if we use the mesh M1 for calculation, the stability limit is much lower being at Ra  3600. As illustrated in Fig. 6, two stable m0 modes are in the same branch due to the same Nu; one is marked by m0U and one by m0D. The m0U family has an upward ¯ow at the centerline, while m0D has a downward ¯ow. These two solutions are stable for Ra < 8000 at least. Beyond this point, we are not able to get a stable solution. The di€erence between our results

and previous ones is believed to be the mesh e€ect. However, even the m0 mode can be stable up to higher Ra, it is still very close to the neutral stability. As a re-sult, it can become unstable easily as an imperfection is introduced. For lower Pr values, the m0 mode becomes unstable at a much lower Ra value, and this can also be taken into account to our lower stability limit.

In Fig. 5(a) or Fig. 6, interestingly, the second pri-mary bifurcation to the m2 family is independent of the m0 mode. The ¯ow patterns at the intersection of the two branches are di€erent; one is m0 and one is m2. The solution of the m2 mode at Ra ˆ 8000 is shown in the b1 of Fig. 5(b). From its isotherms at x3_{ˆ 0:5, there are}

two hot spots (upward ¯ows) and two cold spots (downward ¯ows). The other m2 solution at lower Ra is the solution a1 in Fig. 5(a). This solution branch is di€erent from that proposed by Hardin and Sani [12] for Pr ˆ 6:7. They showed that the intersection was a transcritical bifurcation for Pr ˆ 6:7, where the m0 and m2 modes exchanged stability; the m0 mode became unstable after this point, while the m2 became stable. The m2 mode branching from Racˆ 2470 seems to be

quite stable. Its dynamic response for Ra ˆ 4000 with respect to the tilt disturbance (c ˆ 0:5° sin…pt† for t < 1; c ˆ 0 for t > 1) is shown in Fig. 7(a). Interestingly, through mesh re®nement, Wagner et al. [13] commented that for Pr ˆ 6:7, the appearance of the m2 mode (at Ra ˆ 2800) reported previously [10,11] was due to the use of a too coarse grid, and the m2 mode was not stable. There are other types of m2 modes, as shown in the solutions a2 and a3 in Fig. 5(a), and they can also be m2U and m2D for the solution a2; similar ®nding for Pr ˆ 6:7 was reported by Muller et al. [10]. Also, the solution a3 is also a typical symmetric m2 mode, and its symmetry-breaking bifurcation leads to the solution like a5. The solution a6 also bifurcates from the static line near Ra ˆ 2700. However, these solutions are not stable and hence not so interesting. Furthermore, Wanschura et al. [14] also obtained the m2 family for Pr ˆ 6:7, but

Fig. 7. Dynamic evolution of the basic state to an angular disturbance: (a) an m2 mode at Ra ˆ 4000; (b) an m1 mode at Ra ˆ 3000; the initial and ®nal thermal ®elds are shown in the ®gures.

Fig. 6. A close-up view of Fig. 5 for 2000 < Ra < 4000 and the singularity test function sijalong the static and m1 branches.

they suspected that there was a small range that the m2 mode may be unstable. From the geometric point of view, it is believed that a multi-cell family is favored for smaller aspect ratios. However, for A ˆ 1, it is not easy to judge a stable solution from such a consideration. As will be shown shortly, the m2 modes are hard to survive after a 5° tilt is introduced. Therefore, in practice, they may not be important as well.

The m1 mode is obtained from the third primary bifurcation point in Fig. 5(a) or Fig. 6 is found unstable until Ra ˆ 3250. Wanschura et al. [14] obtained the neutral value at about Ra > 4224. We are able to get a stable m1 solution at Ra ˆ 4000 by both our and Fluent codes. The comparison of the m1 solutions has been shown in Fig. 3. The solution connectivity at the sec-ondary bifurcation is still unknown. The time-dependent evolution of the m1 mode at Ra ˆ 3000 after an angular disturbance is shown in Fig. 7(b). As shown, the solution runs away, and ®nally becomes an m0 mode. Clearly, as shown in Fig. 6, four stable solutions may exist after Ra ˆ 3250. At smaller Ra, the solutions m0 and m2 seem to be more stable.

Other modes, such as the m3 family like solutions b4 and b5 in Fig. 5(b), are found. Again, these modes are unstable and thus not much interesting in practice.

In addition to using the pulse disturbance, the stability of a solution family can be further investigated though symmetry breaking. In practice, the system may not be perfectly aligned with the gravity. Therefore, a little tilt against the gravity is a typical symmetry breaking that leads to imperfect bifurcations. In general, less stable modes can be removed or become unstable, and as a result, the bifurcation structures become

sim-pler and the multiplicity decreases. As shown in Fig. 8 (a), at 5° tilt, the m0 and m2 modes disappear from our bifurcation diagram. At least, we could not obtain their family anymore by using the solution at c ˆ 0 as an initial guess. Indeed, as discussed by Lan et al. [15], the tilt reduces multiplicity and results in a much simpler solution structure. As shown in Fig. 8 (a), the primary mode now becomes m1, which starts from Ra ˆ 0. Typical isotherms and ¯ow patterns are a1 and b1 in Fig. 8(b). In this branch, the ¯ow direction follows the tilt, as shown by the solution a1. On the other hand, the other ¯ow direction against the tilt, such as in the solution b2, is also stable. Without tilt, the solutions b1 and b2 are images to each other; physically, they are the same. Hence, the tilt breaks the symmetry. The reason for b2 being stable is quite simple. As a ¯ow like b2 is developed, the small tilt against the ¯ow is not able to reverse it due to the inertial e€ect. However, at weaker convection below the turning point, the inertial term is not large enough to against the tilt, and as a result, only the b1 family exists.

Tracking the solution family downward from b2, we obtain an unstable family (b3) through a turning point. Such an imperfect bifurcation is generic, which is in general obtained from a pitchfork bifurcation. Similar structures were found in the 2D problem [15]. The other imperfect m3 modes was also obtained, such as the families for b4 or b5, but they are still unstable.

By comparing Figs. 5 and 8, it is clear that the m1 family seems to be more stable for A ˆ 1, while the ax-isymmetric m0 modes and m2 modes are less stable with respect to the symmetry breaking. Nevertheless, multi-plicity still exists, and the ®nal stable mode is indeed

Fig. 8. Calculated results for the single-phase problem (Ls=Lmˆ 0) with a 5° tilt: (a) bifurcation diagram; (b) some thermal and ¯ow

path dependent. Furthermore, the branching families with even modes, such as m0 or m2, become less stable when an asymmetric condition is introduced unless there are other stabilizing mechanisms involved. As will be illustrated shortly, the interface deformation may form a barrier to inhibit the m0 mode from running away. 4.2. Deformable interface (A ˆ 1 and Ls=Lmˆ 1)

As the interface is allowed to deform, the bifurcation structure changes signi®cantly. We take a basic con-duction mode with the same amount of solid and melt, i.e., Ls=Lmˆ 1, as an example for illustration; the heat

¯ux from the bottom is kept the same as before, so that Nu ˆ 1. The bifurcation diagram is shown in Fig. 9(a) and some typical thermal and ¯ow ®elds are shown in Fig. 9(b). As shown, the onset mode is no longer m0, but m1 and the bifurcation becomes subcritical. The primary bifurcation point (Racˆ 2130) is also signi®cantly lower

than the ®xed one (2250 for m0 and 2700 for m1), even though the zone length and the aspect ratio at the onset point are exactly the same as before. The m1 family bi-furcating from the static mode is unstable ®rst (dashed-line), but becomes stable after the turning point (Ra ˆ 2000). Typical stable m1 solutions like a1, b1, and c2 are shown in the ®rst column of Fig. 9(b). As shown, as the ¯ow develops, the interface shape is shaped by the convection signi®cantly. The melt height increases

sig-ni®cantly as well. As a result the solutions move away from the ®xed-interface ones. However, as the melt as-pect ratio A increases, the m1 mode is preferred. The codimension-two point for the ®xed-interface problem is at about A ˆ 1:12 and Ra ˆ 2700 [14], where the m0 and m1 modes may appear at the same time. Therefore, away from the static line, due to the larger melt height, the m1 mode predominates.

At the static line, the subcritical nature is due to the coupling of the ¯ow and the interface. Davis et al. [18] ®rst reported such a phenomena for a thin layer melted from below, and the transition to a hexagon mode was subcritical, i.e., with a jump in the melt-layer thickness. Good agreement between experiments and their weakly nonlinear analysis was obtained. They also found that the subcritical nature predominated for a thick solid layer. We have also obtained a similar conclusion that the subcriticality becomes obvious as the solid amount increases, which will be illustrated shortly. Furthermore, once the m1 mode has been developed, it will need a much lower Ra (below the turning point) to get back to its conductive state. The decrease of the heat transfer rate at the turning point is also dramatic. For a thick solid layer, this quenching point can even be lower than 1708, i.e., the classic onset point.

Beside the ®rst appearing mode, the second primary bifurcation is to an axisymmetric m0 mode, and again, it is subcritical (not very clear). However, due to the

de-Fig. 9. Calculated results for the two-phase problem (A ˆ 1 and Ls=Lmˆ 1): (a) bifurcation diagram; (b) some thermal and ¯ow

formation of the interface caused by the convection, the m0D and m0U modes do not have the same Nu, i.e., having di€erent heat transfer rates across the interface. For the upper branch, it soon becomes stable after the turning point. As shown in the solutions a2, b3, or c3 of Fig. 9(b) (the third column), this mode is m0U; the melt ¯ows upward at the center. Due to the melt convection, the interface becomes concave. On the other hand, the lower branch has a di€erent ¯ow direction, i.e., m0D. As shown in the fourth column of Fig. 9(b), the typical solutions are b4 and c4, and their interfaces are convex. Indeed, right at the static line, there is a transition from the m0D to m0U modes, where the interface concavity also changes. To better illustrate the bifurcation, we represent the bifurcation diagram using a new Nu, de-noted by Nu, by subtracting the original Nu by one (the static mode) and adding a minus sign to it for the m0D mode. As shown in Fig. 10, near the static line, it is clearly a transcritical bifurcation. Beside the dash-line, which is an m0U mode, both m0D and m0U families are stable. Similar transcritical bifurcation was also ob-served in the 2D case (a rectangular system) [15], and the bifurcation point was also due to the change of interface concavity; the turning point also appear in the concave-interface branch. Furthermore, as compared with the m0 mode in Fig. 6, the m0 modes here is more stable, even using the mesh M1. The interface shape, on the other hand, may provides a stabilizing mechanism that con-®nes the m0 ¯ow and hinders it from running away. This will become obvious as the 5° tilt is introduced, which will be discussed shortly.

As illustrated in Fig. 9(a), an m2 mode appears lately behind m0D, but it is unstable; it is stable for the ®xed-interface case. This is also an indication that the

pres-ence of the deformable interface may destabilize the system, easing the onset of the convection. One can also imagine that the interface morphology of the m2 mode is more complicated, and this may discourage its forma-tion. Also, the increasing melt aspect ratio due to con-vection does not favor the m2 family. Typical solutions are shown in the second column of Fig. 9(b) like solu-tions b2 and c1. Other m2 or higher-order modes are also found but they are not stable and thus not inter-esting. Therefore, we do not attempt to go any further. Furthermore, even though the m2 mode is unstable for Ra > 2300, we still have three stable modes, i.e., m1, m0U, and m0D. In reality, which mode being selected is indeed path dependent. Interestingly, these solutions are, at least, stable up to Ra ˆ 8000. The multiplicity exists over such a wide range is also di€erent from that with a rigid interface. Again, it believed that the inter-face shape con®nes the ¯ows and retain them to a higher Ra value. As the Ra number increases, the interface deformation increases as well, and this further inhibits the ¯ow bifurcation. The stability of these modes can be further illustrated through imperfect bifurcations.

An imperfect bifurcation due to 5° tilting is shown in Fig. 11(a). This bifurcation diagram is similar to the ®xed-interface one (Fig. 8). The subcriticality also dis-appears; it exists within a very small tilt angle. Again, the most stable family is still m1 branching from Ra ˆ 0, and it has a ¯ow following the tilt direction. Typical solutions from this branch are illustrated by the solu-tions a1, b1, and c1 in Fig. 11(b). Other two m1 modes are related to the solutions b2 (or c2) and b3 (or c3), respectively. Again, similar to the ®xed-interface ones, one is stable and one is not.

Interestingly, an imperfect m0 mode is still retained, and its upper branch is stable, such as the solutions b4 and c4. As shown in the fourth column of Fig. 11(b), this stable branch is an imperfect m0U mode. The lower branch, which is believed to be related to the m0D mode, is unstable. Typical solutions of the unstable imperfect m0D mode are illustrated in the solutions b5 and c5 of Fig. 11(b). Again, one may compare Figs. 9(a) and 11(a) and ®nd that the m0 family is detached from the static line after the 5° tilting because its perfect mode does not exist. This solution structure here is also somewhat dif-ferent from that of the 2D case [15], where two pairs of convex and concave solutions appear after tilt. Fur-thermore, even with 5° tilt, we still have three stable modes after the turning point (Ra  2800) of the m1 mode.

As mentioned previously, the subcriticality of the primary bifurcation is a€ected by the solid amount. To further illustrate that, we take the thickness ratio Ls=Lm

as a parameter to illustrate the primary bifurcation. The heat ¯ux from the bottom is kept the same by changing the upper thermal boundary condition. As shown in Fig. 12(a), with the increasing solid layer thickness (or solid

Fig. 10. A redrawn transcritical bifurcation diagram from Fig. 9(a) for the m0 modes; Nuis de®ned by subtracting Nu by one and adding a minus for the m0U mode. Two solutions at Ra ˆ 3000 are shown.

amount), the bifurcation of the m1 mode changes indeed from supercritical to subcritical. Solutions at Ra ˆ 2500 are also shown in Fig. 12(b); the m0 mode is also in-cluded for comparison. As shown, when the solid layer thickness is zero, which can be done numerically, the solution is the same as the ®xed-interface one. However, there is a dramatic change of the solution when the solid layer becomes ®nite. The subcriticality occurs at about Ls=Lm 0:5. Further increasing the solid amount, both

the primary and turning points move forward (to the left). However, for Ls=Lm> 4, the change becomes much

smaller. This result is also consistent with the weakly nonlinear analysis for a thin-layer of material melted

from below in an in®nite domain by Davis et al. [18] that the subcriticality to hexagon cells predominates for a thick solid layer.

On the other hand, the e€ect of the solid layer has little e€ect on the bifurcation point of the m0 mode. Therefore, as the solid layer increases, the primary bi-furcation point of the m1 mode moves forward, and to a certain point, at about Ls=Lmˆ 0:3, both modes appear

at the same time, i.e., a co-dimension-two point. At this point, the onset of the convection can be m0 or m1 de-pending on the initial disturbance. Furthermore, the subcriticality nature is particularly interesting in prac-tice. The onset of convection occurs dramatically to a

Fig. 12. The e€ect of solid amount on the primary bifurcation: (a) bifurcation diagram; (b) thermal ®elds and interface shapes at Ra ˆ 2500. The heat ¯ux for all cases is kept the same.

®nite convective mode as the subcritical point is passed. On the other hand, a ®nite convective mode can be quenched to a conductive mode by reducing the driving force to be below the turning point. Such a hystersis may cause a process dicult to control or operate.

5. Conclusions

In this study, we have explored the 3D nonlinear phenomena of a generic two-phase Rayleigh±Benard problem in a cylinder. The fully nonlinear analysis is based upon a ®nite-volume/Newton's method with a solution tracking capability. A singularity test function and dynamic calculations have also been used to locate the bifurcation points and to examine the stability of solution families. The bifurcation diagrams for Pr ˆ 1 and A ˆ 1 have been constructed for both ®xed and free interfaces.

It is found that the solid layer plays a critical role in the two-phase problem, both on the ¯ow onset and the convection modes. As the solid layer diminishes to zero, the free-interface solution is the same as the ®xed-in-terface one. For the ®xed-in®xed-in-terface problem, the ®rst primary bifurcation is either m0U or m0D, the second appearing mode is an m2 mode, and then is an m1 mode. Other unstable or high-order bifurcations are also found, but they seem to be unimportant in practice and are vulnerable to imperfections. As the solid amount increases, the onset of the m1 mode moves forward (Rac

is reduced), while the onset points of the m0 and m2 modes are not changed much. As a result up to a certain solid amount, the ®rst appearing mode becomes m1. As the solid amount is further increased, the primary bi-furcation becomes subcritical. Due to the con®nement of the interface shape, the m0 modes becomes more stable, and its transcritical bifurcation is similar to the 2D one [15], where the transcritical point marks the change of interface concavity. Furthermore, with the interface deformation, as compared with a ®xed-interface one, the m0U mode is also more stable and can be retained even with a 5° tilt.

Although the solution structures obtained here are for Pr ˆ 1, according to the 2D results [15], we do not expect much change along the static line for di€erent Pr numbers. However, away from the static line, the e€ects are not trivial. Nevertheless, the bifurcation structures presented here are quite generic, which may provide useful information for further study.

We have also found multiple stable steady states for both ®xed and free-interface problems and the illus-trated solution structures also provide the basis for further analysis. The imperfect bifurcations due to 5° tilt also help further identify the stability of a solution family. As mentioned, those modes near the neutral stability, which may be quite controversial from di€erent

calculations, can be destroyed easily. Again, the modes that survive under the imperfection are more stable, and may be easier to be observed in experiments. Further-more, in practice, all the stable modes are path depen-dent, and are sensitive to the initial conditions and disturbances. Nevertheless, with the solution structure in mind, it may be easier to better control the system. In addition, other parameters, such as the thermal con-ductivity ratio of the melt and the solid phases and the e€ects of boundary conditions, etc., may be important as well, but they are beyond the scope of this study.

On the other hand for 3D problems, numerical ac-curacy is indeed a€ected by grid resolution. However, the computation e€ort increases signi®cantly as the mesh is re®ned. Therefore, even for the ®xed-interface prob-lem, there are still some disagreements on the solution stability and structures in the literature, particularly, for the bifurcation points with a higher null dimension. On the other hand, because these high-order bifurcations are usually structurely unstable, they are dicult to ®nd experimentally. In other words, these bifurcation may not be so interesting, even though they are helpful in constructing detailed diagrams. Furthermore, they are, unfortunately, usually very sensitive to the mesh and hard to calculate accurately. On the contrary, imperfect structures are more generic and structurely stable and can be easily obtained both experimentally and numer-ically with a much better con®dence. Therefore, the tilt introduced here also help map the key picture of the nonlinear bifurcations.

Finally, solutal e€ects [26,27] and time-dependent modes [14] are not considered in this work. It is believed that similar solution structures due to the solutal e€ects may be found. Furthermore, investigating the 3D time-dependent modes, which are believed to occur at higher Ra numbers, still remains a grand challenge for nu-merical simulation. Again, the grid resolution is a major concern. The control of the bifurcation through rotation or other external forces like magnetic ®elds may be in-teresting as well, and they will be discussed in the near future.

Acknowledgements

This work is supported by the National Science Council of the Republic of China under Grant No. NSC 85-2214-E-008-014.

References

[1] J. Pannu, G. Joglekar, P.A. Rice, Natural convection to cylinders of phase change material used for thermal storage, AIChE Symposium Series (1980) 47±55.

[2] H.W. Hale, R. Viskanta, Solid±liquid phase change heat transfer and interface motion in materials cooled from above and below, Int. J. Heat Mass Transfer 23 (1980) 283±292.

[3] G. Muller, A. Ostrogorsky, Convection in melt growth, in: D.T.J. Hurle (Ed.), Handbook of Crystal Growth, North-Holland, Amsterdam, 1993, pp. 709±819.

[4] D.T.J. Hurle, E. Jakeman, Introduction to the techniques of crystal growth, PhysicoChem. Hydrodyn. 2 (1981) 237± 244.

[5] D.T.J. Hurle, E. Jakeman, A.A. Wheeler, Hydrodynamic stability of the melt during solidi®cation of a binary alloy, Phys. Fluids 23 (1983) 624±627.

[6] S.F. Liang, A. Vidal, A. Acrivos, Buoyancy-driven con-vection in cylinder geometries, J. Fluid Mech. 86 (1969) 239±256.

[7] G.S. Charlson, R.L. Sani, Thermoconvective instability in a bounded cylindrical ¯uid layer, Int. J. Heat Mass Transfer 13 (1970) 1479±1496.

[8] G.S. Charlson, R.L. Sani, Finite amplitude axisymmetric thermoconvective ¯ow in a bounded cylindrical layer of ¯uid, Int. J. Heat Mass Transfer 14 (1971) 2157±2160. [9] Y. Yamaguchi, C.J. Chang, R.A. Rrown, Multiple

buoy-ancy-driven ¯ows in vertical cylinder heated from below, Philos. Trans. R. Soc. Lond. A 312 (1984) 519±552. [10] G. Muller, G. Neumann, W. Weber, Natural convection in

vertical Bridgman con®gurations, J. Crystal Growth 70 (1984) 78±93.

[11] G. Neumann, Three-dimensional numerical simulation of buoyancy-driven convection in vertical cylinders heated from below, J. Fluid Mech. 214 (1990) 559±678.

[12] G.R. Hardin, R.L. Sani, Buoyancy-driven instability in a vertical cylinder: binary ¯uids with Soret e€ect. Part II: weakly non-linear solutions, Int. J. Numer. Methods Fluids 17 (1993) 755±786.

[13] C. Wagner, R.C. Friedrich, R. Narayanan, Comments on the numerical investigation of Rayleigh and Marangoni convection in a cylinder, Phys. Fluids 6 (1994) 1425±1433. [14] M. Wanschura, H.C. Kuhlmann, H.J. Rath, Three-dimen-sional instability of axisymmetric buoyant convection in

cylinders heated from below, J. Fluid Mech. 326 (1996) 399±415.

[15] C.W. Lan, M.C. Liang, M.K. Chen, Stability and bifur-cation analyses of a two-phase Rayleigh Benard problem in a cavity, Phys. Fluids 10 (1998) 1329±1343.

[16] C.J. Chang, R.A. Brown, Natural convection in steady solidi®cation: ®nite element analysis of a two-phase Ray-leigh±Benard problem, J. Comput. Phys. 53 (1984) 1±26. [17] M.C. Liang, C.W. Lan, A ®nite-volume/Newton method

for a two-phase ¯ow problem using primitive variables and collocated grids, J. Comput. Phys. 127 (1996) 330±345. [18] S.H. Davis, U. Muller, C. Dietsche, Pattern selection in

single-component systems coupling Benard convection and solidi®cation, J. Fluid Mech. 144 (1984) 133±151. [19] L.D. Landau, E.M. Lifshitz, Fluid Mechanics, second ed.,

vol. 6, Pergamon Press, Elmsford, 1987, p. 217.

[20] C.W. Lan, M.C. Liang, Multigrid methods for incom-pressible heat ¯ow problems with an unknown interface, J. Comput. Phys. 152 (1999) 55±77.

[21] H.B. Keller, Numerical solution of bifurcation and non-linear eigenvalue problems, in: P.H. Rabinowi (Ed.), Applications of Bifurcation Theory, Academic Press, New York, 1977, pp. 359±384.

[22] Y. Saad, M.H. Schultz, GMRES: a generalized minimal residual algorithm for solving nonsymmetric linear sys-tems, SIAM J. Sci. Stat. Comput. 7 (1986) 856±869. [23] R. Seydel, Numerical computation of branch points in

nonlinear equations, Numer. Math. 33 (1979) 339±352. [24] Fluent UNS 5.0 User Manual, Fluent Inc., 1999. [25] C.W. Lan, M.C. Liang, Three-dimensional simulation of

vertical zone-melting crystal growth: from symmetry breaking to multiple states, J. Crystal Growth 208 (2000) 327±340.

[26] G.B. McFadden, S.R. Coriell, Solutal convection during directional solidi®cation, in: Proceedings of the First National Fluid Dynamics Congress, Cincinnati, OH, 25± 28 July 1988.

[27] M.D. Impey, D.S. Riley, A.A. Wheeler, K.H. Winters, Bifurcation analysis of solutal convection during direc-tional solidi®cation, Phys. Fluids A 3 (1991) 535±550.