- Email: [email protected]
Stability and perturbation analyses for nonlinear convection in a mushy layer and in the presence of viscous heating
Mechanics Research Communications 39 (2012) 18–22
Contents lists available at SciVerse ScienceDirect
Mechanics Research Communications journal homepage: www.elsevier.com/locate/mechrescom
Stability and perturbation analyses for nonlinear convection in a mushy layer and in the presence of viscous heating D.N. Riahi ∗ Department of Mathematics, 1201 W. University Dr., University of Texas-Pan American, Edinburg, TX 78539-2999, USA
a r t i c l e
i n f o
Article history: Received 13 August 2011 Received in revised form 11 October 2011 Available online 19 October 2011 MSC: 76Exx 76Sxx 76R10 Keywords: Convection Mushy layer Convective flow Viscous heating Nonlinear flow
a b s t r a c t This paper studies nonlinear steady convection in a horizontal mushy layer and in the presence of viscous dissipation in the temperature equation, which often is referred to as viscous heating. We consider the appropriate system of equations and the associated boundary conditions for the flow in the mushy layer. Under certain assumptions and conditions, we determine the weakly nonlinear solutions of the resulting system and their stability by using perturbation and stability analyses. We find, in particular, that for a wide range of values of the viscous heating parameter and sufficiently small amplitude of the flow, the only stable convective flow is in the form of subcritical down-hexagons with down-flow at the cells’ centers and up-flow at the cells’ boundaries. This result appears to agree with the available experimental results for the flow pattern near the onset of motion in a mushy layer. © 2011 Elsevier Ltd. All rights reserved.
1. Introduction Since the first single layer model of the flow in the mushy layer developed by Amberg and Homsy (1993), there have been a number of studies of convection in mushy layers that used such type of model including those in Anderson and Worster (1995), Chung and Chen (2000), Riahi (2002, 2003), Roper et al. (2008) and Guba and Worster (2010) in various cases to understand the flow and its planform in such systems. All such studied did not include the effect of a viscous heating term, which is usually very small in magnitude for moderate flow speed, in the temperature equation. However, these studies did not satisfactorily explain the realized planform of the convection in a mushy layer that was observed by Tait et al. (1992). The experimental study by these authors is the only one that has been reported so far for weak convection in a mushy layer, which successfully detected the corresponding flow pattern. These authors cooled and solidified a 28% ammonium chloride–water solution from below very slowly and carefully so that they were able to observe the realized flow pattern in the form of down-hexagonal cells with down-flow at the cells’ centers and up-flow at the cells’ boundaries.
∗ Tel.: +1 956 665 7063; fax: +1 956 665 5091. E-mail addresses: [email protected], [email protected] 0093-6413/$ – see front matter © 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.mechrescom.2011.10.002
In the present study, we take into account the presence of viscous heating ˚ in the temperature equation (Nield, 2007) for the first time in order to determine the preferred flow that takes place at the lowest possible value of the controlling parameter (Rayleigh number R) in the flow systems. Although ˚ is generally very small and is neglected in convective flow system under the Boussinesq approximation (Chandrasekhar, 1961), it turns out that it can become significant under the present modeling formulation and scaling and for sufficiently small amplitude ε of convection near onset of motion. Our main result is the stability of subcritical down-hexagons for sufficiently small amplitude of motion and under a wide range of values of the viscous dissipation parameter, which appears to agree with the experimental observation due to Tait et al. (1992). 2. Governing system We consider a horizontal layer of a binary alloy melt of some constant composition C0 and temperature T∞ , which is cooled from below and is solidified at a constant speed V0 , with the eutectic temperature Te at the position z = 0 held fixed in a frame moving with the solidification speed in the vertical z-direction. Here z-axis is anti-parallel to the gravity vector. Within the layer of melt, there is a mushy layer adjacent to the solidifying surface and of thickness d, where the solid dendrites and the liquid melt coexist (Fig. 1).
D.N. Riahi / Mechanics Research Communications 39 (2012) 18–22
19
large value of the Lewis number. The governing equations (1) are subject to the following boundary conditions: + 1 = w = 0 at z = 0 and
= w = = 0 at z = ı,
(2a-b)
where ı = dV0 / is the dimensionless depth of the mushy layer. Following Amberg and Homsy (1993) and Anderson and Worster (1995) in reducing the model asymptotically, we assume certain rescaling to be given below in the limit of sufficiently small ı Cr =
C , ε ı 1, (x, y, z) = ı(x, y, z), t = ı2 t, R2 = ıR, ı
˜ =0+ u
εR ı
u(x, y, z, t), P = RPB (z) + εRP(x, y, z, t),
D = ıD, = B (z) + ε(x, y, z, t), = B (z) + ε(x, y, z, t), (3) Fig. 1. Mushy layer system.
Here x and y are the horizontal variables along the x- and y-axes in the horizontal plane z = 0. In this paper we focus on the so-called mushy-layer mode (Worster, 1992) at the onset of convection for the flow within the mushy layer only. We consider the equations for momentum, continuity, heat and solute for the mushy region (0 < z < d) in the frame oxyz moving with the solid-mush interface (z = 0) at speed V0 . The governing system is non-dimensionalized by using V0 , /V0 , /V02 , ˇC0 g/V0 , C and T as scales for velocity, length, time, pressure, solute and temperature, respectively. Here is the thermal diffusivity, 0 is a reference (constant) density, ˇ = ˇ* − M˛*, ˛* and ˇ* are the expansion coefficients for the heat and solute, respectively and M is assumed to be constant, C = C0 − Ce , T = TL (C0 ) − Te , and Te is the eutectic temperature. The non-dimensional form of the equations for momentum, continuity, temperature and solute concentration in the mushy layer, which is assumed to be a porous layer and under Darcy’s law (Worster, 1992), are ˜ = −∇ P − Rz, K()u
(1a)
∇ · u˜ = 0,
(1b)
∂ ∂ − ∂t ∂z ∂ ∂ − ∂t ∂z
˜ · ∇ 2 + ˚, ( − St ) + u
(1c)
(1d)
˜ = ux + vy + wz = (1 − )U is the volume flux per unit area, where u ˚ is the non-dimensional form of the viscous heating (Nield, 2007) given by ˚ = DK()[u2 + v2 + w2 ],
K() = 1 + K1 + K2 2 + . . . ,
(4)
where K1 and K2 are constant, and as in Anderson and Worster (1995), K1 = ıKc is assumed, where Kc is a constant of at most order one. For the analysis to be presented in the next section, it is found convenient to use the representation u = GV + E , GV ≡ ∇ × ∇ × (V z), E
˜ · ∇ 2 = 0, [(1 − ) + Cr ] + u
where C is an order one quantity as ı → 0, D is the scaled viscous dissipation parameter, which is assumed to be at most of order one, and the quantities with subscript B are the basic flow variables for the motionless state and are assumed to be a function of the vertical variable only. The small deviation of each dependent variable from its basic quantity is measured by a perturbation amplitude ε and is, in general, depends on the spatial and time variables. The rescaling given in (3) is then used in the governing system (1) and (2). The resulting system of equations and boundary conditions admits a motionless basic state, which is steady and horizontally uniform. The basic state motionless solution is then determined to leading orders in terms of the asymptotic expansions for ı 1, which turns out to be the same as in Anderson and Worster (1995) and will not be repeated here. Since the basic state solution for the solid fraction is found to be small of order ı (Anderson and Worster, 1995), the following expansion for K() will be used later in the governing system
(1e)
is the local solid fraction within the mushy layer, U is the veloc˜ along the ity vector, u and v are the horizontal components of u x- and y-directions, respectively, x and y are unit vectors along the positive x- and y-directions, w is the vertical component of ˜ along the z-direction, is the non-dimensional temperature, u R = ˇCg˘(0)/(V0 ) is the Rayleigh number, ˘(0) is the reference value at = 0 of the permeability ˘() of the porous medium, is the kinematic viscosity, g is acceleration due to gravity, K() ≡ ˘(0)/˘(), D ≡ /[˘(0)Cp T] is the viscous dissipation parameter, Cp is specific heat at constant pressure, St = L/(Cl T) is the Stefan number, Cl is the specific heat per unit volume, L is the latent heat of solidification per unit volume, Cr = (Cs − C0 )/C 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
≡ ∇ × ( z),
(5)
for the divergent-free vector u (Chandrasekhar, 1961), where V and are the poloidal and toroidal functions for u, respectively. Taking the vertical component of the curl of (1a) it can be shown that the toroidal part E of u must vanish. Hence, we set = 0 hereafter. Taking the vertical component of the double curl of (1a), using (1b), (3) and (5) in (1) and (2), we find the following system, which will be analyzed in the next section:
∇ [K(B + ε)2 V ] + 2
∂ ∂z
[GV · ∇ K(B + ε)] − R2 = 0, (6a)
∂ ı∂ − ∂t ∂z
(− + St ) + R
dB dz
2 V − ∇ 2 + εR2 ıDK(
+ε)(GV · GV ) = εR(GV · ∇ ),
∂ ∂ − ∂t ∂z
(−1 + B ) + B + ε −
= εR(GV · ∇ ),
(6b)
C +R ı
dB dz
2 V (6c)
20
D.N. Riahi / Mechanics Research Communications 39 (2012) 18–22
= V = 0 at z = 0,
(6d)
= = V = 0 at z = 1,
(6e)
where 2 ≡ (∂2 /∂x2 + ∂2 /∂y2 ). 3. Analyses and results Since the analyses for the small finite-amplitude steady solutions and their stability are essentially based on the same method of approach provided in Riahi (2003), here we briefly describe the key aspects of this paper in the context of the present problem and refer the readers to Riahi (2003) for additional details of the analyses and the procedure. The analysis to determine the finiteamplitude solutions for the dependent variables is based on the double-series expansions in powers of the two small parameters ı and ε which satisfy the condition given in (3). The appropriate expansions are for V, , and R. In the following descriptions, the coefficient of f ≡ (V, , , R) in the order εm ın of such double-series expansions is designated by fmn . Considering the system (6) to the lowest order in ε, we find the linear problem. At order ε0 ı0 , the linear system yield the following results
(V00 , 00 , 00 ) =
2 + a2 − C
2 + a2 R00 a2
sin( z), − sin( z),
N
[1 + cos( z)]
An Wn ,
(7a)
n=−N
Table 1 Values of F and R11 for D = 3(10−4 ), C = 2, St = 3, K2 = 0, ı = 0.001 and several values of Kc . Kc
F
R11
0 0.2 0.4 0.7
−1.2701 −2.4789 −3.6877 −5.5008
5.1832 3.9744 2.7656 0.9525
which is the expression of the coefficient for R in the order ε0 ı1 at a = ac . Hence, the critical Rayleigh number Rc can be written as Rc = R00c + ıR01c + O(ı2 ).
Next, we analyze the nonlinear problem for the finite-amplitude solutions. The solvability conditions for the nonlinear systems require the following special solutions V00n and 00n of the linear system:
(V00n , 00n ) =
2 + a2 , −1 sin( z)An Wn . R00 a2
2 R00
(7b)
Here a is the horizontal wave number, Wn = exp(i an ·r), i is the pure imaginary number, subscript n takes only non-zero integer values from −N to N, N is a positive integer representing the number of distinct modes, r is the position vector, and the horizontal wave number vectors an satisfy the properties an · z = 0, |an | = a, a−n = −an .
(8)
The coefficients An are constants and satisfy the conditions N
An A∗n = 1,
A∗n = A−n ,
(9a)
n=−N
where the asterisk indicates complex conjugate. We restrict our attention to the simplest types of solutions, which include rolls (N = 1), squares (N = 2), rectangles (N = 2) and hexagons (N = 3), where (9a) reduces to |A1 |2 = . . . = |AN |2 =
1 . 2N
(9b)
Minimizing the expression for R00 with respect to the wave number a, we find R00c = 2 , ac = ,
(10a)
where R00c is the minimum value of R00 achieved at a = ac . We now consider the governing system (6) in the order ε0 ı1 . Multiplying Eq. (6a) by −V00 and (6b) by 00 , adding these two equations, making use of the corresponding boundary conditions and averaging over the mushy layer, we find the condition for the existence of the solutions in this order. This condition yields R01c =
−R00c St , 2C
(10b)
(11)
Consider the governing system (6) at order ε. Multiplying Eqs. (6a) and (6b) by −V00n and 00n , respectively, adding, applying the corresponding boundary conditions and averaging over the mushy layer, we find the solvability condition in this order which yields R10 = 0. However, the solvability condition for the system (6) in the next leading order εı was found to provide an equation for the next leading order coefficient R11 for the contribution of R in the order εı, which is found to be non-zero only for the hexagonal modes. Its expression for hexagons is of the form
2
( 2 + a2 ) = . a2
(10c)
R11 = F(C, St , Kc , K2 ) +
16 2 √ 3 6
D,
(12)
where the expression for F, which depends in general on the parameters C, St , Kc and K2 , has the same form as in the case where the effect of the viscous dissipation in not taken into account. The expression for F is very lengthy and will not be given here, but it was determined in Riahi (2003) in the zero limit of the rotation parameter, and its sign is negative regardless of the values of the parameters C, St , K2 and Kc . Some typical results about the values of F and R11 > 0 for the convective flow in the form of subcritical down-hexagons, where ε, εR11 and R are all negative, are given in Table 1 for ı = 0.001, D = 0.0003 (D ≡ ıD), C = 0.09, St = 3.2, K2 = 0 and for several values of Kc . The values of C and St that are chosen here correspond closely to those values in the laboratory experiments (Tait et al., 1992) using aqueous solutions of ammonium chloride. It can be seen from the results given in Table 1 and in (12) that the values of R11 are in fact positive here for a much wider range of the values of the viscous heating parameter D ≥ 3(10−4 ). However, if the viscous heating effect is not taken into account, then (12) and the results given in Table 1 imply that the values of R11 should be negative implying existence of subcritical up-hexagons, where ε > 0 but both εR11 and R are negative, instead of down-hexagons. Thus, we have shown that simply subcritical down-hexagons can exist instead of subcritical up-hexagon if we take into account the presence of the viscous dissipation for very small values of the viscous dissipation parameter of order 10−4 or higher. Using (3) for the relation between the viscous dissipation parameter D and its rescaled form D, we find from (12) that subcritical down-hexagons can exist if D>
−ıF . 21.511
(13)
Later in this section we explain that our stability results of the finite-amplitude solutions indicated that for sufficiently small
D.N. Riahi / Mechanics Research Communications 39 (2012) 18–22
amplitude of motion the subcritical down-hexagons are the only stable solution in the present problem if such solutions can exist as, for example, in the case provided in Table 1. It is also of interest to note that such a small value of order 10−4 for the viscous heating parameter D in a porous medium can correspond to the value estimated in the past (Chandrasekhar, 1961) for the viscous heating in an ordinary medium, where the viscous dissipation term is based on the velocity gradients terms (Chandrasekhar, 1961), under which viscous heating terms in the temperature equation were neglected for thermal convection problems obeyed by the Boussinesq approximation. Next, the leading order solutions to the ε-system are found. These solutions are functions of independent variables and the dimensionless parameters. We now consider the system (6) in order ε2 . The solvability conditions for the system in the order ε2 yield the expression of the coefficient R20 for R in the order ε2 which together with R11 were used to study the steady secondary solutions. Thus, using the result that R10 is zero and taking into account the coefficients R11 and R20 , which are determined in the present study, we have the expression for R up to the second order in ε in the form R = Rc + εıR11 + ε2 R20 .
(14)
It can be seen from (14) that these coefficients represent leading contributions to the change in R required for obtaining finiteamplitude ε for a nonlinear solution. As stated earlier, we have restricted the analyses to the simplest types of solutions in the form of rolls, squares, rectangles and hexagons. In later references for a solution in the form of rectangles, an angle , which is less than 90◦ , is defined to be the angle between two adjacent wave number vectors of any rectangular cell. For the linear system infinite number of steady solutions is possible, but such linear degeneracy is reduced significantly by the consideration of the nonlinear system of the problem. For the nonlinear case, the preferred flow structure is considered to be the one due to the solutions that correspond to the lowest value of R. Further more, the amplitude |ε| for the preferred solution, which can depend on R20 and R11 , is found to increase with R, a result that corresponds to the observation. To distinguish the physically realizable solution from all possible steady solutions, the stability of the finite-amplitude solutions V, and is investigated by superposing on the finite-amplitude solutions perturbations of infinitesimal amplitude and with the addition of a time dependence of the form exp( t), where is the growth rate. We then obtain the stability system for the disturbances with arbitrary wave number vectors but under the restriction that the magnitude of all such vectors is the same as ac . This system was solved by an expansion similar to that for the finite-amplitude system presented before. The growth rate to the lowest order in ε of the most critical disturbances is found to be zero. In the order ε of the stability system, we found that the corresponding solvability condition yields 10 = 0. Next, we applied the solvability conditions in the orders εı and ε2 of the stability system to determine 11 and 20 . Since 11 and 20 may not in general be zero for a particular solution, we define ∗ = εı 11 + ε2 20 ,
(15)
as the leading order growth rate and combine the solvability conditions in the orders εı and ε2 to obtain the system for * whose very lengthy expression will not be given here. The growth rate * of the disturbances acting on the simple finite-amplitude steady solutions in the form of rolls, squares, rectangles with different angles and hexagons are then determined. Our stability results verified that those bifurcation branches for which the amplitude decreases with increasing R are unstable and,
21
thus, not physically realizable. So the information about the coefficients for R in leading orders in ε given in (14) was found to be useful in the sense that they could indicate preference and stability of particular solutions and vice versa. For the case where the coefficient R11 for R in the order ε is zero, which was found to correspond to steady non-hexagonal type solutions, then we found that R20 > 0 for rolls, rectangles and squares. Thus, such solutions are supercritical, where R > Rc , and so the amplitude of the convective flow is largest, provided the value of R20 is smallest among all such solutions to the nonlinear problem. In addition, R20 is also positive for hexagons. The value of R20 was found to be smallest in the case of rolls. For the case of hexagons, where R11 is non-zero and given by (12), we found that subcritical down-hexagons, where R11 > 0, can exist over a wide range of values of the viscous dissipation and permeability parameters as the results in Table 1 provides an example. However, up-hexagons with up-flow at the cells’ centers and down-flow at the cells’ boundaries are supercritical over the domain where subcritical down-hexagons exist. We also calculated the values of the solid fraction at the center and at any node of a down-hexagonal cell and an up-hexagonal cell for different values of the parameters. We found, in particular that regardless of the values of the parameters, the total solid fraction at the node of the down-hexagonal cell is smaller than the corresponding basic state value, while the result is reversed in the case of the up-hexagonal cell. These results indicate that there is tendency for chimney formation at the nodes of down-hexagons, which is in agreement with the experimental observation due to Tait et al. (1992), while in the case of up-hexagons there is tendency for solid fraction enhancement at the nodes. The growth rate * given by (15) of the disturbances acting on the finite-amplitude solutions has been computed for rolls, squares, rectangles with different angles and hexagons. In all the cases that we investigated for different range of values of the parameters and the small values of |ε|, we found the following stability results: Supercritical rolls are stable only if √ 3ı|R11 | (16a) , |ε| ≥ ε1 , ε1 = H where the expression for H, which is positive, is very length and will not be given here. Down (up)-hexagons are stable if ε2 ≤ |ε| ≤ ε3 , ε3 =
ε2 =
0.5ı|R11 | ε0 , = R20 2
3ı|R11 | , R11 > (<)0, H
(16b)
where R20 refers to this coefficient for the hexagons and ε0 corresponds to the boundary between the subcritical and supercritical domains of hexagons where R = Rc . As can be seen from the expressions for the positive quantities εi (i = 0, 1, 2, 3) given in (9a) and 9(b), the stable hexagons are subcritical for smallest values of |ε| and them become supercritical for larger values of |ε|, and there is a domain of overlap where both supercritical rolls and supercritical hexagons are stable. As an example, we evaluated the quantities εi (i = 0, 1, 2, 3) for C = 0.09, St = 3.2, Kc = 0, K2 = 0 and ı = 0.001 at the viscous dissipation parameter D = (10−4 ). We found R11 = 0.881, ε0 = 0.000040, ε1 = 0.00017, ε2 = 0.000020 and ε3 = 0.00029. Thus, here the subcritical domain for the stable down-hexagons is 0.000020 ≤ |ε| ≤ 0.000040, the supercritical domain for the stable down-hexagons is 0.000040 < |ε| ≤ 0.00029, and the supercritical domain for the stable rolls is |ε| ≥ 0.00017. There is a region of overlap for 0.00017 ≤ |ε| ≤ 0.00029 where both supercritical downhexagons and supercritical rolls are stable and hysteresis effects can exist.
22
D.N. Riahi / Mechanics Research Communications 39 (2012) 18–22
4. Conclusion We studied nonlinear steady convective flow in a mushy layer and in the presence of viscous dissipation. Using perturbation and stability analyses, we first determined the weakly nonlinear steady solutions in the form of rolls, rectangles, squares, down-hexagons and up-hexagons, and then we determined the stability of such solutions with respect to arbitrary disturbances. We found, in particular, that over a wide range of values of the viscous dissipation parameter, the only stable solution is in the form of subcritical down-hexagons for sufficiently small amplitude of the convection, and there is tendency for chimney formation at the nodes of the stable down-hexagons. These results are in agreement with the experimental observation (Tait et al., 1992). References Amberg, G., Homsy, G.M., 1993. Nonlinear analysis of buoyant convection in binary solidification with application to channel formation. J. Fluid Mech. 252, 79–98.
Anderson, D.M., Worster, G.M., 1995. Weakly nonlinear analysis of convection in mushy layers during the solidification of binary alloys. J. Fluid Mech. 302, 307–331. Chandrasekhar, S., 1961. Hydrodynamic and Hydromagnetic Stability. Clarenton Press, Oxford, UK. Chung, C.A., Chen, F., 2000. Onset of plume convection in mushy layers. J. Fluid Mech. 408, 53–82. Guba, P., Worster, M.G., 2010. Interaction between steady and oscillatory convection in mushy layers. J. Fluid Mech. 645, 411–434. Nield, D.A., 2007. The modeling of viscous dissipation in a saturated porous medium. J. Heat Transfer 129, 1459–1460. Riahi, D.N., 2002. On nonlinear convection in mushy layers. Part 1: oscillatory modes of convection. J. Fluid Mech. 467, 331–359. Riahi, D.N., 2003. Nonlinear steady convection in rotating mushy layers. J. Fluid Mech. 485, 279–306. Roper, S.M., Davis, S.H., Voorhees, P.W., 2008. An analysis of convection in a mushy layer with a deformable permeable interface. J. Fluid Mech. 596, 333–352. Tait, S., Jahrling, K., Jaupart, C., 1992. The planform of compositional convection and chimney formation in a mushy layer. Nature 359, 406–408. Worster, M.G., 1992. Instabilities of the liquid and mushy regions during solidification of alloys. J. Fluid Mech. 237, 649–669.
Contents lists available at SciVerse ScienceDirect
Mechanics Research Communications journal homepage: www.elsevier.com/locate/mechrescom
Stability and perturbation analyses for nonlinear convection in a mushy layer and in the presence of viscous heating D.N. Riahi ∗ Department of Mathematics, 1201 W. University Dr., University of Texas-Pan American, Edinburg, TX 78539-2999, USA
a r t i c l e
i n f o
Article history: Received 13 August 2011 Received in revised form 11 October 2011 Available online 19 October 2011 MSC: 76Exx 76Sxx 76R10 Keywords: Convection Mushy layer Convective flow Viscous heating Nonlinear flow
a b s t r a c t This paper studies nonlinear steady convection in a horizontal mushy layer and in the presence of viscous dissipation in the temperature equation, which often is referred to as viscous heating. We consider the appropriate system of equations and the associated boundary conditions for the flow in the mushy layer. Under certain assumptions and conditions, we determine the weakly nonlinear solutions of the resulting system and their stability by using perturbation and stability analyses. We find, in particular, that for a wide range of values of the viscous heating parameter and sufficiently small amplitude of the flow, the only stable convective flow is in the form of subcritical down-hexagons with down-flow at the cells’ centers and up-flow at the cells’ boundaries. This result appears to agree with the available experimental results for the flow pattern near the onset of motion in a mushy layer. © 2011 Elsevier Ltd. All rights reserved.
1. Introduction Since the first single layer model of the flow in the mushy layer developed by Amberg and Homsy (1993), there have been a number of studies of convection in mushy layers that used such type of model including those in Anderson and Worster (1995), Chung and Chen (2000), Riahi (2002, 2003), Roper et al. (2008) and Guba and Worster (2010) in various cases to understand the flow and its planform in such systems. All such studied did not include the effect of a viscous heating term, which is usually very small in magnitude for moderate flow speed, in the temperature equation. However, these studies did not satisfactorily explain the realized planform of the convection in a mushy layer that was observed by Tait et al. (1992). The experimental study by these authors is the only one that has been reported so far for weak convection in a mushy layer, which successfully detected the corresponding flow pattern. These authors cooled and solidified a 28% ammonium chloride–water solution from below very slowly and carefully so that they were able to observe the realized flow pattern in the form of down-hexagonal cells with down-flow at the cells’ centers and up-flow at the cells’ boundaries.
∗ Tel.: +1 956 665 7063; fax: +1 956 665 5091. E-mail addresses: [email protected], [email protected] 0093-6413/$ – see front matter © 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.mechrescom.2011.10.002
In the present study, we take into account the presence of viscous heating ˚ in the temperature equation (Nield, 2007) for the first time in order to determine the preferred flow that takes place at the lowest possible value of the controlling parameter (Rayleigh number R) in the flow systems. Although ˚ is generally very small and is neglected in convective flow system under the Boussinesq approximation (Chandrasekhar, 1961), it turns out that it can become significant under the present modeling formulation and scaling and for sufficiently small amplitude ε of convection near onset of motion. Our main result is the stability of subcritical down-hexagons for sufficiently small amplitude of motion and under a wide range of values of the viscous dissipation parameter, which appears to agree with the experimental observation due to Tait et al. (1992). 2. Governing system We consider a horizontal layer of a binary alloy melt of some constant composition C0 and temperature T∞ , which is cooled from below and is solidified at a constant speed V0 , with the eutectic temperature Te at the position z = 0 held fixed in a frame moving with the solidification speed in the vertical z-direction. Here z-axis is anti-parallel to the gravity vector. Within the layer of melt, there is a mushy layer adjacent to the solidifying surface and of thickness d, where the solid dendrites and the liquid melt coexist (Fig. 1).
D.N. Riahi / Mechanics Research Communications 39 (2012) 18–22
19
large value of the Lewis number. The governing equations (1) are subject to the following boundary conditions: + 1 = w = 0 at z = 0 and
= w = = 0 at z = ı,
(2a-b)
where ı = dV0 / is the dimensionless depth of the mushy layer. Following Amberg and Homsy (1993) and Anderson and Worster (1995) in reducing the model asymptotically, we assume certain rescaling to be given below in the limit of sufficiently small ı Cr =
C , ε ı 1, (x, y, z) = ı(x, y, z), t = ı2 t, R2 = ıR, ı
˜ =0+ u
εR ı
u(x, y, z, t), P = RPB (z) + εRP(x, y, z, t),
D = ıD, = B (z) + ε(x, y, z, t), = B (z) + ε(x, y, z, t), (3) Fig. 1. Mushy layer system.
Here x and y are the horizontal variables along the x- and y-axes in the horizontal plane z = 0. In this paper we focus on the so-called mushy-layer mode (Worster, 1992) at the onset of convection for the flow within the mushy layer only. We consider the equations for momentum, continuity, heat and solute for the mushy region (0 < z < d) in the frame oxyz moving with the solid-mush interface (z = 0) at speed V0 . The governing system is non-dimensionalized by using V0 , /V0 , /V02 , ˇC0 g/V0 , C and T as scales for velocity, length, time, pressure, solute and temperature, respectively. Here is the thermal diffusivity, 0 is a reference (constant) density, ˇ = ˇ* − M˛*, ˛* and ˇ* are the expansion coefficients for the heat and solute, respectively and M is assumed to be constant, C = C0 − Ce , T = TL (C0 ) − Te , and Te is the eutectic temperature. The non-dimensional form of the equations for momentum, continuity, temperature and solute concentration in the mushy layer, which is assumed to be a porous layer and under Darcy’s law (Worster, 1992), are ˜ = −∇ P − Rz, K()u
(1a)
∇ · u˜ = 0,
(1b)
∂ ∂ − ∂t ∂z ∂ ∂ − ∂t ∂z
˜ · ∇ 2 + ˚, ( − St ) + u
(1c)
(1d)
˜ = ux + vy + wz = (1 − )U is the volume flux per unit area, where u ˚ is the non-dimensional form of the viscous heating (Nield, 2007) given by ˚ = DK()[u2 + v2 + w2 ],
K() = 1 + K1 + K2 2 + . . . ,
(4)
where K1 and K2 are constant, and as in Anderson and Worster (1995), K1 = ıKc is assumed, where Kc is a constant of at most order one. For the analysis to be presented in the next section, it is found convenient to use the representation u = GV + E , GV ≡ ∇ × ∇ × (V z), E
˜ · ∇ 2 = 0, [(1 − ) + Cr ] + u
where C is an order one quantity as ı → 0, D is the scaled viscous dissipation parameter, which is assumed to be at most of order one, and the quantities with subscript B are the basic flow variables for the motionless state and are assumed to be a function of the vertical variable only. The small deviation of each dependent variable from its basic quantity is measured by a perturbation amplitude ε and is, in general, depends on the spatial and time variables. The rescaling given in (3) is then used in the governing system (1) and (2). The resulting system of equations and boundary conditions admits a motionless basic state, which is steady and horizontally uniform. The basic state motionless solution is then determined to leading orders in terms of the asymptotic expansions for ı 1, which turns out to be the same as in Anderson and Worster (1995) and will not be repeated here. Since the basic state solution for the solid fraction is found to be small of order ı (Anderson and Worster, 1995), the following expansion for K() will be used later in the governing system
(1e)
is the local solid fraction within the mushy layer, U is the veloc˜ along the ity vector, u and v are the horizontal components of u x- and y-directions, respectively, x and y are unit vectors along the positive x- and y-directions, w is the vertical component of ˜ along the z-direction, is the non-dimensional temperature, u R = ˇCg˘(0)/(V0 ) is the Rayleigh number, ˘(0) is the reference value at = 0 of the permeability ˘() of the porous medium, is the kinematic viscosity, g is acceleration due to gravity, K() ≡ ˘(0)/˘(), D ≡ /[˘(0)Cp T] is the viscous dissipation parameter, Cp is specific heat at constant pressure, St = L/(Cl T) is the Stefan number, Cl is the specific heat per unit volume, L is the latent heat of solidification per unit volume, Cr = (Cs − C0 )/C 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
≡ ∇ × ( z),
(5)
for the divergent-free vector u (Chandrasekhar, 1961), where V and are the poloidal and toroidal functions for u, respectively. Taking the vertical component of the curl of (1a) it can be shown that the toroidal part E of u must vanish. Hence, we set = 0 hereafter. Taking the vertical component of the double curl of (1a), using (1b), (3) and (5) in (1) and (2), we find the following system, which will be analyzed in the next section:
∇ [K(B + ε)2 V ] + 2
∂ ∂z
[GV · ∇ K(B + ε)] − R2 = 0, (6a)
∂ ı∂ − ∂t ∂z
(− + St ) + R
dB dz
2 V − ∇ 2 + εR2 ıDK(
+ε)(GV · GV ) = εR(GV · ∇ ),
∂ ∂ − ∂t ∂z
(−1 + B ) + B + ε −
= εR(GV · ∇ ),
(6b)
C +R ı
dB dz
2 V (6c)
20
D.N. Riahi / Mechanics Research Communications 39 (2012) 18–22
= V = 0 at z = 0,
(6d)
= = V = 0 at z = 1,
(6e)
where 2 ≡ (∂2 /∂x2 + ∂2 /∂y2 ). 3. Analyses and results Since the analyses for the small finite-amplitude steady solutions and their stability are essentially based on the same method of approach provided in Riahi (2003), here we briefly describe the key aspects of this paper in the context of the present problem and refer the readers to Riahi (2003) for additional details of the analyses and the procedure. The analysis to determine the finiteamplitude solutions for the dependent variables is based on the double-series expansions in powers of the two small parameters ı and ε which satisfy the condition given in (3). The appropriate expansions are for V, , and R. In the following descriptions, the coefficient of f ≡ (V, , , R) in the order εm ın of such double-series expansions is designated by fmn . Considering the system (6) to the lowest order in ε, we find the linear problem. At order ε0 ı0 , the linear system yield the following results
(V00 , 00 , 00 ) =
2 + a2 − C
2 + a2 R00 a2
sin( z), − sin( z),
N
[1 + cos( z)]
An Wn ,
(7a)
n=−N
Table 1 Values of F and R11 for D = 3(10−4 ), C = 2, St = 3, K2 = 0, ı = 0.001 and several values of Kc . Kc
F
R11
0 0.2 0.4 0.7
−1.2701 −2.4789 −3.6877 −5.5008
5.1832 3.9744 2.7656 0.9525
which is the expression of the coefficient for R in the order ε0 ı1 at a = ac . Hence, the critical Rayleigh number Rc can be written as Rc = R00c + ıR01c + O(ı2 ).
Next, we analyze the nonlinear problem for the finite-amplitude solutions. The solvability conditions for the nonlinear systems require the following special solutions V00n and 00n of the linear system:
(V00n , 00n ) =
2 + a2 , −1 sin( z)An Wn . R00 a2
2 R00
(7b)
Here a is the horizontal wave number, Wn = exp(i an ·r), i is the pure imaginary number, subscript n takes only non-zero integer values from −N to N, N is a positive integer representing the number of distinct modes, r is the position vector, and the horizontal wave number vectors an satisfy the properties an · z = 0, |an | = a, a−n = −an .
(8)
The coefficients An are constants and satisfy the conditions N
An A∗n = 1,
A∗n = A−n ,
(9a)
n=−N
where the asterisk indicates complex conjugate. We restrict our attention to the simplest types of solutions, which include rolls (N = 1), squares (N = 2), rectangles (N = 2) and hexagons (N = 3), where (9a) reduces to |A1 |2 = . . . = |AN |2 =
1 . 2N
(9b)
Minimizing the expression for R00 with respect to the wave number a, we find R00c = 2 , ac = ,
(10a)
where R00c is the minimum value of R00 achieved at a = ac . We now consider the governing system (6) in the order ε0 ı1 . Multiplying Eq. (6a) by −V00 and (6b) by 00 , adding these two equations, making use of the corresponding boundary conditions and averaging over the mushy layer, we find the condition for the existence of the solutions in this order. This condition yields R01c =
−R00c St , 2C
(10b)
(11)
Consider the governing system (6) at order ε. Multiplying Eqs. (6a) and (6b) by −V00n and 00n , respectively, adding, applying the corresponding boundary conditions and averaging over the mushy layer, we find the solvability condition in this order which yields R10 = 0. However, the solvability condition for the system (6) in the next leading order εı was found to provide an equation for the next leading order coefficient R11 for the contribution of R in the order εı, which is found to be non-zero only for the hexagonal modes. Its expression for hexagons is of the form
2
( 2 + a2 ) = . a2
(10c)
R11 = F(C, St , Kc , K2 ) +
16 2 √ 3 6
D,
(12)
where the expression for F, which depends in general on the parameters C, St , Kc and K2 , has the same form as in the case where the effect of the viscous dissipation in not taken into account. The expression for F is very lengthy and will not be given here, but it was determined in Riahi (2003) in the zero limit of the rotation parameter, and its sign is negative regardless of the values of the parameters C, St , K2 and Kc . Some typical results about the values of F and R11 > 0 for the convective flow in the form of subcritical down-hexagons, where ε, εR11 and R are all negative, are given in Table 1 for ı = 0.001, D = 0.0003 (D ≡ ıD), C = 0.09, St = 3.2, K2 = 0 and for several values of Kc . The values of C and St that are chosen here correspond closely to those values in the laboratory experiments (Tait et al., 1992) using aqueous solutions of ammonium chloride. It can be seen from the results given in Table 1 and in (12) that the values of R11 are in fact positive here for a much wider range of the values of the viscous heating parameter D ≥ 3(10−4 ). However, if the viscous heating effect is not taken into account, then (12) and the results given in Table 1 imply that the values of R11 should be negative implying existence of subcritical up-hexagons, where ε > 0 but both εR11 and R are negative, instead of down-hexagons. Thus, we have shown that simply subcritical down-hexagons can exist instead of subcritical up-hexagon if we take into account the presence of the viscous dissipation for very small values of the viscous dissipation parameter of order 10−4 or higher. Using (3) for the relation between the viscous dissipation parameter D and its rescaled form D, we find from (12) that subcritical down-hexagons can exist if D>
−ıF . 21.511
(13)
Later in this section we explain that our stability results of the finite-amplitude solutions indicated that for sufficiently small
D.N. Riahi / Mechanics Research Communications 39 (2012) 18–22
amplitude of motion the subcritical down-hexagons are the only stable solution in the present problem if such solutions can exist as, for example, in the case provided in Table 1. It is also of interest to note that such a small value of order 10−4 for the viscous heating parameter D in a porous medium can correspond to the value estimated in the past (Chandrasekhar, 1961) for the viscous heating in an ordinary medium, where the viscous dissipation term is based on the velocity gradients terms (Chandrasekhar, 1961), under which viscous heating terms in the temperature equation were neglected for thermal convection problems obeyed by the Boussinesq approximation. Next, the leading order solutions to the ε-system are found. These solutions are functions of independent variables and the dimensionless parameters. We now consider the system (6) in order ε2 . The solvability conditions for the system in the order ε2 yield the expression of the coefficient R20 for R in the order ε2 which together with R11 were used to study the steady secondary solutions. Thus, using the result that R10 is zero and taking into account the coefficients R11 and R20 , which are determined in the present study, we have the expression for R up to the second order in ε in the form R = Rc + εıR11 + ε2 R20 .
(14)
It can be seen from (14) that these coefficients represent leading contributions to the change in R required for obtaining finiteamplitude ε for a nonlinear solution. As stated earlier, we have restricted the analyses to the simplest types of solutions in the form of rolls, squares, rectangles and hexagons. In later references for a solution in the form of rectangles, an angle , which is less than 90◦ , is defined to be the angle between two adjacent wave number vectors of any rectangular cell. For the linear system infinite number of steady solutions is possible, but such linear degeneracy is reduced significantly by the consideration of the nonlinear system of the problem. For the nonlinear case, the preferred flow structure is considered to be the one due to the solutions that correspond to the lowest value of R. Further more, the amplitude |ε| for the preferred solution, which can depend on R20 and R11 , is found to increase with R, a result that corresponds to the observation. To distinguish the physically realizable solution from all possible steady solutions, the stability of the finite-amplitude solutions V, and is investigated by superposing on the finite-amplitude solutions perturbations of infinitesimal amplitude and with the addition of a time dependence of the form exp( t), where is the growth rate. We then obtain the stability system for the disturbances with arbitrary wave number vectors but under the restriction that the magnitude of all such vectors is the same as ac . This system was solved by an expansion similar to that for the finite-amplitude system presented before. The growth rate to the lowest order in ε of the most critical disturbances is found to be zero. In the order ε of the stability system, we found that the corresponding solvability condition yields 10 = 0. Next, we applied the solvability conditions in the orders εı and ε2 of the stability system to determine 11 and 20 . Since 11 and 20 may not in general be zero for a particular solution, we define ∗ = εı 11 + ε2 20 ,
(15)
as the leading order growth rate and combine the solvability conditions in the orders εı and ε2 to obtain the system for * whose very lengthy expression will not be given here. The growth rate * of the disturbances acting on the simple finite-amplitude steady solutions in the form of rolls, squares, rectangles with different angles and hexagons are then determined. Our stability results verified that those bifurcation branches for which the amplitude decreases with increasing R are unstable and,
21
thus, not physically realizable. So the information about the coefficients for R in leading orders in ε given in (14) was found to be useful in the sense that they could indicate preference and stability of particular solutions and vice versa. For the case where the coefficient R11 for R in the order ε is zero, which was found to correspond to steady non-hexagonal type solutions, then we found that R20 > 0 for rolls, rectangles and squares. Thus, such solutions are supercritical, where R > Rc , and so the amplitude of the convective flow is largest, provided the value of R20 is smallest among all such solutions to the nonlinear problem. In addition, R20 is also positive for hexagons. The value of R20 was found to be smallest in the case of rolls. For the case of hexagons, where R11 is non-zero and given by (12), we found that subcritical down-hexagons, where R11 > 0, can exist over a wide range of values of the viscous dissipation and permeability parameters as the results in Table 1 provides an example. However, up-hexagons with up-flow at the cells’ centers and down-flow at the cells’ boundaries are supercritical over the domain where subcritical down-hexagons exist. We also calculated the values of the solid fraction at the center and at any node of a down-hexagonal cell and an up-hexagonal cell for different values of the parameters. We found, in particular that regardless of the values of the parameters, the total solid fraction at the node of the down-hexagonal cell is smaller than the corresponding basic state value, while the result is reversed in the case of the up-hexagonal cell. These results indicate that there is tendency for chimney formation at the nodes of down-hexagons, which is in agreement with the experimental observation due to Tait et al. (1992), while in the case of up-hexagons there is tendency for solid fraction enhancement at the nodes. The growth rate * given by (15) of the disturbances acting on the finite-amplitude solutions has been computed for rolls, squares, rectangles with different angles and hexagons. In all the cases that we investigated for different range of values of the parameters and the small values of |ε|, we found the following stability results: Supercritical rolls are stable only if √ 3ı|R11 | (16a) , |ε| ≥ ε1 , ε1 = H where the expression for H, which is positive, is very length and will not be given here. Down (up)-hexagons are stable if ε2 ≤ |ε| ≤ ε3 , ε3 =
ε2 =
0.5ı|R11 | ε0 , = R20 2
3ı|R11 | , R11 > (<)0, H
(16b)
where R20 refers to this coefficient for the hexagons and ε0 corresponds to the boundary between the subcritical and supercritical domains of hexagons where R = Rc . As can be seen from the expressions for the positive quantities εi (i = 0, 1, 2, 3) given in (9a) and 9(b), the stable hexagons are subcritical for smallest values of |ε| and them become supercritical for larger values of |ε|, and there is a domain of overlap where both supercritical rolls and supercritical hexagons are stable. As an example, we evaluated the quantities εi (i = 0, 1, 2, 3) for C = 0.09, St = 3.2, Kc = 0, K2 = 0 and ı = 0.001 at the viscous dissipation parameter D = (10−4 ). We found R11 = 0.881, ε0 = 0.000040, ε1 = 0.00017, ε2 = 0.000020 and ε3 = 0.00029. Thus, here the subcritical domain for the stable down-hexagons is 0.000020 ≤ |ε| ≤ 0.000040, the supercritical domain for the stable down-hexagons is 0.000040 < |ε| ≤ 0.00029, and the supercritical domain for the stable rolls is |ε| ≥ 0.00017. There is a region of overlap for 0.00017 ≤ |ε| ≤ 0.00029 where both supercritical downhexagons and supercritical rolls are stable and hysteresis effects can exist.
22
D.N. Riahi / Mechanics Research Communications 39 (2012) 18–22
4. Conclusion We studied nonlinear steady convective flow in a mushy layer and in the presence of viscous dissipation. Using perturbation and stability analyses, we first determined the weakly nonlinear steady solutions in the form of rolls, rectangles, squares, down-hexagons and up-hexagons, and then we determined the stability of such solutions with respect to arbitrary disturbances. We found, in particular, that over a wide range of values of the viscous dissipation parameter, the only stable solution is in the form of subcritical down-hexagons for sufficiently small amplitude of the convection, and there is tendency for chimney formation at the nodes of the stable down-hexagons. These results are in agreement with the experimental observation (Tait et al., 1992). References Amberg, G., Homsy, G.M., 1993. Nonlinear analysis of buoyant convection in binary solidification with application to channel formation. J. Fluid Mech. 252, 79–98.
Anderson, D.M., Worster, G.M., 1995. Weakly nonlinear analysis of convection in mushy layers during the solidification of binary alloys. J. Fluid Mech. 302, 307–331. Chandrasekhar, S., 1961. Hydrodynamic and Hydromagnetic Stability. Clarenton Press, Oxford, UK. Chung, C.A., Chen, F., 2000. Onset of plume convection in mushy layers. J. Fluid Mech. 408, 53–82. Guba, P., Worster, M.G., 2010. Interaction between steady and oscillatory convection in mushy layers. J. Fluid Mech. 645, 411–434. Nield, D.A., 2007. The modeling of viscous dissipation in a saturated porous medium. J. Heat Transfer 129, 1459–1460. Riahi, D.N., 2002. On nonlinear convection in mushy layers. Part 1: oscillatory modes of convection. J. Fluid Mech. 467, 331–359. Riahi, D.N., 2003. Nonlinear steady convection in rotating mushy layers. J. Fluid Mech. 485, 279–306. Roper, S.M., Davis, S.H., Voorhees, P.W., 2008. An analysis of convection in a mushy layer with a deformable permeable interface. J. Fluid Mech. 596, 333–352. Tait, S., Jahrling, K., Jaupart, C., 1992. The planform of compositional convection and chimney formation in a mushy layer. Nature 359, 406–408. Worster, M.G., 1992. Instabilities of the liquid and mushy regions during solidification of alloys. J. Fluid Mech. 237, 649–669.