- Email: [email protected]
Gravitational modulation of thermosolutal convection during directional solidification
CRYSTAL GROWTH
Journal of Crystal Growth 129 (1993) 71)—h))
North-Holland
Gravitational modulation of thermosolutal convection during directional solidification B.T. Murray, S.R. Coriell, GB. McFadden, A.A. Wheeler and B.V. Saunders 3~ationa/Institute of .~taodardsand Tcclinoloç’v. Gaithervhuri.’,Marvla,uI 2O~9V,GSA Received 2 July 1992
1)uring directional solidil cation of a binary allos at constant velocity, thermosolutal convection may occur due to the temperature and solute gradients associated with the solidification process. For vertical growth in an ideal furnace (lacking horizontal gradients) a quiescent state is possible. The effect ol a time—periodic vertical gravitational acceleration (or cquivalcntly vibration) on the onset of thermosolutal convection is calculated based on linear stability using Floquet theory. Numerical calculations for the onset of instability have been carried out for a semiconductor alloy with Schniidt number of I)) and Prandt I number of 11.1 with primary emphasis on large modulation frequencies in a microgravity environment k)r which the background gravitational acceleration is negligible. The numerical results demonstrate that there is a significant difference in stability depending on whether a heavier or lighter solute is rejected. For large modulation frequencies, the stability behavior can be described by either the method of averaging or an asymptotic resonant mode analysis.
1. Introduction In a previous paper, we treated the effect of sinusoidal gravity modulation on the onset of solutal convection during vertical directional solidification of a binary alloy [1]; thermal effects were assumed negligible. Here, we extend this work to treat the onset of convective instability due to both thermal and solutal fields using a numerical implementation of Floquet theory. These calculations are computationally intensive; however, for high frequency modulation, we demonstrate that asymptotic methods can be used to simplify the stability analysis [2,3]. The high frequency limit is relevant for a range of processing conditions and a variety of alloy systems. The effect of gravitational modulation on thermosolutal instabilities in an infinite horizontal layer has been studied recently [4,5]. The review article by Davis [6] contains an extensive bibliography of the early work on periodic modulation for convective flows in general, while the reviews by Ostrach [7], Monti et at. [8], and Alexander [9] U))22-l)245/93/$OoJX) s
1993
—
provide a comprehensive list of references more specific to modulated convection in a low-gravity environment. Additional references pertaining to thermosolutal convection and the effects of residual acceleration in directional solidification are contained in Saunders et al. [4]. In addition to being relevant to materials processing in a microgravity environment, periodic modulation is being used as a means of dynamically controlling convection effects for earth based processing [10]. For growth vertically upwards with rejection of the more dense component of the alloy, both the thermal and solutal fields cause the liquid density to decrease with height, i.e., the thermal and solutal fields act together. In contrast, when a lighter solute is rejected the solute field alone would cause the liquid density to increase with height, i.e., the thermal and solutal fields oppose each other. The focus in the present article is on the effect of periodic gravitational acceleration in the absence of steady background acceleration: this is essentially the situation under microgravity conditions. The theory and approximate analysis
Elsevier Science Publishers By. All rights reserved
B. T Murray et a!.
/ Grai’itational modulation of thermosolutal coniection during DS
are also applicable in the presence of a steady background acceleration. When the thermal and solutal fields act together, during half the modulation cycle the liquid density is stably stratified, while for the other half it is unstably stratified. When the two fields oppose one another, the situation is more complicated since due to double-diffusive effects [11] the density stratification is not a meaningful criterion for convective instability. The directional solidification problem is further complicated because the solutal field typically acts over a much smaller length scale than the thermal field. For gravity modulation of a single component fluid there is an important distinction between modes of instability depending on whether the system is stable or unstable in the absence of modulation. If the unmodulated system is unstable there is what is referred to as a fundamental mode of instability; that is, the instability that occurs in the absence of modulation persists if the modulation amplitude is small enough. It is found that the fundamental mode becomes stable for large enough modulation amplitudes. The amount of modulation necessary to stabilize the mode depends on the modulation frequency; for large frequencies, the modulation is unable to stabilize the fundamental mode, and for any fixed modulation amplitude the system is unstable for large enough frequencies. If the unmodulated system is stable, then modulation of sufficient amplitude may cause resonant modes of instability; these resonant modes also occur when the unmodulated system is unstable. The fundamental mode has the same frequency as the modulation frequency (synchronous response), while the resonant modes are either synchronous or subharmonic (response with half the modulation frequency) depending on the modulation frequency. For high frequencies these resonant modes show subharmonic response [1,31. For the two-component system, the situation is more complicated, as was shown for the infinite layer with constant gradients and stress-free boundaries [4]. In addition to synchronous and subharmonic modes of instability, complex conjugate modes can occur with temporal frequency unrelated to the modulation frequency. The fun-
71
damental mode shows either synchronous or complex conjugate response depending on the temporal character of the unmodulated system: the response is synchronous in the fingering regime of the unmodulated system, and there is a complex conjugate response in the diffusive regime of the unmodulated system. Surprisingly, in the two-component case the fundamental mode of instability persists at high enough frequencies and modulation amplitudes, even under conditions for which the unmodulated system is stable. In this case, it has the form of a lobe, since the modulated system must regain stability if the modulation amplitude is small enough. The most dangerous mode at high frequencies can be either the fundamental mode or a resonant mode. In the next section the model used to study the onset of thermosolutal convection in vertical directional solidification is presented. Also, we briefly describe the numerical solution procedures used to solve the resulting Floquet equations when sinusoidal gravitational acceleration is imposed at a given modulation frequency. In the limit of large modulation frequency, we describe two types of approximate analyses which significantly simplify the system of equations governing the stability behavior. We carry out calculations for a model semiconductor alloy with Schmidt number of 10, Prandtl number of 0.1 and distribution coefficient of 0.3, for large modulation frequencies. We calculate the modulation amplitude required for instability for both a heavier and lighter solute. We find that the alloy which rejects a lighter solute is substantially more unstable, with instability occurring by a fundamental mode with synchronous response. For the alloy rejecting a heavier solute, instability occurs by a resonant mode with subharmonic response.
2. Governing equations We consider growth at constant velocity, V, of a binary alloy with the growth direction and the gravitational field in the z-direction. The heat flow, fluid flow, and solute diffusion equations are written in terms of a Cartesian coordinate
B. T. [Vito-i-arci al.
72
/
(7rai ‘itationaf modulation oJ iliermosolu to! ion! cetion during D.S
system (x, y, z) which is attached to the planar crystal—melt interface; the fluid velocity is measured in the laboratory frame. The system is laterally infinite, and in the base state in the moving frame, the temperature and concentration fields are functions of z alone and there is no flow (the density change on solidification is neglected). We use dimensionless variables and measure length, time, fluid velocity, concentra-
tion, and temperature in units of D/V, D/V, V. c~,and G1D/V. respectively. Here. D is the solute diffusion coefficient, c~is the bulk solute concentration, and G1 is the temperature gradient in the liquid at the interface. We consider small perturbations from the base state, and Fourier analyze in the lateral directions, that the x and y dependence of the perturbed quantities is of the form exp(ia~x+ ia .y), where a and ai- are wavenumbers in the x and y directions. Since the x and y directions are equivalent, only the equations. quantity a Theya~ + a~ entersequainto the resulting perturbation ~.
=
lions for the z-component of the fluid velocity, w(z. t), the concentration field, c(z, t), and the temperature fields, T(z, t) and T(z. t), in the fluid and crystal are
~(
Lw
Sc
—
)
Lw
L2w
=
—.
ratio of the crystal and liquid thermal conductivities. The unperturbed gradients are given h~’ 0) C~
=
kk
—
‘T.~° exp[ =
I
1)/k exp[—:], (Pr/Sc) z~j
-
where k is the segregation coefficient. The thermal and solutal Rayleigh numbers are defined as Ra(t) =g(t)aG1(D/V)
a /oK,
Rs(t) =g(t)f3c~(D/V)~/vi), where g(t) is the gravitational acceleration, and and f3 are the thermal and solutal expansion coefficients, respectively. We assume a gravitational acceleration of the form, g(l) g1 + g cos(12t), where g11 ts the background acceleration, g1 is the modulation amplitude, and 12 is the dimensionless modulation frequency, and we write t0t + Ra~cos( lit). Ra(t) Ra I Rs( t ) Rs0 + Rst 1 cos( lit). a
=
= =
(Ia)
Since we have not allowed the interface to deform, it is necessary to relinquish one of the usual solidification boundary conditions that hold at a crystal—melt interface [12]. The dimensionless boundary conditions that we choose to impose at the planar interface. z 0, are w 0, (2a)
(lb)
c
a[Rs(t)c
/
=
+
(Sc/Pr) Ra(
t
) TJ.
=
=
C 1
—
C~ +
wc~ Lc. =
Sc
-
T2
=
(Id)
Here, subscripts indicate and 2/az partial —a2). derivatives The Schmidt the operator L ratio (d of the kinematic viscosity, number, Sc, is the t’, to the solute diffusivity; the Prandtl number, Pr, is the ratio of the kinematic viscosity to the liquid thermal diffusivity, i<; and q is the ratio of the crystal and liquid thermal diffusivities. Since we assume that the heat capacities of the crystal and liquid are identical, q also represents the =
(I
k )c
=
0.
(2h)
1= T, (lc)
q Sc —~—LT.
+
(2c)
~ (2d) In the above formulation, we satistv the conservation relations at the interface, but neglect the relation between interface temperature and concentration. For should the convective instabilities considered here this he a good approximation. The perturbations are required to vanish at the upper boundary, z h~.The temperature perturhations in the solid vanish at Solutions to the above set of equations which have time-periodic coefficients can he obtained using the framework of Floquet theory [13], in which the solutions are represented as the prod=
=
B. T Murray et a!.
/ Gracitational modulation of thermosolutal connection during DS
uct of a temporally periodic function and an exponential function of time. Two different flumerical implementations of Floquet theory were employed to solve the stability problem formulated above and are discussed briefly below. The first approach consists of representing the periodic component of the solutions by a truncated complex Fourier series in time. The solutions for the perturbation quantities are represented by the product of a periodic Fourier series and an exponential term with complex growth rate ~ f(z,
t)
=
ei~ ~
f,n(Z)
eimth,
where f(z, t) represents any of the perturbation quantities in the liquid. Since the temperature in the crystal is only coupled to the temperature in the liquid through the boundary conditions, it is possible to solve for the Fourier components of the crystal temperature analytically. Substitution of the series expansion into the equations in the liquid and boundary conditions yields a set of 32M + 16 coupled two-point boundary value problems in the spatial variable z for the cornplex Fourier coefficients of the periodic solution components. The resulting set of coupled ordinary differential equations subject to the homogeneous boundary conditions yields an eigenvalue problem that is solved as previously described [1] The set of equations contains the complex growth rate if as a parameter. The base state subject to periodic forcing is linearly stable for a given set of parameters if, forrate all udisturbances, 0r of the growth is negative. the the realcalculations, part In setting 0’r 0 allows for the determination of marginal values of the modulation amplitude, ~ and the imaginary part of the growth rate, u 1, for fixed values of the remaining parameters. For many cases, the imaginary part of ci is equal to zero (synchronous response) or 12/2 (subharmonic response). Restricting ci to these values, the number of unknowns is reduced to 16M + 8. The second approach employed consists of approximating the spatial behavior of the disturbance quantities via the pseudospectral technique =
73
in the physical domain as previously described [I]. The approach corresponds to expanding the solutions in terms of a truncated series of Chebyshev polynomials I~(s),
f(z
N
t)
=
~ f~(t) 7~(s), =
where s
=
(2z/h1,)
1 in the liquid region. The temperature in the crystal is calculated numerically in the same manner in this approach. The pseudospectral discretization requires that the expansions satisfy the governing equations at specific collocation points for the Chebyshev polynomials. When implemented in the physical domain the unknowns are the solution values at the collocation points. The spatial differential operators in the governing partial differential equations are replaced by discrete matrix operators. As a result, the governing set of partial differential equations and boundary conditions becomes a set of coupled ordinary differential-algebraic equations in time for the unknown solution values at the collocation points. In this second solution approach, Floquet theory is implemented by constructing a fundamental solution matrix. The columns of this matrix are linearly independent calculated solutions for the unknowns at the end of one forcing period. The eigenvalues of this matrix are the Floquet multipliers from which the complex growth rate if is obtained. The utility of the second approach is that o- itself is the eigenvalue, so that its magnitude is obtained for fixed values of the remaining parameters, which simplifies the search for marginal and gives moremodes; information regarding thevalues proximity of other in particular, this approach allows the calculation of complex conjugate modes which may occur in the thermosolutal problem [4]. —
3. High frequency analysis In this section we discuss two methods of approximation which are valid in the limit of large modulation frequency. For many alloys and processing conditions, this is a relevant limit. The basic methods of analysis have been previously
74
B. 1. Murray ci 0/. / (,ra, Oaiiona/ inodu/anoti of ilierinosolutal tout 1’ctioui iIurin~’1)5
described for solutal convection [3] and the extension for thermosolutal convection is straightforward. Hence, we simply give the resulting cxtensions of these approximations to thermosolutal convection. In the method of averaging, the modulated Rayleigh numbers, Ra1 and Rst I) are proportional to the modulation frequency Ii in the limit of large frequency. The solution variables are decomposed into an averaged part and rapidly varying part, where the integral of the rapidly varying components over one modulation period (synchronous behavior) vanishes. Using this decomposition and averaging the governing equations produces a set of equations for the averaged quantities with coefficients that no longer depend on time. The stability of these equations can now he treated by a standard normal mode analysis, where the time-dependence is of the form exp(ut). The averaged equations for thermosolutal convection are -
(ciL~ L~ ) Sc
=
a2 Rs0~+
—
L2~
:
R~11T
.
~ ~
2
—~
Pr
~(
2
.
+ -
T~).
~-~-
Pr
‘
-
3a) (IC
~: +
—
ciT
—
T
=
+
(3h) (3c)
=
-
-
1 =
L~, Sc —LT. Pr
.t~/ +
~2
Sc
Sc
—,
iT
.
(3d)
L
+
‘:
(1
k)~=0.
—
~1’~
(41-i)
0.
=
(4e)
where A
Pr --
=
~
Pr
+
~—
.
+
~——---
Sc
a2
+ r
~ Sc
At the tipper boundary. conditions are
=
Ii
.
the houndar~
..
~
=
=
-
(4d)
=
By employing the averaging technique, the stahility problem for the effect of sinusoidal niorlulation in the limit of high frequency is reduced from the system (1) with time-periodic coefficients to (3) which depend on alone. In the averaged analysis. the order of the system of ordinary differential equations is increased by two over the usual set of unmodulated equations: the resulting computational times arc greatly reduced over the full set of Floquet equations. Another type of asymptotic analysis was used to describe resonant mode behavior in the limit of large frequency for the single component case [3]. We now extend this analysis to the doublediffusive directional solidification case. We appl~ this analysis in the absence of a background grayitational acceleration, and set Rs~ atid Ra0 equal to zero. For the analysis, we assume the modulated Rayleigh numbers scale as 112 and define 1/
=
Rs1t~/1l2. II
=
Ra~t/l22.
Pr
-
We introduce the new variables lit. I and ~ lIC. In the limit of large 12. (1) reduces to =
where the overhar denotes averaged quantities and the unknown function ~ arises from the analysis and represents the contribution of the rapidly varying components in the averaged equations. The scaled modulation amplitudes are given by .1. Ra1 >/li and Rs1 1/11, The averaged temperature in the crystal can he determined analytically and the resulting boundary conditions on the averaged quantities in the liquid at the interface are =
— — —
~:
=
— —
— —
0
( a)
=
I2T.
=
I ~~1.tt; C~ +
~ +
=
—tu cos =
~
=
=
Sc!! Pr T+ II~
(1. 0
.
(~~i)
(71-i) (Se) -
-
In this approximate analysis, the spatial order of the governing equations has been reduced trom eight to two, only conditions are and that the ii 0 at remaining 0 and Iiboundary =
=
B. T. Murray et a!.
/
Gravitational modulation of thermoso!utal convection during DS
The set of partial differential equations (5) is separable, and we introduce the functions 0(T) and Z(z) such that w(z, r) ê(z, T) T(z, r)
= =
OT(r)Z(z), —0(r)Z(z)c~°~, 0(r)Z(z)T~ ~
=
(~.ccos T)O 2~
Z~.
--
a2
—
ji
=
Ratt0 0. Thus, the system is stable in the absence of modulation and we investigate the modulation amplitude at which instability occurs. For typical growth velocities, modulation frequencies that are greater than 1 Hz correspond to very large dimensionless frequencies. For example, for a growth velocity, V= iO~cm/s, and liquid diffusion coefficient, D l0~cm2/s, a dimension=
=
which results in the following two ordinary differenttal equations 0~+
75
0,
(6a)
s
—TI. T~°~ + Hc~°~ z Pr Z
=
o. (6b)
The equation for 0 is the Mathieu equation which describes the temporal behavior. As the magnitude of j.c increases, the first transition to instability occurs as a subharmonic mode at ~ 0.454. The spatial dependence is described by the function Z which requires solution of a two-point boundary value problem since Z(0) Z(h~) 0. This is an eigenvalue problem and nontrivial solutions are obtained only for specific values of the parameters. Here we choose H~as the eigenvalue parameter and numerically calculate the lowest value for fixed values of the remaining parameters. In the case where Prh,/Sc ~z 1, the unperturbed temperature gradient, T..1tt~, is approximately constant and the eigenvalue problem can be reduced to an algebraic equation involving Bessel functions, The differential eigenvalue problems that arise in the two approximate analyses described above are solved numerically using the approach previously described [12]. =
=
=
4. Results and discussion We carry out a series of calculations for a typical semiconductor alloy with Sc 10, Pr 0.1, and k 0.3. We consider conditions appropriate to processing in a microgravity environment and assume that the steady background gravitational acceleration can be neglected, i.e., we set Rst0~ =
=
=
=
less frequency, 12 i0~is equal to a dimensional frequency of 1.6 Hz. For most of the calculations we set the dimensionless height, hf= 10, and the dimensionless wavenumber, a 0.5. In principle, we should vary the wavenumber to determine the minimum modulation amplitude, but in order to reduce the required computations we only do this for selected cases. Unless otherwise stated, the values of the dimensionless parameters given here will be used in all the calculations. For growth vertically upwards with rejection of a lighter solute (temperature field stabilizing; solute field destabilizing), both the thermal and solutal Rayleigh numbers are defined to be positive. We first consider an alloy rejecting a heavier solute for which during half of the modulation cycle both the temperature and concentration field are stabilizing, while during the other half of the cycle both fields are destabilizing; thus Rs°~ and Ra(i) have different signs. In fig. 1, for Li 1000, we plot the modulated solutal Rayleigh number as a function of the modulated thermal Rayleigh number calculated from the full set of equations and boundary conditions, eqs. (1) and (2). The system is stable below the curve and unstable above. If there were no interaction between the thermal and solutal fields, the stability boundary would be comprised of a horizontal line and a vertical line. The shape of the stability boundary in fig. 1 indicates only a weak interaction between the thermal and solutal fields; the system is unstable for Rs~° greater than 3>< i0~ or Ra~°greater than 700. For this case, the instability is a subharmonic resonant mode which oscillates with half the modulation frequency. We will demonstrate shortly that at this frequency the instability is well characterized by the asymptotic resonant mode analysis. We next consider an alloy for which the temperature field is stabilizing when the solute field =
.
=
=
70
B. 71 Murray ci a!.
/
(iratitational moi/ulatii,n of t/icrmosolutal
Thermosolutal Convection with Vertical g—Jitt.er Sc = 10 k =1°>=0.0 0.3 hi Pr == 10 0.1 Rs(°>=0.0 Ra — . —~
I
—
=
~
coot etilon iluriuig
1)5
tion between the two fields is necessary for occurrence the synchronous synchronous mode mode. can As we show below,of the be scribed by the method of averaging at this
quency. As discussed previously. the method of averaging and the resonant mode analysis assume that the modulation amplitude scales as Ii and 122. respectively, for sufficiently large Ii. In fig.3. we plot Rst I asafunction of Ii for Rst it I OORat I) and the other parameters the same as in figs. I and 2. The points ieprcsent the results Ironi the numerical calculations, while the solid line corre— sponds to Rst I) b12 and the dashed line to Rst 1 1)7.122 where the constants I~and b-. were calculated from the method of averaging and resonant mode analysis, respectively. It is clear that the large frequency asymptotic behavior holds with only slight deviations at the lowest frequency 12 102. In view of this, in fig. 4 we show the full numerical and the numerical asymptotic results for Rs1 it as a function of Rs1 i 1/ Ra1 for 12 =
Stable
the will defre-
±
=
=
111111
I t~
~
to
Fig. I. The modulated solutal Rayleigh number as a function of the modulated thermal Rayleigh number at the onset of suhharmonic instability in the absence of background modulation. The temperature and solutal fields act together,
is destabilizing over half the cycle and vice versa during the other half of the cycle; thus Rsw and Ra~’~ have the same sign. Fig. 2 is a corresponding plot to fig. 1. except that RaW is now positive. For sufficiently small Rst or for sufficiently small Ra~1.subharmonic resonant modes are obtained which are in agreement with the results shown in fig. I. However, away from these limiting values, the instability occurs with synchronous response at significantly lower modulation amplitude than in fig. I. For this synchronous mode, the system is unstable for a range of Ra(i) for Rsw greater than approximately l0~and for a range of RsW if Ra~°is greater than approximately 20. Comparing these results with those in fig. 1, it is interesting to note that instability occurs ft~rmodulation amplitudes two orders of magnitude smaller when the temperature and solute fields oppose each other. This can be attributed to a double-diffusive coupling, since each component taken individually would lead to suhharmonic resonant instahility. While the shape of the stability boundary for the synchronous mode exhibits a weak interaction between thermal and solutal fields, the interac-.
=
.
-
=
Thermosolutal Convection with Vertical g—Jitter Sc = 10 k = 0.3 hi = 10 Rs~°1= 0.0 Rat0~=0.0 Pr = 0.1 0 = 1000 a = 0.5
-
/ /
—
2
~
‘~-~
/
/ ~-
--
I
I flTflf
—
Stable I
rrrrrrr
Ill
Il/UI
liii
I rTI
Ra(i) Fig. 2. The modulated solutal Rayleigh number as a function of the modulated thermal Rayleigh number at the oiiset 0! instability in the absence of background modulation. ihe ..
.
-
temperature and solufal fields oppose each other. Fhe solid line indicates synchronous response and the (lashed line subharnionie response.
B. T Murray ci a!.
/ Gravitational modulation of’thermo,so!uta! connection during DS
l0~.The results from the full equations are mdicated by the data points, while the results for the synchronous modes from the method of averaging are indicated by the solid line and the results for the subharmonic mode from the resonant mode analysis are given by the dashed line. The asymptotic results are in excellent agreement with the full numerical results at this frequency. Using the averaged equations, we have also minimized RsW with respect to wavenumber, a, and the results are shown by the lower dashed curve. This curve deviates only slightly from the a 0.5 curve with the largest difference at lower values of t° Rs°~/Ra~° Wefor notea that quantity is a constant giventhealloy and Rs°~/Ra processing conditions and independent of the amplitude of the gravitational modulation, g The results presented in figs. 1—4 were primarily for a fixed spatial wavenumber, a In figs. 5 and 6, we plot the modulated solutal Rayleigh number as a function of wavenumber for Rst°/Ra°~ ±100. In fig. 5 the thermal and solutal fields act together (Rs~i)/Ra(t) 100) and the resonant instability is subharmonic. The
77
Thermosolutal Convection with Vertical g—Jitter Sc = 10 k = 0.3 h 1 = 10 Rs(°)=0.0 Ra(°)=0.0 Pr = 0.1 0 1000 a 0.5
—
=
=
0.5.
=
—
=
Thermosolutal Convection with Vertical g—Jitter Sc = 10 k = 0.3 h 1 = 10 Rs~°~= 0.0 Ra(°> 0/~Ra(’)~ 0.0 Pr = 100 = 0.1 a = 0.5 Rs(
‘~ ‘
111111/
111111111
id~
FF1111111
11111115
111111111
II
id
if~ 1O~ IO~ Rst’t/IRa0~ Fig. 4. The modulated solutal Rayleigh number as a function of the ratio of the modulated solutal and thermal Rayleigh numbers at the onset of instability. The points represent the full numerical calculations with circles corresponding to Rst i)/Ra(i) positive, triangles to Rs~° /Ra’° negative: solid circles correspond to synchronous modes and open symbols to suhharmonic modes. The solid and dashed curves represent the method of averaging and the resonant mode asymptotic analyses, respectively, for a = 0.5; the lower dashed curve represents the minimum value with respect to wavenumber from the method of averaging.
minimum Rs~1occurs for a wavenumber a 1.5; the value of RsW for a 0.5 is 1.4 times this minimum value. The structure of the curve for larger wavenumbers suggests that a second minimum might occur for somewhat different values of Rs~/Ra~.The dashed curve is calculated from the asymptotic resonant mode analysis. The approximate analysis yields good agreement at small wavenumbers, hut underestimates the prediction from the full numerical calculation for larger wavenumbers. A similar behavior was found for the onset of solutal convection in ref. [31, where the limitations of the resonant mode analysis have been discussed. In fig. 6 the thermal and solutal fields oppose each other (Rs°~/Ra°~ 100) and the instability is synchronous. The mini=
111111111
11111111
10 0 Fig. 3. The modulated solutal Rayleigh number as a function of modulation frequency at the onset of instability. The solid line (synchronous response) has slope unity and the dashed line (suhharmonic response) has slope two,
=
mum Rs~~ occurs for a wavenumber a 0.5 in this case. The dashed curve is calculated from the
70
B. ‘/1 Murray en a!.
/
(inn ‘itationaf modulation
Thermosolutal Convection with Vertical g—Jitter Sc10 k0.3 h =10 10~ 0.0 Ra~°~ 0.0 = 0.1 Rs 0 = 1000 Rs(’)/EaO) = —100
o °
0~
of t/ienInovo/uia/
co,i/i’itntln
t/uu’in,ç’ 1)5
niethod of averaging and is in ti)od agreement with the full numerical calculations with increas.
.
.
ing deviation as the wavenumber increases. How-
ever, in this case, the approximate analysis accurately predicts the minimum value. In space experiments the background gravitational acceleration is very small but not zero. We have carried out a few calculations with a nonzero background acceleration. For Rsttt/Ra°~ 101) and Rst°t 1, 0. 1, the values of Rst I) :tt the onset of instahility are 1735, 1605, 1464, respeclively. As Rst°~increases from zero, the stability of the unmodulated state decreases and thus the modulation amplitude required t’or instability decreases. For Rs~t/Rat I) 100 and Rst 0. 1. there is less than 0.01 ~ change in the modulation amplitude: resonant modes of insta=
~
=
0
H
—
=
0
l’ll’l’
0.0
2.0
4.0
8.0
8.0
10.0
a Fig, 5. The nnodulated solutal Rayleigh number asa function of the wavenuniher at the Onset of suhharmonic instability. The temperature and solutal fields act together. The unstable region is shaded; the dashed curve is based on the asymptotic
—
~:
—
. .
hility are insensitive to the background acceleration level. Based on this limited number of calculations, the case of small background acceleration is reasonably well approximated by setting the background acceleration to zero.
analysis tor resonais t mode behavior.
In our previous paper on solutal convection [I] we neglected thermal effects. For the properties and processing conditions considered here, this is Thermosolutal Convection with Vertical g—Jitter Sc = 10 k = 0.3 hi = 10 Rs~°~ 0.0 Rat0>= 0.0 Pr = 0.1 (4 = 1000 Rs(’>/Ra0) = 100
=
a good approximation when Rat is less than about 600 (see fig. I) for the case when the thermal and solutal fields act together. hut Rat It must he less than approximately It) when the fields oppose each other (see fig. 2), In the ab-
~:~I~1:;l:t5.the resonant mode analwhere ~i ((.454 is the first transition point of the Mathieu equation and r~ is the first zero of the Bessel function. For a 0.5. Sc 10. A ((.3. =
‘~1i’
=
=
=
and 12 1000, the above tormula gives RsC 2.14 X l0~.which is in reasonable agreement with the calculations (see figs. I and 2). From eq. (6) with =
‘‘‘‘ia’’’ ~ a Fig. f. The modulated solutal Rayleigh numter as a function °
~
a’
‘
‘
~.s
of the wavenumber at the onset of synchronous instability. The temperature and solutal fields oppose each other. The unstable region is shaded; the dashed curve is based on the asymptotic analysis using the method of averaging.
t =
the assumption that the temperature gradient is constant, it is possible to obtain a similar expression for the purely thermal problem. namely R’ Ci
i
— —
/5
Pr 12
Sc2
1
+
~-
a/i
-
B. T. Murray et a!.
/ (irai’itational modulation of thermoso!utal convection
during DS
79
For a 0.5, Sc 10, Pr 0.1, h 10, and LI 1001), the formula yields Rate 630, which again agrees with the results of figs. 1 and 2. The appearance of Sc in the above formula for the pure thermal case is a result of our choice of a length scale based on D/V. We recall that the resonant mode analysis is a good approximation for small values of wavenumber, a, but it significantly underestimates the minimum value of the modulation amplitude with respect to wavenumher. All the calculations discussed above were for the value h 1 10. The above formula for the purely thermal problem indicates that the dependence on h( for the resonant modes is fairly weak. A few numerical calculations using the full equations have also shown a weak dependence of the modulation amplitude on h~.We have investigated the dependence of the fundamental mode on h1. using both the full equations and the method of averaging. The dependence is very weak; for example, for R5W/RaW 100 the modulation amplitude for instability calculated from the full equations decreases by 0.4% as h, increases from 10 to 40. For the averaged equations, no further change occurs in the range investigated (Ii, 100).
oppose each other a fundamental, synchronous mode with Rs°~ a 12 occurs as well. The problem of thermosolutal convection in directional solidification with gravity modulation depends on a large number of parameters. We have only investigated a small subset of these parameters here. We have shown that for high frequency modulation the method of averaging provides an accurate description of the fundamental mode which can be the most dangerous mode when the temperature and solutal fields oppose each other. The asymptotic analysis of
5. (.onclusions
similar the fundamental modes manner describedto here in that they synchronous can be the most dangerous modes at high modulation frequency. In addition to the fundamental oscillatory mode, resonant modes occur and can exhibit strong resonance with the modulation frequency.
=
=
=
=
=
=
=
=
In summary, we have studied the stability of thermosolutal convection in a model of vertical directional solidification subject to purely sinusoidal gravitational modulation using a numerical implementation of Floquet theory and two types of approximate analyses for large modulation frequencies. The calculations show that the presence of both the temperature and solute fields can yield a significant coupling which causes instability at lower modulation amplitude than is required for the fields taken individually. This coupling is more significant when the two fields oppose each other than when they act together. The calculations show that when the fields act together, a subharmonic, resonant mode of instability occurs with Rs(i) a ~2 When the fields
resonant modes yields an accurate prediction of the modulation amplitude for instability at low wavenumbers. These approximate analyses significantly reduce the computational effort and should allow for investigation of a wider range of parameters. Although we have concentrated on the situation with zero background gravitational acceleration, the analysis is not limited to this case. For constant gravitational acceleration, the onset of thermosolutal convection can occur by an oscillatory mode. The computational methods presented here can be used to study the effect of sinusoidal modulation on these fundamental oscillatory modes. For stress-free boundaries and constant temperature and solute gradients [4], these fundamental oscillatory modes behave in a
Acknowledgement This work was conducted with the support of the Microgravity Science and Applications Division of the National Aeronautics and Space Administration. R ~ e erences [11B.T. Murray, SR.
Coriell and GB. McFadden, J. Crystal Growth 110 (1991) 713.
B. T Murray
Of)
[21G.Z. 131 [4] [5] [6]
17] 181
Ct
al.
/ Grai’itational
modulation of thcrmo,so!utal coni’ccnon during DS
Gershuni and F.M. Zhukhovitskii. Convective Stahility of lneompressihle Fluids (Keter, Jerusalem. 1976) ch. 8 A.A. Wheeler, GB. McFadden, B.T. Murray and SR. Coriell, Phys. Fluids A 3 (1991) 2847. By. Saunders, B.T. Murray. GB. McFadden. SR. Coriell and A.A. Wheeler, Phys. Fluids A 4(1992)1176. G. Terrones and CF. Chen, J. Fluid Mech., in press. S.H. Davis. Ann. Rev. Fluid Mech. 8 (1976) 57. 5. Ostrach, Ann. Rev. Fluid Mech. 14(1982) 313. R. Monti. J.-J. Favier and D. Langhein. in: Fluid Sciences and Material Science in Space, A European Perspective, Ed. H.U. Walter (Springer. Berlin 1987) p. 637.
191 J.1.D. Alexander. Microgravity Sci. Technol. .7! 109)1) 52. [101 R. Caram, M. Banan and W.R. Wilco~<,J. Crystal Growth 114 (1991) 249. Ill] iS. Turner. Buoyancy Effects in Fluids (Cambridge tiniversity Press, Cambridge, 1973). [12] SR. Coriell, MR. Cordes, Wi. Boeftinger and R.F Sekerka, J. Crystal Growth 49 (198ff) 13. [13] V.A. Yakubovich and V.M. Starzhinskii, Linear [)iff’erenhal Equations with Periodic Coefficients (Wiles’, Ne~ York. 1975).
Journal of Crystal Growth 129 (1993) 71)—h))
North-Holland
Gravitational modulation of thermosolutal convection during directional solidification B.T. Murray, S.R. Coriell, GB. McFadden, A.A. Wheeler and B.V. Saunders 3~ationa/Institute of .~taodardsand Tcclinoloç’v. Gaithervhuri.’,Marvla,uI 2O~9V,GSA Received 2 July 1992
1)uring directional solidil cation of a binary allos at constant velocity, thermosolutal convection may occur due to the temperature and solute gradients associated with the solidification process. For vertical growth in an ideal furnace (lacking horizontal gradients) a quiescent state is possible. The effect ol a time—periodic vertical gravitational acceleration (or cquivalcntly vibration) on the onset of thermosolutal convection is calculated based on linear stability using Floquet theory. Numerical calculations for the onset of instability have been carried out for a semiconductor alloy with Schniidt number of I)) and Prandt I number of 11.1 with primary emphasis on large modulation frequencies in a microgravity environment k)r which the background gravitational acceleration is negligible. The numerical results demonstrate that there is a significant difference in stability depending on whether a heavier or lighter solute is rejected. For large modulation frequencies, the stability behavior can be described by either the method of averaging or an asymptotic resonant mode analysis.
1. Introduction In a previous paper, we treated the effect of sinusoidal gravity modulation on the onset of solutal convection during vertical directional solidification of a binary alloy [1]; thermal effects were assumed negligible. Here, we extend this work to treat the onset of convective instability due to both thermal and solutal fields using a numerical implementation of Floquet theory. These calculations are computationally intensive; however, for high frequency modulation, we demonstrate that asymptotic methods can be used to simplify the stability analysis [2,3]. The high frequency limit is relevant for a range of processing conditions and a variety of alloy systems. The effect of gravitational modulation on thermosolutal instabilities in an infinite horizontal layer has been studied recently [4,5]. The review article by Davis [6] contains an extensive bibliography of the early work on periodic modulation for convective flows in general, while the reviews by Ostrach [7], Monti et at. [8], and Alexander [9] U))22-l)245/93/$OoJX) s
1993
—
provide a comprehensive list of references more specific to modulated convection in a low-gravity environment. Additional references pertaining to thermosolutal convection and the effects of residual acceleration in directional solidification are contained in Saunders et al. [4]. In addition to being relevant to materials processing in a microgravity environment, periodic modulation is being used as a means of dynamically controlling convection effects for earth based processing [10]. For growth vertically upwards with rejection of the more dense component of the alloy, both the thermal and solutal fields cause the liquid density to decrease with height, i.e., the thermal and solutal fields act together. In contrast, when a lighter solute is rejected the solute field alone would cause the liquid density to increase with height, i.e., the thermal and solutal fields oppose each other. The focus in the present article is on the effect of periodic gravitational acceleration in the absence of steady background acceleration: this is essentially the situation under microgravity conditions. The theory and approximate analysis
Elsevier Science Publishers By. All rights reserved
B. T Murray et a!.
/ Grai’itational modulation of thermosolutal coniection during DS
are also applicable in the presence of a steady background acceleration. When the thermal and solutal fields act together, during half the modulation cycle the liquid density is stably stratified, while for the other half it is unstably stratified. When the two fields oppose one another, the situation is more complicated since due to double-diffusive effects [11] the density stratification is not a meaningful criterion for convective instability. The directional solidification problem is further complicated because the solutal field typically acts over a much smaller length scale than the thermal field. For gravity modulation of a single component fluid there is an important distinction between modes of instability depending on whether the system is stable or unstable in the absence of modulation. If the unmodulated system is unstable there is what is referred to as a fundamental mode of instability; that is, the instability that occurs in the absence of modulation persists if the modulation amplitude is small enough. It is found that the fundamental mode becomes stable for large enough modulation amplitudes. The amount of modulation necessary to stabilize the mode depends on the modulation frequency; for large frequencies, the modulation is unable to stabilize the fundamental mode, and for any fixed modulation amplitude the system is unstable for large enough frequencies. If the unmodulated system is stable, then modulation of sufficient amplitude may cause resonant modes of instability; these resonant modes also occur when the unmodulated system is unstable. The fundamental mode has the same frequency as the modulation frequency (synchronous response), while the resonant modes are either synchronous or subharmonic (response with half the modulation frequency) depending on the modulation frequency. For high frequencies these resonant modes show subharmonic response [1,31. For the two-component system, the situation is more complicated, as was shown for the infinite layer with constant gradients and stress-free boundaries [4]. In addition to synchronous and subharmonic modes of instability, complex conjugate modes can occur with temporal frequency unrelated to the modulation frequency. The fun-
71
damental mode shows either synchronous or complex conjugate response depending on the temporal character of the unmodulated system: the response is synchronous in the fingering regime of the unmodulated system, and there is a complex conjugate response in the diffusive regime of the unmodulated system. Surprisingly, in the two-component case the fundamental mode of instability persists at high enough frequencies and modulation amplitudes, even under conditions for which the unmodulated system is stable. In this case, it has the form of a lobe, since the modulated system must regain stability if the modulation amplitude is small enough. The most dangerous mode at high frequencies can be either the fundamental mode or a resonant mode. In the next section the model used to study the onset of thermosolutal convection in vertical directional solidification is presented. Also, we briefly describe the numerical solution procedures used to solve the resulting Floquet equations when sinusoidal gravitational acceleration is imposed at a given modulation frequency. In the limit of large modulation frequency, we describe two types of approximate analyses which significantly simplify the system of equations governing the stability behavior. We carry out calculations for a model semiconductor alloy with Schmidt number of 10, Prandtl number of 0.1 and distribution coefficient of 0.3, for large modulation frequencies. We calculate the modulation amplitude required for instability for both a heavier and lighter solute. We find that the alloy which rejects a lighter solute is substantially more unstable, with instability occurring by a fundamental mode with synchronous response. For the alloy rejecting a heavier solute, instability occurs by a resonant mode with subharmonic response.
2. Governing equations We consider growth at constant velocity, V, of a binary alloy with the growth direction and the gravitational field in the z-direction. The heat flow, fluid flow, and solute diffusion equations are written in terms of a Cartesian coordinate
B. T. [Vito-i-arci al.
72
/
(7rai ‘itationaf modulation oJ iliermosolu to! ion! cetion during D.S
system (x, y, z) which is attached to the planar crystal—melt interface; the fluid velocity is measured in the laboratory frame. The system is laterally infinite, and in the base state in the moving frame, the temperature and concentration fields are functions of z alone and there is no flow (the density change on solidification is neglected). We use dimensionless variables and measure length, time, fluid velocity, concentra-
tion, and temperature in units of D/V, D/V, V. c~,and G1D/V. respectively. Here. D is the solute diffusion coefficient, c~is the bulk solute concentration, and G1 is the temperature gradient in the liquid at the interface. We consider small perturbations from the base state, and Fourier analyze in the lateral directions, that the x and y dependence of the perturbed quantities is of the form exp(ia~x+ ia .y), where a and ai- are wavenumbers in the x and y directions. Since the x and y directions are equivalent, only the equations. quantity a Theya~ + a~ entersequainto the resulting perturbation ~.
=
lions for the z-component of the fluid velocity, w(z. t), the concentration field, c(z, t), and the temperature fields, T(z, t) and T(z. t), in the fluid and crystal are
~(
Lw
Sc
—
)
Lw
L2w
=
—.
ratio of the crystal and liquid thermal conductivities. The unperturbed gradients are given h~’ 0) C~
=
kk
—
‘T.~° exp[ =
I
1)/k exp[—:], (Pr/Sc) z~j
-
where k is the segregation coefficient. The thermal and solutal Rayleigh numbers are defined as Ra(t) =g(t)aG1(D/V)
a /oK,
Rs(t) =g(t)f3c~(D/V)~/vi), where g(t) is the gravitational acceleration, and and f3 are the thermal and solutal expansion coefficients, respectively. We assume a gravitational acceleration of the form, g(l) g1 + g cos(12t), where g11 ts the background acceleration, g1 is the modulation amplitude, and 12 is the dimensionless modulation frequency, and we write t0t + Ra~cos( lit). Ra(t) Ra I Rs( t ) Rs0 + Rst 1 cos( lit). a
=
= =
(Ia)
Since we have not allowed the interface to deform, it is necessary to relinquish one of the usual solidification boundary conditions that hold at a crystal—melt interface [12]. The dimensionless boundary conditions that we choose to impose at the planar interface. z 0, are w 0, (2a)
(lb)
c
a[Rs(t)c
/
=
+
(Sc/Pr) Ra(
t
) TJ.
=
=
C 1
—
C~ +
wc~ Lc. =
Sc
-
T2
=
(Id)
Here, subscripts indicate and 2/az partial —a2). derivatives The Schmidt the operator L ratio (d of the kinematic viscosity, number, Sc, is the t’, to the solute diffusivity; the Prandtl number, Pr, is the ratio of the kinematic viscosity to the liquid thermal diffusivity, i<; and q is the ratio of the crystal and liquid thermal diffusivities. Since we assume that the heat capacities of the crystal and liquid are identical, q also represents the =
(I
k )c
=
0.
(2h)
1= T, (lc)
q Sc —~—LT.
+
(2c)
~ (2d) In the above formulation, we satistv the conservation relations at the interface, but neglect the relation between interface temperature and concentration. For should the convective instabilities considered here this he a good approximation. The perturbations are required to vanish at the upper boundary, z h~.The temperature perturhations in the solid vanish at Solutions to the above set of equations which have time-periodic coefficients can he obtained using the framework of Floquet theory [13], in which the solutions are represented as the prod=
=
B. T Murray et a!.
/ Gracitational modulation of thermosolutal connection during DS
uct of a temporally periodic function and an exponential function of time. Two different flumerical implementations of Floquet theory were employed to solve the stability problem formulated above and are discussed briefly below. The first approach consists of representing the periodic component of the solutions by a truncated complex Fourier series in time. The solutions for the perturbation quantities are represented by the product of a periodic Fourier series and an exponential term with complex growth rate ~ f(z,
t)
=
ei~ ~
f,n(Z)
eimth,
where f(z, t) represents any of the perturbation quantities in the liquid. Since the temperature in the crystal is only coupled to the temperature in the liquid through the boundary conditions, it is possible to solve for the Fourier components of the crystal temperature analytically. Substitution of the series expansion into the equations in the liquid and boundary conditions yields a set of 32M + 16 coupled two-point boundary value problems in the spatial variable z for the cornplex Fourier coefficients of the periodic solution components. The resulting set of coupled ordinary differential equations subject to the homogeneous boundary conditions yields an eigenvalue problem that is solved as previously described [1] The set of equations contains the complex growth rate if as a parameter. The base state subject to periodic forcing is linearly stable for a given set of parameters if, forrate all udisturbances, 0r of the growth is negative. the the realcalculations, part In setting 0’r 0 allows for the determination of marginal values of the modulation amplitude, ~ and the imaginary part of the growth rate, u 1, for fixed values of the remaining parameters. For many cases, the imaginary part of ci is equal to zero (synchronous response) or 12/2 (subharmonic response). Restricting ci to these values, the number of unknowns is reduced to 16M + 8. The second approach employed consists of approximating the spatial behavior of the disturbance quantities via the pseudospectral technique =
73
in the physical domain as previously described [I]. The approach corresponds to expanding the solutions in terms of a truncated series of Chebyshev polynomials I~(s),
f(z
N
t)
=
~ f~(t) 7~(s), =
where s
=
(2z/h1,)
1 in the liquid region. The temperature in the crystal is calculated numerically in the same manner in this approach. The pseudospectral discretization requires that the expansions satisfy the governing equations at specific collocation points for the Chebyshev polynomials. When implemented in the physical domain the unknowns are the solution values at the collocation points. The spatial differential operators in the governing partial differential equations are replaced by discrete matrix operators. As a result, the governing set of partial differential equations and boundary conditions becomes a set of coupled ordinary differential-algebraic equations in time for the unknown solution values at the collocation points. In this second solution approach, Floquet theory is implemented by constructing a fundamental solution matrix. The columns of this matrix are linearly independent calculated solutions for the unknowns at the end of one forcing period. The eigenvalues of this matrix are the Floquet multipliers from which the complex growth rate if is obtained. The utility of the second approach is that o- itself is the eigenvalue, so that its magnitude is obtained for fixed values of the remaining parameters, which simplifies the search for marginal and gives moremodes; information regarding thevalues proximity of other in particular, this approach allows the calculation of complex conjugate modes which may occur in the thermosolutal problem [4]. —
3. High frequency analysis In this section we discuss two methods of approximation which are valid in the limit of large modulation frequency. For many alloys and processing conditions, this is a relevant limit. The basic methods of analysis have been previously
74
B. 1. Murray ci 0/. / (,ra, Oaiiona/ inodu/anoti of ilierinosolutal tout 1’ctioui iIurin~’1)5
described for solutal convection [3] and the extension for thermosolutal convection is straightforward. Hence, we simply give the resulting cxtensions of these approximations to thermosolutal convection. In the method of averaging, the modulated Rayleigh numbers, Ra1 and Rst I) are proportional to the modulation frequency Ii in the limit of large frequency. The solution variables are decomposed into an averaged part and rapidly varying part, where the integral of the rapidly varying components over one modulation period (synchronous behavior) vanishes. Using this decomposition and averaging the governing equations produces a set of equations for the averaged quantities with coefficients that no longer depend on time. The stability of these equations can now he treated by a standard normal mode analysis, where the time-dependence is of the form exp(ut). The averaged equations for thermosolutal convection are -
(ciL~ L~ ) Sc
=
a2 Rs0~+
—
L2~
:
R~11T
.
~ ~
2
—~
Pr
~(
2
.
+ -
T~).
~-~-
Pr
‘
-
3a) (IC
~: +
—
ciT
—
T
=
+
(3h) (3c)
=
-
-
1 =
L~, Sc —LT. Pr
.t~/ +
~2
Sc
Sc
—,
iT
.
(3d)
L
+
‘:
(1
k)~=0.
—
~1’~
(41-i)
0.
=
(4e)
where A
Pr --
=
~
Pr
+
~—
.
+
~——---
Sc
a2
+ r
~ Sc
At the tipper boundary. conditions are
=
Ii
.
the houndar~
..
~
=
=
-
(4d)
=
By employing the averaging technique, the stahility problem for the effect of sinusoidal niorlulation in the limit of high frequency is reduced from the system (1) with time-periodic coefficients to (3) which depend on alone. In the averaged analysis. the order of the system of ordinary differential equations is increased by two over the usual set of unmodulated equations: the resulting computational times arc greatly reduced over the full set of Floquet equations. Another type of asymptotic analysis was used to describe resonant mode behavior in the limit of large frequency for the single component case [3]. We now extend this analysis to the doublediffusive directional solidification case. We appl~ this analysis in the absence of a background grayitational acceleration, and set Rs~ atid Ra0 equal to zero. For the analysis, we assume the modulated Rayleigh numbers scale as 112 and define 1/
=
Rs1t~/1l2. II
=
Ra~t/l22.
Pr
-
We introduce the new variables lit. I and ~ lIC. In the limit of large 12. (1) reduces to =
where the overhar denotes averaged quantities and the unknown function ~ arises from the analysis and represents the contribution of the rapidly varying components in the averaged equations. The scaled modulation amplitudes are given by .1. Ra1 >/li and Rs1 1/11, The averaged temperature in the crystal can he determined analytically and the resulting boundary conditions on the averaged quantities in the liquid at the interface are =
— — —
~:
=
— —
— —
0
( a)
=
I2T.
=
I ~~1.tt; C~ +
~ +
=
—tu cos =
~
=
=
Sc!! Pr T+ II~
(1. 0
.
(~~i)
(71-i) (Se) -
-
In this approximate analysis, the spatial order of the governing equations has been reduced trom eight to two, only conditions are and that the ii 0 at remaining 0 and Iiboundary =
=
B. T. Murray et a!.
/
Gravitational modulation of thermoso!utal convection during DS
The set of partial differential equations (5) is separable, and we introduce the functions 0(T) and Z(z) such that w(z, r) ê(z, T) T(z, r)
= =
OT(r)Z(z), —0(r)Z(z)c~°~, 0(r)Z(z)T~ ~
=
(~.ccos T)O 2~
Z~.
--
a2
—
ji
=
Ratt0 0. Thus, the system is stable in the absence of modulation and we investigate the modulation amplitude at which instability occurs. For typical growth velocities, modulation frequencies that are greater than 1 Hz correspond to very large dimensionless frequencies. For example, for a growth velocity, V= iO~cm/s, and liquid diffusion coefficient, D l0~cm2/s, a dimension=
=
which results in the following two ordinary differenttal equations 0~+
75
0,
(6a)
s
—TI. T~°~ + Hc~°~ z Pr Z
=
o. (6b)
The equation for 0 is the Mathieu equation which describes the temporal behavior. As the magnitude of j.c increases, the first transition to instability occurs as a subharmonic mode at ~ 0.454. The spatial dependence is described by the function Z which requires solution of a two-point boundary value problem since Z(0) Z(h~) 0. This is an eigenvalue problem and nontrivial solutions are obtained only for specific values of the parameters. Here we choose H~as the eigenvalue parameter and numerically calculate the lowest value for fixed values of the remaining parameters. In the case where Prh,/Sc ~z 1, the unperturbed temperature gradient, T..1tt~, is approximately constant and the eigenvalue problem can be reduced to an algebraic equation involving Bessel functions, The differential eigenvalue problems that arise in the two approximate analyses described above are solved numerically using the approach previously described [12]. =
=
=
4. Results and discussion We carry out a series of calculations for a typical semiconductor alloy with Sc 10, Pr 0.1, and k 0.3. We consider conditions appropriate to processing in a microgravity environment and assume that the steady background gravitational acceleration can be neglected, i.e., we set Rst0~ =
=
=
=
less frequency, 12 i0~is equal to a dimensional frequency of 1.6 Hz. For most of the calculations we set the dimensionless height, hf= 10, and the dimensionless wavenumber, a 0.5. In principle, we should vary the wavenumber to determine the minimum modulation amplitude, but in order to reduce the required computations we only do this for selected cases. Unless otherwise stated, the values of the dimensionless parameters given here will be used in all the calculations. For growth vertically upwards with rejection of a lighter solute (temperature field stabilizing; solute field destabilizing), both the thermal and solutal Rayleigh numbers are defined to be positive. We first consider an alloy rejecting a heavier solute for which during half of the modulation cycle both the temperature and concentration field are stabilizing, while during the other half of the cycle both fields are destabilizing; thus Rs°~ and Ra(i) have different signs. In fig. 1, for Li 1000, we plot the modulated solutal Rayleigh number as a function of the modulated thermal Rayleigh number calculated from the full set of equations and boundary conditions, eqs. (1) and (2). The system is stable below the curve and unstable above. If there were no interaction between the thermal and solutal fields, the stability boundary would be comprised of a horizontal line and a vertical line. The shape of the stability boundary in fig. 1 indicates only a weak interaction between the thermal and solutal fields; the system is unstable for Rs~° greater than 3>< i0~ or Ra~°greater than 700. For this case, the instability is a subharmonic resonant mode which oscillates with half the modulation frequency. We will demonstrate shortly that at this frequency the instability is well characterized by the asymptotic resonant mode analysis. We next consider an alloy for which the temperature field is stabilizing when the solute field =
.
=
=
70
B. 71 Murray ci a!.
/
(iratitational moi/ulatii,n of t/icrmosolutal
Thermosolutal Convection with Vertical g—Jitt.er Sc = 10 k =1°>=0.0 0.3 hi Pr == 10 0.1 Rs(°>=0.0 Ra — . —~
I
—
=
~
coot etilon iluriuig
1)5
tion between the two fields is necessary for occurrence the synchronous synchronous mode mode. can As we show below,of the be scribed by the method of averaging at this
quency. As discussed previously. the method of averaging and the resonant mode analysis assume that the modulation amplitude scales as Ii and 122. respectively, for sufficiently large Ii. In fig.3. we plot Rst I asafunction of Ii for Rst it I OORat I) and the other parameters the same as in figs. I and 2. The points ieprcsent the results Ironi the numerical calculations, while the solid line corre— sponds to Rst I) b12 and the dashed line to Rst 1 1)7.122 where the constants I~and b-. were calculated from the method of averaging and resonant mode analysis, respectively. It is clear that the large frequency asymptotic behavior holds with only slight deviations at the lowest frequency 12 102. In view of this, in fig. 4 we show the full numerical and the numerical asymptotic results for Rs1 it as a function of Rs1 i 1/ Ra1 for 12 =
Stable
the will defre-
±
=
=
111111
I t~
~
to
Fig. I. The modulated solutal Rayleigh number as a function of the modulated thermal Rayleigh number at the onset of suhharmonic instability in the absence of background modulation. The temperature and solutal fields act together,
is destabilizing over half the cycle and vice versa during the other half of the cycle; thus Rsw and Ra~’~ have the same sign. Fig. 2 is a corresponding plot to fig. 1. except that RaW is now positive. For sufficiently small Rst or for sufficiently small Ra~1.subharmonic resonant modes are obtained which are in agreement with the results shown in fig. I. However, away from these limiting values, the instability occurs with synchronous response at significantly lower modulation amplitude than in fig. I. For this synchronous mode, the system is unstable for a range of Ra(i) for Rsw greater than approximately l0~and for a range of RsW if Ra~°is greater than approximately 20. Comparing these results with those in fig. 1, it is interesting to note that instability occurs ft~rmodulation amplitudes two orders of magnitude smaller when the temperature and solute fields oppose each other. This can be attributed to a double-diffusive coupling, since each component taken individually would lead to suhharmonic resonant instahility. While the shape of the stability boundary for the synchronous mode exhibits a weak interaction between thermal and solutal fields, the interac-.
=
.
-
=
Thermosolutal Convection with Vertical g—Jitter Sc = 10 k = 0.3 hi = 10 Rs~°1= 0.0 Rat0~=0.0 Pr = 0.1 0 = 1000 a = 0.5
-
/ /
—
2
~
‘~-~
/
/ ~-
--
I
I flTflf
—
Stable I
rrrrrrr
Ill
Il/UI
liii
I rTI
Ra(i) Fig. 2. The modulated solutal Rayleigh number as a function of the modulated thermal Rayleigh number at the oiiset 0! instability in the absence of background modulation. ihe ..
.
-
temperature and solufal fields oppose each other. Fhe solid line indicates synchronous response and the (lashed line subharnionie response.
B. T Murray ci a!.
/ Gravitational modulation of’thermo,so!uta! connection during DS
l0~.The results from the full equations are mdicated by the data points, while the results for the synchronous modes from the method of averaging are indicated by the solid line and the results for the subharmonic mode from the resonant mode analysis are given by the dashed line. The asymptotic results are in excellent agreement with the full numerical results at this frequency. Using the averaged equations, we have also minimized RsW with respect to wavenumber, a, and the results are shown by the lower dashed curve. This curve deviates only slightly from the a 0.5 curve with the largest difference at lower values of t° Rs°~/Ra~° Wefor notea that quantity is a constant giventhealloy and Rs°~/Ra processing conditions and independent of the amplitude of the gravitational modulation, g The results presented in figs. 1—4 were primarily for a fixed spatial wavenumber, a In figs. 5 and 6, we plot the modulated solutal Rayleigh number as a function of wavenumber for Rst°/Ra°~ ±100. In fig. 5 the thermal and solutal fields act together (Rs~i)/Ra(t) 100) and the resonant instability is subharmonic. The
77
Thermosolutal Convection with Vertical g—Jitter Sc = 10 k = 0.3 h 1 = 10 Rs(°)=0.0 Ra(°)=0.0 Pr = 0.1 0 1000 a 0.5
—
=
=
0.5.
=
—
=
Thermosolutal Convection with Vertical g—Jitter Sc = 10 k = 0.3 h 1 = 10 Rs~°~= 0.0 Ra(°> 0/~Ra(’)~ 0.0 Pr = 100 = 0.1 a = 0.5 Rs(
‘~ ‘
111111/
111111111
id~
FF1111111
11111115
111111111
II
id
if~ 1O~ IO~ Rst’t/IRa0~ Fig. 4. The modulated solutal Rayleigh number as a function of the ratio of the modulated solutal and thermal Rayleigh numbers at the onset of instability. The points represent the full numerical calculations with circles corresponding to Rst i)/Ra(i) positive, triangles to Rs~° /Ra’° negative: solid circles correspond to synchronous modes and open symbols to suhharmonic modes. The solid and dashed curves represent the method of averaging and the resonant mode asymptotic analyses, respectively, for a = 0.5; the lower dashed curve represents the minimum value with respect to wavenumber from the method of averaging.
minimum Rs~1occurs for a wavenumber a 1.5; the value of RsW for a 0.5 is 1.4 times this minimum value. The structure of the curve for larger wavenumbers suggests that a second minimum might occur for somewhat different values of Rs~/Ra~.The dashed curve is calculated from the asymptotic resonant mode analysis. The approximate analysis yields good agreement at small wavenumbers, hut underestimates the prediction from the full numerical calculation for larger wavenumbers. A similar behavior was found for the onset of solutal convection in ref. [31, where the limitations of the resonant mode analysis have been discussed. In fig. 6 the thermal and solutal fields oppose each other (Rs°~/Ra°~ 100) and the instability is synchronous. The mini=
111111111
11111111
10 0 Fig. 3. The modulated solutal Rayleigh number as a function of modulation frequency at the onset of instability. The solid line (synchronous response) has slope unity and the dashed line (suhharmonic response) has slope two,
=
mum Rs~~ occurs for a wavenumber a 0.5 in this case. The dashed curve is calculated from the
70
B. ‘/1 Murray en a!.
/
(inn ‘itationaf modulation
Thermosolutal Convection with Vertical g—Jitter Sc10 k0.3 h =10 10~ 0.0 Ra~°~ 0.0 = 0.1 Rs 0 = 1000 Rs(’)/EaO) = —100
o °
0~
of t/ienInovo/uia/
co,i/i’itntln
t/uu’in,ç’ 1)5
niethod of averaging and is in ti)od agreement with the full numerical calculations with increas.
.
.
ing deviation as the wavenumber increases. How-
ever, in this case, the approximate analysis accurately predicts the minimum value. In space experiments the background gravitational acceleration is very small but not zero. We have carried out a few calculations with a nonzero background acceleration. For Rsttt/Ra°~ 101) and Rst°t 1, 0. 1, the values of Rst I) :tt the onset of instahility are 1735, 1605, 1464, respeclively. As Rst°~increases from zero, the stability of the unmodulated state decreases and thus the modulation amplitude required t’or instability decreases. For Rs~t/Rat I) 100 and Rst 0. 1. there is less than 0.01 ~ change in the modulation amplitude: resonant modes of insta=
~
=
0
H
—
=
0
l’ll’l’
0.0
2.0
4.0
8.0
8.0
10.0
a Fig, 5. The nnodulated solutal Rayleigh number asa function of the wavenuniher at the Onset of suhharmonic instability. The temperature and solutal fields act together. The unstable region is shaded; the dashed curve is based on the asymptotic
—
~:
—
. .
hility are insensitive to the background acceleration level. Based on this limited number of calculations, the case of small background acceleration is reasonably well approximated by setting the background acceleration to zero.
analysis tor resonais t mode behavior.
In our previous paper on solutal convection [I] we neglected thermal effects. For the properties and processing conditions considered here, this is Thermosolutal Convection with Vertical g—Jitter Sc = 10 k = 0.3 hi = 10 Rs~°~ 0.0 Rat0>= 0.0 Pr = 0.1 (4 = 1000 Rs(’>/Ra0) = 100
=
a good approximation when Rat is less than about 600 (see fig. I) for the case when the thermal and solutal fields act together. hut Rat It must he less than approximately It) when the fields oppose each other (see fig. 2), In the ab-
~:~I~1:;l:t5.the resonant mode analwhere ~i ((.454 is the first transition point of the Mathieu equation and r~ is the first zero of the Bessel function. For a 0.5. Sc 10. A ((.3. =
‘~1i’
=
=
=
and 12 1000, the above tormula gives RsC 2.14 X l0~.which is in reasonable agreement with the calculations (see figs. I and 2). From eq. (6) with =
‘‘‘‘ia’’’ ~ a Fig. f. The modulated solutal Rayleigh numter as a function °
~
a’
‘
‘
~.s
of the wavenumber at the onset of synchronous instability. The temperature and solutal fields oppose each other. The unstable region is shaded; the dashed curve is based on the asymptotic analysis using the method of averaging.
t =
the assumption that the temperature gradient is constant, it is possible to obtain a similar expression for the purely thermal problem. namely R’ Ci
i
— —
/5
Pr 12
Sc2
1
+
~-
a/i
-
B. T. Murray et a!.
/ (irai’itational modulation of thermoso!utal convection
during DS
79
For a 0.5, Sc 10, Pr 0.1, h 10, and LI 1001), the formula yields Rate 630, which again agrees with the results of figs. 1 and 2. The appearance of Sc in the above formula for the pure thermal case is a result of our choice of a length scale based on D/V. We recall that the resonant mode analysis is a good approximation for small values of wavenumber, a, but it significantly underestimates the minimum value of the modulation amplitude with respect to wavenumher. All the calculations discussed above were for the value h 1 10. The above formula for the purely thermal problem indicates that the dependence on h( for the resonant modes is fairly weak. A few numerical calculations using the full equations have also shown a weak dependence of the modulation amplitude on h~.We have investigated the dependence of the fundamental mode on h1. using both the full equations and the method of averaging. The dependence is very weak; for example, for R5W/RaW 100 the modulation amplitude for instability calculated from the full equations decreases by 0.4% as h, increases from 10 to 40. For the averaged equations, no further change occurs in the range investigated (Ii, 100).
oppose each other a fundamental, synchronous mode with Rs°~ a 12 occurs as well. The problem of thermosolutal convection in directional solidification with gravity modulation depends on a large number of parameters. We have only investigated a small subset of these parameters here. We have shown that for high frequency modulation the method of averaging provides an accurate description of the fundamental mode which can be the most dangerous mode when the temperature and solutal fields oppose each other. The asymptotic analysis of
5. (.onclusions
similar the fundamental modes manner describedto here in that they synchronous can be the most dangerous modes at high modulation frequency. In addition to the fundamental oscillatory mode, resonant modes occur and can exhibit strong resonance with the modulation frequency.
=
=
=
=
=
=
=
=
In summary, we have studied the stability of thermosolutal convection in a model of vertical directional solidification subject to purely sinusoidal gravitational modulation using a numerical implementation of Floquet theory and two types of approximate analyses for large modulation frequencies. The calculations show that the presence of both the temperature and solute fields can yield a significant coupling which causes instability at lower modulation amplitude than is required for the fields taken individually. This coupling is more significant when the two fields oppose each other than when they act together. The calculations show that when the fields act together, a subharmonic, resonant mode of instability occurs with Rs(i) a ~2 When the fields
resonant modes yields an accurate prediction of the modulation amplitude for instability at low wavenumbers. These approximate analyses significantly reduce the computational effort and should allow for investigation of a wider range of parameters. Although we have concentrated on the situation with zero background gravitational acceleration, the analysis is not limited to this case. For constant gravitational acceleration, the onset of thermosolutal convection can occur by an oscillatory mode. The computational methods presented here can be used to study the effect of sinusoidal modulation on these fundamental oscillatory modes. For stress-free boundaries and constant temperature and solute gradients [4], these fundamental oscillatory modes behave in a
Acknowledgement This work was conducted with the support of the Microgravity Science and Applications Division of the National Aeronautics and Space Administration. R ~ e erences [11B.T. Murray, SR.
Coriell and GB. McFadden, J. Crystal Growth 110 (1991) 713.
B. T Murray
Of)
[21G.Z. 131 [4] [5] [6]
17] 181
Ct
al.
/ Grai’itational
modulation of thcrmo,so!utal coni’ccnon during DS
Gershuni and F.M. Zhukhovitskii. Convective Stahility of lneompressihle Fluids (Keter, Jerusalem. 1976) ch. 8 A.A. Wheeler, GB. McFadden, B.T. Murray and SR. Coriell, Phys. Fluids A 3 (1991) 2847. By. Saunders, B.T. Murray. GB. McFadden. SR. Coriell and A.A. Wheeler, Phys. Fluids A 4(1992)1176. G. Terrones and CF. Chen, J. Fluid Mech., in press. S.H. Davis. Ann. Rev. Fluid Mech. 8 (1976) 57. 5. Ostrach, Ann. Rev. Fluid Mech. 14(1982) 313. R. Monti. J.-J. Favier and D. Langhein. in: Fluid Sciences and Material Science in Space, A European Perspective, Ed. H.U. Walter (Springer. Berlin 1987) p. 637.
191 J.1.D. Alexander. Microgravity Sci. Technol. .7! 109)1) 52. [101 R. Caram, M. Banan and W.R. Wilco~<,J. Crystal Growth 114 (1991) 249. Ill] iS. Turner. Buoyancy Effects in Fluids (Cambridge tiniversity Press, Cambridge, 1973). [12] SR. Coriell, MR. Cordes, Wi. Boeftinger and R.F Sekerka, J. Crystal Growth 49 (198ff) 13. [13] V.A. Yakubovich and V.M. Starzhinskii, Linear [)iff’erenhal Equations with Periodic Coefficients (Wiles’, Ne~ York. 1975).