Numerical insights into magnetic dynamo action ... - repository.tudelft.nl

Numerical insights into magnetic dynamo action in a turbulent regime Sasa ˇ Kenjeresˇ 1 and Kemal Hanjalic´ 2 Department of Multi Scale Physics and JM Burgerscentre for Fluid Dynamics, Delft University of Technology, Lorentzweg 1, 2628 CJ Delft, The Netherlands E-mail: [email protected] New Journal of Physics 9 (2007) 306

Received 16 February 2007 Published 31 August 2007 Online at http://www.njp.org/ doi:10.1088/1367-2630/9/8/306

We report on hybrid numerical simulations of a turbulent magnetic dynamo. The simulated set-up mimics the Riga dynamo experiment characterized by Re ≈ 3.5 × 106 and 15 Rem 20 (Gailitis et al 2000 Phys. Rev. Lett. 84, 4365–8). The simulations were performed by a simultaneous fully coupled solution of the transient Reynolds-averaged Navier–Stokes (T-RANS) equations for the fluid velocity and turbulence field, and the direct numerical solution (DNS) of the magnetic induction equations. This fully integrated hybrid T-RANS/DNS approach, applied in the finite-volume numerical framework with a multi-blockstructured nonorthogonal geometry-fitted computational mesh, reproduced the mechanism of self-generation of a magnetic field in close accordance with the experimental records. In addition to the numerical confirmation of the Riga findings, the numerical simulations provided detailed insights into the temporal and spatial dynamics of flow, turbulence and electromagnetic fields and their reorganization due to mutual interactions, revealing the full four-dimensional picture of a dynamo action in the turbulent regime under realistic working conditions.

Abstract.

1

Author to whom any correspondence should be addressed. Present affiliation: Marie Curie Chair, Department of Mechanics and Aeronautics, University of Rome, ‘La Sapienza’, Italy.

2

New Journal of Physics 9 (2007) 306 1367-2630/07/010306+29$30.00

PII: S1367-2630(07)44493-7 © IOP Publishing Ltd and Deutsche Physikalische Gesellschaft

2

DEUTSCHE PHYSIKALISCHE GESELLSCHAFT

Contents

1. Introduction 2. Numerical model 2.1. Equations and subscale model of turbulence . . . . . . . . . . . . . . . . . . . 2.2. The finite-volume discretization of equations . . . . . . . . . . . . . . . . . . . 2.3. Specification of boundary conditions . . . . . . . . . . . . . . . . . . . . . . . 3. Results and discussion of the numerical simulations 4. Conclusions and future prospects Acknowledgments References

2 4 6 9 10 11 18 27 27

1. Introduction

