On oscillatory modes of nonlinear compositional convection in mushy layers

Nonlinear Analysis: Real World Applications 10 (2009) 209–226 www.elsevier.com/locate/nonrwa

On oscillatory modes of nonlinear compositional convection in mushy layers D.N. Riahi ∗ Department of Mathematics, 1201 West University Drive, University of Texas-Pan American, Edinburg, TX 78541-2999, United States Received 8 March 2007; accepted 27 August 2007

Abstract We consider the problem of nonlinear compositional convection in horizontal mushy layers during the solidification of binary alloys. Under a near-eutectic approximation and the limit of large far-field temperature, we determine a number of weakly nonlinear oscillatory solutions, and the stability of these solutions with respect to arbitrary disturbances is then investigated. The present investigation is an extension of the problem of the oscillatory modes of convection, which was investigated by [D.N. Riahi, On nonlinear convection in mushy layers. Part 1. Oscillatory modes of convection, J. Fluid Mech. 467 (2002) 331–359] in the absence of the main permeability parameter K 1 and very recently by [P. Guba, M.G. Worster, Nonlinear oscillatory convection in mushy layers, J. Fluid Mech. 553 (2006) 419–443] for two-dimensional cases in the presence of K 1 , to include the effects of K 1 , over a range of values of the other parameters, and for both two- and three-dimensional motion. It was found, in particular, that the results reported in [D.N. Riahi, On nonlinear convection in mushy layers. Part 1. Oscillatory modes of convection, J. Fluid Mech. 467 (2002) 331–359; P. Guba, M.G. Worster, Nonlinear oscillatory convection in mushy layers, J. Fluid Mech. 553 (2006) 419–443] are recovered if K 1 is zero or sufficiently small. In such cases two-dimensional solutions in the form of simple-travelling rolls are mostly the only stable and preferred solutions. However, as in the more realistic cases, if K 1 is not sufficiently small, then such solutions are replaced by preferred and stable three-dimensional solutions, which are mostly simple-travelling waves in the form of rectangles, squares or hexagons. c 2007 Elsevier Ltd. All rights reserved.

Keywords: Compositional convection; Mushy layer; Oscillatory convection; Buoyant flow; Nonlinear convection

1. Introduction Compositional convection in horizontal mushy layers during alloy solidification has been the subject of a number of studies [8,15,1–3,13,7,9,12,10]. The present study considers the problem of oscillatory modes of nonlinear compositional convection in mushy layers. A main motivation of the present and most of the above studies has been to understand the qualitative features of the modes of convection in the mushy layer, the onset of plume convection and the formation of chimneys and chimney convection, which originates from plume convection and can lead to highly undesirable imperfections called freckles in the alloy crystal. A subsequent goal will then be to develop ways to effectively control the undesirable plume and chimney convection. ∗ Tel.: +1 217 333 0679; fax: +1 217 244 5707.

E-mail address: [email protected]. c 2007 Elsevier Ltd. All rights reserved. 1468-1218/$ - see front matter doi:10.1016/j.nonrwa.2007.08.028

210

D.N. Riahi / Nonlinear Analysis: Real World Applications 10 (2009) 209–226

The theoretical studies of the compositional convection in mushy layers [1–3,7,9,12] were mainly based on an analytical model that was developed first by Amberg and Homsy [1]. Amberg and Homsy [1] made the simplifying assumptions that both the thickness of the layer and the departure from the eutectic point being small and the mushy layer being isolated from the overlying liquid layer. However, the qualitative features of the results based on this model did not appear to be changed when the impermeable upper boundary of the mushy layer is replaced by a permeable one [7]. Later [12], hereafter referred to as R02, extended the linear work of Anderson and Worster [3] to some extent by studying the problem of nonlinear two- and three-dimensional oscillatory compositional convection in horizontal mushy layers during the solidification of binary alloys. He analysed the oscillatory modes of convection in particular range of the parameter values where the critical value Rc of the scaled Rayleigh number R for the onset of oscillatory convection is distinctly lower than the critical value of the Rayleigh number for the onset of steady convection. His results indicated preference of supercritical simple-travelling rolls over most of the studied range of the parameter values, while supercritical standing rolls could be preferred only over a rather small range of the parameter values. His detailed nonlinear study of the oscillatory modes of convection in mushy layers complemented to some extent the previous nonlinear studies of the stationary convection in mushy layers [1,2]. However, due to the complexity of the three-dimensional oscillatory problem, R02 considered a simplifying assumption by assuming that the main permeability problem K 1 is of order the perturbation amplitude ε (ε  1). As a result of this restriction, the effects of K 1 did not enter the weakly nonlinear results reported in R02. Very recently, using the [1]’s model, Guba and Worster [10] studied analytically the weakly nonlinear two-dimensional oscillatory convection in mushy layers for order one values of the parameter K 1 . They found, in particular, that depending on the parameter values, supercritical travelling or standing rolls can be stable, but travelling rolls were likely to be selected for weak sensitivity of the permeability on the solid fraction. These results were in agreement with those in R02 in the limit of sufficiently small values of K 1 . In the present paper we extend both of the work in R02 and [10] by including, respectively, the parameter K 1 and three-dimensional effects in the analysis and found interesting results. In particular, we found that in contrast to the results reported in R02, which were effectively valid for zero value of K 1 , and depending on the parameter values three-dimensional supercritical solutions in the form of oscillatory rectangles, squares or hexagons can become the preferred and stable form of the compositional convection in a mushy layer, while two-dimensional solutions reported in [10] can become unstable. The preferred solutions that we referred to earlier in the previous paragraphs were based on the assumption of the type adopted before by Busse [5] that the solution, which exists at the lowest value of the control parameter R, is the physically preferred solution. This assumption can follow from the stability results [4,11]. In the present paper we follow [5] and, in addition, carry out stability analysis of the finite-amplitude solutions to determine the preferred and stable oscillatory solutions. The following two Sections 2 and 3 deal with the governing system and the finite amplitude and stability analyses. The results and discussion for the solutions and stability of the oscillatory modes are presented in Section 4, which is followed by conclusion in Section 5. 2. Governing system We consider a binary alloy melt that is cooled from below and is solidified at a constant speed V0 . Following [1] and [2], we consider a horizontal mushy layer of thickness d adjacent and above the solidification front to be physically isolated from the overlying liquid and underlying solid zones. The overlying liquid is assumed to have a composition ˜ C0 > Ce and temperature T∞ > TL (C0 ) far above the mushy layer, where Ce is the eutectic composition, TL (C) is the liquidus temperature of the alloy and C˜ is the composition. It is then assumed that the horizontal mushy layer is bounded from above and below by rigid and isothermal boundaries. We consider the solidification system in a moving frame of reference o x˜ y˜ z˜ , whose origin lies on the solidification front and translating at the speed V0 with the solidification front in the positive z˜ -direction. The reader is referred to R02 for the motivation and justification in using the Amberg and Homsy [1] type of model for the present study. The equations for Darcy-momentum, continuity, heat and solute for the flow in the mushy layer in the already described moving frame are non-dimensionalized by using V0 , k/V0 , k/V02 , β1Cρgk/V0 , 1C and 1T as scales for velocity, length, time, pressure, solute and temperature, respectively. Here k is the thermal diffusivity, ρ is a reference

D.N. Riahi / Nonlinear Analysis: Real World Applications 10 (2009) 209–226

211

(constant) density, β = β ∗ −Γ α ∗ , where α ∗ and β ∗ are the expansion coefficients for the heat and solute, respectively, and Γ is the slope of the liquidus, which is assumed to be a constant, 1C = C0 –Ce , 1T = TL (C0 ) − Te and Te is the eutectic temperature. The non-dimensional form of the equations for Darcy-momentum, continuity, temperature and solute concentration in the mushy layer are then ˜ u˜ = −∇ P˜ − R˜ θ˜z , K (φ)

(1a)

∇.u˜ = 0,

(1b)

˜ + u.∇ (∂/∂ t˜ − ∂/∂ z˜ )(θ˜ − St φ) ˜ θ˜ = ∇ 2 θ˜ ,

(1c)

