Three-dimensional stability of a rotating MHD flow

Summary.The increase of quality of semi-conductors grown under rotating magnetic fields has recently motivated many numerical studies aiming at predicting the flow induced by the associated electromagnetic body forces in a finite length container filled with a liquid metal. The mean feature of the basic flow is now well known but the hydrodynamic stability of the solution is the subject of many researches which are conducted either experimentally or numerically. From these stability studies, it appears that the value of the critical magnetic Taylor number, Tm, which is deduced is a decreasing function of the aspect ratio H/R of the cylinder. Nevertheless, most of the numerical studies have assumed an axial symmetry of the flow. This problem is reconsidered here with a three-dimensional code. Surprisingly, it is found that no significant azimuthal modes develop for Taylor numbers up to twice the critical value. Even small, the non-axisymmetric modes are found to grow faster near the Ekman layers than in the core flow.


Introduction
The principle of the crystal growth technique is to solidify a molten semiconductor and an additional dopant which are placed in a cylindrical volume submitted to a temperature gradient. Many different configurations are used. In the vertical Bridgman technique, the fluid is bounded by the solid walls of the growth ampoule, whereas in the Floating-zone technique, the vertical lateral boundaries of the liquid crystal are free surfaces. The always increasing size of crystal growth installations is responsible for more important buoyancy or Marangoni forces which lead to unsteady flow and temperature field in the melt. This has been a motivation for considering the use of a magnetic field to damp these fluctuations which are detrimental to the solidified crystal. At first, many attempts have been made to use constant magnetic fields which are known to brake all kinds of motion. The results are positive in several situations but the pattern of the magnetic field distribution has to be carefully determined. Recently, a rotating magnetic field (RMF) has been experimentally tested and its effect seems promising as it is detailed below. One of the first experimental evidence of the ability of a RMF to damp temperature and flow fluctuations was given by Dold and Benz [1] who have built a test cell for temperature measurements in a cylindrical cell submitted to a vertical temperature gradient. When no RMF is applied, a Rayleigh number of 1.1Á10 6 produces temperature fluctuations of approximately 3 C around a mean value equal to 720 C. For the same Rayleigh number, a 20 Hz RMF with an intensity as low as 1.63 mT reduces the temperature fluctuations by a factor 10. Moreover, the frequency of the fluctuations was observed to be higher when a RMF is applied [1]. The solidification of a Germanium doped gallium sample has confirmed that low intensity dopant striations are obtained when a RMF is applied and that the width of these striations is also smaller. The same positive influence of a RMF is reported by Dold and Benz [2] when applied to the Floating-zone configuration. In this case, the Marangoni flow generated by the shear stress along the free surface seems to generate an intense velocity field, and a slightly stronger magnetic field intensity (but still reasonable) is necessary to cancel the flow and temperature fluctuations (B=7.5 mT with a frequency of 50 Hz).
The exact physical reason which explains the positive effect of a RMF has never been clearly clarified. One can conjecture that the bulk rotating motion given to the core flow prevents, or at least delays, the triggering of hydrodynamic instabilities in the liquid. One is not yet able to predict the exact magnetic field intensity which will stabilize the flow in a real industrial situation, but, at least, it has been known for a long time that a too strong rotation will generate centrifugal instabilities similar to the well known Taylor vortices. This has been the motivation for many recent numerical works on the stability study of a rotating MHD flow driven by a RMF.
Kaiser and Benz [3] have performed a two-dimensional computation of the flow induced by the RFM. As mentioned earlier by Davidson et al. [4], they have shown that the flow consists of a mean azimuthal movement which induces a meridional recirculation owing to the weakening of the centrifugal forces near the top and bottom wall of the container which contains the liquid metal. Still using a two-dimensional computer code, Marty et al. [5], have expressed the dependence of the critical Taylor number on the aspect ratio h ¼ H=R of the container. In this study it has been shown that an increase of the Taylor number is responsible for the occurrence of unsteady vortices which are converced by the mean meridional flow. A recent stability study of the axially symmetric solution (Grants and Gerbeth [6]) has shown that several linearly unstable solutions can exist close to the stable basic state. The triggering of these instabilities has been studied as a function of the amplitude of the perturbations. Soon after this work, the same authors [7] have also performed a three-dimensional linear stability study of the flow. All the variables have been expanded into Chebyshev polynomials in the meridional plane whereas a Fourier transform has been used for describing the azimuthal variations. Their conclusion is that non-axisymmetric modes first become unstable. When the aspect ratio h is increased, they have observed that the wave number of the most unstable modes decrease from 3 to 2 and from 2 to 1 for h ¼ 1:31 and 2.47, respectively. It is worth mentioning a clever experimental work done by the same authors [8], which confirms some of the results of their previous theoretical analysis. Recently, Walker et al. [9] have presented an extension of the work done by Grants and Gerbeth [7] by adding a rotation of the container either in the same direction as that of the magnetic field or in the opposite direction. Again, a spectral decomposition has been used in the meridional plane while a Fourier transform has been employed in the azimuthal direction. Table 1 presents the value of the critical value of the Taylor number as a function of the aspect ratio h of the cavity. As can be seen, a certain disagreement between these results exists. It can be explained by the difficulty of determining the precise value of the critical Taylor number with a numerical tool. As a matter of fact, the instability can take a long time to grow and it needs a considerable computational effort to obtain visible flow oscillations. These discrepancies between different authors have been a motivation to perform a fully numerical computation of the flow.
Consequently, this study presents the stability results obtained with a three-dimensional numerical computation of the problem. Our approach is a Direct Numerical Simulation (DNS) and differs from Grants et al. (2002) [7], who linearized the Navier-Stokes equations and then performed a stability study of the remaining terms.
We consider a cylinder with height H and R filled with an incompressible liquid metal of kinematic viscosity m, density q and electrical conductivity r ( Fig. 1) The aspect ratio is defined as h ¼ H=R. The whole cylinder is submitted to a rotating magnetic field with an angular frequency x and intensity B. The frequency of the RMF is expected to be small enough to neglect any magnetic field distortion inside the liquid metal (i.e., the magnetic field distribution is the same as it would be in vacuum). Except for numerical tests which will be described later for which additional temperature terms have been added in the vertical momentum equation, the dimensionless flow equations write: where the characteristic scales are R for length, U Ã ¼ BR ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi rx=ð2qÞ p for the velocity, R=U Ã for time and rxRB 2 =2 for the azimuthal Lorentz forces F h which have the following analytical expression [5]: where J 1 is the Bessel function of first kind and where the k k are such as J 0 1 ðk k Þ¼0.  The magnetic Taylor number is defined as: T m ¼ rxB 2 R 4 =ð2qc 2 Þ: it is the square of a Reynolds number built with the typical velocity U Ã ¼ BR ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi rx=ð2qÞ p resulting from a balance between inertial forces in the fluid and electromagnetic forces. The time-averaging of the radial and axial forces F r and F z is equal to zero.
The finite-difference technique has been used in the ðr; zÞ plane and a spectral method has been considered in the azimuthal direction h. A pressure-velocity formulation associated to a staggered ''Marker and Cell mesh'' [10] has been chosen. A second-order fractional step method has been used for the time differencing. The precision of a spectral method is known to be OðN Àa Þ where a is a function of the variable derivability. If this one can be derivated infinitely, a becomes much greater than one and the precision of the calculation increases. Using the Fourier decomposition constitutes a natural approach for the azimuthal discretization. Each variable u(pressure, velocity, temperature and forces) was then written as follows: where k is the azimuthal wave number,û k the k-th Fourier coefficient in the ''spectral space'' of the variable u and S the number of spectral modes. The Gauss-Labatto formula has been used to complete this discretization: The convective terms have been evaluated by usinĝ u j v j :e Àikh j with À S=2 k ðS=2 À 1Þ: Then, the following dimensionless continuity and Navier-Stokes equations have to be solved in the ''spectral space'' for each mode k: In the above equations, the forceF ¼ðF r ; F h ; F z Þ can be considered as a source term. The r and z derivatives have been discretized with a classical second-order finite-difference technique. A second-order fractional step method associated to a Runge-Kutta 3 (RK3)/ Crank-Nicolson scheme has been used for the time differencing. A grid size of the order of 100Â200 has been used in the meridional plane and 16 to 32 modes have been used in the azimuthal direction. The time step is adapted to each case but the value of 0.001 is a mean value.