Many events and phenomena in nature and technology are governed by complex interactions between fluid flow, turbulence, heat transfer and electromagnetic fields. Well-known examples are the origin and self-sustenance of magnetic fields in planets, stars and galaxies, but these interactions are also encountered in deep-space exploration, space craft-propulsion systems, satellite navigation, processes in metallurgy, semiconductor crystal growth, arc-welding, magnetic drug targeting, etc., Moffatt [1, 2], Rüdiger and Hollerbach [3], Wielebinsky and Krause [4], Glatzmaier and Roberts [5], Roberts and Glatzmaier [6], Busse [7], Glatzmaier [8]. The immense variety of scales characterizing these phenomena (ranging from stars to semiconductor crystals) and the intensity of the interactions (weak or strong coupling) combined with the still unsolved questions related to the origin of turbulence—make this subject particularly challenging. In this paper, we focus on the mechanism of partial conversion of the mechanical energy of a flowing electrically conductive medium into the magnetic energy, referenced as a magnetic dynamo. The magnetic dynamo takes place when stretching of a magnetic field overcomes its resistive damping. To provide favourable conditions for a magnetic dynamo, the critical threshold of working parameters must be reached. Since the magnetic dynamo always involves interactions between the fluid-flow and electromagnetic fields, the two sets of non-dimensional parameters are needed. For underlying fluid flow and turbulence intensity this is the hydrodynamic Reynolds number defined as Re = UD/ν representing the ratio between the inertial and viscous forces. The second characteristic parameter is the ratio of inertial forces and electromagnetic damping, defined as the magnetic Reynolds number, Rem = UD/λ, where U, D, ν, λ (= 1/µ0 σ) are the characteristic velocity, length scale, kinematic and magnetic viscosity, respectively, and µ0 and σ are the magnetic permeability and electric conductivity of a working fluid. The basic prerequisite for a magnetic dynamo requires that Rem > 1. In order to experimentally reach favourable conditions for a magnetic dynamo, relatively large length scales and high velocities are needed even for fluids with the highest electrical conductivity (e.g. liquid sodium, λ ≈ 0.1 m2 s−1 ). In such conditions, the fluid flow is highly turbulent, so obviously the turbulence plays an important role in the fluid-flow—electromagnetic-field interactions. It is easy to imagine that experimental studies face many practical limitations—from large dimensions of set-ups to potentially hazardous working fluids (sodium). The first successful realization of a magnetic New Journal of Physics 9 (2007) 306 (http://www.njp.org/)

3


dynamo in a laboratory took place in 1999 when experimentalists in Riga (Gailitis et al [9]–[15]) and in Karlsruhe (Stieglitz and Müller [16], Müller and Stieglitz [17], Müller et al [18, 19]) independently detected the self-generation of a magnetic field. Very recently, experimental groups from Lyon/Paris/Cadarache (VKS2 dynamo) also reported the successful recordings of a self-generated magnetic field, Bourgoin et al [20], Monchaux et al [21]. The most realistic experimental studies have been reported in Berhanu et al [22] where the magnetic fields reversals have been observed for the very first time. All these successful experiments used constrained working conditions with imposed mechanical forcing. As such, their primary goal was not to mimic exactly the geodynamo, but rather to provide experimental proof of a dynamo action with self-generation and selfsustenance of a magnetic field. These experiments opened a route for the new generation of lessconstrained dynamos that will mimic more closely the actual Earth-like scaled conditions, as currently investigated in the Maryland dynamo, Peffley et al [23] or the Grenoble dynamo, Nataf et al [24]. Despite their indisputable success, the past experimental studies mentioned above focused mainly on the time recording of the magnetic field at selected locations. Such single point measurements cannot fully reveal the detailed spatial and temporal dynamics of a selfgenerated magnetic field and its interactions with the underlying turbulent flow. Detailed measurements of fluid flow and turbulence are difficult to perform in opaque working fluids such as liquid metals because of their inaccessibility to optical/laser measurements such as laser doppler anemometry (LDA) or particle image velocimetry (PIV). At present, the only way to provide detailed three-dimensional information on complex physics of timedependent fluid flow, turbulence and electromagnetic field interactions is numerical simulation. Of course, before exploiting the numerical results to gain deeper physical insights into specific phenomena, the computer simulations must be verified and proven to be numerically accurate and the mathematical models used to close the equation set to be based on solid physical foundation. In his recent review, Glatzmaier [8] addressed advances in computational simulations of convection and magnetic field generation in the geodynamo applications. The relevance of practically all (up to now) performed simulations has been questioned since none are able to operate in realistic (typical Earth-like) regimes due to limitations in numerical resolutions. Furthermore, Glatzmaier rightly argued that one of the weakest points in the present generation of geodynamo models lies in the overly crude approximation of turbulent transport. Practically all geodynamo models have been based on simple linear stress-strain models and the assumption that the turbulent diffusivity is constant in space and time even though this is in its essence a time and spatially-dependent tensor quantity. Improvement of these models has been considered as a part of the ‘grand challenge’ that will lead to better understanding of the geodynamo physical phenomena under realistic conditions. Our motives in the present research are to contribute to better modelling of the turbulent transport in presence of the Lorentz force and of the mutual interactions between the velocity and electromagnetic fields in a turbulent regime, as well as to perform numerical simulations of the realistic working conditions. We focused on mimicking the Riga dynamo experimental set-up which provided the first ever experimental proof of the dynamo action in a turbulent regime, but the approach can be applied to any other configuration. Our approach is based on solving ensemble averaged transport equations where the single-point subscale turbulence closure is

New Journal of Physics 9 (2007) 306 (http://www.njp.org/)

4


employed for the unresolved turbulence contributions, Kenjereˇs and Hanjalić [25]. This is the first application of such an approach to turbulent magnetic dynamos.3

2. Numerical model

In our recent publication, Kenjereˇs et al [27] we presented the first attempt to simulate the Riga dynamo experimental set-up under realistic working conditions. There we applied a segregated numerical approach for coupling the solvers for simulations of the fluid-flow and electromagnetic fields. The electromagnetic fields were computed using the Adams–Bashforth finite-difference solver while the velocity and turbulence fields were modelled by a finite-volume Navier– Stokes solver using a Reynolds-averaged Navier–Stokes (RANS) turbulence model. Compared to the previously proposed simple one-dimensional back-reaction model used by Gailitis et al [12, 14], the segregated coupling resulted in improved predictions of the growth rates of the selfexcited magnetic field, while the predicted frequencies were just slightly different. In addition to the improved growth rates, the major advancement on the previous one-dimensional backreaction model was in predicting the radial reorganizations of the axial velocity profiles. Despite these improvements, the segregated solvers coupling also showed some deficiencies. The twoway interactions were performed by essentially feeding into the fluid solver the estimated magnitudes of the Lorentz force obtained from the electromagnetic kinematic solver. The intensity of the Lorentz force was evaluated from the experimentally determined Ohmic losses (from the measured increase in the power consumption) for different propeller rotation rates. The recalculated velocity fields are then submitted back to the electromagnetic solver where these differences in the velocity fields are analyzed. Essentially, in such an approach, only a singleiteration in exchanging the mutual coupling of the variables is performed. On the other hand, from the numerical point of view, the exchange of fields between two conceptually quite different numerical solvers (Maxwell’s solver based on a finite difference discretization for rectangular geometries and a Navier–Stokes solver based on the finite-volume approach for a nonorthogonal numerical mesh) required the introduction of a multi-dimensional interpolation/extrapolating techniques in order to preserve the smoothness of the fields with the divergence-free conditions. In order to perform fully coupled fluid flow/electromagnetic interactions without introducing any presumptions of the Ohmic losses, we decided to proceed with the development of a fully integrated finite-volume based Navier–Stokes/Maxwell solver for three-dimensional nonorthogonal geometries. The ability of the finite-volume-based numerical solver to preserve conservation of all variables is considered as the main advantage compared to the finite-difference or the finite-element based numerical solvers. These features are of high importance when dealing with multi-physics multi-scale phenomena with an imperative on preserving the divergencefree fields. Since the Riga dynamo experimental set-up is characterized by a high value of the hydrodynamic Reynolds number (Re = UD/ν ≈ 3.5 × 106 ) and a relatively low value of the magnetic Reynolds number (Rem ≈ 18), before pursuing such simulations of the fully coupled phenomena, the crucial prerequisite is to have reliable and accurate simulations of highly turbulent flow regimes in complex geometries. The full DNS, resolving all turbulence scales in space (ranging from the largest one defined by characteristic geometrical dimensions to viscous dissipative ones) and time (with a time step smaller than the typical lifetime of the smallest 3

An excerpt of the hybrid T-RANS/direct numerical solutions (DNS) results presented here can be found in Kenjereˇs and Hanjalić [26].


5


eddies before being destroyed by viscosity) is presently out of reach of any computer since it requires a grid refinement proportional to Re3/4 per coordinate direction. The second option is to numerically resolve large eddies and model the effect of all turbulence scales smaller than the underlying numerical mesh (that serves as a spatial filter) using a subgrid-scale turbulence model. Evidently, in such a large eddy simulation (LES) approach, the numerical grid resolution is less demanding than in the DNS method. The problem with the LES approach is in the extreme demand on grid refinement in the near-wall flow regions in order to resolve the thin wall boundary layers, requiring an almost identical grid resolution as for DNS (grid resolution savings are only in the regions further away from a wall), Pope [28]. The current LES of the magnetohydrodynamics (MHD) turbulence have been mainly performed for highly simplified situations such as for initially homogenous decaying turbulence subjected to an external magnetic field for an incompressible (Agullo et al [29]; Müller and Carati [30]; Knaepen and Moin [31]; Haugen and Brandenburg [32]) or a compressible (Chernyshov et al [33]) fluid. In order to overcome the current limitations, the numerical schemes based on enhanced constant viscosity or hyper-viscosity/hyper-diffusivity have been introduced. With an increase in the hyper-viscosity and/or hyper-diffusivity, the dissipative subrange of flow scales is significantly reduced (and so is, consequently, the numerical grid resolution), Biskamp and Müller [34]. Unfortunately, such a simple approach does not account for the local turbulence effects, which strongly influence the velocity field and its interactions with the electromagnetic field. In addition to this deficiency, a replacement of the normal diffusivity with a hyper-diffusivity causes the saturation of the dynamo-generated magnetic fields to occur at significantly higher magnetic Re numbers when compared to the referent normal cases, as shown in Brandenburg and Sarson [35]. The more advanced subgrid models have been introduced recently in the works of Buffett [36] and Matsushima [37] where the scale-similarity based models for heat flux are investigated, but again, for highly simplified geometries (a periodic box). A novel approach has been published by Ponty et al [38] where a mixed numerical scheme is proposed. The magnetic induction was solved on a fully resolved numerical mesh and an LES was applied for the velocity fields. In their more recent study Ponty et al [39] extended the previous method with a Lagrangian-averaged model for velocity fields. Although a relatively simple configuration with a Taylor–Green vortex was simulated and the kinematic Reynolds number was relatively low (Re 1.5 × 104 ), potentials of such a mixed numerical approach were clearly demonstrated. In this study, we propose a hybrid approach involving simultaneous solving of the twoway coupled fluid flow and electromagnetic interactions. For the fluid flow variables (velocity, pressure and turbulence) we introduce the transient RANS (T-RANS) method, whereas a fully resolving (DNS) approach was used for the electromagnetic variables. The justification of a such an approach is motivated by a huge disparity of characteristic length and timescales for velocity and electromagnetic fields, respectively. For working conditions of the Riga dynamo set-up it can 1/4 be estimated that the typical magnetic diffusive length scale, defined as ηB = λ3 /ε , is much 1/4 larger than the typical viscous (Kolmogorov) velocity scale, ηu = ν3 /ε since the magnetic Prandtl number estimated from the liquid sodium properties gives Prm = ν/λ = (ηu /ηB )4/3 ≈ 6.5 × 10−6 . In order to mimic the realistic experimental conditions, the fully developed steady RANS solutions are first performed. In this stage, the magnetic induction equations were not solved. If the velocity fields and turbulence levels are properly captured, such calculated fields will represent a proper basis for possible capturing of the self-generation of a magnetic field New Journal of Physics 9 (2007) 306 (http://www.njp.org/)

6


when critical parameters (saturation levels) are reached. When fully convergent solutions are obtained, the magnetic induction equation is activated. For each time step a series of iterations are performed until the convergent fields are finally obtained. In contrast to the previous segregated numerical simulations, this iterative procedure involves the simultaneous solution of both the momentum and magnetic induction equations with implicitly updated (the most recent) fields. This procedure is advanced in time reproducing the growth of a self-generated magnetic field. The generated magnetic field creates the Lorentz force, which feeds back into the momentum equation and, ultimately, the saturation magnetic regime is achieved. As it will be shown below, the T-RANS approach applied proved to be capable of capturing large-scale instabilities and the consequent dynamo action manifested in magnetic and velocity field oscillatory growth and saturation. However, it should not be forgotten that T-RANS excludes the possibility for recovering any spectral information about the stochastic turbulence. Its effect on the fluid flow and dynamo is fully accounted for by subscale models of the ensemble-averaged turbulence quantities. 2.1. Equations and subscale model of turbulence The governing fluid momentum and magnetic induction equations that describe the two-way coupled fluid flow/electromagnetic interactions can be written as: 1 ∂ ∂Uˆ i ∂Uˆ j 1 ∂Pˆ ∂Bˆ i ∂Bˆ j ∂Uˆ i , (1) + + − Bˆ j = Bˆ j ν − Uˆ i Uˆ j − ∂t ∂xj ∂xj ∂xi ρ ∂xi ρµ0 ∂xj ∂xi

ˆ × Bˆ FL = Jˆ × Bˆ = 1/µ0 ∇ × B ∂Bˆ i ∂ = ∂t ∂xj

1 ∂Bˆ i − Bˆ i Uˆ j + Uˆ i Bˆ j . µ0 σ ∂xj

(2)

The equation set is closed with the divergency free conditions ∂Uˆ j /∂xj = 0, ∂Bˆ j /∂xj = 0 for the velocity and magnetic fields, respectively. It can be seen that these equations are directly interconnected through the Lorentz force in the momentum equation and the convective term in the magnetic induction equation. Since the fully resolving approach will be applied for the magnetic induction, the original form of these equations for the three coordinate directions are discretized. For the momentum equations, in order to obtain the evolution equations for the ensemble-mean variables, the Reynolds decomposition is introduced. All instantaneous variables are represented as a sum of the ensemble averaged and fluctuating contributions, i.e. Uˆ i = Ui + ui = Ui + ui , Pˆ = P + p = P + p, Bˆ i = Bi + bi = Bi + bi , etc. The ensembleaveraged momentum equations can then be written as:

u ∂ 1 ∂P ∂Bi ∂Bj 1 ∂Ui ∂Uj ∂Ui b = + + − Bj ν − Ui Uj − τij − τij − Bj , (3) ∂t ∂xj ∂xj ∂xi ρ ∂xi ρµ0 ∂xj ∂xi where τiju = ui uj , τijb = 1/ρµ0 bi bj are the ensemble averaged second moments of the velocity and magnetic field fluctuations, respectively.4 In order to close the system of the 4

An ensemble averaging of equation (2) would lead to the appearance of a cross-correlation—the turbulent magnetic flux ui bj , but in view of the large-scale separation between the velocity and magnetic field fluctuations this term can be neglected, permitting thus to treat equation (2) in its original form in the DNS fashion. New Journal of Physics 9 (2007) 306 (http://www.njp.org/)

7


ensemble-averaged equations, additional relations (turbulence closure models) are needed to account for the subscale turbulence contributions. The most natural way to obtain these correlations is to derive the full transport equations for ui uj , bi bj , ui bj , etc. Such a system of equations becomes very complex with many terms that should again be approximated (modelled) involving even higher order correlations (triple moments, correlations involving fluctuating pressure, derivatives of the fluctuating velocity and magnetic field, and others). The difficulty is not only in the number of equations but also in our present inability to model and to evaluate such correlations since that process requires simultaneous measurement of the fluctuating velocity and electromagnetic fields. Such details can only be extracted from DNS. Unfortunately, at present, such DNS studies of the induction and the dynamo action at low Prm are still in a rudimentary phase and oriented primarily towards providing magnetic and kinematic spectra, Ponty et al [39, 40], Mininni et al [41]. An interesting approach is followed in Yoshizawa and Hamba [42] where an ensemble-averaged model of MHD turbulence is proposed based on the two-scale directinteraction approximation (TSDIA) statistical theory. The final form of the model, refined in Yoshizawa [42, 43] and Hamba [44, 45], is the eddy-viscosity based model with four equations for the bulk statistical quantities. These include transport equations for the turbulent kinetic energy (TKE) (k), its dissipation rate (ε), the turbulent residual helicity (H) and the turbulent crosshelicity (W ).5 The main difficulty in this approach is the estimation of the model coefficients due to the already mentioned lack of the detailed DNS studies. Recently, Hamba [46] proposed two different sets of coefficients and performed the two-dimensional simulations of an MHD dynamo in a rotating spherical shell. Despite the fact that a further optimization of the model terms and coefficients needs to be done, this approach demonstrated the potentials of the ensemble-averaged approach in simulating the dynamo phenomena. In our earlier research we developed several RANS models at different levels—including a full-second moment closure (Kenjereˇs et al [46]) and a simplified eddy-viscosity based model (Kenjereˇs and Hanjalić [47, 48]; Hanjalić and Kenjereˇs [49]–[52])—for turbulent flows subjected to an external magnetic field. These closures were developed for the subscale τiju turbulent stresses. In the second-moment closure (SMC), the full transport equations for all components of the turbulent stress are solved. This enables the capturing of the mechanism of the fluctuation distributions among the particular velocity components so that the effects of turbulence anisotropy can be properly modelled. The turbulence anisotropy is of crucial importance for accurate predictions of flows in the proximity of solid walls and when subjected to external body forces. For the Riga dynamo experimental set-up with a strongly turbulent helical motion in the inner cylinder (high Reynolds number), the near-wall effects are of secondary importance, provided the wall shear stress and its effect on flow and turbulence properties is well accounted for. Additionally, since it is our first attempt to perform fully coupled two-way magnetic dynamo simulations in realistic working conditions (geometry, flow, turbulence and electromagnetic parameters), we will proceed with a simplified eddy viscosity model. Now, instead of solving the full set of the second-moment equations, the transport equations for the turbulence kinetic energy (k = 0.5ui ui ) and its dissipation rate (ε = 2ν(∂ui /∂xj )2 ) are solved: ∂ ∂Ui ∂k ∂k t = + Dk − Uj k − τiju − ε − Pkb , (4) ν ∂t ∂xj ∂xj ∂xj 5

These variables are defined as: H = b · j − u · ∇ × u , W = u · b .


8


Table 1. Specification of the model coefficients. Cε1

Cε2

CL

σk

σε

Cµ

1.44

1.92

0.025

1.0

1.3

0.09

∂ ∂ε = ∂t ∂xj

∂Ui ε ε2 ∂ε t − Cε2 − Pεb , + Dε − Uj ε − Cε1 τiju ν ∂xj ∂xj k k

(5)

where 1 νt ∂k Dkt = − ui p − uj k = , ρ σk ∂xj 2 ∂p ∂u ν ∂u νt ∂ ε j i Dεt = −2 , − ν uj = ρ ∂xk ∂xk ∂xk σε ∂xj σ 2 k ε σ 2 b , Pεb = Pkb , Pk = B0 k exp −CL B0 ρ ρ ε k

(6)

are turbulent diffusion (Dkt , Dεt ) modelled by a simple gradient diffusion hypothesis, whereas the ‘magnetic’ production/destruction terms (Pkb , Pεb ) are modelled by introducing the locally determined turbulence parameters and the timescale of the magnetic damping (Kenjereˇs and Hanjalić [47]), respectively. Here B0 = Bi2 is the intensity of the magnetic field. The kinematic turbulent stresses are then evaluated as: k2 2 ∂Ui ∂Uj u . (7) + , νt = Cµ τij = kδij − νt 3 ∂xj ∂xi ε The specification of the model coefficients is given in table 1. Note that with this model, both the direct (through the momentum equations) and the indirect (through the MHD terms in the turbulence equations, Pkb , Pεb ) effects of the Lorentz force on the fluid flow and turbulence are introduced. This model was extensively validated in a series of generic situations including fully developed turbulent channel flows (Kenjereˇs et al [46], Kenjereˇs and Hanjalić [47]), turbulent thermal convection in enclosures (Kenjereˇs and Hanjalić [25, 48, 53], Hanjalić and Kenjereˇs [49, 50, 52]), unsteady electromagnetically driven recirculating melt flows (Schwarze and Obermeier [53]) and control of the crossing-shock turbulent layer in magnetogasdynamic applications (Gaitonde and Miller [55]). For all the above mentioned applications, the fluid flow was subjected entirely or locally to external magnetic fields of different orientations and strength, i.e. a wide range of flow regimes (weakly or highly turbulent flows with or without heat transfer), and intensities of electromagnetic interactions (weak or strong interactions) are simulated. This variety of different flow situations and applications gives reasonable confidence in the validity of the above introduced model for simulating turbulent flows subjected to the Lorentz force. The remaining magnetic turbulent stress τijb is negligible for the particular Riga dynamo set-up. The more general derivation of the magnetic turbulent contributions can be performed by following a similar approach as for the SMC level for velocity fluctuations, Kenjereˇs et al [46], or a simplified four-equation eddy-viscosity approach can be considered as suggested by Hamba [45]. This is left for future work. New Journal of Physics 9 (2007) 306 (http://www.njp.org/)

9


2.2. The finite-volume discretization of equations It is recalled that the finite-volume approach, used here for the simultaneous solving of fluid flow, turbulence and electromagnetic fields, is based on the pre-integration of the field equations over an elementary control volume, thus ensuring the conservation of all variables over each grid cell. The transport equations for mass, momentum, magnetic induction, TKE and its dissipation rate, equations (2)–(5), can thus all be written in a general form as t+ t t+ t ∂ ∂ dV dt = − Uj nj dS + S dV dt, (8) ∂xj t t V ∂t S V where the volume integral of diffusive and convective terms (first term on the rhs) is replaced by a sum of the surface integrals over the cell-faces of a control volume. Since the neighbouring control volumes share identical cell-faces, strong conservation properties of the transport variables are preserved. The control volumes are in forms of hexahedrals. All transport variables are located in the geometrical centre of such control volumes (collocated variable arrangement). The Cartesian vector and tensor components are used for representing a general nonorthogonal numerical mesh. The grid-nonorthogonality is included through the local curvilinear coordinate system j representing the cofactor of ∂xi /∂xj coordinate transformation, (xj ) with the Jacobian (J) and βm ∂ 1 ∂ j ≈ · · · β . This together with the fact that the product J x1 x2 x3 represents the m j ∂xi J ∂x exact total volume ( V ) of the control cell around the central position (P), makes it possible that all these coordinate transformation can be easily represented in the Cartesian coordinate system. In order to prevent a decoupling between the velocity and pressure fields (checkerboard pressure oscillations) the Rhie–Chow interpolation is used in the pressure-correction equation. The corrected velocity and pressure fields are iteratively calculated by the semi-implicit method for pressure-linked equations (SIMPLE) algorithm. The numerical accuracy of the entire discretized system is of the second-order. The time integration is performed by the fully implicit secondorder scheme based on three consecutive time steps. The diffusive terms (first part of the first term on the rhs in equation (8)) are discretized by the second-order central differencing scheme. The remaining convective terms (second part of the first term on the rhs in equation (8)) are calculated by the monotonicity preserving total variation diminishing scheme with the UMIST limiter. Under the assumption that the source terms (S ) are uniformly distributed over the whole control volume, they are simply integrated as: P S dV = S,P V. (9) S = V

By putting all parts together, the final form of the discretized system suitable for an iterative solution procedure can be written as: 1 AP new P old = Asur sur + S,P + AP P −1 , (10) αφ αφ sur where the old P represents the value from the previous iteration and α is an under-relaxation parameter. The sum of the discretized set of coefficients, AP = sur Asur includes both the diffusive and convective contributions from all control volumes surrounding (sur ) the central New Journal of Physics 9 (2007) 306 (http://www.njp.org/)

10


node (P). The linearized system of equations is then solved using Stone’s strongly implicit procedure (SIP) based on an incomplete lower-upper (LU)-factorization. The finite-volume solver presented here can be applied for any three-dimensional nonorthogonal geometry distributed over many multi-blocks (locally structured). The local grid refinements can be activated in the pre-specified blocks. The solver can be executed in serial or parallel mode. The latter involves the domain decomposition based message passing interface (MPI) directives. In practical terms, this domain decomposition is efficiently distributed over a separate multi-block that can be part of different domains, e.g. which belong to flow, walls or external electromagnetic layers. 2.3. Specification of boundary conditions The boundary conditions for all fluid flow and turbulence parameters along the solid walls are imposed through the set of wall-functions. This is the only viable alternative for high Reynolds number flows in complex domains because a full numerical resolving of boundary layers along the solid walls at such high Re numbers would require modification of all equations for the viscous and non-viscous wall-blocking effects and the application of a much finer computational grid, which is at present beyond the reach of present simulations. We note, however, that the viscous wall effects play a secondary role in the Riga dynamo set-up since the highly turbulent swirling flow pattern inside the inner cylinder determines to a large degree the flow features in the whole rig. Nevertheless, in order to properly capture the high gradients of velocity in the near-wall regions, the numerical mesh is always clustered in these regions. The numerical mesh is so designed to ensure that the first near-wall grid point is always located at the non-dimensional distance (normalized with the inner-wall scales) in the range 40 r + (x+ , z+ ) 100. Instead of simulating the flow through the rotating propeller, for which a complex moving mesh would be needed, we mimic the propeller by imposing both the axial and tangential velocity projections in the inner cylinder at the propeller exit plane. Since the direct velocity measurements for the realistic working conditions are not available, we scaled-up velocity components measured in the 1:2 scaled-down water model of the Riga dynamo set-up, Gailitis et al [9, 10]. These initial velocity distributions are allowed to be dynamically self-adapting under the influence of the self-generated Lorentz force. This enabled us to perform a single simulation for a fully coupled system, instead of a series of simulations for different values of the rotation rates. This resulted in significant reductions of the computational time. Since the experimental data for the fully coupled phenomena are reported for different propeller rotation rates, the rotational rates for particular time instants are estimated from the instantaneous velocity components at the propeller location. In contrast to momentum and turbulence parameters equations, the magnetic induction equation is solved in the inner and outer cylinders, in the dividing walls and in the surrounding ring of the sodium at rest. Additional difficulty appeared in the specification of the nonlocal boundary conditions for the electromagnetic variables. In the previous work of Gailitis et al [14], where the two-dimensional finite-difference based kinematic solver was tested for a simplified Riga dynamo set-up, the Laplace equation was solved in the vacuum exterior and matching conditions are imposed at the outer walls of the surrounding ring of sodium at rest. Such an approach for three-dimensional geometries significantly reduce the numerical effectiveness of the solver. This is the main reason why we followed a simplified approach by which we imposed the vertical magnetic field boundary condition at the outer wall of the surrounding ring of sodium at rest in New Journal of Physics 9 (2007) 306 (http://www.njp.org/)

11


accordance with Brandenburg et al [56, 57] and Hamba [45]. This boundary condition allows an escape of the self-generated magnetic field from the discretized domain while still keeping reasonable critical dynamo thresholds (a difference of 20% between the theoretically determined and numerically estimated critical thresholds in the Riga dynamo set-up, Gailitis et al [14]). In this study we performed additional extensive validation of this boundary condition in the kinematic mode (without back-reaction of the self-generated magnetic field on the fluid flow) and proved that the critical dynamo threshold is reached for Rem > 15. This value is in good agreement with the theoretically estimated (based on studies of the convective instabilities of the Ponomarenko dynamo, Gailitis et al [11, 14]) critical threshold of Recm = 17.7. In the work of Zuniga et al [58] an assessment of different electromagnetic boundary conditions on the onset of a magnetic dynamo action is performed. Both the Riga dynamo and the Karlsruhe-dynamo experimental set-ups were analyzed in detail. It was concluded that an inclusion of a stagnant surrounding layer always resulted in reducing the critical threshold, what appeared to be very convenient from the practical point of a view for the realization of a magnetic field self-excitation. Obviously, as long as this layer is included in the numerical simulations, simplified boundary conditions for electromagnetic variables (the vertical field conditions) will be a reasonable first approximation. Finally, it is concluded that for the first series of fully-coupled two-way simulations targeting numerical reproduction of the self-generation of a magnetic field in the Riga dynamo set-up, the vertical magnetic field condition was found sufficient. For future studies a novel approach based on a combined integro-differential formulation and mixed boundary element/finite volume approach can be used, Iskakov et al [59].

3. Results and discussion of the numerical simulations

The sketch of the Riga dynamo set-up is shown in figure 1 (right). It consists of the inner cylinder where a strong helical motion is generated by a propeller (driven by two electric motors with a power up to 100 kW each), the outer-cylinder with a back-flow and surrounding ring filled with sodium at rest. Although this set-up was not designed to actually mimic any real geodynamo situations, many similarities in its conceptual design can be identified in a simplified mechanism of the multi-columnar convective vortex patterns inside of Earth, Busse [7], as shown in figure 1 (left). In our simulations, the geometry of the Riga dynamo set-up is divided in structured multi-block sub-domains consisting of hexahedral control volumes, figure 2. Because the two-way coupled simulations were to be performed, the entire configuration was discretized including—next to the domains occupied with moving fluid (inner and outer cylinders), also both walls between the compartments and the surrounding ring of sodium at rest. Then different multi-blocks (indicated by different colours in figure 2) are grouped together in specific zones in such a way that an optimal load balancing is achieved when the numerical code is executed in the parallel mode. In total, 88 structured multi-blocks are used with approx. 4 × 106 cell control volumes. The numerical mesh is refined in the near-wall regions and in the regions where the flow streamlines exhibit strong curvature. This is necessary in order to properly resolve the velocity gradients (that in turn return accurately the turbulence energy production) in these regions and to ensure that the non-dimensional wall-distance of the wall-nearest to the cell centres lies in the range (40 < x+ , z+ < 100) considered as optimal for the applications of wall-functions.


12


Figure 1. Left panel: sketch of the columnar vortex pattern of convection

inside Earth according to Busse [7]. Right panel: Riga dynamo experimental set-up, Gailitis et al [9, 10]: 1—propeller, 2—inner cylinder with helical flow, 3—outer passage with back flow, 4—surrounding ring of sodium at rest, 5—thermal insulation. Despite its relative simplicity, the Riga dynamo set-up exhibits a few quite challenging flow features that should be properly captured by a fluid flow solver and its RANS subscale turbulence model. Since the dynamics of the imposed helical fluid flow in the inner cylinder is of crucial importance for reaching critical conditions for magnetic field self-generation, its accurate reproduction was one of the simulations’ primary goals. The second challenge is a strong, 180◦ bend in the flow direction from the inner to the outer cylinder with possible flow separation and subsequent re-attachment along the inner walls after the bend. In order to compare predictive performances of the introduced numerical model for the fluid flow and turbulence, a steady RANS simulation of a pure hydrodynamical flow in a 1:2 scaled-down experimental set-up with water at Re ≈ 106 were first performed. The comparison of the experimental and numerical results for the axial and tangential velocity components at different locations in the inner cylinder are shown in figure 3. It can be seen that the eddy-viscosity RANS subscale closure model reproduces well the experimental profiles at all locations. The characteristic peaks of the velocity components—in the centre of the inner cylinder for the axial and in the proximity of the wall for the tangential velocity, are well captured. This satisfactory reproduction of the water experiment may sound surprising in view of the fact that the applied and similar linear eddy-viscosity models have been notoriously known as inadequate for reproducing accurately confined swirling flows, because of their tendency to return a solid-body-rotation irrespective of the manner in which the swirl was generated, Jakirlić et al [61], Kitoh [62], Speziale et al [63]. However, in the present experiment, the helical flow is designed on purpose to produce tangential velocity profiles similar to the solid-body-rotation in order to achieve the most suitable conditions for the self-excitation under limited power resources, Gailitis et al [9]–[13]. This explains why the applied eddy-viscosity RANS models New Journal of Physics 9 (2007) 306 (http://www.njp.org/)

13


Figure 2. Multi-block structure of the numerical grid used for simulations of

the Riga dynamo experimental set-up. Please note that the magnetic induction equation is solved in the entire domain while fluid flow equations are solved only in the inner cylinder and outer loop parts. Different colours indicate specific zones grouped together in order to provide optimal load balancing for parallel execution of simulations. (EVM) returned results in close agreement with the experiments. In order to double-check these simulations, we also applied a SMC model. Because an SMC implies the solution of separate transport equations for each component of the turbulent Reynolds stress (τiju ), it has been proved to reproduce much better the rotating and swirling flows than the eddy viscosity models, Jakirlić et al [61], Speziale [64]. Figure 3 shows that the SMC model returned some improvements with respect to the intensity of the decay of the axial velocity intensity along the central axis (r = 0) for the central and lower locations. On the other hand, considering all additional computational requirements and numerical stiffness associated with the application of an SMC approach, we decided to keep the EVM, which proved to provide sufficiently accurate velocity and turbulence parameters that can be used for the full-scale two-way coupled simulations. It is worth mentioning that in order to achieve good agreement with experiments, it was necessary to use the higher order bounded convective schemes total variation diminishing (TVD) in transport equations for momentum and turbulent parameters. With upwind based convective schemes significant under-predictions of the tangential velocity profiles are observed. After completing the simulations of the 1:2 scaled-down set-up with water, the RANS simulations were performed for the full scale experimental set-up with sodium as a working fluid at Re ≈ 3.5 × 106 . When a fully convergent steady state is achieved, the magnetic induction is activated and we proceeded with the simultaneous coupled solutions. The weak seed magnetic field (BO ≈ 10−5 T) satisfying the divergency free condition is initiated at the first time step. The time step was specified to be 1/100 of the typical experimentally recorded period of around 1 s New Journal of Physics 9 (2007) 306 (http://www.njp.org/)

14

DEUTSCHE PHYSIKALISCHE GESELLSCHAFT 10 Experiment EVM SMC

10

Experiment EVM SMC

5

5 Upper

Upper 0

Uy

Ut

5

5

Central

Central 0

5

5 Lower 0

0

0.025

0.05

Lower 0 0

r

0.025

0.05

r

Figure 3. Comparison of experimental (LDA) (experiments performed at the

Dresden University of Technology, Germany, Christen et al [60]) and numerically obtained (EVM, SMC) velocity profiles (left panel: axial, right panel: tangential) at characteristic vertical locations in the inner cylinder of the 1:2 scaled-down water model of the Riga dynamo set-up (Uy , Ut in (m s−1 ) and r in (m)). of the axial magnetic field component. This choice is based on the fact that we do not solve the full turbulence spectrum but use a single-point RANS closure that recognizes only a single time scale expressed as τ = k/ε. It is, however, the system time scale (period of the magnetic oscillations) and not the turbulence scale that is relevant for choosing the time step. The results of numerically obtained frequencies of the axial magnetic field as a function of the propeller rotation rates, presented in figure 4, show good agreement with the experimental record. In the kinematic regime there is no back-reaction of the self-generated magnetic field onto the fluid flow, and an exponential growth or decay of the seed magnetic field takes place. After this initial growth of the self-generated magnetic field in the kinematic regime, the Lorentz force grows and soon becomes strong enough to significantly affect the underlying fluid flow and turbulence. Since an increase in the Lorentz force imposes a stronger braking of fluid motion, it essentially reduces the source of self-generation of the magnetic field. Eventually, in the saturation regime a balance is reached, characterized by a zero-growth rate of magnetic field intensity. A snapshot of the velocity and magnetic field obtained by such a two-way coupled simulation when the saturation regime is reached is shown in figure 5. In order to portray spatial distribution of the velocity field, a series of massless particle trajectories are connected in space (stream-traces) and coloured with the magnitude of the axial velocity component, figure 5 (left). The two distinct regions can be observed— the core of the inner cylinder with highly helical flow and the outer-cylinder with a weak tangential component. The magnetic flux lines coloured with the magnitude of the axial magnetic field component show characteristic double-helix spatial New Journal of Physics 9 (2007) 306 (http://www.njp.org/)

15 Frequency (1–1 s)


1.5

Experimental kinematic Experimental saturation Simulations

1.0 Kinematic

0.5 1500

Rotation rate

Saturation

2000

(1min–1)

Figure 4. The rotation rate dependency of the self-induced magnetic field

frequency in the Riga dynamo set-up for kinematic (one-way) and saturation (two-way coupled) regimes. Experimental data from Gailitis et al [14].

Figure 5. Typical swirling flow pattern inside of the Riga dynamo set-up: streamlines coloured by axial velocity (in (m s−1 )) (left), and visualization of the self-induced magnetic field with magnetic flux lines coloured by the axial magnetic field (in (T)) magnitude (right)—a snapshot.

distribution, figure 5 (right). While the flow pattern is confined in the inner and outer cylinders, the self-generated magnetic field is allowed to leave the discretized domain, as indicated by the open magnetic flux lines (for this particular snapshot in the upper part of the set-up). The dynamic response of the axial and tangential velocity profiles at specific locations in the upper and lower segments of the set-up are shown in figure 6. The axial velocity profiles are significantly reduced at all locations. It can be seen that after initial adjustments, for the time instants after cycle 20, only mild changes can be observed.6 Since the axial Lorentz force is

Note that the cycle time instants actually correspond to the following time instants: cycle 5 = 1 s, cycle 10 = 2 s, cycle 15 = 2.85 s, cycle 20 = 3.5 s, cycle 30 = 5.5 s, cycle 40 = 7.5 s, cycle 50 = 9.5 s, respectively. 6


16


Cycle 5 Cycle 20 Cycle 30 Cycle 40 Cycle 50

14

13

Uy

12

Lower

Upper

11

10 –0.125 –0.100 –0.075 – 0.05 –0.025

0

0.025

0.050

0.075

0.100

0.125

0.075

0.100

0.125

r 10 Cycle 5 Cycle 20 Cycle 30 Cycle 40 Cycle 50

7.5

Ut

5.0

Lower

Upper

2.5

0 – 0.125 – 0.100 – 0.075 – 0.050 – 0.025

0

0.025

0.050

r Figure 6. The dynamic response of the axial and tangential velocity profiles

at upper (y = −1 m) and lower (y = 1 m) vertical positions (Uy , Ut in (m s−1 ), r in (m)).


17


strongest in the central part of the inner cylinder, the axial velocity is strongly dampened there, figure 6 (top). The tangential velocity profiles in the upper part show almost no changes because of the strong propeller forcing at this location. In the lower part, there is a significant reduction of the tangential component, figure 6 (bottom). The dynamic response of the TKE profiles at the same locations as for the velocity fields discussed above shows a selective damping or enhancement of the velocity fluctuations, figure 7. It is noted that even for the strongest damping in a fully developed saturation regime, the profiles of the TKE show finite and in some flow regions significant values indicating a persistence of the turbulent character of the underlying fluid flow. In the inner cylinder, the profiles indicate that the Lorentz force suppresses the velocity fluctuations for all selected locations, figure 7 (top). In contrast to the inner cylinder, locations along the outer cylinder show both temporal and spatial dependency of suppression/enhancement intensity caused by the self-generated magnetic field and resulting Lorentz force. The upper part of the outer cylinder shows significantly higher average levels of the TKE compared with its lower part, figure 7 (bottom). This is caused by a strong curvature of streamlines (180◦ flow redirection from inner to outer cylinder) that produces a strong generation of the TKE. The temporal dependency shows a non-monotonic behaviour—the peak value shows an initial increase (cycle 20), then reduction (cycle 30) and finally, oscillatory behaviour is established. Similar temporal dependency is observed for profiles in the upper part of the outer cylinder. The main difference is that the finally established levels of the TKE clearly show an increase compared to its initial levels, figure 7 (bottom). As intermediate conclusions we would like to recall the significance of the proper modelling of turbulence in the Riga-dynamo experiment. Note that smooth streamlines shown in figure 5 do not indicate a laminar flow but just an instantaneous snapshot of the ensemble-averaged turbulent field. The flow fluctuations are not neglected—these subscale contributions are calculated from full transport differential equations for the TKE (k) and its dissipation rate (ε). The T-RANS approach presented here does not capture the fine flow structure oscillations directly—but only the large-scale (relatively slow) unsteadiness that does not interact directly with turbulence. The laminar solutions of pure hydrodynamical flow for such a high value of Re will produce erroneous solutions if the fully resolving numerical mesh is not applied (in which case we would have a DNS), and consequently, the structure of the underlying fluid flow will not be able to produce critical conditions for possible self-excitation of a magnetic field. The resulting self-generated and self-sustained magnetic field profiles are shown in figure 8. Here, the symbols represent the experimentally recorded maximum values at the particular radial locations for different propeller rotation rates —ranging from 1900 to 2500 rpm, Gailitis et al [15]. All profiles are rescaled with the maximum value along the radial direction. It can be seen that for both the axial and radial magnetic field component, the experiments show a slight deviation in the profile symmetry. In contrast to the experimental profiles, the numerical results show perfect symmetry indicating full convergence of the magnetic induction equation. Again, we plotted the time evolution for the numerical data, showing that after cycle 30 a fully saturated regime is achieved. It can be seen that very good agreement between the experimental and numerical recordings is achieved for the axial magnetic field component, figure 8 (top). The results for the radial magnetic field component show a larger scatter, but again, agreement can be considered as good, figure 8 (bottom). In contrast to the zero axial magnetic field at the outer boundary of the surrounding sodium at rest, the radial component correctly shows a final value. This boundary peak value is slightly under-predicted as compared with measurements, but shows clearly that the boundary conditions for the vertical magnetic field is properly imposed. New Journal of Physics 9 (2007) 306 (http://www.njp.org/)

18


In contrast to the monotonic time evolution of the axial magnetic field profiles, the radial profiles show an interesting variation in the radial location of the peak value. The peak value travels from the centre of the outer cylinder for earlier time instants (cycle 10–cycle 15) and finally settles at the middle radial distance in the inner cylinder, figure 8 (bottom). The time dependence of the axial magnetic field (By ) in the saturation regime, at the upper and lower probe, are shown in figure 9. The amplitudes of the magnetic oscillations are very well captured by numerical simulations (especially at the upper probe), while the frequency is slightly under-predicted. The entire time-series of the axial magnetic field self-generation are shown in figure 10. Here the monitoring points are distributed along the vertical lines that are located at different positions, as specified in table 2. The evolutions show highest levels of amplification in the middle part of the set-up (MON13). In addition to this vertical dependency, a strong radial variation of the recorded signals is clearly visible with the strongest amplification in proximity of the wall dividing the inner from the outer cylinder (in accordance with figure 8 (top)). A similar oscillatory behaviour is observed for the axial velocity profiles, figure 11. Starting from the statistically steady and convergent solutions, as the self-generated magnetic field grows, the Lorentz force increases and begins to influence the underlying fluid flow. In the fully saturated regime (for time instants >3 s) periods with intensive disturbances are observed. Such oscillatory behaviour of the velocity field illustrates that fully coupled two-way interactions between the fluid flow and electromagnetic field are simulated. It is interesting to note the selective response of the TKE in the inner and outer cylinder, figure 12. While the TKE is highly suppressed in the inner cylinder at all locations (figure 12 (top)), it is significantly enhanced in the outer cylinder (figure 12 (bottom)). The vertical snapshots of the magnetic fluxes and contours of the axial Lorentz force are shown in figure 13. In addition to the rotational magnetic field patterns observed in the previous figures, a significant magnetic field reorganization in the vertical direction can be observed here. The Lorentz force contours show dramatic non-uniform braking of fluid motion in both the inner and outer cylinder. Animations reveal strong outbursts of the time-dependent Lorentz force in this plane, indicating clearly the intensity of interactions between the fluid flow and self-generated magnetic field. The horizontal snapshots of the magnetic flux lines distributions with contours of the axial magnetic field are shown in figure 14. The strongest stretching of the magnetic flux lines corresponds to regions with maximal helicity (the inner cylinder). In contrast to the measurements, the numerical simulations can provide full threedimensional spatial distributions of the self-generated magnetic field. In figure 15 the isosurfaces of the axial magnetic field of opposite signs (red and blue) are shown together with stream-traces of the instantaneous velocity field (grey tubes). The strongest self-amplification takes place in the lower part of the set-up and then it moves upwards. This vertical upward shift is observed in experiments too and explained in terms of a simultaneous reduction of both the axial and tangential velocity components in the lower part of the inner cylinder, as shown in figure 6. Animations reveal a pattern of a slowly rotating asymmetric magnetic field, with a frequency of approx. 1 Hz. This is significantly different from the typical frequency of the propeller that is approx. 30 Hz. 4. Conclusions and future prospects

This paper reports on two-way-coupled hybrid numerical simulations of the full-scale Riga dynamo experiment, in which the experimental configuration, fluid properties and externally New Journal of Physics 9 (2007) 306 (http://www.njp.org/)

19


5

Cycle 5 Cycle 20 Cycle 30 Cycle 40 Cycle 50

4

3

Lower

TKE

Upper

2

1

0 –0.125

–0.1

–0.075

–0.05

–0.025

0

0.025

0.05

0.075

0.1

0.125

r Cycle 5 Cycle 20 Cycle 30 Cycle 40 Cycle 50

2.25

TKE

1.5

0.75 Upper

Lower

0

–0.075

–0.050

–0.025

0

0.025

0.050

0.075

r(+ –)0.125 Figure 7. The dynamic response of the TKE profiles at upper (y = −1 m) and lower (y = 1 m) vertical positions (TKE = 0.5(uu + vv + ww) in (m2 s−2 )).


20


1.50 Cycle 10 Cycle 12 Cycle 15 Cycle 20 Cycle 30 Cycle 40 Cycle 50

1.25

By/Bymax

1.0

0.75

0.50

0.25

0 –0.4

–0.2

0

0.2

0.4

r 1.50 Cycle 10 Cycle 12 Cycle 15 Cycle 20 Cycle 30 Cycle 40 Cycle 50

1.25

Br/Brmax

1.0

0.75

0.50

0.25

0 –0.4

–0.2

0

0.2

0.4

r Figure 8. Radial dependency of the non-dimensional magnetic field distributions:

above, axial; below, radial magnetic field component. The symbols represent experimental values obtained at the different rotational rates of the propeller, ranging from 1900 to 2500 rpm. The lines represent numerically obtained values for different time instants. Experimental data from Gailitis et al [15].


21


Upper probe Lower probe

0.075 Max: 0.054 0.050

Max: 0.044

0.025

By 0

–0.025

–0.050 6

7

8

9

Time Upper probe Lower probe

0.075 Max: 0.055

Max: 0.048

0.050

0.025

By 0

–0.025

–0.050 6

7

8

9

Time Figure 9. Time evolutions of the axial magnetic field component at the upper and lower probe locations with indications of maximum amplitudes: above, experimental recordings for the 2100 rpm propeller working regime; below, numerical simulations under identical conditions (By in [T ], time in (s)). Experimental data from Gailitis et al [15].


22 MON: By


0.1 0 –0.1

MON: By

0 0.1

MON: By

1

2

3

4

5

6

7

8

9

10

0 –0.1

0.1

1

2

3

4

5

6

7

8

9

10

0 –0.1

0.1

1

2

3

4

5

6

7

8

9

10

0 –0.1 1

2

11

12

11

12

MON06 MON07 MON08 MON09 MON10

Location: rest

0

12


Location: outer

0

11 MON16 MON17 MON18 MON19 MON20

Location: inner 2

0

MON: By


Location: inner 1

3

4

5

6

7

8

9

10

11

12

Time Figure 10. Time evolution of the axial magnetic field component (By ) at characteristic locations distributed along vertical lines in an inner cylinder (inner 1 and inner 2), outer passage (outer) and surrounding ring of sodium (rest), respectively. Table 2. Specification of the monitoring point locations. MON01–MON05 (x = 0.16; z = 0; y = −1, −0.5, 0, 0.5, 1) m MON06–MON10 (x = 0.25; z = 0; y = −1, −0.5, 0, 0.5, 1) m MON11–MON15 (x = 0.1; z = 0; y = −1, −0.5, 0, 0.5, 1) m MON16–MON20 (x = 0.05; z = 0; y = −1, −0.5, 0, 0.5, 1) m

imposed parameters have all been closely replicated, thus reproducing realistic working conditions with Re = 3.5 × 106 , Rem = 18. The simulations are performed by solving the timedependent (‘transient’) transport equations for the ensemble average fluid flow and turbulence variables (T-RANS) simultaneously with the DNS of the magnetic induction equations. The T-RANS equations have been closed using an earlier developed eddy-viscosity model for the subscale turbulence contributions in which the effects of the fluctuating Lorentz force have New Journal of Physics 9 (2007) 306 (http://www.njp.org/)

23


MON: Uy

12

Location: inner 1

10 0

MON: Uy


11

–4

2

1

4

3

5

6

7

8

9

10

–5 –6 1

2

12


Location: outer

0

11

3

4

5

6

7

8

9

10

11

12

Time Figure 11. Time evolutions of axial velocity components (Uy ) at particular

MON: TKE

monitoring locations. 4 2 0 0

MON: TKE


Location: inner 1

1

2

3

4

5

6

7

8

9

10

11

12

3 MON01 MON02 MON03 MON04 MON05

2 Location: outer

1 0 0

1

2

3

4

5

6

7

8

9

10

11

12

Time

Figure 12. Time evolutions of the TKE (TKE =k = 0.5u2i in (m 2 s−2 )) at

monitoring locations. been included. The fluid-flow and electromagnetic equations are solved using a fully integrated finite-volume Navier–Stokes/Maxwell solver for three-dimensional nonorthogonal geometries. The fluid-flow solver was first tested by simulating a realistic 1:2 scaled-down water replica of the actual Riga experiment for which a detailed experimental database exists. The magnetic induction solver was separately tested in the kinematic mode (one-way coupling, the velocity fields are imposed and kept constant) and the results compared with experimentally determined frequencies of the self-generated magnetic field as a function of the propeller rotation rates. In both cases, good agreement with the experimental results is obtained. The fully coupled (two-way) simulations closely reproduced the experimentally recorded self-excitation and the subsequent saturated self-sustenance of the generated magnetic field. New Journal of Physics 9 (2007) 306 (http://www.njp.org/)

24


Figure 13. The magnetic flux lines with the contours of the axial Lorentz force

component (in (N m−3 )) in the central vertical (x–y)-plane (z = 0 m) for different time instants: cycle 12, cycle 15, cycle 20, cycle 40, respectively.

Detailed comparisons revealed that the simulated features of the self-excitation and sustenance, including frequencies, amplitudes and spatial distributions of magnetic fields are all in good agreement with available experimental data. In addition to the numerical verification of the main experimental findings—the histogram of the magnetic field self-excitation, saturation and selfsustenance, the simulations provided data on full time and space dynamics of the fluid velocity, turbulence statistics and magnetic field, from which practically any information and deductions could be extracted—of course within the framework of the ensemble-averaged time dynamics. A small selection of such snapshots shown in figures 5 and 12–14 illustrate the point by revealing the self-organization of the vortical and magnetic structures and the double-helix as the basic dynamo mechanism in the configuration considered. The simulations presented here were primarily intended to demonstrate that we have succeeded in replicating numerically the action of a hydro-magnetic dynamo in a turbulent regime without any external excitation or artificial inputs. The configuration and working conditions of the Riga dynamo were chosen for two reasons: because of the availability of all information about the experiments and results, needed to mimic the real experiment, and because the setup appeared to be quite challenging and thus suitable to demonstrate that the method could be applied to any other complex dynamo configuration. For the reasons outlined above, we have confined our discussion to detecting only the basic flow and magnetic field dynamics. Of course, much more information about the actual dynamo physics can be deduced, but we will leave these for future work. Prior to a more extensive analysis of the physics details, further refinements of both the numerical and modelling aspects are desirable to clarify the small discrepancies between the measured and simulated frequencies. The first task should be to use a higher order model for turbulent stresses τiju = ui uj , such as the full SMC proposed by Kenjereˇs et al [46]. This should make it possible to properly capture the turbulence anisotropy, which is known to be enhanced by preferential damping of some of the turbulent stress components by the Lorentz force, depending on the local configuration of the magnetic field. Whilst improved New Journal of Physics 9 (2007) 306 (http://www.njp.org/)

25


Figure 14. The magnetic flux lines superimposed on the axial magnetic field

component By (contours in (T)) in the central horizontal plane (y = −0.5 m) for different time instants: cycle 10, cycle 12, cycle 15, cycle 20, cycle 30, cycle 40, respectively.


26


Figure 15. Time evolution of the axial magnetic field (By = 0.05 (red), −0.05

(blue) in (T)) with streak-lines (grey tubes): cycle 10, cycle 12, cycle 15, cycle 20, cycle 30, cycle 40, respectively.


27


turbulence modelling will probably have little effect on the overall dynamics, it is expected that the potential of the SMC in better capturing of the stress anisotropy and of the effect of the streamline curvature and swirl should improve the overall quality of the simulations, especially the frequency of the self-generated magnetic field, figure 9. Such simulations should also provide better understanding of the redistributive mechanism of the energy of velocity fluctuations among the particular stress components. The next step can be to further sensitize the T-RANS approach to higher-frequency instabilities caused by the generated Lorentz force, by applying a hybrid RANS/LES approach. Our recent test of a seamless hybrid method in simulation of a high Rayleigh number turbulent thermal convection showed that such an approach could serve as an intermediate route between LES and T-RANS techniques, suitable for high Re and Ra situations, Kenjereˇs and Hanjalić [65]. A major future step will be to move towards the less-constrained turbulent magnetic dynamos, such as the Maryland dynamo, Shew and Lathrop [66]. In that case, effects of both the thermal convection and rotation in addition to the Lorentz force should be included in the turbulence subscale closure. In this context, an additional but important aspect should be the development of closure models for the turbulent magnetic stress (τijb = bi bj ) and the magnetic flux (τijub = ui bj ), since in such configurations the Rem is significantly larger (Rem ≈ 600). Finally, numerical treatment can be improved by applying the more refined boundary conditions for the electromagnetic variables. A novel approach based on a combination of boundary element and finite volume approach of Iskakov et al [59] could be an interesting path to follow. Acknowledgments

The results reported in this paper emerged from the research project entitled: ‘MAGDYN: Magnetic Field Dynamos—Laboratory Studies based on the Riga Dynamo Facility’, sponsored by the European Commission contract HPRI-CT-2001-50027. In the framework of this research project, Dr F Stefani, Dr G Gerbeth and Professor A Gailitis are gratefully acknowledged for numerous discussions and for providing the experimental data. Mr S Renaudier is acknowledged for his involvement in generation of the numerical grid. The research of Dr S Kenjereˇs has been possible by a fellowship of the Royal Netherlands Academy of Arts and Sciences (KNAW). The supercomputing facilities have been provided by the National Computing Facilities Foundation (NCF) and sponsored by the Netherlands Organization for Scientific Research (NWO). References [1] Moffatt H K 1978 Magnetic Field Generation in Electrically Conducting Fluids (Cambridge: Cambridge University Press) [2] Moffatt H K 1989 Nature 341 285 [3] Rüdiger G and Hollerbach The Magnetic Universe: Geophysical and Astrophysical Dynamo Theory (Wienhem: Wiley) [4] Wielebinski R and Krause F 1993 Astron. Astrophys. Rev. 4 449 [5] Glatzmaier G A and Roberts P H 1995 Nature 337 203 [6] Roberts P H and Glatzmaier G A 2000 Rev. Mod. Phys. 72 1081 New Journal of Physics 9 (2007) 306 (http://www.njp.org/)

28 [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33] [34] [35] [36] [37] [38] [39] [40] [41] [42] [43] [44] [45] [46] [47] [48] [49] [50] [51]


Busse F H 2000 Annu. Rev. Fluid Mech. 32 383 Glatzmaier G A 2002 Ann. Rev. Earth Planet. Sci. 30 237 Gailitis A et al 2000 Phys. Rev. Lett. 84 4365 Gailitis A et al 2001 Phys. Rev. Lett. 86 3024 Gailitis A, Lielausis O, Platacis E, Gerbeth G and Stefani F 2002 Rev. Mod. Phys. 74 973 Gailitis A, Lielausis O, Platacis E, Gerbeth G and Stefani F 2002 Magnetohydrodynamics 38 5 Gaialitis A, Lielausis O, Platacis E, Gerberth G and Stefani F 2003 Surv. Geophys. 24 247 Gailitis A, Lielausis O, Platacis E, Gerbeth G and Stefani F 2004 Phys. Plasmas 11 2838 Gailitis A, Lielausis O, Platacis E, Gerbeth G and Stefani F 2006 Proc. 15th Riga and 6th PAMIR Conf. on Fundamental and Applied MHD (Riga, Latvia) p 15 Stieglitz R and Müller U 2001 Phys. Fluids 13 561 Müller U and Stieglitz R 2002 Nonlinear Processes Geophys. 9 165 Müller U, Stieglitz R and Horanyi S 2004 J. Fluid Mech. 498 31 Müller U, Stieglitz R and Horanyi S 2006 J. Fluid Mech. 552 419 Bourgoin M, Volk R, Plihon N, Augier P, Odier P and Pinton J F 2006 New J. Phys. 8 329 Monchaux R et al 2007 Phys. Rev. Lett. 98 044502 Berhanu M et al 2007 Europhys. Lett. 77 59001 Peffley N L, Cawthorne A B and Lathrop D P 2000 Phys. Rev. E 61 5287 Nataf H C, Allboussiere T, Brito D, Cardin P, Gagniere N, Jault D, Masson J P and Schmitt D 2006 Geophys. Astrophys. Fluid Dyn. 100 281 Kenjereˇs S and Hanjalić K 1999 Int. J. Heat Fluid Flow 20 329 Kenjereˇs S and Hanjalić K 2007 Phys. Rev. Lett. 98 104501 Kenjereˇs S, Hanjalić K, Renaudier S, Stefani F, Gerbeth G and Gailitis A 2006 Phys. Plasmas 13 122308 Pope S B 2004 New J. Phys. 6 35 Agullo O, Müller W C, Knaepen B and Carati D 2001 Phys. Plasmas 8 3502 Müller W C and Carati D 2002 Phys. Plasmas 9 824 Knaepen B and Moin P 2004 Phys. Fluids 16 1255 Haugen N E and Brandeburg A 2006 Phys. Fluids 18 075106 Chernyshov A A, Karelsky K V and Petrosyan A S 2006 Phys. Plasmas 13 104501 Biskamp D and Müller W C 2000 Phys. Plasmas 7 4889 Brandenburg A and Sarson G R 2002 Phys. Rev. Lett. 88 055003 Buffett B A 2003 Geophys. J. Int. 153 753 Matsushima M 2005 Phys. Earth Planet. Inter. 153 74 Ponty Y, Politano H and Pinton J F 2004 Phys. Rev. Lett 92 144503 Ponty Y, Mininni P D, Montgomery D C, Pinton J F, Politano H and Pouquet A 2005 Phys. Rev. Lett. 94 164502 Mininni P D, Ponty Y, Montgomery D C, Pinton J F, Politano H and Pouquet A 2005 Astrophys. J. 626 853 Yoshizawa A and Hamba F 1988 Phys. Fluids 31 2276 Yoshizawa A 1990 Phys. Fluids B 2 3064 Yoshizawa A 1998 Hydrodynamic and Magnetohydrodynamic Turbulent Flows: Modelling and Statistical Theory (Dordrecht: Kluwer) Hamba F 1992 Phys. Fluids A 4 441 Hamba F 2004 Phys. Plasmas 11 5316 Kenjereˇs S, Hanjalić K and Bal D 2004 Phys. Fluids 16 1229 Kenjereˇs S and Hanjalić K 2000 Int. J. Heat Fluid Flow 21 329 Kenjereˇs S and Hanjalić K 2004 Int. J. Heat Fluid Flow 25 559 Hanjalić K and Kenjereˇs S 2000 J. Turbul. 1 1 Hanjalić K and Kenjereˇs S 2001 Flow Turbul. Combust. 66 427 Hanjalić K and Kenjereˇs S 2002 Ann. NY Acad. Sci. 972 19


29 [52] [53] [54] [55] [56] [57] [58] [59] [60] [61] [62] [63] [64] [65] [66]


Hanjalić K and Kenjereˇs S 2006 J. Appl. Mech—Trans. ASME 73 430 Kenjereˇs S and Hanjalić K 2002 Phys. Rev. E 66 036307 Schwarze R and Obermeier F 2004 Modelling Simul. Mater. Sci. Eng. 12 985 Gaitonde D V and Miller J H 2005 AIAA J. 43 1064 Brandenburg A, Nordlund A, Stein R F and Torkelsson U 1995 Astrophys. J. 446 741 Brandenburg A, Jennings R L, Nordlund A, Rieutord M, Stein R F and Tuominen I 1996 J. Fluid Mech. 306 325 Avalos-Zuniga R, Plunian F and Gailitis A 2003 Phys. Rev. E 68 066307 Iskakov A B, Descombes S and Dormy E 2004 J. Comput. Phys. 197 540 Christen M, Hänel H and Will G 1998 Beiträge zu Fluidenergiemaschienen 4 ed W H Faragallah and G Grabow (Sulzbach: Faragallah-Verlag und Bildarchiv) Jakirlić S, Hanjalić K and Tropea C 2002 AIAA J. 40 1984 Kitoh O 1991 J. Fluid Mech. 225 445 Speziale C G, Younis B A and Berger S A 2000 J. Fluid Mech. 407 1 Speziale C G, Sarkar S and Gatski T B 1991 J. Fluid Mech. 227 245 Kenjereˇs S and Hanjalić K 2006 Int. J. Heat Fluid Flow 27 800 Shew W L and Lathrop D P 2005 Phys. Earth Planet. Inter. 153 136