˜ θ˜ + Cr φ] ˜ + u.∇ ˜ θ˜ = 0, (∂/∂ t˜ − ∂/∂ z˜ )[(1 − φ)

(1d)

where u˜ = ux ˜ + vy ˜ + wz ˜ is the volume flux per unit area, u˜ and v˜ are the horizontal components of u˜ in the x˜ and y˜ -directions, respectively, x and y are unit vectors in the positive x˜ and y˜ -directions, w˜ is the vertical component of u˜ in the z˜ -direction, z is a unit vector in the positive z˜ -direction, P˜ is the modified pressure, θ˜ is the non-dimensional form of either composition or temperature as shown in R02, t˜ is the time variable, φ˜ is the local solid fraction, ˜ of the R˜ = β1CgΠ (0)/(V0 ν) is the Rayleigh number, Π (0) is reference value at φ˜ = 0 of the permeability Π (φ) ˜ ˜ porous medium, ν is the kinematic viscosity, g is acceleration due to gravity, K (φ) ≡ Π (0)/Π (φ), St = L/(C L 1T ) is the Stefan number, C L is the specific heat per unit volume, L is the latent heat of solidification per unit volume, Cr = (Cs − C0 )/1C is a concentration ratio, and Cs is the composition of the solid phase forming the dendrites. Eq. (1d) is based on the limit of sufficiently large value of the Lewis number k/ks , where ks is the solute diffusivity. The above equations are subject to the following boundary conditions: θ˜ + 1 = w˜ = 0

at z˜ = 0,

(1e)

θ˜ = w˜ = φ˜ = 0

at z˜ = δ,

(1f)

where δ = dV0 /k is a growth Peclet number representing the dimensionless depth of the layer. Next, we consider the following rescaling in the limit of sufficiently small δ: Cr = C/δ,

St = S/δ,

(x, ˜ y˜ , z˜ ) = δ(x, y, z),

ε  δ  1, t˜ = δ t,

˜ R = δ R,

2

2

˜ = (θ B , φ B , 0, PB ) + ε[θ (x, y, z, t), φ(x, y, z, t), (R/δ)u(x, y, z, t), RP(x, y, z, t)], ˜ u, (θ˜ , φ, ˜ P)

(2a) (2b) (2c)

where C and S are the corresponding rescaled parameters, ε is small perturbation amplitude, and the quantities with subscript ‘B’ are the basic flow variables for the motionless state, which are a function of z only and are given below in terms of asymptotic expansions for small δ θ B = (z − 1) + δ(z − z 2 )G/2 + · · · ,

(3a)

φ B = δ(1 − z)/C + δ [−(1 − z) /C + (z − z)G/(2C)] + · · · ,

(3b)

PB = P0 + R[(z − z 2 /2) + δ(z 2 /2 − z 3 /3)G/2 + · · ·],

(3c)

2

2

2

2

˜ is considered where G ≡ 1 + S/C and P0 is a constant. Since φ is small, the following expansion for K (φ) ˜ = 1 + K 1 φ˜ + K 2 φ˜ 2 + · · · , K (φ)

(4)

where K 1 and K 2 are constants. For the analysis presented in the next section, it was found to be convenient to use the general representation u = Ω V + EΨ , Ω ≡ ∇ × ∇ × z, E ≡ ∇ × z,

(5)

for the divergent-free vector field u [6], where V and Ψ are the poloidal and toroidal functions for u, respectively. By taking the vertical component of the curl of (1a), it can be shown that the toroidal part EΨ of u must vanish. Taking the vertical component of the double curl of (1a) and using (5) and (1b) in (1), we find the final version of the governing system ∇ 2 [K (φ B + εφ)∆2 V ] + (∂/∂z)[Ω V.∇ K (φ B + εφ)] − R∆2 θ = 0,

(6a)

212

D.N. Riahi / Nonlinear Analysis: Real World Applications 10 (2009) 209–226

(∂/∂t − δ∂/∂z)(−θ + Sφ/δ) + R(dθ B /dz)∆2 V + ∇ 2 θ = ε RΩ V.∇θ, (∂/∂t − δ∂/∂z)[(−1 + φ B )θ + θ B φ + εφθ − Cφ/δ] + R(dθ B /dz)∆2 V = ε RΩ V.∇θ,

(6b) (6c)

θ = V = 0 at z = 0, θ = V = φ = 0 at z = 1,

(6d) (6e)

where ∆2 ≡ ∂ 2 /∂ x 2 + ∂ 2 /∂ y 2 . 3. Analysis Similar to the analyses presented in R02, we first carry out a weakly nonlinear analysis, based on a double-series expansions in powers of δ and ε, to determine the finite-amplitude oscillatory solutions and then investigate the stability of such solutions. Because of the similarity of the present analysis with the corresponding ones in R02, we make the description of the analysis here as short as possible and at the same time make the presentation somewhat self-explanatory to the reader. As in the finite-amplitude analysis carried out in R02, we first make a formal asymptotic expansion in ε and then at each order in ε make a formal asymptotic expansion in δ. Since we investigate oscillatory modes of convection, the following expansions are for the dependent variables of the perturbation system (1), R and the frequency ω for the oscillatory modes of convection: (V, θ, φ, R, ω) = [(V00 + δV01 + · · ·), (θ00 + δθ01 + · · ·), (φ00 + δφ01 + · · ·), (R00 + δ R01 + · · ·), (ω00 + δω01 + · · ·)] + ε[(V10 + δV11 + · · ·), (θ10 + δθ11 + · · ·), (φ10 + δφ11 + · · ·), (R10 + δ R11 + · · ·), (ω10 + δω11 + · · ·)] + ε 2 [(V20 + δV21 + · · ·), (θ20 + δθ21 + · · ·), (φ2(−1) /δ + φ20 + δφ21 + · · ·), (R20 + δ R21 + · · ·), (ω20 + δω21 + · · ·)] + · · · .

(7)

The analyses for the oscillatory modes were already done in R02. Hence, no details will be provided here and, instead, the main procedure and results of the analysis are given briefly for the oscillatory modes. We consider first the linear problem. At order 1/δ we find ω00 = 0. At order δ 0 we find N X

V00 = {(π 2 + a 2 )/[R00 (a)2 G]} sin(π z)

+ − − (A+ m ηm + Am ηm ),

(8a)

m=−N

θ00 = − sin(π z)

N X

+ − − (A+ m ηm + Am ηm ),

(8b)

m=−N N X

φ00 =

+ ∗ − − [ f m (z)A+ m ηm + f m (z)Am ηm ],

(8c)

m=−N

where ± ηm ≡ exp[i(am .r ± Sm ω01 t)],

R00 = [(π + a ) /(a G)]] 2

2 2

2

f m (z) = {−2π /[CG(π

(8d)

0.5

,

2 − ω01 )]}{(iω01 Sm /π ) sin(π z) + cos(π z) + exp[iω01 Sm (z

3

2

for m > 0

and

(8e) − 1)]}

(8f)

and −1 for m < 0. (8g) √ Here i is the pure imaginary number (i ≡ − 1), subscript ‘m’ takes only non-zero integer values from −N to N , N is a positive integer, r is the position vector, and the horizontal wave number vectors am satisfy the properties Sm = 1

am .z = 0,

|am | = a,

a−m = −am .

(9)

D.N. Riahi / Nonlinear Analysis: Real World Applications 10 (2009) 209–226

213

− The coefficients A+ m and Am satisfy the conditions N X

2 − 2 (|A+ m | + |Am | ) = 2,

± A±∗ m = A−m ,

(10)

m=−N

where the asterisk indicates the complex conjugate. Minimizing the expression for R00 , given in (8e), with respect to the wave number a, we find √ R00c = 2π/ G, ac = π. (11) It turns out that for all the values of the parameters where the oscillatory modes are calculated in the present study, the critical value Rc of R at the onset of motion for the oscillatory modes (R02) is smaller than the corresponding one for the steady modes. Next, the nonlinear problem for the governing system (6) at order ε is analysed. At order ε/δ, we find ω10 = 0. At order ε the system (6a)–(6e) can be reduced to the form given by (A.1) in the Appendix. Following the same reasoning given in R02, the solvability conditions for the nonlinear system (A.1) yield R10 = 0, while those for the system at order εδ yield R11 = 0. The types of solutions that can be realized or be preferred by the nonlinear system are described here briefly as follows. The simple-travelling form of oscillatory solutions can be in the form of either right-travelling mode, where the phase velocity of the mode is in the direction of the component of the mode’s wave vector along r, or in the form of left-travelling mode where the phase velocity of the mode is in the direction opposite to that of the component of the mode’s wave vector along r. Since identical nonlinear results are obtained in the present study for either leftor right-travelling modes, the analysis in this paper for the case of simple-travelling modes is presented only in the form of right-travelling modes. For the right-travelling solutions in the so-called ‘semi-regular’ case, in which scalar product between any one of the a-vectors and its two neighboring a-vectors assumes the constant values α1 and α2 ([4]; R02), we have + − 2 − 2 |A+ 1 | = · · · = |A N | = 0, |A1 | = · · · = |A N | = 1/N .

