# Three-dimensional stability of a rotating MHD flow

Thierry Tomasino, Philippe Marty

► To cite this version:

Thierry Tomasino, Philippe Marty. Three-dimensional stability of a rotating MHD flow. Acta Mechanica, 2004, 172 (3-4), pp.135-150. 10.1007/s00707-004-0161-7. hal-00265715

HAL Id: hal-00265715

https://hal.science/hal-00265715v1

Submitted on 24 Mar 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Three-dimensional stability of a rotating MHD flow

T. Tomasino and P. Marty, Grenoble cedex, France

**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, $T_m$, 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.

## 1 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 \cdot 10^6$ produces temperature fluctuations of approximately $3^\circ C$ around a mean value equal to $720^\circ 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 converged 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.
<table><caption>**Table 1.** Critical Taylor number $T_{mc}$ as a function of the aspect ratio $h = H/R$ according to different authors</caption><thead><tr><th>h = H/R</th><th>1</th><th>2</th><th>4</th><th>6</th><th>10</th><th>12</th></tr></thead><tbody><tr><td>Kaiser and Benz<br>(1998) [3]</td><td></td><td>0.54·10<sup>5</sup></td><td>11250</td><td></td><td>5·10<sup>3</sup></td><td></td></tr><tr><td>Marty et al.<br>(1999) [5]</td><td>0.76–0.88·10<sup>6</sup></td><td>1.7–1.8·10<sup>5</sup></td><td>2.6–3.4·10<sup>4</sup></td><td>1.1–1.8·10<sup>4</sup></td><td>4–7·10<sup>3</sup></td><td></td></tr><tr><td>Grants and Gerbeth<br>(2001) [6]</td><td></td><td>1.635·10<sup>5</sup></td><td>2.96·10<sup>4</sup></td><td></td><td></td><td></td></tr><tr><td>Present work</td><td>1.85·10<sup>6</sup></td><td>1.7·10<sup>5</sup></td><td>3.05·10<sup>4</sup></td><td>1.46·10<sup>4</sup></td><td></td><td>5.75·10<sup>3</sup></td></tr></tbody></table>

## 2 Numerical formulation

We consider a cylinder with height $H$ and $R$ filled with an incompressible liquid metal of kinematic viscosity $\nu$, density $\rho$ and electrical conductivity $\sigma$ (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 $\omega$ 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:

$$ \mathrm{div}\vec{V} = 0, \qquad (1.1) $$

$$ \frac{D\vec{V}}{Dt} = -\vec{\nabla}p + \frac{1}{T_{\mathrm{in}}^{1/2}} \nabla^2 \vec{V} + \vec{F}, \qquad (1.2) $$

where the characteristic scales are $R$ for length, $U^* = BR\sqrt{\sigma\omega/(2\rho)}$ for the velocity, $R/U^*$ for time and $\sigma\omega RB^2/2$ for the azimuthal Lorentz forces $F_0$ which have the following analytical expression [5]:

$$ F_0 = r - \sum_{k=0}^{\infty} \frac{2J_1(\lambda_k r)}{(\lambda_k^2 - 1)J_1(\lambda_k)} \left\{ \frac{\sinh(\lambda_k z) + \sinh(\lambda_k (h-z))}{\sinh(\lambda_k h)} \right\}, \qquad (2) $$

where $J_1$ is the Bessel function of first kind and where the $\lambda_k$ are such as $J'_1(\lambda_k) = 0$.

<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Geometry under consideration showing the probing points where the velocity signals have been recorded</figcaption></figure>
The magnetic Taylor number is defined as: $T_m = \sigma\omega B^2 R^4 / (2\rho\gamma^2)$; it is the square of a Reynolds number built with the typical velocity $U^* = BR\sqrt{\sigma\omega/(2\rho)}$ 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 $\theta$. 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^{-\alpha})$ where $\alpha$ is a function of the variable derivability. If this one can be derivated infinitely, $\alpha$ 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 $\varphi(\text{pressure, velocity, temperature and forces})$ was then written as follows:

$$ \varphi(r, z, \theta, t) = \sum_{k=-S/2}^{S/2-1} \hat{\varphi}^k(r, z, t) e^{ik\theta}, \quad (3) $$