Validation tests
The importance of the accuracy of the numerical model has been shown in [6]. This is the reason why a significant effort has been done for validating our numerical model. The meridional discretization in the ðr; zÞ plane was first validated by comparing the velocity profiles for various magnetic Taylor numbers to those obtained with the 2-D code used in [5]. Another validation has been done with the well-known axisymmetric vortex breakdown phenomena. This phenomenon appears when one end wall of a cylindrical container which contains a fluid is driven into rotation at a constant angular velocity X. For particular Reynolds numbers ðRe ¼ XR 2 =cÞ and aspect ratios ðH=RÞ, the flow is purely axisymmetric and one vortex or more appear along the central axis. This causes the vertical velocity along the z axis to change its sign. This is visible in Fig. 2a  results from [11] and with the numerical computations from [12] Fig. 2b. The transient development of the vortices near the Oz-axis is also in fine agreement with Escudier (1984) experiments: this is reported in Fig. 3. This first test validates the ðr; zÞ plane discretization and the time-discretization scheme of our code. The second validation test is the Rayleigh-Be´nard problem which has allowed us to validate the azimuthal discretization and the spectral method effectiveness. This test has also confirmed that the code was able to describe the transition from 2-D to 3-D flows.
The energy equation has been added and discretized similarly to the Navier-Stokes equations. A buoyancy term using the Boussinesq approximation has been incorporated in the vertical momentum equation. In such a situation, the energy and the Navier-Stokes equations are strongly coupled.
The bottom end wall was then maintained at a temperature greater than the top wall whereas the vertical walls are thermally insulated. It is well known that a convection flow appears when the temperature difference DT ¼ T þ À T À between the end walls exceeds a certain threshold.  same time, the number of spectral modes has been changed from S ¼ 8 to S ¼ 32. The critical value of the Rayleigh number ðRa c Þ at which the flow becomes 3-D (i.e., for k=1) has been measured (Fig. 4). A good agreement has been found since our critical value Ra c =3900 is only 4% higher than numerical results obtained by Verzicco and Camussi [13].
For higher values of Ra, when the flow becomes turbulent, the asymptotic law Nu ¼ 0:155Ra 0:265 has been found for the Nusselt number between the two cylinder end walls. This result is in correct agreement with the experimental results of Cioni et al. [14] obtained on mercury: Nu ¼ 0:140Ra 0:260 , and with the numerical results of Verzicco and Camussi [13]. Nu ¼ 0:165Ra 0:255 .