(12a)

For the case of regular solutions, we have α1 = α2 . For the standing wave solutions, we have ± 2 2 |A± 1 | = · · · = |A N | = 1/(2N ).

(12b)

For the oscillatory solutions in the form of general-travelling waves of the types introduced in R02 and in the regular or semi-regular case, we have + 0.5 A+ 1 = · · · = A N = [(0.5 − b)/N ] ,

− 0.5 A− 1 = · · · = A N = [(0.5 + b)/N ] .

(12c)

The constant b in (12c), which is restricted in the range |b| < 0.5

(12d)

(R02), is the parameter for the general-travelling wave, whose specific value in the range (12d) provides particular general-travelling wave solutions. The system of equations and boundary conditions at order ε are already given by (A.1) in the Appendix. The expressions for the solutions V10 , θ10 and φ10 of this system, which will be needed for the analysis of the governing system at O(ε2 ), are given below (V10 , θ10 ) =

N X

(o)

(o)

+ + − + − + [(Blp , Elp )(Al+ A+ p ηl η p + Al A p ηl η p )

l, p=−N (o)∗ + (Blp ,

(o)∗

+ − − − − − Elp )(Al+ A− p ηl η p + Al A p ηl η p )],

(13a)

214

D.N. Riahi / Nonlinear Analysis: Real World Applications 10 (2009) 209–226 N X

φ10 =

(11)

(12)

(21)

(22)

+ + + − + − − + − + − − − − [ flp Al+ A+ p ηl η p + f lp Al A p ηl η p + f lp Al A p ηl η p + f lp Al A p ηl η p ]

l, p=−N N X

+ − − − + + + − − − ( f m+ A+ m ηm + f m Am ηm + gm Am ηm + gm Am ηm ),

+

(13b)

m=−N (o)

(o)

(i j)

+ and g − , which are generally where the expressions for the coefficients Blp , Elp , flp (i, j = 0, 1, 2), f m+ , f m− , gm m functions of z, are given by (A.3) and (A.4) in the Appendix. The solvability conditions for the system at order ε 2 /δ yield ω20 = 0 and trivial (zero) solutions follow for the dependent variables. At order ε 2 the system (6a)–(6e) are reduced to the form given by (A.4) in the Appendix. The solvability conditions for this system yield a system of equations, which are given by (A.5) in the Appendix. The ± η± η± i, which differ from zero only if system (A.5) contain integral expressions of the form hηn± ηm l p

an + am + al + a p = 0

(14a)

±Sn ± Sm ± Sl ± S p = 0

(14b)

and

hold. Here an angular bracket indicates a total average over the layer and in time. Using the conditions (14a) and (14b) in (A.5), we find the following simplified equations: N X √ 2 − 2 (o1) + 2 + 2 (o2) − − + R20 Gπ(|A+ | + |A | ) = {[Tnm |An | |Am | + Tnm An Am A−n A+ −m δ(Sn + Sm ) n n m=−N (o3) − 2 + 2 (o4) + − − + Tnm |An | |Am | + Tnm An Am A−n A+ −m δ(Sn − Sm ) (o5) − + + (o6) + 2 − 2 + Tnm An Am A−n A− −m δ(Sm − Sn ) + Tnm |An | |Am | (o7) + + − (o8) − 2 − 2 An Am A−n A− + Tnm −m δ(Sn + Sm ) + Tnm |An | |Am | ]},

(n = −N , . . . , −1, 1, . . . , N ),

(15a)

where δ(a) = 1

for a = 0

and

0 for a 6= 0,

(15b)

(oi)

and the expressions for Tnm (i = 1, . . . , 8) are given, respectively, by (A.6a)–(A.6h) in the Appendix. The (oi) expressions for Tnm turn out to satisfy the symmetry conditions of the form (oi) (oi) Tnm = Tmn

(i = 1, . . . , 8).

(15c)

Since the simplest types of solutions are often observed in the applications, we restrict our attention to the simplest types of oscillatory solutions, which are either regular or semi-regular. Simple form of oscillatory solutions correspond to the cases of two-dimensional rolls (N = 1), squares (N = 2) and hexagons (N = 3), while those for a semiregular solution correspond to the cases with different values of angle γ of oscillatory rectangles (N = 2). Here, γ (0 < γ < 90◦ ) or 180◦ − γ is the angle between any two adjacent wave number vectors of a rectangular cell for a particular solution in the form of rectangles. To distinguish the physically realizable finite-amplitude solutions among all the possible oscillatory solutions, the stability of V , θ, φ with respect to arbitrary three-dimensional disturbances Vd , θd , φd is investigated. The timedependent disturbances can be assumed in the form (Vd , θd , φd ) = [V 0 (x, y, z, t), θ 0 (x, y, z, t), φ 0 (x, y, z, t)] exp(σ t),

(16)

where σ is the growth rate of the disturbances. When the governing equations and the boundary conditions of the form (6a)–(6e) for the finite-amplitude flow are subtracted from the corresponding equations and boundary conditions for the total dependent variables for the finite-amplitude flow and the disturbance quantities, and the resulting system is

D.N. Riahi / Nonlinear Analysis: Real World Applications 10 (2009) 209–226

215

linearized with respect to the disturbance quantities, we obtain the stability system, which is given by (A.7a)–(A.7e) in the Appendix. When the expansion (7) is used in (A.7a)–(A.7e), it becomes evident that the stability system can be solved by a similar expansion 0 0 0 0 0 0 (V 0 , θ 0 , φ 0 , ω0 , σ ) = [(V00 + δV01 + · · ·), (θ00 + δθ01 + · · ·), (φ00 + δφ01 + · · ·), 0 0 0 0 (ω00 + δω01 + · · ·), (σ00 + δσ01 + · · ·)] + ε[(V10 + δV11 + · · ·), 0 0 0 0 0 0 (θ10 + δθ11 + · · ·), (φ10 + δφ11 + · · ·), (ω10 + δω11 + · · ·), (σ10 + δσ11 + · · ·)] 0 0 0 0 0 + ε 2 [(V20 + δV21 + · · ·), (θ20 + δθ21 + · · ·), (φ02(−1) /δ + φ20 + δ φ021 + · · ·), 0 0 (ω20 + δω21 + · · ·), (σ20 + δσ21 + · · ·)] + · · · ,

(17)

where the expansion for φ 0 is singular at order ε 2 as δ → 0, but again as in the case R02, the O(1/δ) term is needed in the stability analysis since the O(ε2 ) stability problem is found to be forced by a term of order 1/δ in the solute equation for the disturbances. Following the assumptions and procedures described in R02, we find the most critical disturbances, which have the maximum growth rate, are characterized by σ0 = 0, where σ0 = σ00 + δσ01 + · · · . Then the linear and nonlinear solutions for the dependent variables of the disturbances at order δ 0 and ε are found, where the lengthy solvability conditions for the disturbance system at order ε yield σ10 = 0. Next, we applied the very lengthy solvability conditions for (A.7a)–(A.7e) at order ε 2 , which will not be given here, to determine σ20 . The leading order growth rate σ ∗ = εσ10 + ε 2 σ20 of the disturbances acting on the finite-amplitude oscillatory motion can then be determined from these systems following the method of approach due to [4], which is now a standard stability procedure. 4. Results and discussion 4.1. Linear problem The linear system with its eigenvalue problem for oscillatory modes, which led to the results (8)–(11), are in general, functions of the physical parameters C,S and K 1 . Although the scaled Stefan number S and the scaled compositional ratio C were represented earlier partly in terms of the composite parameters G and G t for the simplicity in the formulation, we are interested here to present and discuss the results for the oscillatory case in terms of the physical parameters. It should also be noted that in typical experiments with ammonium chloride–water solutions, like those of [14], the values of the scaled concentration ratio and the scaled Stefan number are about C ≈ 20δ and S ≈ 5δ, so that for a sufficiently thin mushy layer, the results presented in this paper, which are for the range of values 0.0355 ≤ C < 0.1 and S = 0.25C, could at least be qualitatively relevant. The linear results presented in [3] and in R02 (see also online supplementary materials to R02) were mainly in terms of the composite parameters G and G t . Here we briefly present the linear results in terms of the physical parameters instead. The period of oscillation of the solution does not depend on K 1 . For S = 0.25C, the period of oscillation decreases with C. It has a higher rate of increase with C at lower values of C. The dependence of Rc on the parameters implies that the linear system is stabilizing as K 1 increases or as S or C decreases for fixed relation between S and C. Since S represents a measure of the latent heat relative to the heat content and C represents a measure of the difference between the characteristic compositions of the solid and liquid phases and the compositional variation of the liquid, then the linear system is stabilized as C increases for a given S, or as S decreases, for a given C. Due to (11) and since S = 0.25C in the present calculation, the effect of the Stefan number dominates over that of the composition ratio and the flow is destabilizing as C increases. The stabilizing effect on the linear system when K 1 increases, is consistent with the physical role played by K 1 since the permeability of the mushy layer decreases with increasing K 1 . 4.2. Nonlinear problem Important quantity due to the nonlinear oscillatory effects is the coefficients R20 , which is calculated in the present study. Since R10 = 0, then it can be seen from the expansions (7) that R20 represents leading contributions to the