where $k$ is the azimuthal wave number, $\hat{\varphi}^k$ the $k$-th Fourier coefficient in the “spectral space” of the variable $\varphi$ and $S$ the number of spectral modes. The Gauss-Labatto formula has been used to complete this discretization:

$$ \theta_j = \frac{2\pi}{S} \quad j = 0, 1, \dots, S-1. \quad (4) $$

The convective terms have been evaluated by using

$$ \dot{w}^k = \frac{1}{S} \sum_{j=0}^{S-1} u_j v_j e^{-ik\theta_j} \quad \text{with} \quad -S/2 \le k \le (S/2 - 1). $$

Then, the following dimensionless continuity and Navier-Stokes equations have to be solved in the “spectral space” for each mode $k$:

$$ \frac{\partial \hat{V}_r}{\partial r} + ik\hat{V}_0 + r \frac{\partial \hat{V}_z}{\partial z} = 0, \quad (5) $$

$$ \begin{aligned} & \frac{\partial \hat{V}_r}{\partial t} + \frac{1}{r} \frac{\partial}{\partial r} \left( \frac{r \hat{V}_r^2}{r} \right) + ik \frac{V_0 \hat{V}_r}{r} + \frac{\partial}{\partial z} \left( V_r \hat{V}_z \right) - \frac{V_0^2}{r} \\ & = \hat{F}_r - \frac{\partial \hat{P}}{\partial r} + \frac{1}{T_m^{1/2}} \left\{ \frac{\partial}{\partial r} \left( \frac{1}{r} \frac{\partial(r \hat{V}_r)}{\partial r} \right) - \frac{k^2}{r^2} \hat{V}_r + \frac{\partial^2 \hat{V}_r}{\partial z^2} - \frac{2}{r^2} ik \hat{V}_0 \right\}, \end{aligned} \quad (6) $$

$$ \begin{aligned} & \frac{\partial \hat{V}_0}{\partial t} + \frac{1}{r} \frac{\partial}{\partial r} \left( r V_0 \hat{V}_r \right) + \frac{V_0 \hat{V}_r}{r} + \frac{ik}{r} \hat{V}_0^2 + \frac{\partial}{\partial z} \left( V_0 \hat{V}_z \right) \\ & = \hat{F}_0 - \frac{ik}{r} \hat{P} + \frac{1}{T_m^{1/2}} \left\{ \frac{\partial}{\partial r} \left( \frac{1}{r} \frac{\partial(r \hat{V}_0)}{\partial r} \right) - \frac{k^2}{r^2} \hat{V}_0 + \frac{\partial^2 \hat{V}_0}{\partial z^2} + \frac{2}{r^2} ik \hat{V}_r \right\}, \end{aligned} \quad (7) $$