Results
The flow stability has been studied by starting the time-dependent computation at t ¼ 0 with a fluid at rest and for a given value of T m and of the aspect ratio h. Two different cases have been considered: near the critical stability threshold, i.e., for T m % T mc and for magnetic Taylor numbers approximately equal to twice the critical value, i.e., for T m % 2T mc : The values of h which have been studied vary from h ¼ 1 to h ¼ 12. The grid size is 100Â200 in the ðr; zÞ plane with S=16 modes in the azimuthal direction. A time step of the order of 0.001 is used.

Results for T m % T mc
When the magnetic Taylor number T m is smaller than the critical value T mc , the time history of each velocity component at a given point of the flow shows a progressive evolution of the amplitude which finally keeps a constant value after the spin-up time has been reached. For h ¼ 2 and T m = 160000, Fig. 5a shows the time evolution of the axial velocity v z in r ¼ 0:5 and z ¼ 1, i.e., at mid-height of the cylinder. One can see that the flow has become almost steady for t larger than approximately 50. For higher values of T m , an unstable oscillatory regime can be observed: Fig. 5b and c shows this behavior for T m =170000 and 180000, respectively. For h ¼ 2 the value T mc =170000 can be obtained. When h is varied from 1 to 12, the value of the critical Taylor number which has been obtained is reported in the last row of Table 1. These results are in good agreement with previous results obtained with 2-D, axially symmetric codes [5], [6], but differ from those from [3]. The most striking result from this work is that for each value of h which has been considered, absolutely no azimuthal variation has been noticed, i.e., no spectral mode except k ¼ 0 is sufficiently large to qualify the flow as non-axisymmetric. This means that the vortices which are created stay axially symmetric. A frequency analysis of the time history of the velocity components shows that the signal is not rigorously a sinusoidal function of time. This has been shown through the recording of three different probes located in the ðr; zÞ plane in (0.5; h/2), (0.75 ; 0.05h) and (0.75 ; 0.95h), respectively (see Fig. 1b). The first probe is then located in the mid-plane whereas the two others are close to the end walls and situated near the Ekman layers.
-For h=1 and T m =1.9Á10 6 , the mid-plane probe reveals that v h has a 90 shift with v r and v z . This is shown in Fig. 6a. A Fourier analysis shows that the signal consists of a main frequency, say f 0 , which is such as f 0 =0.5, together with the frequency 3 2 f 0 . The amplitude of the 100Â16Â100 nodes and Dt ¼ 0:005: a modes k ¼ 0 and k ¼ 1, b modes k ¼ 2 and k ¼ 3; the energy of the modes k ¼ 1, 2 and 3 has been scaled with that of the mode k ¼ 0 latter frequency is only 10% of that of f 0 . The end walls probe signal is shown in Fig. 6b where we see that the almost purely sinusoidal nature is lost. A frequency analysis confirms that, in addition to the frequencies f 0 and 3 2 f 0 , a frequency equal to 2f 0 has appeared with an amplitude equal to 30% of that of f 0 . This would seem to indicate that the place where instabilities develop with the faster growth rate is near the Ekman layers close to the endwalls. A careful observation of the meridional streamlines over one period (which has the value 1=f 0 ) does not allow to see any noticeable change (Fig. 7).
-Giving a higher aspect ratio h ¼ 2 to the cylinder together with T m ¼ 1:8 Á 10 5 , which is a few percent larger than the critical threshold, also shows that the core probe located in ðr ¼ 0:5; z ¼ h=2Þ displays an almost purely sinusoidal variation with the frequency f 0 =0.66. The velocity v r is shifted 90 with respect to v h and v z which are in opposite phase (Fig. 8a). The tendency to the period-doubling phenomenon which has been observed for h ¼ 1 on the end-wall probes is confirmed in Fig. 8b. Nevertheless the frequency analysis on these end probes shows that only f 0 and 2f 0 are present in the signal. Following the streamlines evolution over one period f À1 0 =1.5 shows that vortices are periodically generated along the central axis r ¼ 0 (Fig. 9). This phenomenon is reminiscent of the vortex breakdown problem which has been discussed earlier. These vortices are quite different from the so-called Taylor vortices which have been reported in most of the previous studies on this problem. Finally, it is interesting to compare the value of the frequency f 0 which has been observed with that predicted by the recent analytical stability study published in [6]. For h ¼ 2 and