216

D.N. Riahi / Nonlinear Analysis: Real World Applications 10 (2009) 209–226

change in R required to obtain finite amplitude ε for a nonlinear solution. In terms of these coefficients the amplitude of convection is of order |ε| = [(R − Rc )/R20 ]0.5 .

(18)

It should also be noted that for a given value of R > Rc , which corresponds to the supercritical convection state, the expression (18) for |ε| is maximum, among all the solutions of the nonlinear problem, only if R20 has the smallest positive value for particular solution, and in this sense such particular solution is referred to as the preferred solution in this paper. The sign of R20 determines whether the oscillatory solution exists for values of R above or below Rc . For R20 < 0, which corresponds to the subcritical convection state, R < Rc and |ε| decreases with increasing R, which corresponds to an unstable state as our stability analysis actually indicated. In the present problem the coefficient R20 is due to the nonlinear convective terms in the temperature equation and the nonlinear interactions between the flow velocity and the non-uniform and nonlinear permeability associated with the perturbation to the basic state solid fraction. 4.2.1. Oscillatory hexagons We begin the presentation of the results and the corresponding discussion for the cases of oscillatory solutions in the form of hexagons. Using (15a), the coefficients R20 for the solutions in the form of simple-travelling hexagons, standing hexagons and general-travelling hexagons, were computed for S = 0.25C and various values of C, K 1 and K 2 . Since the expression for R20 contains a number of definite integrals, simple numerical integrations of type of an extended trapezoidal rule were carried out to determine the numerical values of R20 with sufficient accuracy. The sign (o) of R20 determines whether the oscillatory solution exists for values of R above or below Rc . Here we are interested to study the preferred oscillatory hexagonal type of solutions, which, as to be presented later in this section, turn out to be stable, for particular range of values of the parameters, and correspond to the relatively lowest values of R. These types of oscillatory solutions should correspond to the smallest values of R20 > 0. We calculated the values of R20 for the oscillatory solutions for S = 0.25C but for different values of C, K 1 and K 2 , where the experimentally relevant range of S and C referred to earlier in this section could be covered. In all the calculations that we carried out, we found that simple-travelling hexagons have mostly smaller values of R20 as compared with other types of hexagonal solutions. The value of R20 increases with K 2 , and for C ≥ 0.1 the value of positive R20 is generally non-monotonic with respect to C, often decreases with either increasing K 1 or decreasing C and eventually R20 becomes negative. The non-monotonic property of R20 with respect to C is found to increase sharply for the value of C beyond 0.1. When the coefficient R20 just changes sign for particular values of the parameters, there is neighboring points in parameter space where such coefficient is positive but has very small magnitude and, thus, the corresponding solution is often preferred. Our generated data for R20 for the above stated three types of solutions and for different parameter values indicate that the effect of K 2 is generally stabilizing in the sense that the value of R20 increases with K 2 . Depending on the parameter values and relative to the solutions in the form of standing hexagons and general-travelling hexagons, the solution in the form of simple-travelling hexagons is mostly preferred. In addition, the solution in the form of the simple-travelling hexagons is found to become the preferred and stable solution among all the other two- and three-dimensional solutions that have been investigated over particular range of the parameter values. For particular values of the parameters, such evidence can be seen from Table 1 which presents the preferred and stable oscillatory solutions for S = 0.25C, K 2 = 0 and for selected values of K 1 and C. The results presented in the Table 1 indicate that sensitivity of the permeability with respect to the local solid fraction favors three-dimensional types of solution. About the effects of S, C and K 1 on the preferred hexagons, some results are presented in Fig. 1 for R20 of the simpletravelling hexagonal solution for S = 0.25C, K 2 = 0 and for several values of K 1 . As can be seen from this figure, the effect of K 1 is destabilizing, and the stabilizing effect of C dominates over the destabilizing effect of S. Our other generated data for wider range of values for C indicates that the permeability parameter K 1 has a non-monotonic effect on the solution in the sense that it is destabilizing in the most of the range for C and stabilizing elsewhere. The present investigation is restricted to the range K 1 < 0.1 and K 2 < 0.1 since for value of each of these parameters beyond this range, the order of magnitude of the coefficient R20 increases rapidly with either K 1 or K 2 becoming much larger than unity by about at least 3 orders of magnitude, so that the modeling assumption of the type (7), which assumes that the coefficients, such as |R20 |, in those double expansions in powers of ε and δ be of order unity can no longer be justified.

217

D.N. Riahi / Nonlinear Analysis: Real World Applications 10 (2009) 209–226

Fig. 1. R20 versus C(S = 0.25C) for the oscillatory solution in the form of simple-travelling hexagons. Here the graphs are for K 2 = 0 and for three different values of K 1 = 0.045 (black solid line with circle symbol), 0.065 (red dashed-dotted line with square symbol) and 0.085 (blue dotted line with triangle symbol). Table 1 Preferred supercritical oscillatory solutions for S = 0.25C, K 2 = 0 and different values of C and K 1

0.0355 0.0392 0.0450 0.0521 0.0610 0.0717 0.0845 0.0995 0.1166 0.1357 0.1565 0.1786 0.2013 0.2241 0.2461

0.000

0.025

0.045

0.065

0.085

4 4 4 4 4 4 4 4 4 4 4 4 4 4 4

7d 1 1 1 1 1 7d 7d 1 7d 1 10 10 7a 7d

12(.4) 7b 7b 7d 7d 1 1 1 1 8d 1 1 1 7a 9a (.4)

3(.2) 3(.2) 9d (.4) 10 7d 7d 3(.4) 1 1 1 3(.2) 1 1 9a (.3) 8a

9d (.4) 9d (.4) 7c 7c 7d 7d 3(.3) 1 1 10 3(.1) 3(.3) 7b 9a (.1) 9a (.1)

Here the simple-travelling, standing and general-travelling solutions for hexagons are designated, respectively, by the numbers 1, 2 and 3; for rolls are designated, respectively, by the numbers 4, 5 and 6; for rectangles are designated, respectively, by the numbers 7, 8 and 9; and for squares are designated, respectively, by the numbers 10, 11 and 12. The first column and the first row present, respectively, the values of C and K 1 . The value of the parameter b for the general-travelling solution is given in parentheses. The subscripts a, b, c and d represent, respectively, the angles γ = 15◦ , 30◦ , 50◦ and 70◦ .

We also examined the vertical distribution of solid fraction at different time and location in the horizontal direction for the preferred oscillatory hexagons. Our generated data, over most of the periodic domain in time, at centers and at the nodes of such solution indicated useful information about magnitude and the sign of the perturbation to solid fraction at the nodes and at the centers of the cells. Some typical results are presented in Figs. 2 and 3 for the vertical distribution of the basic state (solid line) and total solid fraction at both a center (dashed line) and a node (dotted line) of a preferred standing hexagonal solution. In these calculations δ = 0.15, S = 0.25C, C = 0.1565, K 1 = 0.085, K 2 = 0.0, R20 = 12.46, ω01 = 5.0 and the value ε is chosen to be the maximum value of |ε| = 0.001 beyond which the solid fraction becomes negative. This is based on the physical ground that the value of the perturbation to the solid fraction cannot be such that total solid fraction becomes negative. It is seen from the Fig. 2, which is drawn at the beginning of the period of oscillation 2π/ω01 , that there is mostly tendency for chimney formation at the node and solid dendrite formation at the center of the cell. Fig. 3, which is drawn at a time about half of the period of

218

D.N. Riahi / Nonlinear Analysis: Real World Applications 10 (2009) 209–226