$$ \begin{aligned} & \frac{\partial \hat{V}_z}{\partial t} + \frac{1}{r} \frac{\partial}{\partial r} \left( r V_z \hat{V}_r \right) + \frac{ik}{r} \left( V_0 \hat{V}_z \right) + \frac{\partial}{\partial z} \left( \hat{V}_z^2 \right) \\ & = \hat{F}_z - \frac{\partial \hat{P}}{\partial z} + \frac{1}{T_m^{1/2}} \left\{ \frac{1}{r} \frac{\partial}{\partial r} \left( r \frac{\partial(\hat{V}_z)}{\partial r} \right) - \frac{k^2}{r^2} \hat{V}_z + \frac{\partial^2 \hat{V}_z}{\partial z^2} \right\}. \end{aligned} \quad (8) $$
In the above equations, the force $\vec{F} = (F_r, F_\theta, 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.

### 3 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 $\Omega$. For particular Reynolds numbers ($Re = \Omega R^2/\gamma$) 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 which is in excellent agreement with the experimental

<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Vertical profiles of the axial velocity $-V_z$ at $r=0$ for a vortex breakdown induced when $Re=1854$ and $H/R=2$; a our results, b Sorensen et al. [12] results for different meshes (the curve indicated by the arrow corresponds to the finest grid)</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>Fig. 3. Time evolution of the streamlines in the vortex breakdown problem for $Re=1854$ and $H/R=2$; grid size: 70×140 nodes; time step: 0.01; the time scale which has been used here is the rotation period of the lid, $\Omega^{-1}$</figcaption></figure>

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-Bé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 $\Delta T = T^+ - T^-$ between the end walls exceeds a certain threshold. The dimensionless numbers of this Rayleigh-Bénard problem are the Prandtl number and the Rayleigh number built with the characteristic length $H$ and the characteristic velocity $V = \sqrt{g\beta H \Delta T}$. They write as follows: $Ra = \frac{\beta g H^3 \Delta T}{\gamma^\alpha}$ and $Pr = \frac{\gamma}{\alpha}$, where $\alpha$ is the thermal diffusivity of the fluid.

An aspect ratio $h=H/R=2$ and a low Prandtl number, $Pr = 0.022$, equal to that of mercury, have been considered for different values of the Rayleigh number. The meridional grid has been refined from 33 × 65 for small values of $Ra$ to 65 × 129 for $Ra \approx 10^6$. In the
<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Time evolution of the non-axisymmetric mode $k=1$ of the “azimuthal energy” during the Rayleigh-Bénard transient regime for different Rayleigh numbers ($Pr=0.022$ and $h=2$): **a** Verzicco et al. [13] results with, from the lower to the upper curve, $Ra = 3600, Ra = 3700, Ra = 3750, Ra = 3800, Ra = 4000$; **b** our results: $Ra = 3800, Ra = 3900, Ra = 3950, Ra = 4000, Ra = 4100$</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Time evolution of the axial velocity $v_z$ in ($r=0.5$; $z=1$) for $h=2$ and **a** $T_m = 1.6 \cdot 10^5$, **b** $1.7 \cdot 10^5$, and **c** $1.8 \cdot 10^5$</figcaption></figure>

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
<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Time evolution of the different velocity components for $h = 1$ and $T_m = 1.9 \cdot 10^6$: **a** mid-plane probe in ($r = 0.5$; $z = h/2$), and **b** upper end wall probe in ($r = 0.75$; $z = 0.95 h$)</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Fig. 7. Streamline pattern for $h = 1$ and $T_m = 1.9 \cdot 10^6$</figcaption></figure>

<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 8. Time evolution of the different velocity components for $h = 2$ and $T_m = 1.8 \cdot 10^5$: **a** mid-plane probe in $(r = 0.5; z = h/2)$, and **b** upper end wall probe in $(r = 0.75; z = 0.95 h)$</figcaption></figure>

$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}$.
<figure><img src="image_11.png" /><figcaption>Fig. 9. Time evolution of the meridional streamlines for $T_m = 180000$ and $h = 2$</figcaption></figure>

# 4 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 \approx T_{mc}$ and for magnetic Taylor numbers approximately equal to twice the critical value, i.e., for $T_m \approx 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.

## 4.1 Results for $T_m \approx 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
<figure><img src="image_12.png" /><figcaption>**Fig. 11.** Time evolution of the modes $k=0,1,2$ and $3$ of the azimuthal velocity ($v_0$) at mid-height ($z=0.5h$, $r=0.75$) for $h=1$, $T_m=3.6 \cdot 10^6$, $100 \times 16 \times 100$ nodes and $\Delta t=0.005$: **a** modes $k=0$ and $k=1$, **b** modes $k=2$ et $k=3$; the energy of the modes $k=1,2$ and $3$ has been scaled with that of the mode $k=0$</figcaption></figure>

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 end-walls. 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 \cdot 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_0$ 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_0^{-1}=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
<figure><img src="image_13.png" /><figcaption>**Fig. 12.** Time evolution of the modes $k = 0, 1, 2$ and $3$ of the azimuthal velocity ($v_0$) at mid-height ($z = 0.5h$, $r = 0.75$) for $h = 6$, $T_m = 2.5 \cdot 10^4$, $55 \times 16 \times 300$ nodes and $\Delta t = 0.001$: **a** modes $k=0$ and $k=1$, **b** modes $k=2$ et $k=3$; the energy of the modes $k=1, 2$ and $3$ has been scaled with that of the mode $k=0$</figcaption></figure>