Results from T m % 2T mc
In order to know how long the axial symmetry will survive when the Taylor number is increased further, numerical simulations have been done for T m % 2T mc and h=1, 2 and 6.  Fig. 10 shows the time evolution of the modes k=0,1,2 and 3 of the azimuthal velocity in one of the two near-wall probes. It is shown that the non-axisymmetric modes are almost equal and limited to 0.1% of the mode k=0. They start to grow near t=20 which corresponds to the spin-up time of this problem. An analysis of the mid-plane probe displays similar results except that the velocity signals look less turbulent and are 10 times smaller than near the walls (Fig. 11). -For h=2 and T m ¼ 3:0 Á 10 5 ¼ 1:8T mc , the modes k=1,2 and 3 near the end walls are roughly 10 times greater than in the mid-plane. The mode k=2 is now the dominant nonaxisymmetric mode and reaches 10% of the mode m=0 in the Ekman layers.  Fig. 13. Time evolution of the modes k ¼ 0; 1; 2 and 3 of the azimuthal velocity ðv 0 Þ near the lower end wall (z ¼ 0:5h, r ¼ 0:75) for h ¼ 6, T m ¼ 2:5Á10 4 ,55Â16Â300 nodes and Dt ¼ 0:001 : a modes k ¼ 0 and k ¼ 1, b modes k ¼ 2 and k ¼ 3; the energy of the modes k ¼ 1, 2 and 3 has been scaled with that of the mode k ¼ 0 -For h=6 and T m ¼ 2:5 Á 10 4 ¼ 1:7T mc , the same remarks can be made, i.e., the mode k=1 has the faster growth rate and its intensity is roughly 10 times larger near the end walls than in the core (see Figs. 12 and 13).
For each one of these three values of h, the Fourier spectrum of the velocity signal is no more discrete as it was near the critical threshold T mc : it is more continuous and limited by the dimensionless frequency f $ 1.

Conclusion
This work shows that the axial symmetry of a rotating flow driven by a RMF is kept far beyond the critical threshold at which unsteady vortices appear. Near the threshold, i.e., for T m % T mc , the velocity fluctuations are almost sinusoidal with a main frequency f 0 which compares well with recent analytical predictions [6]. The intensity of the secondary higher frequencies is found higher near the end walls than in the core near the mid-plane. For h ¼ 2, a periodic growth of vortices is noticed near the flow axis.
For values of T m twice as large as the critical value, the non-axisymmetric modes remain smaller than the mode k=0 but the mode 2 can reach 10% of the mode 0 for h ¼ 2.
Whatever the value of T m , the growth rate of the instabilities is found larger near the end walls. This seems to indicate that the Ekman layers can be the primary source of instability, at least for the aspect ratios which have been considered.