Fig. 2. Solid fraction versus z for the oscillatory solution in the form of a preferred standing hexagons with ε = 0.001, C = 0.1565, S = 0.25C, K 1 = 0.085 and K 2 = 0. Here the black solid line with circle symbol, red dashed-dotted line with square symbol and blue dotted line with ˜ ˜ triangle symbol present the basic solid fraction φ B , φa(x = y = t = 0) and φa(x = 4/3, y = t = 0), respectively.

Fig. 3. The same as the Fig. 2 but for half-period of oscillation later.

oscillation later, shows mostly an opposite picture to that of the Fig. 2. That is, the Fig. 3 shows that there is mostly tendency for channel formation at the center and solid dendrite formation at the node of the same cell. This result may indicate a beneficial effect of the oscillatory mode to at least reduce the tendency for channel formation in the cell. This result is a consequence of the fact that oscillatory hexagons can be preferred only if K 1 is not too small (Table 1), and it indicates that the dependence of the permeability with respect to the local solid fraction may favor tendency for reducing chimney formation within the mushy layer. Similar to the discussion provided in R02 for the solid fraction due to the standing modes, the chimneys and the compositional strips are in the vertical direction since the phase speed of these modes is zero, but the vertical extent of the chimneys can vary depending on the variation of the solid fraction perturbation with respect to time. However, in the case of a travelling wave type of solution, the chimneys and the compositional strips in a travelling wave state can be inclined because of the non-zero values of the phase speed of such a wave relative to the uniform upward speed of the mushy layer ([3]; R02). 4.2.2. Oscillatory rolls Here again R20 is the important coefficient for the nonlinear aspects of the rolls, and the solutions can be in the form of simple-travelling rolls, standing rolls or general-travelling rolls. The coefficient R20 was computed for various

D.N. Riahi / Nonlinear Analysis: Real World Applications 10 (2009) 209–226

219

Fig. 4. R20 versus C (S = 0.25C) for the oscillatory solutions in the form of simple-travelling rolls (black solid line with circle symbol), standing rolls (red dashed-dotted line with square symbol) and general-travelling rolls with b = 0.3 (blue dotted line with triangle symbol). The graphs are for K 1 = 0.065 and K 2 = 0.05.

values of the parameters. As was the case for all the oscillatory solutions, the effect of K 2 was found to be stabilizing in the sense that R20 increases with K 2 . The effect of K 1 on R20 was found to be mostly stabilizing in the sense that this coefficient mostly increases with increasing K 1 . In addition, it was found that the rate of increase of R20 with respect to K 2 is much larger than the rate at which R20 increases with K 1 . For the least stabilizing case where K 2 = 0, the main results about the coefficient R20 are as follows. For not too large value of C the coefficient R20 is generally positive and, thus, the flow due to rolls is supercritical. For sufficiently small value of K 1 , the preferred solution is that due to simple-travelling rolls (Table 1). Our generated data for the value of R20 indicated that depending on the parameter values, the value of R20 for simple travelling rolls or standing rolls can be smaller or larger than that for general-travelling rolls. Some typical results for the oscillatory rolls in the form of simple-travelling, standing and a type of general-travelling waves (b = 0.3) are given in Fig. 4. It can be seen from this figure that the destabilizing effect of S dominates over the stabilizing effect of C for the lower range of values of C, while the stabilizing effect of C dominates over the destabilizing effect of S for the upper range of values for C. In addition, the simple-travelling rolls are preferred over the other two types of rolls in the lower range of values of C, general-travelling rolls become preferred over the other two rolls-types in a small intermediate range of values of C, and standing rolls are preferred over the other two types of rolls in an upper range of values of C. Comparing values of R20 for oscillatory rolls and hexagons for different values of the parameters that we have determined so far, we find that oscillatory hexagons and particularly simple-travelling hexagons are preferred mostly over the oscillatory rolls if K 1 is not too small (Table 1). As to be discussed later in this section, all the subcritical oscillatory solutions with R20 < 0 are found to be unstable and, therefore, such solutions are not preferred. 4.2.3. Oscillatory rectangles Here we consider oscillatory rectangular solutions, which include those sets of classes of solutions with γ 6= 90◦ and those oscillatory squares with γ = 90◦ . Again, as in the case of oscillatory rolls and hexagons, the important coefficient for the nonlinear effects is R20 for the solutions in the form of simple-travelling rectangles, standing rectangles, general-travelling rectangles, simple-travelling squares, standing squares and general-travelling squares. The coefficient R20 for each of these solutions was computed for various values of the parameters. For C beyond some value, the effect of K 1 on the values of this coefficient for such solutions was found to be non-monotonic in the sense that the effect of this parameter can be either stabilizing or destabilizing depending on the value of K 1 , so that the values of this coefficient can increase and then decrease, or vice versa, with increasing K 1 . For the least stabilizing case where K 2 = 0, the main results about R20 are as follows. The coefficient R20 is mostly positive and, thus, the flows due to the rectangular solutions are mostly supercritical. The effect of C on the value of R20 is generally non-monotonic for C > 0.1 in the sense that C has stabilizing or destabilizing effect in different range of

220

D.N. Riahi / Nonlinear Analysis: Real World Applications 10 (2009) 209–226

Fig. 5. R20 versus C (S = 0.25C) for the oscillatory solution in the form of simple-travelling squares for K 2 = 0 and three values of K 1 = 0.025 (black solid line with circle symbol), 0.045 (red dashed-dotted line with square symbol) and 0.065 (blue dotted line with triangle symbol).

Fig. 6. R20 versus C (S = 0.25C) for the oscillatory solutions in the form of simple-travelling rectangles with three values of the angle γ = 70◦ (black solid line with circle symbol), 30◦ (red dashed-dotted line with square symbol) and 15◦ (blue dotted line with triangle symbol). The graphs are for K 1 = 0.045 and K 2 = 0.

values for C > 0.1, and such stabilizing or destabilizing effect can be quite sharp close to certain values of C. The simple-travelling rectangles are mostly preferred as compared to other types of oscillatory rectangles. Some results about the dependence of R20 for the solutions in the form of simple-travelling squares on C and K 1 are given in Fig. 5, where R20 versus C is plotted for S = 0.25C, K 2 = 0 and for several values of K 1 . It can be seen from this figure that the stabilizing effect of C dominates over the destabilizing effect S and K 1 is destabilizing over the considered range of values of C. The parameter K 1 is destabilizing over most of the investigated range of values of C, except near the upper range of values of C where K 1 is stabilizing. Fig. 6 presents R20 of simple-travelling rectangles versus C for S = 0.25C, K 1 = 0.045, K 2 = 0.0 and for three different values of the angle γ . It can be seen from this figure that γ is destabilizing, and the stabilizing effect of C dominates over the destabilizing effect of S if γ is not too small. Also subcritical domain seems to be diminished with decreasing γ . For the rectangular solution with γ = 70◦ , where a subcritical domain is shown in the figure, there are values of C at which R20 can be sufficiently small and positive leading to the preference of such solution. Comparing the values of R20 for the least stable case and for different two- and three-dimensional oscillatory solutions that we investigated so far, we find that over most of the values of the parameters, simple-travelling rectangles and simple-

D.N. Riahi / Nonlinear Analysis: Real World Applications 10 (2009) 209–226

221

Fig. 7. R20 versus C (S = 0.25C) for the oscillatory solutions in the form of simple-travelling rolls (black solid line with circle symbol), simpletravelling rectangles with γ = 70◦ (red dashed-dotted line with square symbol) and simple-travelling hexagons (blue dotted line with triangle symbol). The graphs are for K 1 = 0.025 and K 2 = 0.