$T_m = 1.6 \cdot 10^5$, these authors find a dimensionless pulsation $\omega = 2\pi f$ equal to 1659. Using our physical scales changes this value into $\omega = \frac{1659}{T_m^{1/2}} = 4.147$. This value is remarkably close to the value $\omega = 2\pi f_0 = 4.14$ which is deduced from our value of $f_0$.

## 4.2 Results from $T_m \approx 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 \approx 2T_{mc}$ and $h=1, 2$ and $6$.
<figure><img src="image_14.png" /><figcaption>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 \cdot 10^4, 55 \times 16 \times 300$ nodes and $\Delta t = 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$</figcaption></figure>

- For $h=1$ and $T_m = 3.6 \cdot 10^6 = 1.9T_{mc}$, 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 \cdot 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 non-axisymmetric mode and reaches 10% of the mode $m=0$ in the Ekman layers.
- For $h=6$ and $T_m = 2.5 \cdot 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 \sim 1$.

## 5 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 \approx 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.

## Acknowledgement

This work has been supported by the French-German CNRS-DFG program. Dr. Martin-Witkowski at Urbana-Champaign is thanked for giving constructive remarks on the first manuscript. Fruitful discussions with Dr. Richard Pasquetti in Nice and with Dr. Gunter Gerbeth in Dresden are acknowledged.

## References

[1] Dold, P., Benz, K. W.: Modification of fluid flow and heat transport in vertical Bridgman configuration by rotating magnetic fields. Crys. Res. Technol. 32(1), 51–60 (1997).

[2] Dold, P., Benz, K. W.: Rotating magnetic fields: fluid flow and crystal growth applications. Progress in Crystal Growth and Characterization of Materials 38(1–4), 7–38 (1999).

[3] Kaiser, T., Benz, K. W.: Taylor vortex instabilities induced by a rotating magnetic field : a numerical approach. Phys. Fluids 10(5), 1070–1109 (1998).

[4] Davidson, P. A., Hunt, J. C. R.: Swirling recirculating flow in a liquid metal column generated by a rotating magnetic field. J. Fluid Mech. 185, 67–106 (1987).

[5] Marty, P., Martin Witkowski, L., Trombetta, P., Tomasino, T., Garandet, J. P.: On the stability of rotating MHD flows. In: Transfer phenomena in magnetohydrodynamic and electroconducting flows (Alemany, A., Marty, Ph., Thibault, J. P., eds.). Kluwer Academic Publishers 1999.

[6] Grants, I., Gerbeth, G.: Stability of axially symmetric flow driven by a rotating magnetic field in a cylindrical cavity. J. Fluid Mech. 431, 407–426 (2001).

[7] Grants, I., Gerbeth, G.: Linear three-dimensional instability of a magnetically driven rotating flow. J. Fluid Mech. 463, 229–239 (2002).
[8] Grants, I., Gerbeth, G.: Experimental study of non-normal nonlinear transition to turbulence in a rotating magnetic field driven flow. Phys. Fluids 15(10), 2803–2809 (2003).

[9] Walker, J. S., Martin Witkowski, L.: Linear stability analysis for a rotating cylinder with a rotating magnetic field. Phys. Fluids 16(7), 2294–2299 (2004).

[10] Harlow, F. H., Welsh, J. E.: Numerical calculation of time dependent viscous incompressible flow of fluid with free surface. Phys. Fluids 8, 2182–2189 (1965).

[11] Escudier, M. P.: Observations of the flow produced in a cylindrical container by a rotating end wall. Exp. Fluids 2, 189–196 (1984).

[12] Sorensen, J. N., Phuoc Loc, T.: High-order axisymmetric Navier-Stokes code: description and evaluation of boundary conditions. Int. J. Num. Meth. Fluids 9, 1517–1537 (1989).

[13] Verzicco, R., Camussi, R.: Transitional regimes of low-Prandtl thermal convection in a cylindrical cell. Phys. Fluids 9, 1287–1295 (1997).

[14] Cioni, S., Ciliberto, S., Sommeria, J.: Strongly turbulent Rayleigh-Bénard convection in mercury: comparison with results at moderate Prandtl number. J. Fluid Mech. 335, 111–140 (1997).

**Authors' address:** T. Tomasino and P. Marty, LEG1, BP53, 38041 Grenoble cedex, France
(E-mail: philippe.marty@hmg.inpg.fr)