travelling hexagons are mostly preferred over other types of solutions. Some typical results are provided in Fig. 7, which presents R20 for simple-travelling solutions in the form of rolls, rectangles with γ = 70◦ and hexagons versus C for S = 0.25C, K 1 = 0.025 and K 2 = 0.0. It can be seen from this figure that hexagons are the preferred pattern over most of the investigated range of values in C, while rectangles replace hexagons only for small intermediate range of values in C or if C is sufficiently small. Our further generated data indicates that the values of C near which the destabilizing effect of S or the stabilizing effect of C can be strong within small intervals in C, are different for different types of flow pattern. 4.3. Stability of finite-amplitude oscillatory solutions Following standard stability procedure [4], the systems for the growth rate σ ∗ of the disturbances acting on the finite-amplitude oscillatory solutions have been simplified, and the expression for σ ∗ has been computed for different types of solutions. In all the cases that have been investigated only supercritical oscillatory solutions in the form of rolls, rectangles and hexagons are found to be possibly stable in particular range of the values for the non-dimensional parameters. For sufficiently small K 1 , supercritical simple-travelling rolls are stable over most of the domain of the values of C, while supercritical standing rolls are stable over a small domain of relatively large but not too large values of C. The region of stable standing rolls widen somewhat with the increasing values of K 2 . For not too small value of K 1 and depending on the values of the parameters, supercritical simple-travelling hexagons or supercritical simple-travelling rectangles can mostly be stable and be preferred over the rest of the detected solutions. There are small regions in the parameter space where the other three-dimensional supercritical solutions in the form of standing hexagons, generaltravelling hexagons, standing rectangles or general-travelling rectangles can be stable and be preferred over other types of solution. No subcritical solution was found to be stable. 4.4. Comparison with the work and results in [10] Very recently weakly nonlinear two-dimensional oscillatory convection in horizontal mushy layers was investigated analytically by Guba and Worster [10] using the [1] type of model in the limit of large Stefan number. Guba and Worster [10] studied two-dimensional oscillatory modes in the form of only simple-travelling or standing rolls. Depending on the parameter values, travelling or standing rolls were found to be supercritically stable. The authors concluded that for weak sensitivity of permeability on the solid fraction, travelling rolls were more likely to be selected. These results are all in agreement with the corresponding ones in R02. Their other result that decrease in permeability promoted subcritical state was already speculated in R02 and is similar to that due to Amberg and Homsy [1] for steady rolls.

222

D.N. Riahi / Nonlinear Analysis: Real World Applications 10 (2009) 209–226

Guba and Worster [10] considered order one values of K 1 for simple-travelling and standing rolls, but they did not consider general-travelling rolls. They also determined stability regions in (S, C)-space for both two-dimensional oscillatory and steady cases and found region of unstable rolls. Their two-dimensional analysis could not determine the preferred solutions in the region where their investigated oscillatory rolls were subcritically unstable. We now compare the results due to Guba and Worster [10] to those in this paper for the oscillatory modes. In the present study of the oscillatory modes we investigated both two- and three-dimensional states with the effects of K 1 included, determined the corresponding two- and three-dimensional solutions and their stability with respect to threedimensional disturbances. For particular parameter values, we found those two-dimensional oscillatory solutions that were predicted stable in [10], were, in fact, unstable to three-dimensional oscillatory solutions (Table 1) for the same parameter values. The other results reported in [10] for sufficiently small values of K 1 , were generally in agreement with those in the present study and in qualitative agreement with those in R02. 5. Conclusion We investigated the problem of nonlinear convection, due to oscillatory modes, in horizontal mushy layers during the solidification of binary alloys. We carried out both two- and three-dimensional finite amplitude and stability investigations for different values of the parameters including the permeability parameters K 1 and K 2 . We found that, if K 1 is not too small and depending on the values of the other parameters, three-dimensional oscillatory solutions mostly in the form of supercritical simple-travelling rectangles or supercritical simple-travelling hexagons can be the stable and preferred form of the flow solutions. No subcritical solution was found to be stable. Relatively strong sensitivity of the permeability with respect to the local solid fraction was found to favor preference of threedimensional types of solution. Our results about the vertical distribution of the solid fraction for the hexagonal solutions indicated possible beneficial effect of the oscillatory modes in reducing the tendency for chimney formation within the mushy layer. Appendix The system of equations and boundary conditions at order ε are given below ∆2 (−∇ 2 V10 + R00 θ10 + R10 θ00 ) = K 1 {∇ 2 (φ00 ∆2 V00 ) + (∂/∂z)[(∂ 2 V00 /∂ x∂z)(∂φ00 /∂ x) + (∂ 2 V00 /∂ y∂z)(∂φ00 /∂ y) + π 2 V00 (∂φ00 /∂z)], ∇ 2 θ10 + G(R00 ∆2 V10 + R10 ∆2 V00 + R00 Ω V00 .∇θ00 ) = 0, S[(∂/∂t1 ) − (∂/∂z)]φ10 + ∇ 2 θ10 + R00 ∆2 V10 = −R10 ∆2 V00 − (Sω11 /ω01 )∂φ00 /∂t1 + R00 Ω V00 .∇θ00 , V10 = θ10 = 0 at z = 0, V10 = θ10 = φ10 = 0 at z = 1, where t1 = ωt/ω01 }. (o)

(o)

(0)

(1)

(A.1)

The coefficient functions Blp and Elp , which are introduced in (13a) and (13b), have the following expressions: (o)

(o)

(0)

(1)

(2)

(2)

(Blp , Elp ) = (Blp , Elp ) + (Blp , Elp ) cos(2π z) + (Blp , Elp ) sin(2π z) (3)

(3)

(4)

(4)

+ (Blp , Elp ) cos(π z) exp[iω01 S p (z − 1)] + (Blp , Elp ) sin(π z) exp[iω01 S p (z − 1)]. (A.2a) Here (0) (0) (1) (1) (2) 2√ 2√ 2√ (Blp , Elp ) = (1, alp G)C1 exp(rlp z) + (1, alp G)C2 exp(−rlp z) + (1, −alp G)C3 exp(rlp z) (2) (2) 2√ + (1, −alp G)C4 exp(−rlp z) + (1, −GR00 )L lp for ψlp 6= ±1,

where √ (1) 2 2 ψlp = al .a p /π 2 , alp ≡ 2π 2 (1 + ψlp ), rlp ≡ (alp + G R00 alp )0.5 , √ √ (2) (2) 2 2 2 rlp ≡ (alp − G R00 alp )0.5 , L lp ≡ iω01 K 1 S p π 3 /[CG G(π 2 − ω01 )(alp − GR200 )], (1)

(1)

(3)

(1)

(1)

C1 ≡ [Mlp exp(−rlp ) − Mlp ]/[exp(rlp ) − exp(−rlp )], (1)

(1)

(3)

(1)

(1)

C2 ≡ [Mlp exp(rlp ) − Mlp ]/[exp(rlp ) − exp(−rlp )],

(A.2b)

223

D.N. Riahi / Nonlinear Analysis: Real World Applications 10 (2009) 209–226 (2)

(2)

(4)

(2)

(2)

C3 ≡ [Mlp exp(−rlp ) − Mlp ]/[exp(rlp ) − exp(−rlp )], (2)

(2)

(4)

(2)

(2)

C4 ≡ [Mlp exp(rlp ) − Mlp ]/[exp(rlp ) − exp(−rlp )], √ (1) (2) 2 (2) (3) (4) (3) (1) (3) Mlp ≡ 0.5[ G R00 L lp /alp − L lp + L lp + L lp ], L lp ≡ −Blp − Blp exp(−iω01 S p ), √ (4) (1) (3) (2) (2) 2 (2) (3) (4) 2 L lp ≡ [−Elp − Elp exp(−iω01 S p )]/alp , Mlp ≡ 0.5[− G R00 L lp /alp − L lp + L lp − L lp ], √ (3) (2) 2 (2) (5) (6) (5) (1) (3) Mlp ≡ 0.5[ G R00 L lp /alp − L lp + L lp + L lp ], L lp ≡ −Blp + Blp , √ (6) (1) (3) (4) (2) 2 (2) (5) (6) 2 L lp ≡ (−Elp + Elp )/alp , Mlp ≡ 0.5[− G R00 L lp /alp − L lp + L lp − L lp ]; (A.2c) (0)

(6)

(4)

(4)

(1)

(3)

(1)

(6)

(4)

2 2 2 Blp = R00 [alp (L lp − L lp )z 3 /6 + alp L lp z 2 /2] + L lp z 2 /2 + {−L lp − 0.5L lp − R00 [alp (L lp − L lp )/6 (4)

(5)

(3)

2 + alp L lp /2] + L lp }z + L lp (0) Elp

=

(6) 2 alp (L lp

−

(4) L lp )z

for ψlp = −1,

2 (4) + alp L lp

for ψlp = −1,

(A.2d)

where √ (1) 2 L lp ≡ −iω01 π 3 S p K 1 /[CG G(π 2 − ω01 )]; (0) Blp (0) Elp

= =

(A.2e)

(1) (1) (1) C5 exp(rlp z) + C6 exp(−rlp z) + C7 z + C8 + 0.25L lp z for ψlp = 1, √ √ (1) (1) (1) 2π G[C5 exp(rlp z) + C6 exp(−rlp z) − C7 z − C8 ] − G L lp (0.5π z 2

+ 0.25/π )

for ψlp = 1,

(A.2f)

where (8)

(7)

(1)

(1)

(1)

C5 ≡ [L lp − L lp exp(−rlp )]/[exp(rlp ) − exp(−rlp )], √ (7) (3) (4) 2 (1) L lp ≡ 0.5L lp + L lp alp /(4π G) + L lp /(16π 2 ), √ √ √ (8) (5) (6) 2 (1) L lp ≡ 0.5L lp + L lp alp /(4π G) − L lp [1/8 − π G/4 − G/(8π )], (7)

(1)

(8)

(1)

(1)

C6 ≡ [L lp exp(rlp ) − L lp ]/[exp(rlp ) − exp(−rlp )], (1)

(1)

(1)

(5)

(7)

(3)

C7 ≡ −C5 exp(rlp ) − C6 exp(−rlp ) − C8 − 0.25L lp + L lp , C8 ≡ −L lp + L lp ; √ (1) 3 2 2 2 2 2 2 Blp = {−2iω01 S p K 1 π /[CG G(π − ω01 )]}/{alp + 4π − GR00 alp /(alp + 4π 2 )}, (1)

(A.2g)

(1)

2 2 Elp = −GR00 alp Blp /(alp + 4π 2 );

(A.2h) √ (2) 2 2 2 2 Blp = {π 2 (1 − ψlp ) − 2K 1 π 4 (alp + 4π 2 )/[R00 CG G(π 2 − ω01 )]}/{GR00 alp − (alp + 4π 2 )2 /R00 }, √ (2) (2) 2 2 Elp = 2K 1 π 4 /[R00 CG G(π 2 − ω01 )] − (alp + 4π 2 )Blp /R00 ; (A.2i) √ (3) (4) (3) 2 2 3 Elp = Elp = 0 for ψlp = −1, Blp = −iω01 S p K 1 π 3 (ω01 − 5π 2 )/[CG G(π 2 − ω01 ) ] √ (4) 2 2 3 for ψlp = −1, Blp = K 1 π 4 (3π 2 + ω01 )/[CG G(π 2 − ω01 ) ] for ψlp = −1; (A.2j) (3)

(4)

(3)

(4)

(Elp , Elp ) = [(C9 C13 − C11 C15 ), (C10 C15 − C9 C12 )]/(C10 C13 − C11 C12 )

for ψlp 6= −1,

(3)

2 2 (Blp , Blp ) = {Elp [(−alp − π 2 − ω01 ), (−2iω01 S p π )] (4)

2 2 2 + Elp [(2iω01 S p π ), (−alp − π 2 − ω01 )]}/(GR00 alp )

for ψlp 6= −1,

(A.2k)

√ 2 2 C14 ≡ −3K 1 π 4 alp /[CG G(π 2 − ω01 )].

(A.2l)

where √ 2 2 C9 ≡ iω01 S p K 1 π 3 alp /[CG G(π 2 − ω01 )], 2 2 2 2 2 C10 ≡ C13 ≡ −alp R00 + [(alp + π 2 + ω01 ) + 4ω01 π ]/GR00 , 2 2 C11 ≡ −4iω01 S p π(π 2 + ω01 + alp )/GR00 ≡ −C12 ,

224

D.N. Riahi / Nonlinear Analysis: Real World Applications 10 (2009) 209–226 (i j)

+ , g − and f The coefficient functions f m+ , f m− , gm m lp expressions:

(i, j = 0, 1, 2), which are given in (13), have the following

2 2 2 f m+ ≡ f m = {−2π 3 ω11 /[CG(π 2 − ω01 ) ]}{iSm [(π 2 + ω01 )/π ] sin(π z) + 2ω01 cos(π z) 2 2 + [2ω01 + iSm (z − 1)(π − ω01 )] exp[iω01 Sm (z − 1)]}, f m− = f m∗ ; √ + 2 gm ≡ gm = {−R10 π/[C G(π 2 − ω01 )]}{π cos(π z) + iω01 Sm sin(π z) + π exp[iω01 Sm (z − 1)]},

(A.3a)

− ∗ gm = gm ;

(A.3b)

(i j) flp

=

(i j) (i j) (i j) (i j) (i j) −hHlp exp(−h lp z)i1 exp(−h lp z) + hHlp exp(h lp z)iz

(i j) exp(−h lp z)

(i, j = 0, 1, 2),

(A.3c)

where Z h f iz ≡ (00)

(h lp , (10)

(h lp , (01) (h lp , (02) (h lp , (12) (h lp , (22)

z

Z f dz,

0 (00) Hlp ) (10) Hlp ) (01) Hlp ) (02) Hlp ) (12) Hlp ) (22)

≡ ≡ ≡

1

f dz,

h f i1 ≡

0 2 (s) {0, [R00 alp Blp − π 2 sin(2π z)(1 − ψlp )/G]/C}, (00) (20) (20) (00) {−iω01 Sl , Hlp }, (h lp , Hlp ) ≡ {iω01 Sl , Hlp }, 2 (o) {−iω01 S p , [R00 alp Blp − π 2 sin(2π z)(1 − ψlp )/G]/C}, (01)∗

≡ {h lp

(01)∗

, Hlp

(11)

(01)∗

≡ {−iω01 (Sl − S p ), Hlp (11)∗

(h lp , Hlp ) ≡ {h lp

(01)∗

, Hlp

(11)

(01)

(h lp , Hlp ) ≡ {−iω01 (Sl + S p ), Hlp },

},

},

(21)

(21)

(12)∗

(h lp , Hlp ) ≡ {h lp

(01)

, Hlp }, (A.3d)

}.

The system of equations and boundary conditions at order ε 2 is {∆2 (−∇ 2 V20 + R00 θ20 + R20 θ00 ) = ∂/∂z[K 1 (Ω V00 .∇φ10 + Ω V10 .∇φ00 ) 2 2 + K 2 Ω V00 .∇(φ00 )] + ∇ 2 [K 1 (φ10 ∆2 V00 + φ00 ∆2 V10 ) + K 2 φ00 ∆2 V00 ], ∇ 2 θ20 + G∆2 (R00 V20 + R20 V00 ) = GR00 (Ω V00 .∇θ10 + Ω V10 .∇θ00 ), S(∂/∂t1 − ∂/∂z)φ20 + ∆2 (R00 V20 + R20 V00 ) + ∇ 2 θ20 = −(Sω21 /ω01 )(∂/∂t1 )φ00 − (Sω11 /ω01 )(∂/∂t1 )φ10 + R00 (Ω V00 .∇θ10 + Ω V10 .∇θ00 ), θ20 = V20 = 0 at z = 0, θ20 = V20 = φ20 = 0 at z = 1}.

(A.4)

The solvability conditions for the system (A.4) are reduced to the following equations: X √ (1) 2 − 2 + + + + + + + − + + + − + + + R20 Gπ(|A+ {[Elpm (A+ n | + |An | ) = n Am Al A p hηn ηm ηl η p i + An Am Al A p hηn ηm ηl η p i) l, p,m (2)

− + + + − + + − − + + − − + + + Elpm (A+ n Am Al A p hηn ηm ηl η p i + An Am Al A p hηn ηm ηl η p i) (3)

+ + − + + + − − + + − − + + − + Elpm (A+ n Am Al A p hηn ηm ηl η p i + An Am Al A p hηn ηm ηl η p i) (4)

− + − + − + − − − + − − − + − + Elpm (A+ n Am Al A p hηn ηm ηl η p i + An Am Al A p hηn ηm ηl η p i) (5)

+ − + + + − + − + − + − + − + + Elpm (A+ n Am Al A p hηn ηm ηl η p i + An Am Al A p hηn ηm ηl η p i) (6)

− − + + − − + − − − + − − − + + Elpm (A+ n Am Al A p hηn ηm ηl η p i + An Am Al A p hηn ηm ηl η p i) (7)

+ − − + + − − − + − − − + − − + Elpm (A+ n Am Al A p hηn ηm ηl η p i + An Am Al A p hηn ηm ηl η p i) (8)

− − − + − − − − − − − − − − − + Elpm (A+ n Am Al A p hηn ηm ηl η p i + An Am Al A p hηn ηm ηl η p i)]},

(n = −N , . . . , −1, 1, . . . , N ).

(A.5a)

225

D.N. Riahi / Nonlinear Analysis: Real World Applications 10 (2009) 209–226 (i)

Here the coefficient Elpm (i = 1, . . . , 8) given in (A.5a) have the following expressions: (1)

(11)

(11)

Elpm = K 1 π 2 (ψml + ψmp )hcos2 (π z) flp i − 0.5K 1 πhsin(2π z)(d f lp /dz)i √ √ 2 (o) (o) hcos(π z)Blp (d f m /dz)i + K 1 Gπ 2 (ψml + ψmp )hcos(π z)Blp f m i − K 1 Galp + 2K 2 π 2 ψlp hcos2 (π z) f m f p i − K 2 πhsin(2π z) f m (d f p /dz)i √ (o) (11) 2 2 )(Blp f m )i − K 1 hsin(π z)(d 2 /dz 2 − alpm )[sin(π z) flp ]i − (K 1 G/π )hsin(π z)(d 2 /dz 2 − alpm √ (o) 2 − K 2 hsin(π z)(d 2 /dz 2 − alpm )[sin(π z) fl f p ]i + 0.5 Gπ 2 R00 (ψml + ψmp )hsin(2π z)Elp i √ (o) (o) − Gπ R00 hsin2 (π z)(dE lp /dz)i − GR00 π 2 (ψml + ψmp )hsin2 (π z)(dB lp /dz)i (o)

2 + 0.5GR00 πalp hsin(2π z)Blp i,

2 alpm ≡ 2π 2 (1.5 + ψlp + ψlm + ψ pm );

(A.5b)

(2)

the expression for Elpm has the same form as the one given in (A.5b), provided f m is replaced by f m∗ ; the expression (3)

(11)

(o)

(o)

for Elpm has the same form as the one given in (A.5b), provided flp , Blp , Elp and f p are replaced, respectively, (12)

(o)∗

(o)∗

by flp , Blp , Elp

(4)

(3)

and f p∗ ; the expression for Elpm is the same as the one for Elpm , provided f m is replaced by f m∗ ;

(5)

(1)

(11)

the expression for Elpm is the same as the one for Elpm , provided flp fl∗ ; the expression for

(6) Elpm

is the same as the one for and f p∗ ; the expression for The expressions for the

(21)

and fl are replaced, respectively by flp

(1)

(1)

(1)

(1)

(1)

(1)

(o1) Tnm = E −m,−n,m δnm + (E −m,−n,m + E −n,−m,m )δn,−m + (E −m.−n,m + E −n,−m,m + E m,−m,−n ) × (1 − δnm )(1 − δn,−m ), (2)

(2)

(2)

(2)

(2)

(3)

(3)

(5)

(3)

(5)

=

(o5) Tnm =

(5) E −m,−n,m δnm

(4) + (E −m,−n,m

+

(6) E −n,−m,m )δn,−m

(4) + (E −m,−n,m

+

(3) E −n,−m,m )δn,−m

(5) + (E −m,−n,m

+

(6) E −n,−m,

+

(3) E −n,−m,m

+

× (1 − δnm )(1 − δn,−m ), (5) + (E −m,−n,m (6)

+

(3) E m,−m,−n )

(A.6e) (4)

(6)

(4)

(7)

(o6) Tnm = E −m,−n,m δnm + (E −m,−n,m + E −n,−m,m )δn,−m + (E −m,−n,m + E −n,−m,m + E m,−m,−n ) × (1 − δnm )(1 − δn,−m ), (7)

(7)

(7)

(7)

(7)

(8)

(8)

(8)

(A.6f)

(4)

(o7) Tnm = E −m,−n,m δnm + (E −m,−n,m + E −n,−m,m )δn,−m + (E −m,−n,m + E −n,−m,m + E m,−m,−n ) × (1 − δnm )(1 − δn,−m ), (8)

(A.6c)

(6) E m,−m,−n )

(A.6d)

× (1 − δnm )(1 − δn,−m ), (6)

(A.6b)

(2)

(o3) Tnm = E −m,−n,m δnm + (E −m,−n,m + E −n,−m,m )δn,−m + (E −m,−n,m + E −n,−m,m + E m,−m,−n ) × (1 − δnm )(1 − δn,−m ), (4) E −m,−n,m δnm

(A.6a)

(5)

(o2) Tnm = E −m,−n,m δnm + (E −m,−n,m + E −n,−m,m )δn,−m + (E −m,−n,m + E −n,−m,m + E m,−m,−n ) × (1 − δnm )(1 − δn,−m ),

(o4) Tnm

and

(5) (7) is the same as the one for Elpm , provided f m is replaced by f m∗ ; the expression for Elpm (5) (21) (o) (o) (22) (o)∗ (o)∗ Elpm , provided flp , Blp , Elp and f p are replaced, respectively, by flp , Blp , Elp (8) (7) Elpm is the same as the one for Elpm , provided f m is replaced by f m∗ . (oi) coefficients Tnm (i = 1, . . . , 8), which were introduced in (15a), are given below

(8)

(A.6g)

(8)

(o8) Tnm = E −m,−n,m δnm + (E −m,−n,m + E −n,−m,m )δn,−m + (E −m,−n,m + E −n,−m,m + E m,−m,−n ) × (1 − δnm )(1 − δn,−m ),

(A.6h)

where δnm = 1

for n = m

and

0 for n 6= m.

The stability system is given below ˜ (φc)∆2 V + K (φ)∆ ˜ 2 V 0 ] + (∂/∂z){εΩ V.∇[φ 0 (d/dφ)K ˜ (φ)] ˜ ∇ 2 [εφ 0 (d/dφ)K

(A.6i)

226

D.N. Riahi / Nonlinear Analysis: Real World Applications 10 (2009) 209–226

˜ − R∆2 θ 0 = 0, + Ω V 0 .∇ K (φ)}

(A.7a)

(∂/∂t + σ − δ∂/∂z)(−θ + Sφ /δ) + R(dθ B /dz)∆2 V + ∇ θ = ε R(Ω V.∇θ + Ω V .∇θ ), (∂/∂t + σ − δ∂/∂z)[(−1 + φ B )θ 0 + θ B φ 0 + εφθ 0 + εφ 0 θ − Cφ 0 /δ] + R(dθ B /dz)∆2 V 0 0

0

= ε R(Ω V.∇θ 0 + Ω V 0 .∇θ ), V = θ0 = 0 0

at z = 0,

V 0 = θ 0 = φ0 = 0

at z = 1.

0

2 0

0

0

(A.7b) (A.7c) (A.7d) (A.7e)

References [1] G. Amberg, G.M. Homsy, Nonlinear analysis of buoyant convection in binary solidification with application to channel formation, J. Fluid Mech. 252 (1993) 79–98. [2] D.M. Anderson, M.G. Worster, Weakly nonlinear analysis of convection in mushy layers during the solidification of binary alloys, J. Fluid Mech. 302 (1995) 307–331. [3] D.M. Anderson, M.G. Worster, A new oscillatory instability in a mushy layer during the solidification of binary alloys, J. Fluid Mech. 307 (1996) 245–267. [4] F.H. Busse, The stability of finite amplitude convection and its relation to an extremum principal, J. Fluid Mech. 30 (1967) 625–649. [5] F.H. Busse, Patterns of convection in spherical shells, J. Fluid Mech. 72 (1975) 67–85. [6] S. Chandrasekhar, Hydrodynamic and Hydromagnetic Stability, Clarendon, 1961. [7] C.A. Chung, F. Chen, Onset of plume convection in mushy layers, J. Fluid Mech. 408 (2000) 53–82. [8] A.C. Fowler, The formation of freckles in binary alloys, IMA J. Appl. Math. 35 (1985) 159–174. [9] P. Guba, On the finite amplitude steady convection in rotating mushy layers, J. Fluid Mech. 437 (2001) 337–365. [10] P. Guba, M.G. Worster, Nonlinear oscillatory convection in mushy layers, J. Fluid Mech. 553 (2006) 419–443. [11] D.N. Riahi, Nonlinear convection in a porous layer with finite conducting boundaries, J. Fluid Mech. 129 (1983) 153–171. [12] D.N. Riahi, On nonlinear convection in mushy layers. Part 1. Oscillatory modes of convection, J. Fluid Mech. 467 (2002) 331–359. [13] T.P. Schulze, M.G. Worster, Weak convection, liquid inclusions and the formation of chimneys in mushy layers, J. Fluid Mech. 388 (1999) 197–215. [14] S. Tait, K. Jahrling, C. Jaupart, The planform of compositional convection and chimney formation in a mushy layer, Nature 359 (1992) 406–408. [15] M.G. Worster, Instabilities of liquid and mushy regions during solidification of alloys, J. Fluid Mech. 237 (1992) 649–669.