Direct numerical simulation of quasi-two-dimensional MHD turbulent shear flows
DDirect numerical simulation of quasi-two-dimensional MHD turbulent shear flows
Long Chen , Alban Poth´erat , Ming-Jiu Ni and Ren´e Moreau School of Engineering Science, University of Chinese Academy of Sciences,Beijing 101408, China, [email protected] Fluid and Complex Systems Research Centre,Coventry University, Coventry CV15FB, United Kingdom Universit´e de Grenoble, Laboratoire SIMAP, GroupeEPM, BP 75, 38402 Saint Martin d’H`eres, FranceFebruary 4, 2021
Abstract
High-resolution direct numerical simulations (DNS) have been performed to study the tur-bulent shear flow of an electrically conducting fluid in a cylindrical container. The flow isdriven by an annular azimuthal Lorentz force induced by the interaction between the radialelectric currents ( I ) injected through a large number of small electrodes placed in the bottomwall and the magnetic field ( B ) imposed in the axial direction. All the numerical parameters,including the geometry of the container, the value of the external electric currents and thestrength of the magnetic fields, are set to be in line with the experiment performed by Mes-sadek & Moreau (2002) (J. Fluid Mech. , 137-159). Firstly, maintaining the Hartmannlayers in the laminar regime, three dimensional simulations are carried out to reproduce someof the experimental observations, such as the global angular momentum and the velocity pro-files. In this regime, the variation laws of the wall shear stresses, the energy spectra and thevisualizations of the flow structures near the side wall indicate the presence of separation orturbulence within the side wall layers, even though the current injection electrodes are farfrom the side wall. Furthermore, a parametric analysis of the flow also reveals that the Ek-man recirculations have significant influence on the vortex size, the free shear layer, and theglobal dissipation. Second, we recover the scaling laws of the cutoff scale that separate thelarge quasi-two-dimensional scales from the small three-dimensional ones (Sommeria & Moreau(1982) J. Fluid Mech. , 507-518), and thus establish their validity in sheared MHD tur-bulence. Furthermore, we find that three-componentality are and the three-dimensionalityappear concurrently and that both the two-dimensional cutoff frequency and the mean energyassociated to the axial component of velocity scale with N t , respectively as 0 . N . t and0 . N − . t . In this paper, we apply the three-dimensional (3D) direct numerical simulations (DNS) to studythe flow of an electrically conducting incompressible fluid in a cylindrical container, popularlyknown as MATUR (MAgnetohydrodynamics TURbulence), designed to investigate the quasi-two-dimensional (Q2D) turbulence (Alboussi`ere et al. , 1999) in the presence of external magnetic fields(MHD). MATUR is shown in Fig. 1, where the magnetic field is applied along the axial direction.1 a r X i v : . [ phy s i c s . f l u - dyn ] F e b . Chen et al. DNS of MATUR For such a laboratory scale configuration, the magnetic Reynolds number is much smaller thanunity (here, R m = µ m σU L ≈ . µ m denotes the magnetic permeability of vacuum, σ is the electrical conductivity, U = 0 . L = 0 .
01 m are the typical characteristic velocityand the characteristic length scale, respectively), then the induced magnetic field b is much smallerthan the imposed one B ( b ∼ R m B (cid:28) B ), and the Lorentz force acting on the flow is obtainedas F Lorentz = j × B , where j denotes the electric current density (Roberts, 1967). If the magneticfield is strong enough, velocity variations along the field lines are damped by the Lorentz force andthe flow tends to be Q2D.This tendency was observed in several experiments. Kolesnikov & Tsinober (1974) and David-son (1997) stated that this evolution towards a quasi-two-dimensional regime was a consequenceof the invariance of the angular momentum component parallel to the magnetic field, when itsperpendicular components decay exponentially ( ∼ e σB t/ρ ). Eckert et al. (2001) conducted anexperiment on MHD turbulence in a sodium channel flow exposed to a transverse magnetic field,and the measured turbulence intensity and energy spectra were found to exhibit a spectral slopevarying with the magnetic interaction parameter ( N = σB L/ρU ) from a k − / law for N ≤ − N (cid:39) et al. (1999)and Poth´erat et al. (2000) carried out experimental and theoretical studies on the Q2D turbulentshear flows, where the transport of a scalar quantity and the free surface effect were considered,using the MATUR equipment. Messadek & Moreau (2002) further provided experiment data onthe MHD turbulent shear flows in a wide range of Ha and Re , and highlighted the importantrole of the Hartmann layers where the Joule effect and viscosity dissipate most of the kineticenergy. Recently, Stelzer et al. (2015 b ) built a new experimental device called ZUCCHINI (ZUrichCylindrical CHannel INstability Investigation), which, as MATUR, featured a free shear layer atthe edge of inner disk electrode. Combining it with finite element simulations (based on a 2Daxisymmetric model), they studied the instabilities of the free shear layer and identified severalflow regimes characterized by the nature of the instabilities of the Kelvin-Helmholtz type (Stelzer et al. , 2015 a ). Based on the FLOWCUBE platform, a more homogeneous type of turbulencebetween Hartmann walls was produced from the destabilisation of vortex arrays (Klein & Poth´erat(2010); Poth´erat & Klein (2014); Baker et al. (2018)). These authors focused on the transitionbetween 3D and Q2D turbulence. In particular, the cut-off length scale (cid:98) l c ⊥ ( ∼ N / t , where N t = σB aρU (cid:16) aL (cid:17) (1)is the true interaction parameter) first theorised by Sommeria & Moreau (1982), that separate 3Dfrom Q2D fluctuations were obtained experimentally, as well as the evidences of inverse and directenergy cascades in 3D magnetohydrodynamic turbulence.However, a major disadvantage of experimental approaches is that the liquid metal used fortheir high electrical conductivity is non-transparent. Although the velocity fields can be measuredby ultrasonic Doppler velocimetry or potential probe techniques, more complete information, e.gthe distribution of the flow fields and the electromagnetic quantities, are rather difficult to ob-tain. Therefore, numerical simulations, which can complement the experimental measurements,have been developed recently to study MHD turbulence. Taking advantage of the Q2D propertyof the MHD flows in case of high N and Ha , several simplified effective 2D models have been . Chen et al. DNS of MATUR developed by averaging the full Navier-Stokes equations along the direction of the magnetic fields.The advantages of using these 2D models are evident, not only to save the costs compared to afull 3D numerical approach, but also to provide accurate results where 3D numerical solutionscannot fully resolve the boundary layer in case of high Ha . Sommeria & Moreau (1982) deriveda two-dimensional model (denoted as SM82 hereafter) based on the simple exponential profile ofHartmann layers. It gave good results in the flow regime where inertia is small but failed to de-scribe flows where strong rotation induces secondary flows, such as Ekman pumping. The 2D modeldeveloped by Poth´erat et al. (2000) (denoted as PSM hereafter), accounting for some 3D effects,gave more accurate prediction in the Q2D flows. With PSM, both of the velocity profiles and theglobal angular momentum measurements from MATUR (Poth´erat et al. , 2005) were reproduced,and it was proved that the local and global Ekman recirculations altered the shape of the flowsignificantly as well as the global dissipation. However, both the SM82 model and the PSM modelbreak down if the Hartmann layer becomes turbulent, where the flow may still remain quasi-two-dimensional but the boundary layer friction is altered. Poth´erat & Schweitzer (2011) established analternative shallow water model specifically for this case, and recovered experimentally measuredvelocity profiles and global momentum in this regime.In mainly azimuthal flows such as in a toroidal containers, the dynamics of the side wall layer andthe free shear layer near the injected electrodes on the flow is complex because of rotation effects.Even when the Hartmann layers are stable, significant flow alterations may occur, including non-trivial 3D effects (Tabeling & Chabrerie, 1981), which could not be observed easily in experimentsor with any Q2D model. It has been proven that turbulence may remain localized in a layer nearthe outer cylinder wall prior to transition happening in Hartmann layers as indicated by Zhao& Zikanov (2012), who conducted a series of 3D DNS of MHD turbulence flows in a toroidalduct. In addition, for cases with lower values of Hartmann number and higher values of Reynoldsnumber, three-dimensionality would be more pronounced, even within the Hartmann − B ¨ odewadt layers, which have been studied theoretically by Davidson & Poth´erat (2002) and Moresco &Alboussi`ere (2003). All of these discoveries encourage us to perform 3D DNS on the flows inMATUR configuration (corresponding to the realistic experiment of Messadek & Moreau (2002)).Besides reproducing the results obtained in the experiments, theories and Q2D simulations, wefocus on answering the following questions in the present work.1. Does the separation or turbulence emerge within the side wall layer when the electrodes arefar away from the side wall and while the Hartmann layer remains laminar?2. What causes the angular momentum dissipation in regimes where the Hartmann layer islaminar?3. How much and what type of three-dimensionality subsists in sheared MHD turbulence athigh Ha ? In particular,(a) How much energy subsists in the secondary flow?(b) Is there a cutoff lengthscale between quasi-2D and 3D lengthscales in sheared turbulencetoo?However, two factors restrict the investigated range of Reynolds number and Hartmann numberin the present DNS studies. One is the computing resource, because very fine grids are requiredto capture the small-scale turbulent structure and to resolve the thin Hartmann boundary layers.Another is the lack of robust computational schemes capable of dealing with nonlinear unsteadyhigh- Ha flows. In particular, when non-orthogonal grids are used, extra non-orthogonal correctionschemes are required. By applying Large eddy simulations (Kobayashi, 2006, 2008), Re could besomewhat increased, but the resolution requirements for the Hartmann layers remained essentially . Chen et al. DNS of MATUR the same as those in DNS, since no reliably accurate wall function models were known for thecase of turbulent flows. Here, the problem of inadequate computational resources is overcome byemploying massively parallel computing. As for the numerical method, we apply the finite volumemethod based on the consistent and conservative scheme developed by Ni et al. (2007), which canbe used to accurately simulate MHD flows at a high Hartmann number. Therefore, the Hartmannnumber is allowed to vary from 55 to 792 (magnetic fields change from 0 . §
2, a short description of the physical model underlyingthis work and the flow conditions in the MATUR cell are given. Particular attention is given tothe modifications dealing with the electrical conductive wall and the injected current density. Thenumerical algorithm and the detailed computational grid study are also presented in this part. Themain properties of the MHD turbulence are described and discussed in §
3, including the generalaspect of the flow, the secondary flows, the properties of the free shear layer and side wall layer,the global angular momentum as well as three dimensionality. Finally, we offer concluding remarksin § The physical model of MATUR is shown in Fig.1. It is a cylindrical container (radius (cid:101) r = 0 . a = 0 .
01 m, with (cid:101) distinguishing the dimensional quantity from their dimensionlesscounterpart), in which the bottom and the upper walls are electrically insulating while the verticalwalls are conducting. Electric currents are injected at the bottom of the container through alarge number of point-electrodes spread along a concentric ring parallel with the vertical wall. Asindicated by Messadek & Moreau (2002), a continuous electrode ring would induce a strong localdamping in the flows, and thus, a series of discrete point-electrodes are positioned to reduce thisunwanted effect. The concentric circles are located at (cid:101) r e = 0 .
054 m or (cid:101) r e = 0 .
093 m, respectively.In the present study, only (cid:101) r e = 0 .
054 m is considered. The container is filled with mercury, andexposed to a constant homogeneous vertical magnetic field parallel to the axis of MATUR. Theinjected currents leave the fluid through the vertical wall, to induce radial electric current that giverise to an azimuthal force on the fluid in the annulus between the electrode circle and the outerwall. The material properties of the fluid at room temperature, such as the mass density ρ , thekinematic viscosity ν and the electrical conductivity σ , are assumed constant ( ρ = 1 . × kg/m , ν = 1 . × − m /s and σ = 1 . × S/m). An external homogeneous magneticfield of amplitude B is applied along the axial direction. At a low magnetic Reynolds number, thefull system of the induction equation and the Navier-Stokes equations for an incompressible fluidcan be approximated to the first order O ( Rm ). Thus, the non-dimensional magnetohydrodynamicequations governing the flow can be written as (Roberts, 1967) ∇ · v = 0 , (2) ∂ v ∂t + ( v · ∇ ) v = −∇ p + 1 Re ∆ v + N ( j × e z ) , (3) j = −∇ ϕ + v × e z , (4) ∇ · j = 0 . (5) . Chen et al. DNS of MATUR Figure 1: Sketch of the experimental set up. A typical electric circuit including one of the point-electrodes mounted flush at the bottom Hartmann layer is represented.where the variables j , ϕ , v , p denote the current density, the electric potential, the velocity andthe pressure, respectively. Here, lengths are scaled by a , the velocity by a scale U to be specifiedshortly, and the current by σBU . The typical scales for the other variables are as follows: ρU for the pressure, U Ba for the electrical potential, a/U for time. The Hartmann number Ha andthe interaction parameter N are defined as Ha = Ba (cid:114) σρν , N = σB aρU , (6)and the Reynolds number is given as Re = Ha /N . In the present work, all these non-dimensionalnumbers are based on the thickness of the container a . For this experiment, an approximateazimuthal velocity can be derived from the theory in Sommeria & Moreau (1982). Indeed, Poth´erat et al. (2000) derived an approximate expression for the z -averaged azimuthal velocity in the inviscidlaminar, axisymmetric Q2D regime, using a Dirac delta function centred at the electrodes r = r e to describe the injected current j W , where the integral is equal to the total injected current I : j W = I/ πr e δ ( r − r e ). (cid:26) U SM82 θ ( r ) = I πr √ σρν , r e < r < r ,U SM82 θ ( r ) = 0 , r ≤ r e . (7)and U SM82 r = 0. Based on this, we choose the velocity scale as U = I πr √ σρν . Finally, velocityfluctuations are defined as v (cid:48) = v − (cid:104) v (cid:105) .The boundary conditions for v and ϕ are as follows. For the velocity, we apply standard no-slipconditions at all walls. As for the electric potential, we impose perfectly conducting side walls andperfectly insulating Hartmann walls except for the locations where the electrodes are located, i.e. :at the top wall, v = 0 , ∂ z ϕ = 0 at z = 1 (top wall) , (8)at the surface of point electrodes located at the bottom wall: v = 0 , ∂ z ϕ = IA re σ BU at z = 0 (9) . Chen et al. DNS of MATUR at bottom wall, outside point electrode surfaces: v = 0 , ∂ z ϕ = 0 at z = 0 (10)at the lateral wall (BC1), v = 0 , ϕ = 0 at r = r (11)Here, 128 points-electrodes (0 .
001 m diameter and 0 . r e at the bottom wall. The surface of each point-electrode contains 102cells with those on the periphery cut by the electrode edge. The area A re , which appears in theboundary condition of (9), is the total area of the small electrodes. Note that we do not need anextra forcing term in (3) to drive the flow because the injected current from the bottom electrodesinteracts with the magnetic field, and generates Lorentz force that drives the flow. The currentcirculate inside the mercury, as shown in figure 1, and leaves the domain at the vertical wall.Moreover, because of the perfectly conducting properties of the wall, the electrical potential acrossthe wall can be regarded as a constant, and set to zero in (11). The direct numerical simulations of the governing equations are performed based on the finitevolume approach. For the pressure − velocity coupling, a second-order temporal accurate pressure-correction algorithm has been used. Based on a consistent and conservative scheme (Ni et al. , 2007),the electrical potential Poisson equation is then solved to obtain ϕ . The detailed process withineach time step could be split into: a) obtain a predicted velocity by solving the momentum equationwith pressure from the previous iteration; b) calculate the predicted velocity fluxes which are usedas the source term of the pressure difference Poisson equation, which will be solved to obtain thepressure difference; then apply the updated pressure difference to update the velocity and pressure;c) solve the Poisson equation for the electric potential to get the electrical potential, which is usedto calculate the current density fluxes on cell faces with the consistent and conservative scheme.Then the current density at each cell centre is reconstructed through a conservative interpolation j = ∇ · ( jr ) with r the position vector; d) calculate the Lorentz force F Lorentz = j × e z at the cellcentre based on the reconstructed current density, which is used as the source term of momentumequation of next time step. Step (b) is iterated 3 times before solving the electrical potentialPoisson equation.The central scheme is applied for all convective-term approximations. All inviscid terms and thepressure gradient are approximated with a second-order accuracy. A second order implicit Eulermethod is used for time integration. In order to guarantee a robust solution for unsteady flows andmake the temporal cutoff frequency match the spatial cutoff frequency, the present simulations arerun with a constant time step which satisfies the Courant-Friedrichs-Lewy condition.A Preconditioned Bi-Conjugate Gradient solver (PBiCG) applicable to asymmetric matriceshas been used for the solution of the velocity-pressure coupling equation, together with a Diagonal-Incomplete LU (DILU) decomposition for preconditioning. The pre-conditioned conjugate gradient(PCG) iterative solver with a diagonal incomplete Cholesky (DIC) pre-conditioner, which dealswith symmetric matrices, has been applied for the solution of pressure and electric potential equa-tions. Note that for all the iterative solutions, velocity, pressure and electric potential, a constantconvergence criterion of 10 − is used.In order to verify the accuracy of our numerical code, the classic Hunt’s flow has been simulatedfor comparison with the analytical solution (Hunt 1965). Herein, the parameters are set to Re =100, Ha = 5000, and the conductance ratio of C w = σ w t w /σ f L f is set to 0.01. As illustratedin figure 2, the velocity matches well with Hunt’s analytical result, especially within the thinboundary layer from z = 0 . z = 1. Moreover, the calculated pressure gradient matches well . Chen et al. DNS of MATUR Figure 2: The typical ”M-shape” velocity profile of Hunt’s case (a) and the comparison of thenumerical result with Hunt’s analytical solution (b) for Ha = 5000.Figure 3: The radial distribution of the mean azimuthal velocity (a) and the rms of velocityfluctuations (b) for BC1 and BC2 at Ha = 66 and Re = 15972. . Chen et al. DNS of MATUR with the analytical pressure gradient, with a relative discrepancy lower than 0 .
08% based on anerror estimation of | ( ∇ p ) Anal − ( ∇ p ) Numr ( ∇ p ) Anal | .In addition, it should be noted that the boundary condition of a perfectly conducting sidewall which we used in all simulations ((11), denoted as BC1) is not entirely consistent with thereal experiment where the total current was imposed through the lateral wall. The experimentalconditions would be represented by replacing (11) by the electric boundary condition v = 0 , ∂ r ϕ = − IA side σ BU at r = r , (12)denoted as BC2, and where A side is the surface area of the side wall. For this reason, we conductedtwo further validation steps for Ha = 66 and Re = 15972. Firstly, we compared the total currentthrough the side wall with the total injected current I based on BC1. The relative error of 0 . .
8% and 3%, respectively.This indicates that the choice of either BC1 or BC2 does not have any significant impact on theresulting flow field and that both are compatible with good conservation of charge.As a further validation, the numerical solutions of the average azimuthal velocity (cid:104) U θ (cid:105) θ andthe time and space average of the angular momentum, L lam , are also compared with the resultsproduced by the Q2D model, as shown in Table.1. Here, U i , ( i = θ, r, z ) represents individualvelocity components averaged over time in a quasi-steady state, (cid:104)·(cid:105) i , ( i = θ, r, z ) represents theaverage along the specific direction (Hartmann layers are excluded), and (cid:104) L lam (cid:105) t = 1 T (cid:90) T (cid:104) L lam (cid:105) V dt, (cid:104) L lam (cid:105) V = 1 V Ω (cid:90) rv θ ( r ) d Ω , (13)where V Ω = 2 π ˜ r is the non-dimensional volume of the computational domain, (cid:104)·(cid:105) V is the volumeaverage and v i is the instantaneous velocity. According to the theory of Sommeria & Moreau(1982), an approximate global angular momentum can be derived (Poth´erat et al. , 2000), underthe assumption of axisymmetry, i.e. L SM = I πU a √ ρνσ (1 − r e r ) . (14)Across the range of considered parameters, the maximum relative errors of (cid:104) v θ (cid:105) θ,t /U SM82 θ ( r )and (cid:104) L lam (cid:105) t /L SM are less than 2 . Due to the localisation of the Lorentz force within r ∈ [5 . , r = 5 .
4. In order to capture more precise flowinformation, highly refined grid resolution are required in both of the free shear layer and theouter wall side layer. Note that the unstructured computational grids are made of hexahedra andprisms in the present study, and the grid details in case of Ha = 792 and Re = 15972 are usedfor illustration. In the radial direction, N r = 864 grid points are generated, 30 ( resp.
25) of whichare devoted to the side wall ( resp. free shear) layer located at r = 11 ( r = 5 . γ r = 1 . r = 11 . Chen et al. DNS of MATUR N z N Ha N Sh τ Sh τ Ha (cid:104) v θ (cid:105) θ,t /U SM82 θ ( r, z ) (cid:104) L lam (cid:105) t /L SM G G G G G N Ha , N Sh denote the grid points within the Hartmann layer andside wall layer, respectively. τ Sh , τ Ha , (cid:104) L lam (cid:105) t /L SM and (cid:104) v θ (cid:105) θ,t /U SM θ ( r, z ) at r = 9 . z = 0 . Ha = 792 , Re = 15972.and r = 5 . r min ≈ . r max ≈ . N θ = 4096 uniformly spaced gridpoints and the axial direction uses N z = 300 non-uniformly spaced grid points. Note that to fullyresolve the Hartmann layer along the axial direction, 25 uniformly spaced grid points are devotedto each layer and a smooth transition is set toward the core region where a coarser grid resolutionis sufficient. Grid points are spread according to a geometric sequence of ratio γ z = 1 . z ≈ Ha − and z ≈ − Ha − . The point nearest to the Hartmann walls is located ∆ z min ≈ × − away from them, and the largest step is ∆ z max ≈ . × − .In order to assess the quality of the grids, wall coordinates are introduced r + = Re Shτ r, z + = Re Haτ z, (15)where Re Shτ = (cid:112) Reτ Sh , Re Haτ = (cid:112) Reτ Ha . (16) τ Sh = 1 T A Sh (cid:90) (cid:90) ( − ∂v θ ∂r ) (cid:12)(cid:12)(cid:12) r =1 dtds, τ Ha = 12 T A Ha (cid:90) (cid:90) ( ∂v θ ∂z (cid:12)(cid:12)(cid:12) z =0 − ∂v θ ∂z (cid:12)(cid:12)(cid:12) z =1 ) dtds. (17)Here, τ Sh , τ Ha denote the associated dimensionless forms of the mean stress at the side walland Hartmann wall, respectively. T is the time interval for average, A Sh = 22 π and A Ha = 121 π .Hence, the value of the smallest wall-normal grid step in the ∆ r + units varies from 0.08 (seeTable.1). The respective variation in the ∆ z + units is from 0.18 (see Table 1). Moreover, thesimulated results show that the highest velocity occurs at r ≈ .
04, where the azimuthal grid step∆ θ ≈ .
011 is sufficiently small. Meanwhile, in order to ensure that the full range of the dissipativescales is resolved, the smallest turbulent scales ( l min ⊥ and l min z ) predicted by Poth´erat & Dymkou(2010) are used to evaluate the grids quality. The adopted grids indicate ∆ max ⊥ l min ⊥ (cid:39) . , ∆ max (cid:107) l min z (cid:39) . Re t = U L S L ν ).Here, S L is the size of the large scales, which is evaluated from the profiles of RMS of relativeazimuthal velocity fluctuations (see figure 8(a) of Poth´erat & Schweitzer (2011)), and U L is thevelocity of the large scales, which is calculated according to U L = ( (cid:82) v (cid:48) i v (cid:48) i dv ) / , i = ( θ, r, z ), where v (cid:48) i denotes the velocity fluctuations. The sizes of smallest scales according to Poth´erat & Dymkou(2010) for all cases are listed in Table 2.In addition, the grid independence studies are also conducted on the numerical case of Ha = 792and Re = 15972. Note that not only the grid sizes, but also the grid points along the z -direction,within the Hartmann layer and within the side wall layer are tested. The time averaged wall stress τ Sh , τ Ha , the time averaged angular momentum (cid:104) L lam (cid:105) t and the time-space average azimuthalvelocity (cid:104) v θ (cid:105) θ,t , at r = 9 . z = 0 . Ha = 792 , Re = 15972, (7) and (14) are not strictly valid since the flow is not axisymmetricbut remain sufficiently accurate to roughly assess the accuracy of the simulations. . Chen et al. DNS of MATUR Figure 4: Comparison of (cid:104) v θ (cid:105) θ,z between experiment (symbols), present numeric (solid lines), 2Dnumeric based on PSM model(dashed dot lines) and U SM82 θ ( r ) (dashed lines) for cases at Ha = 792, Re = 15972 ( a ) , and Ha = 792, Re = 31944 ( b ). The dashed line depicts an algebraic law r − (predicted by Eqn.(7)).All the simulations are stopped when the total angular momentum of the flow is statisticallysteady, i.e. after 3 t Ha , where t Ha = ˜ t Ha / ( a/U ) and ˜ t Ha = a / ( νHa ) (18)denote the non-dimensional and dimensional Hartmann damping times, respectively. Average andRMS quantities are then evaluated over a time interval of 4 t Ha , and the computed values of theseparameters are compared to evaluate the grid resolution within the Hartmann layer and side walllayer. Evidently, the results are very close to each other, even with the worst spatial resolution,as shown in Table.1. Firstly, the reliability and the accuracy of the DNS results are confirmed bya difference of less than 2 .
1% between the present results and the solutions predicted by (7) and(14). Moreover, grid-independent solutions are also achieved with the grids under consideration.For example, less than 3 .
2% difference is found for all the predicted values on the coarsest gridand the finest grid, and this discrepancy is even further reduced when the two finest grids arecompared. Hence, one can conclude that grid G Ha = 792 and Re = 15972. The mesh is, however, further refined when 3D effects becomesignificant, e.g. in the case of Ha = 55 and Re = 15972.For different numerical cases, the dimensionless time steps used in the computations are alteredaccording to the parameters applied in simulation, such that ∆ t = 5 . × − for ( Ha, Re ) =(264 , . × − for ( Ha, Re ) = (264 , . × − for all other cases. The valuesare determined by the limits of numerical stability, which highly depends on the viscous term andthe convective term of the momentum equations at high Re . In addition, higher Re or Ha demandsmaller time steps because of the higher azimuthal velocity and the thinner Hartmann layers. Ha (cid:29) and N (cid:29) As far as the authors know, it is the very first attempt to reproduce the MATUR experimentby performing 3D direct numerical simulations, and hence as a necessary validation procedure, a . Chen et al. DNS of MATUR detailed comparison with the available experimental results need to be carried out. In addition,since the PSM model can deal with the cases when Ha (cid:29) N t = N ( (cid:101) r /a ) (cid:29)
1, somenumerical cases falling into this space are also investigated for validation. However, note thatthis model becomes imprecise when either of the parameters, Ha or N v become of the order of 1,( N v = N ( (cid:101) r /a ) is the interaction parameter based on the horizontal scale), which, again, stressesthe importance of conducting 3D direct numerical simulations.The predictions of the mean azimuthal velocities from different approaches, denoted by (cid:104) v θ (cid:105) θ,z ,are plotted in figure 4. Clearly, a good agreement is found between the experimental results,the numerical data and the laminar theoretical prediction. In the outer region r e < r < r , themaximum relative discrepancy between the results of DNS and experiment is less than 6 .
7% (8 . Re = 15972 ( Re = 31944) based on an error estimation of (cid:107)(cid:104) v DNS θ (cid:105) t −(cid:104) v exp θ (cid:105) t (cid:107) (cid:107)(cid:104) v exp θ (cid:105) t (cid:107) . Themaximum relative discrepancy between the results of DNS and the PSM model is less than 1 . . Re = 15972 ( Re = 31944) based on an error estimation of (cid:107)(cid:104) v DNS θ (cid:105) t −(cid:104) v PSM θ (cid:105) t (cid:107) (cid:107)(cid:104) v PSM θ (cid:105) t (cid:107) .In particular, the velocities exhibit the characteristic feature that they increase sharply across thefree shear layer ( r = r e ) due to the current injection. Accordingly, the shear layer separates theflow into an outer and inner regions. Moreover, from both of the experimental and the numericaldata, the downward trend of the azimuthal velocity in regions between the injected electrodesand the vertical wall follow the expected scaling law of (cid:104) v θ (cid:105) θ,z ∼ r − , as predicted by (7), andwhich reflects the geometrical spreading of the radial forcing current in the Hartmann layer, i.e. j r ∼ (4 πr ) − (see figure 5(c)).A typical instantaneous distribution of v θ obtained numerically over a radial cross-section θ = 0for moderate forcing current is shown in figure 5(a). Under a strong magnetic field, the velocitygradient along the magnetic field lines is remarkably damped, except in the Hartmann layers wherean exponential profile subsists. The detailed velocity distribution in the Hartmann layer is shownin figure 5(b), and a good agreement is observed between the present numerical results and theexact solution, given as v θ = v core θ (1 − exp( Haz )), with v core θ indicating the azimuthal velocity in thecore flow. It also demonstrates that the thickness of the Hartmann layer at r = 7 . r = 9 is thesame. Besides, figure 5(c) reveals that the vertical profiles of radial current density within the topand bottom Hartmann layers collapse with each other, implying that the electric current intensity I injected at the electrodes divides in two equal parts between the two symmetric Hartmann layers.In addition, the the radial current density is much higher near the Hartmann wall, so the Jouledissipation mainly takes place in the thin Hartmann layers, in line with the laminar Hartmannlayer theory. Therefore, the annular fluid domain located between the selected circular electrodesand the cathode (5 . ≤ r ≤
11) is driven in the azimuthal direction by the Lorentz force, whilethe central fluid domain ( r < .
4) is entrained by friction within the free shear layer. Interestingly,the current density at the upper wall stands a little lower than at the bottom wall, showing thatdespite the excellent agreement between PSM and experimental data, the flow is ever so slightlythree-dimensional.In the following part, the evolution of the large structures and the spectral analysis are dis-cussed. Subsequently, we study the secondary flow induced by Ekman pumping. We also investigatethe characteristics of the free shear layer and the side wall layer before presenting the turbulentstatistics, global angular momentum and three-dimensionality.
The different cases investigated are listed in Table 2, where the dimensionless parameter R (= Re/Ha ) represents the Reynolds number scaled on the thickness of the Hartmann layer. Accordingto the experiments of Moresco & Alboussiere (2004), the flow within the Hartmann layer becomesturbulent when R ≥ . Chen et al. DNS of MATUR Figure 5: ( a ) The distribution of the instantaneous azimuthal velocity and current streamlines onplane θ = 0 at Ha = 264, Re = 15972. ( b ) The vertical profiles of the instantaneous azimuthalvelocity within the bottom Hartmann layer along r = 7 . r = 9. ( c ) the distribution of radialcurrent density both across top Hartmann layer and bottom Hartmann layer along r = 7 . r = 9, where the origin of the distribution along r = 9 is shifted to n = 0 .
02, where n is thewall normal coordinate i.e. for the bottom wall, n = z , and for the top wall n = 1 − z and v coreθ ( r ) = (cid:104) v θ ( r, θ, z = 0 . , t ) (cid:104) θ,t . Here v θ is averaged in time and θ . Ha
264 528 792 792 264 264 132 110 99 80 66 55 Re N v R l min ⊥ ( × − ) 3.95 3.05 3.10 2.52 2.91 2.45 2.88 2.84 2.81 2.77 2.71 2.67 l minz ( × − ) 1.92 1.31 1.45 1.34 1.22 1.12 1.10 1.08 1.01 0.96 0.92 0.89Table 2: Non-dimensional parameters in cases calculated numerically. . Chen et al. DNS of MATUR Therefore, only cases with
R <
380 are considered in this paper. Furthermore, five of the relevantinteraction parameters scaled on the horizontal length N v are at the order of unity, aiming to studythe three dimensionality of the flows.In the calculations, the electrical current is injected at t = 0 when the fluid is at rest andremains constant during the whole simulation, following the actual experimental procedure. Fordifferent electrical current intensities, the flow goes through a sequence of evolution and reachesdifferent equilibrium and quasi-equilibrium states presented on figure 6. We shall now give anoverall view of the numerical results, while more details of the evolution and local quantities willbe reported later.For R ≤ ≤ R ≤ R ≥ z axis in near-solid body rotation. They mainlyremain localized near the free shear layer once generated there, as shown in figure 6(a). Thus, thevelocity fluctuations in the inner region and the annular outer region are of much lower intensity,and the azimuthal velocity contours reveal that the wall side layer and the Hartmann layer arestable. Besides the low value of R , part of the reason for the stability of the side layer is that theselarge vorticity structures remain distant from it, and little interaction between them takes place.However, the thickness of the side wall layer is still smaller than the scaling for a straight duct( Ha − / ), due to the recirculations induced by Ekman pumping, a point we will analyse in detailin § R ≥
121 small scale three-dimensional turbulence appears in the side layer. The onset ofthree-dimensional turbulence within R (cid:46)
121 is consistent with the value of 138 reported by Zhao &Zikanov (2012), albeit a little lower. The difference in curvature of the external wall ( a/ ˜ r = 1 / a/ ˜ r = 1 / R ≥ .
2, the size of the large structures increases, leading to highly turbulent fluctuationsin both the inner and the outer regions. Accordingly, long azimuthal vorticity streaks of relativelyhigh intensity exist in the outer region that induce instabilities within the wall side layer. Con-currently, the wall side boundary layer separates from the wall, suggesting that separation resultsfrom the fluctuations in the outer region, as observed previously observed by Poth´erat & Schweitzer(2011). For Ha = 55 , Re = 15972 ( R = 290), Ekman recirculations are strong, the Hartmann layeris relatively thick so the centripetal radial velocity in the Hartmann layer is relatively strong (seefigure 12). Thus, in the vicinity of a free shear layer, figure 6(f) conveys that the azimuthal velocitynear the Hartmann layer is higher than in the core. Within the side wall layer, one can observethe dramatic variation of the azimuthal velocity and vertical velocity (see figure 7(a)) along themagnetic field lines besides the separation of the layer.Finally, in the examples shown here, three dimensional effects are only noticeable within theshear layer for R ≥
121 ( Ha =132, Re = 15972, this is in fact better seen from the analysis of thevertical velocity in the side layer in section 3.5). By contrast, for R (cid:39)
290 ( Ha = 55, Re = 15972),weak three-dimensionality, where flow patterns are topologically identical but less intense near thetop wall, exists outside the shear layer (see figure 7(b), (c)). The presence of three-dimensionality, . Chen et al. DNS of MATUR Figure 6: Typical snapshots of equilibrium or quasi-equilibrium states obtained from numericalsimulations. Contours of the magnitude of the vorticity on plane z = 0 . θ = 0 (right column). ( a ), ( b ) Ha = 792 , Re = 15972.( c ), ( d ) Ha = 264 , Re = 31944. ( e ), ( f ) Ha = 55 , Re = 15972. The velocity are normalised with U . . Chen et al. DNS of MATUR Figure 7: Snapshots of velocity contours for Ha = 55, Re = 15972. ( a ) Distributions of axialvelocity v z in the plane θ = 0 (near the side wall). ( b ) Distribution of vorticity in the plane nearthe top Hartmann wall at z = 0 .
8. ( c ) Distribution of vorticity in the plane near the bottomHartmann wall at z = 0 .
2. Note that the contour levels for vorticity and v z are chosen so as toenhance the visibility of turbulent flow structures.however, is controlled by the true interaction parameter at the scale of the considered structure.A consequence is that inertia-induced three-dimensionality is expected to appear in a parallellayers of thickness δ || ∼ aHa − / when the local turnover time δ || /U becomes smaller than thetwo-dimensionalisation time at that scale ρa /σB δ || , i.e. when the Reynolds number based onthe parallel layer thickness R || = U δ || /ν exceeds unity. By contrast, since the separation of thewall-side layer is induced by the tail of large two-dimensional structures, which is mostly quasi-twodimensional, it can be expected to be controlled by R . For flows without boundary layer separation, the parameters Ha = 132, Re = 15972 are chosento illustrate the typical evolution process. According to the variation of the azimuthal velocitysignals, as shown in figure 8(a), three different stages can be distinguished before the flow reachesits final quasi-equilibrium state. The acceleration of the fluid in a laminar regime corresponds tostage 1. After a short time, the laminar shear layer at r = 5 . . ≤ r ≤
11 is driven in rotation by the Lorentz force.Poth´erat et al. (2005) estimated the stability threshold for this layer as
Re/ √ Ha < .
5, implyingthat the circular free shear layer becomes linearly unstable when the Reynolds number based onits thickness exceeds the threshold of 2.5. For all injected current intensities considered here, thiscritical value on the azimuthal velocity is reached very quickly. Then the circular free shear issubject to a Kelvin-Helmholtz instability, which breaks it up into small vortices ( t/t Ha = 0 . ω z , is presented in figure 9. These vortices merge into larger structures very soon aftertheir inception ( t/t Ha = 1 . , . , .
261 9(b-d)). They become distorted because of the shear.As shown in figure 8(a), the velocity at two representative locations keeps increasing during thisstage. The amplitude of velocity oscillation at r = 8 . r = 9 .
6, theturbulent structures exert a weaker influence in the region far from the electrodes.Stage 3 corresponds to the quasi-equilibrium state of the flow, which, at this point, cyclesthrough a recurring sequence. First, large structures progressively loose intensity and elongatealong the free layer ( t/t Ha = 5 . t/t Ha = 5 . t/t Ha = 5 . . Chen et al. DNS of MATUR Figure 8: ( a )Typical instantaneous azimuthal velocity signals vs. time at (6 . , , .
5) and(9 . , , .
5) at Ha = 132 , Re = 15972. The solid lines denote the numerical results and thedashed lines denote the theoretical value derived by Messadek & Moreau (2002), i.e. U SM82 θ ( r, t ) = I πrU √ σρν (1 − exp( − tt Ha )). ( b ) Corresponding power density spectra at point (5 . , , .
5) withinthe free shear layer.Figure 9: Evolution of the flow with time at Ha = 132, Re = 15972. The distribution of the axialvortex structures, ω z , on plane z = 0 . . Chen et al. DNS of MATUR Figure 10: ( a )The average side wall shear stress signals variation with time in case of Ha = 55, Re = 15972 and Ha = 132, Re = 15972. ( b )-( f ) Contours of the magnitude of the vorticity onplane z = 0 . Ha = 55, Re = 15972.( t/t Ha = 5 . f .
26) shown in figure 8(b).As R is increased, the most striking feature is the appearance of separation and turbulencewithin the side wall layer. Figure 10(a) shows the evolution of the mean shear stress on the sidewall τ Sh , with a sudden increase of τ Sh between t and t for Ha = 55 , Re = 15972. In this interval,the increase of τ Sh can be ascribed to the random turbulent fluctuations, seen more in detail insection 3.5. From t (figure 10(c)), vorticity streaks attached to the vortices generated at the freeshear layer start reaching out to the outer side layer, and incur local variations of its thickness.At t (figure 10(d)), these variations have become severe to the point of incurring boundary layerseparation. The decrease of τ Sh at t shown in figure 10(a), confirms the occurrence of separationat the wall side layer. This is consistent with the visualizations of the vortex structures shown infigure 10(d). For Ha = 132 and Re = 15972, by contrast, the evolution of τ Sh is rather smoothand no brutal change in shear stress is observed, indicating the absence of separation. Moreover,the evolution of case at Ha = 55, Re = 15972 in the later stage ( t , t , t ) is similar to that ofcase at Ha = 132, Re = 15972, which goes through the same recurring sequence. The power density spectra of several typical cases are analysed in this section shown in figure 11. Tocalculate of the average, we employ 60 probes, uniformly distributed along the angular direction,and take the average of the measured signal as v ( t ) = Σ v i ( t ). Spectra are obtained usingWelch’s averaged periodogram method, while a Hamming window was applied to each overlappingsegment of data. Additionally, spatial energy spectra can be deducted from these spectra by takingadvantage of the large average azimuthal velocity and using Taylor’s hypothesis, 2 πf = (cid:104) U θ (cid:105) θ k . . Chen et al. DNS of MATUR f t Ha E ( f ) f t Ha E ( f ) f tr f t Ha E ( f ) f tr
3 5/3 R f t r t H a k
150 200 2501020304050 012 f tr , DNSk,DNSf tr , 0.16Rk=1 Figure 11: Frequency spectra obtained from velocity time-series. For reference, the power laws f − . , f − / and f − , denoted with dash-dot line, are also shown. Solid lines correspond to spectrafiltered using Bezier spline (red), time-series of azimuthal velocity time-series at point ( r = 6 . z = 0 .
5, green (within the free shear layer)), axial velocity at point ( r = 6 . z = 0 .
5, red), axialvelocity at point ( r =10.98 (b) 10.95 (c) , z=0.5, blue (within the side wall layer)). azimuthalvelocity at point ( r =10.95, z = 0 . a ) Ha = 528, Re = 15972 ( R = 30 .
3) ( b ) Ha = 264, Re = 31944 ( R = 121), ( c ) Ha = 55, Re = 15972 ( R = 290 . d ) The separation frequency ofthe spectra slops f − / and f − and the forcing scales for different R . . Chen et al. DNS of MATUR For R = 30 .
3, the spectrum of v θ ( t ) exhibits an expected strong peak at a fundamental fre-quency corresponding to the passage of the large structures through the measuring probes. Theother noticeable peaks represent the harmonic and sub-harmonic frequencies, as shown in figure11(a). The scaling of the spectrum in the inertial region obeys scales as f − . . According to Eckert et al. (2001), the spectral exponent in duct flow turbulence for N ≈ . R = 121, the peak of the base frequency in the spectra of azimuthal velocity fluctuationscan still be observed (see figure 11(b)), and thus the large structures continue to rotate aroundthe axis. The spectrum in the inertial region can be separated into two parts, with a transitionfrequency f tr t Ha (cid:39)
20. For f tr t Ha ≤
20, the power spectral density scales as f − / , while f tr t Ha >
20, is scales as f − . By applying the Taylor’s hypothesis, we find the corresponding azimuthalwavenumber of k tr (cid:39) .
10. This value is in accordance with the results of Messadek & Moreau(2002), i.e. ˜ k tr (cid:39) − . The authors claimed that the split spectrum arises as the result of weakJoule dissipation. The spectrum of v z ( t ) within the free shear layer exhibits axial flow that is ofmuch lower energy than the other two components, confirming that the turbulence is dominatedby its 2D horizontal components.For R = 290 .
4, the spectra exhibit different features (see figure 11(c)). In the region of thefree shear layer, the f − (high frequencies) and the f − / (low frequencies) power laws are stilldistinguishable on the azimuthal velocity. The transition frequency, f tr t Ha (cid:39)
46, is higher butthe corresponding non-dimensional wavenumber almost remains the same, k tr (cid:39) .
09. Within theside wall layer, by contrast, the power density spectra of azimuthal velocity (orange line) and axialvelocity (blue line) in the inertial region, exhibit a scaling of the form f − / .The apparent constance of k tr prompts us to compare the transition frequency and the forcingscale for cases of R ≥ R follows f tr t Ha (cid:39) . R . Since in all cases, k (cid:39)
1, the separation between the two slopesmay simply result from the forcing geometry. In this case, the separation between k − / and k − may reflect the usual 2D split between an inverse energy cascade and a direct enstrophy cascade,as already found in MHD flows by Sommeria (1986). In this case, k tr may be interpreted as theforcing scale, at which the mean flow transfers energy to the mean flow (Alexakis & Biferale, 2018). Poth´erat et al. (2000) showed that the recirculations induced by Ekman pumping significantlyinfluence the flow. They can be identified by looking at the variations of the velocity profilein figure 6(e). Although the PSM model already provides a good understanding on the Ekmanpumping effect, the detailed flow information across the fluid layer and the structures on the planeparallel to the magnetic fields remain unexplored. This is one of the motivations for carrying out3D DNS.The streamlines of the average fluid velocity shown in figure 12(b) reveal that the Ekmanpumping induces a centripetal flow ( i.e. away from the side wall) in the Hartmann layers and acentrifugal flow ( i.e. towards the side wall) in the core. In addition, the streamlines are almostsymmetric with respect to the centre plane z = 0 . . Chen et al. DNS of MATUR Figure 12: The distribution of mean radial velocity ( a ), and mean axial velocity, mean streamlines( b ) on plane θ = 0 at Ha = 80 , Re = 15972.in figure 12(b). This mechanism explains that recirculations are stronger within the free shear layerand the side wall layer. Ekman pumping incurs a net centrifugal transport of angular momentumas the velocity is smaller in the Hartmann layer. This has two consequences: a squeeze of theside wall layer and an increased dissipation at the side wall layer when these recirculations areimportant.Poth´erat et al. (2000) have derived the analytical expression for the vertical velocity at theinterface between the Hartmann layer and the core, u − z = − λHaN ∇ ⊥ · [( u −⊥ · ∇ ⊥ ) u −⊥ ] , (19)Here, u − denotes the velocity at the edge of the Hartmann layer and the subscript ⊥ denotes thevector projection in the direction perpendicular to the magnetic field. Under the assumptions ofPSM, any vertical component of velocity is associated to recirculations, whether local or global.Therefore, to assess the limits of the PSM approach, we compare the energies associated to thevertical velocity component obtained with DNS at Re = 15972 to the energy obtained from (19),distinguishing energies E mz and E (cid:48) z associated to the average flow and to the fluctuations. TheDNS results show that < E mz > t ∼ Ha − . and < E (cid:48) z > t ∼ Ha − . (see figure 13(c), bearing inmind that Re is kept constant in these scalings). Values of E z obtained by integrating (19) fromthe results of PSM simulations follow the same scaling even for values of N as low as N (cid:39) . Ha = 55) where the model is expected to break down. This indicates that PSM predicts theglobal recirculations associated to the mean flow very accurately. The energy associated to thefluctuations predicted by PSM, by contrast, only matches DNS precisely for N (cid:38) Ha (cid:38) . E (cid:48) z considerably. The origin of the discrepancy can be foundin the radial profiles of vertical velocity fluctuations (figure 13(a), figure 13(b)): the discrepancybetween the profiles obtained with PSM and DNS is exclusively concentrated in the free shear layerand the wall side layer. More specifically, while this discrepancy grows continuously as Ha decreasesbut remains moderate in the shear layer, the two profiles brutally depart from one another in thewall side layer for Ha = 55, when small scale turbulent fluctuations appear. Since the profiles of (cid:104) u (cid:48) z (cid:105) / obtained with either PSM or DNS match everywhere else, even at for N (cid:46) .
19, it canbe concluded that PSM remains robust at predicting both global and local recirculation down to N (cid:46)
1, but breaks down when small scale turbulent fluctuations not driven by Ekman pumpingappear.Figure 13(d) shows the detailed evolution of E z . One can see that after the injection of theelectrical current density at t/t Ha = 0, part of the energy is converted to the kinetic energy alongthe magnetic field lines, which quickly results in a maximum E z . After that, E z tends to decreaseuntil a constant time-averaged value is reached. E z reaches a lower constant value when strongermagnetic fields are imposed for flows initialised in turbulent states. Since the entire secondary flow . Chen et al. DNS of MATUR Ha < E m z > t , < E ’ z > t
150 300 450 60010
Figure 13: Comparison of the rms (root mean square) of the fluctuations of the vertical velocity( u rms z ( r, θ, z )) along the radial direction near the side wall layer with Re = 15972 ( a ) near thecurrent injection position r = r e ( b ) near the side wall. Values were averaged along θ and z . ( c )Average part E mz = (cid:82) (cid:104) u z (cid:105) t dv and fluctuating part E (cid:48) z = (cid:82) ( u (cid:48) z ) dv of the energy in the z − velocitycomponent ( d ) Evolution of E z , where E z = E mz + E (cid:48) z . Re is fixed at 15972 for all 4 figures. . Chen et al. DNS of MATUR Figure 14: ( a ) Radial profiles of (cid:104) U θ (cid:105) θ,z . ( b ) Log-log plots of the free turbulent shear layerthickness δ f and the large scale vortex size δ vortex .transits through thin parallel layers, some residual axial flows are always observed within theseboundary layers. Therefore, E z always stabilises at a non-zero constant value in all the numericalcases, even though it is very small at Ha = 528 ( E z < − ). Interestingly, the walls have oppositeeffects on the energy in the third component, depending on whether it is driven by recirculationsor turbulence: Here the residual value of E z being mostly due to Ekman pumping, it is driven byfriction at the Hartmann walls. On the other hand, when turbulence freely decays in the presenceof solid Hartmann walls, the energy in the third component associated to random fluctuationsvanishes (in the sense that E z /E → t → ∞ , (Poth´erat & Kornet, 2015)). In unboundedor periodic domains, by contrast turbulence decays to a state where E z /E = 1 / One of the main purposes of the experimental study conducted by Messadek & Moreau (2002) wasto investigate the thickness of the free shear layer when turbulence is well established. Figure 14(a)shows the radial distribution of the mean azimuthal velocity, from which the thickness of the freeshear layer could be estimated. As mentioned in section 3.4, the free shear layer develops quicklyto an unstable state as soon as the current is injected. Therefore, a fraction of the momentum isconveyed from the annulus to the inner region and the boundary layer thickness increases visibly.The resulting entrainment of the fluid in the inner domain is characterized by a lower maximumvalue of U θ compared to that predicted by the laminar theory. For the annulus region, the (non-dimensional) value of (cid:104) u θ (cid:105) θ,z,t increases with Ha and decreases with Re due to different Ekmanpumping effects. To be more specific, at moderate Ha , either increasing Re or decreasing Ha enhances the recirculations, and thus the energy dissipation is expected to grow. Conversely, forthe inner region, the radial transfer of the momentum associated with recirculations set the fluid inrotation. Thus, the value of (cid:104) u θ (cid:105) θ,z,t in that region increases when recirculations become stronger.In other words, the thickness of the free shear layer increases when the Ekman pumping effectbecomes stronger, which is opposite to their effect on the laminar side layer ( Ha − / ).We now use the mean azimuthal velocity profiles within a wide parameter spaces of { Re, Ha } to determine the thickness of the shear layer δ f , which is defined as, δ f = ∆ U θ ( dU/dr ) max . (20) . Chen et al. DNS of MATUR ReHa δ S W × (ReHa ) × (ReHa ) Figure 15: Variations of the wall side layer thickness δ SW against parameter C = ReHa / (Poth´erat et al. , 2000) for Re = 15972 (red), Re = 31944 (green) and Re = 4792 (black). Fullcircles: boundary layer in laminar state, full squares: attached boundary layer with 3D turbulence,hollow squares: separated boundary layer.Here, ∆ U θ = U θmax − U θmin , U θmin = 0 and U θmax is the mean velocity value at the intersectionof the maximum slope and the mean velocity profile predicted by the laminar theory (7). Assuggested by Messadek & Moreau (2002), the layer thickness δ f depends on both Re and Ha according to the relation δ f = C f ( R ) /n . (21)As shown in figure 14(b), the best fit from our data yields C f (cid:39) .
44 and n (cid:39) . . Ha, Re ) = (792 , δ f ≈ .
64, agrees wellwith the experiment one, δ f ≈ .
61. This scaling is very different from the theoretical laminar one( δ f = Ha − / , independent of Re ) and reflects the role of two-dimensional inertia, measured by R ,in determining the thickness of the free shear layer. In addition, the size of large scales estimatedfrom the velocity fluctuations (see Poth´erat et al. (2005)) is also shown. Here, R L represents theReynolds number based on the velocity and the size of the large vortices and he size of the largevortices is estimated from the profiles of rms of azimuthal velocity fluctuations. As shown in figure14(b), their size follows a very similar scaling to the size of the boundary layer δ v (cid:39) . R . L witha maximum relative discrepancy to that law lower than 5 . C = ReHa − / , which (Poth´erat et al. , 2000) identified as the governing parameterwhen secondary flows dominate. The thickness δ SW was defined as the distance from the wall tothe point where the velocity magnitude reaches 90% of U SW θ , where U SW θ is the velocity value atthe intersection of the minimum slope and the mean azimuthal velocity profile. The results arereported on Figure 15. For low values of C , recirculations are weak and δ SW is expected to scaleas δ SW ∼ Ha − / . As C increases, recirculations become more prominent and δ SW approaches thetheoretical scaling of δ SW ∼ C − . The picture changes slightly before (in the sense of increasing C )three-dimensional turbulence appears in the boundary layer: the thickness suddenly increases tosettle on a larger scaling characterised by δ SW (cid:39) . ReHa − / ) − . . More simulations would . Chen et al. DNS of MATUR Figure 16: Variations of the global average angular momentum (cid:104) L lam (cid:105) t /L SM versus R . Regionwith R ≥
380 (Moresco & Alboussiere, 2004) is defined as regime where turbulence formed withinthe Hartmann layer. Region with 145 . ≤ R ≤ . R <
121 is corresponding to flows without separation of sidelayer.be needed to confirm that C remains the relevant parameter in this regime (though the continuedprominence of the secondary flows would suggest this may be the case) and to confirm the expo-nent of Re in this scaling. Interestingly, boundary layer separation has little visible impact on thescaling of δ SW , most likely because of the small ratio of the surface where it occurs to the totalsurface of the wall. For a first estimate the the global angular momentum, we note that most of the viscous and Jouledissipation takes place within the Hartmann layer for Q2D flows under a strong magnetic field.Since the most intense part of the flow occurs in the outer annulus where the driving force acts,the contribution of the inner region to the total angular momentum can be neglected to derive thesimple expression (14) from the theory of Sommeria & Moreau (1982) for the elementary case of asteady inertialess flow. Note that in this case, the Hartmann layers remain laminar and inertialess.This implies that the angular momentum varies linearly with I and is independent of B .The values of L lam obtained from the present numerical results (circle open symbols) are plottedin figure 16, along with the values of the angular momentum measured in MATUR (Messadek &Moreau, 2002) and predicted with PSM model (triangular open symbols). All the data reportedin this figure is normalised by the value L SM predicted with the theoretical expression (14), thedashed line corresponds to the theoretical prediction derived by Poth´erat & Schweitzer (2011) forturbulent Hartmann layers and the square symbols denote the experimental data. As shown infigure 16, the experimental values, the results of PSM model and our numerical values collapsewell into a single curve, which can be divided into three different zones demarcated by changesin slope around R (cid:39)
121 and R (cid:39) R < R ≥ . Chen et al. DNS of MATUR Ha τ H a τ Sh
800 160010 τ Ha ,DNS τ Sh ,DNS τ Sh ,PSM τ Sh ,laminar τ Ha ,laminar Figure 17: Space and time-averaged wall friction over the entire range of Ha investigated andfor ( Re = 15972): circle full symbols denote wall stress at Hartmann wall and circle openedsymbols denote wall stress on the side wall For comparison, theories for laminar flows (dashed anddashed-dot line) and results of the PSM model (gradient opened symbols) are shown too.(2011), which supposes that the Hartmann layers are turbulent. When 121 ≤ R < et al. (2000), who modifiedthe equation for L lam from the SM82 theory to account for the dissipation in the side layers, dL lam dt = F − S ν − L lam t Ha , (22)where F denotes the global electric forcing and S ν ( ∼ τ Sh ) denotes the viscous dissipation at theside wall layer. At a small forcing, the corresponding viscous effect on the angular momentumis negligible in comparison with the Hartmann friction. Thus, the DNS and PSM results areconsistent with the values predicted with the theoretical expression (14) when R < ≤ R < R ∈ [145 . , . τ Ha , τ Sh ) with Ha areshown in figure 17. The values of wall stress are different on the Hartmann walls and side walls.As shown in figure 17, all the shear stresses on the Hartmann walls collapse well into a singlecurve (laminar solution, τ Ha ∼ Ha/Re ), even when Ha = 55, which indicates that the Hartmannlayers remain laminar for all the cases considered in this work. However, for the side wall shearstresses, the values gradually depart from the laminar solution as Ha decreases ( Ha ≤ Ha (1056, or 1320, or 1584). When Ha ≥ Ha ≤ Ha . Chen et al. DNS of MATUR Figure 18: ( a ) Evolution of the global angular momentum with time. ( b ) Transient time obtainednumerically after switching on the forcing on a fluid at rest versus the non-dimensional parameter N / Ha / .decreases, which results in an increased squeezing of the side wall layer. Therefore, the valuesof τ Sh obtained from PSM simulations are higher than those predicted by the scaling law for astraight laminar Shercliff layer. However, DNS values of τ Sh are still significantly higher than theones from the PSM for Ha ≤ R ≤ .
4. When 121 ≤ R ≤ .
4, the discrepancy between the numericalresults and the theoretical results is associated with either the strong squeezing of the side walllayer or turbulence within side wall layer.Finally, we compare the time variations of the mean angular momentum for different Re and Ha , shown in figure 18(a). For the same Ha , one can see that the amplitude of the oscillations ofthe angular momentum in the quasi-steady-state increase with Re , due to the increasingly turbulentnature of the flow. Furthermore, the transient time ( t qs ) for the system to reach the quasi-steadystate from the fluid being at rest decreases with Re (we estimate this time by measuring the slopeof the L lam ( t ) curve near equilibrium in a log-log diagram, as in Poth´erat et al. (2005)). Thistoo is associated with global dissipation, which is significantly increased by the strong Ekmanpumping at higher Re , shortening the transient time. For a given Re , the variation trend with Ha is opposite, reflecting the damping of recirculations by the Lorentz force. This effect is quantified inFigure 18(b), which shows the variation of the transient time with the combined non-dimensionalparameter N / Ha / = Ha / /Re / , with a scaling t qs /t Ha (cid:39) . N / Ha / ) . , with amaximum relative error between the fitting curve and the numerical data of 9 . t qs /t Ha is governed by the same parameter as predicted by PSM, the scaling exponent of 0.467 stands muchlower than the value of 1, indicated in Poth´erat et al. (2005). Since a lower exponent is indicativeof a higher level of dissipation, this discrepancy may be attributed to the dissipation incurred by3D turbulence in the side layers that the PSM model cannot account for. . Chen et al. DNS of MATUR Figure 19: The results of case with Ha = 55 , Re = 15972. ( a ) Radial profile of correlation C (cid:48) . ( b )The azimuthal velocity fluctuation signals of six different points in the top and bottom Hartmannlayer variation with time. For easy identification, values of u θ at r = 6 .
85 ( r = 4 .
5) is displayed as u θ + 0 .
03 ( u θ − . c ) The azimuthal velocity fluctuations signals of three different points inthe side wall layer variation with time, along r = 10 . According to Poth´erat & Klein (2014) and Baker et al. (2015), the dimensionality of a struc-ture in a channel of gap a is determined by the ratio l z /a . Here l z is the momentum diffusionlength along B by the Lorentz force. This diffusion process takes place in a typical diffusion time τ D ( l ⊥ ) = τ j ( l z /l ⊥ ) , where τ j = ρ/σB is the Joule dissipation time. Hence, l z can be estimatedas (Sommeria & Moreau, 1982) l z ( l ⊥ ) = l ⊥ (cid:112) N ( l ⊥ ) , (23)where N ( l ⊥ , u ( l ⊥ )) = σB l ⊥ /ρu ( l ⊥ ), is a scale-dependent interaction parameter and u ( l ⊥ ) is thevelocity associated to a fluid structure of size of l ⊥ . l z /a (cid:29) a , and that the considered structure isconsequently quasi-two dimensional.However, the separation of the side wall may produce complex 3D structures in the case of Ha =55 , Re = 15972. Therefore, the correlation between V T ( r, θ, t ) and V B ( r, θ, t ) measured at locationswithin either Hartmann layers exactly aligned with the magnetic field lines (i.e. respectively at z = z (with z = 0 . < /Ha ) and z = 1 − z , but at the same coordinates ( r, θ )) is used toassess the three dimensionality of the flow (Klein & Poth´erat, 2010). Here, V T ( r, θ, t ) and V B ( r, θ, t )represent the azimuthal velocity fluctuations. The correlation function is defined as c (cid:48) ( r, θ ) = T i (cid:80) t =0 V T ( r, θ, t ) V B ( r, θ, t ) T i (cid:80) t =0 V B ( r, θ, t ) . (24)where T i is the duration of the recorded signals. Considering the symmetry of the problem, weshall analyse the radial dependence of this correlation through C (cid:48) ( r ) = (cid:104) c (cid:48) ( r, θ ) (cid:105) θ . For Ha =55 , Re = 15972, the correlation profile of C (cid:48) ( r ) in figure 19(a) indicates that the flow is veryclose to quasi-2D in most regions, where correlation C (cid:48) is almost unity. This is also demonstratedon figure 19(b), which shows that the instantaneous azimuthal velocity signals away from theshear layers near the top and bottom Hartmann layers are near identical for the same position( r, θ ). At r = r e , by contrast, the signals slightly differ, so the correlation C (cid:48) is slightly lower . Chen et al. DNS of MATUR f t Ha P S D f c t Ha f t Ha P S D f c t Ha f t Ha P S D bottomtop Figure 20: Power density spectra calculated from instantaneous azimuthal velocity signals acquiredat locations near the free shear layer inside the top (red line) and bottom (green line) Hartmannlayers for different magnetic interaction parameters. ( a ) Ha = 66 , Re = 15972. ( b ) Ha = 264 , Re =31944. ( c ) Ha = 264 , Re = 15972.than unity (0.832). At this location, the signals are mostly identical except for a slightly higheramplitude of the bottom signal and short burst where the signals are weakly correlated. Thedifference in amplitude reflects weak three-dimensionality, as defined by Klein & Poth´erat (2010),in the sense that flow near the top and bottom are identical in topology but differ in intensity.The short bursts, by contrast, indicate strong three-dimensionality where topologies are no moreidentical. The bursts correspond to the passage of coherent structures. Hence, the overall pictureis that while the free shear layer itself and the larger structures responsible for the lower frequencyoscillations are only weakly three-dimensional, the smaller coherent structures that navigate alongit can exhibit strong three-dimensionality.Figure 19(c) confirms that strongly three-dimensionality emerge within the side wall layer,where the signals recorded at all three monitored depths differ noticeably. This is consistent withthe distribution of azimuthal velocity within the side layer shown in figure 6(f), where the flowson different transverse plane are not topologically equivalent any more, because of the presence ofsmall scale three-dimensional turbulence there.To further characterise the scale-dependence of three-dimensionality near the electrodes, weanalyse the spectra obtained from the velocity signals near the top and bottom Hartmann layers,as shown in figure 20. For Ha = 66 , Re = 15972 and Ha = 264 , Re = 31944, pairs of energy spectraobtained near top and bottom Hartmann walls reveal that the higher frequencies carry significantlyless energy in the vicinity of the top wall than that near the bottom wall, which is a clear evidenceof three-dimensionality. By contrast, lower frequencies almost carry the same amount of energy.According to the theory of Sommeria & Moreau (1982) and the experiments of Baker et al. (2018),a cutoff scale k c (corresponding to f c here) separating the quasi-two-dimensional structures fromthe three-dimensional structures can be identified. Here, the velocity signals, averaged over theazimuthal direction, were taken within the free shear layer, so as to minimise the influence of theside walls, i.e. r = r e . For the details about the calculation of f c , the readers are referred toPoth´erat & Klein (2014). When N t decreases, f c decreases as three dimensionality contaminateslarger and larger scales, as shown in figure 20(a,b). For Ha = 264 and Re = 15972, fluctuationswithin the free shear layer are quasi-two dimensional over the entire spectral range (see figure20(c)). The variations of the azimuthally averaged cut-off frequency (cid:104) f c (cid:105) θ ( N t ) are represented infigure 21(a), and reveal that the variations of (cid:104) f c (cid:105) θ across all investigated cases collapse onto asingle curve, and therefore that (cid:104) f c (cid:105) θ is determined by N t with a scaling law (cid:104) f c (cid:105) θ (cid:39) . N . t . (25)This general law gives a clear estimate for the minimum frequency of vortices that are affected by3D inertial effects. Additionally, we obtain the minimum transverse wavenumber of 3D vortices, . Chen et al. DNS of MATUR Figure 21: ( a ) The azimuthal averages of f c , separating Q2D large structures from the small 3Dones, normalised by U /a . ( b ) Radial profiles of azimuthal averages of f c , when Ha = 55 , Re =15972 and Ha = 66 , Re = 15972. (cid:104) k c (cid:105) θ (cid:39) . N . t by applying Taylor’s hypothesis, taking advantage of the strong azimuthalflow component. This law is the first numerical confirmation of the original theoretical law givenby Sommeria & Moreau (1982), k c ∼ N / t , following Baker et al. (2018)’s recent experimentalconfirmation. It is also interesting to note that while this law applies to homogeneous and shearedturbulence alike, the corresponding scaling law for the cutoff frequency (25) differs from the powerlaw with exponent ∼ / (cid:104) f c (cid:105) θ along r ,shown in figure 21(b). (cid:104) f c (cid:105) θ clearly drops at the locus of the free shear layer and the side walllayer. On the other hand, it remains much higher outside of these regions. Hence, structures arequasi-two-dimensional over a greater range of scales, down to smaller ones outside the shear layersand three-dimensional turbulence appears concentrated in the side layer and to a slightly lowerextent, to the free shear layer. This suggests that three-dimensionality arises out of direct energytransfer from the mean shear flow to scales at the scale of the shear layers ( δ ) or lower. This issupported by the argument that the high level of energy imparted to them by the shear makesthe turnover time at this scale ( ∼ δ/ (cid:104) U θ ( r ) (cid:105) , with r = r e or r = R for the free shear and sidelayers respectively) significantly smaller than the two-dimensionalisation time ( ∼ ( σB /ρ )( a/δ ) ).Indeed, the ratio of former to the latter expresses as N ( δ/a ) , which is significantly smaller thanunity in both cases.Conversely, away from the shear layers, the mean flow does not inject energy into the smallscales. Since the recent experiments on MHD turbulence without a strong mean flow of Baker et al. (2018) suggest that even in the presence of moderate three-dimensionality, the energy cascadesupscale, the flow is dominated by larger vortices for which two-dimensionalisation is more efficient.Unlike in these experiments, however, the presence of vorticity streaks in the wake of these vorticessuggests that an additional transfer mechanism less favourable to large scales may be at play inthe outer region of MATUR. . Chen et al. DNS of MATUR Figure 22: The distribution of the average relative turbulent intensity (cid:104) u (cid:48) z u (cid:48) z (cid:105) t,z,θ with (cid:104) f c (cid:105) θ , when Re = 15972, for the different radial values of (cid:104) f c (cid:105) θ shown on figure 21(b). To understand the occurrence of the third component velocity fluctuations not driven by globalEkman pumping, we show instantaneous distributions of axial velocities in the cross-section θ =0 (see figure 7(a)). One can see that the turbulent fluctuations are localized within the layernear the side wall. Moreover, this three-dimensionality contaminates the entire height of thevessel, from the bottom to the top of the Hartmann layers. As conveyed in figure 7(b,c), thedistributions of instantaneous vertical velocity (including intensity and topological structure) aredifferent on the plane near the top wall ( z = 0 .
8) and near the bottom wall ( z = 0 . (cid:104) u (cid:48) z (cid:105) t,z,θ and the cutoff frequency, which provides a measure of three-dimensionalityacross the turbulent spectrum. These quantities for several cases are plotted on figure 22. Thecollapse of the data into a single curve shows that the degree of the two-dimensionality of theenergy spectrum is tightly linked to the amount of energy in the third component. (cid:104) u (cid:48) z (cid:105) t,z,θ is therefore solely determined by the true interaction parameter N t and follows a simple powerlaw, i.e. (cid:104) u (cid:48) z (cid:105) t,z,θ = 2 × − (cid:104) f c (cid:105) − . θ . The maximum relative error between the fitting curveand the numerical data is lower than 2 . f c (cid:39) . N . t , and (cid:104) u (cid:48) z (cid:105) t,z,θ (cid:39) . × − N − . t . It isnoteworthy that this scaling implies a different dependence on Ha (cid:104) u (cid:48) z (cid:105) t,z,θ ∼ Ha − . than thatfor E (cid:48) z = πR a (cid:104)(cid:104) u (cid:48) z (cid:105) t,z,θ (cid:105) r ∼ Ha − . . Since the two quantities only differ though radial averaging,the difference can be understood by noticing that if fluctuations are confined to a thin radial regionof thickness δ D then E (cid:48) z ∼ πR a ( δ D /R ) (cid:104) u (cid:48) z (cid:105) t,z,θ,δ D , where the subscript δ D indicates averagingover that region. Under this assumption, it would follow that δ D /R ∼ Ha − . (finding the scalingwith Re would require extra simulations over a wider range of its values). The corresponding regionis much thinner than the free shear layer and the wall side layer (see section 3.6). However theradial profiles of vertical velocity fluctuations (Fig. 13) suggest that a sharp peak of verticalvelocity fluctuations develops within the later, that would explain this scaling. This suggests thatthe dominant contribution to the three-component turbulence in the regimes explored in this paperarises out of a very thin region of thickness within the wall side layer. . Chen et al. DNS of MATUR The present study reports three-dimensional direct numerical simulations of electrically drivenMHD turbulent shear flows in the MATUR experiment (Messadek & Moreau, 2002). The numericalresults, which are obtained in a configuration where current is injected far from side walls (at aradius of r e = 5 .
4) and in regimes where the Hartmann layers remain laminar (18 . ≤ R ≤ . R = Re/Ha was varied: from R (cid:39) R ≥ .
2, that layer separates from the wall, most likely underthe influence fluctuations induced by the large vortices in the outer region. The thickness δ SW ofthe wall side layers follows a sequence following these changes of regime: in the laminar regime, itconverges to the asymptotic scaling δ SW ∼ ( ReHa / ) − theoretically predicted by Poth´erat et al. (2000), but becomes much thicker with a scaling of δ SW ∼ ( ReHa − / ) − . following the onset ofthree-dimensional turbulence.In addition,the energy spectra exhibit a significant dependence on R : for 18 . ≤ R ≤ . E ( f ) ∼ f − power law. For R ≥ E ( f ) ∼ f − / and high frequency ranges E ( f ) ∼ f − , similar to the split between inverseenergy and direct enstrophy cascades found in quasi-2D MHD flows by Sommeria (1986). Thedynamics of the free shear layer seen in the experiments was also recovered as DNS showed thatthe thickness of the free shear layer varies nearly as the vortex size does (scaling as R / . and R / . L , respectively).Detailed analysis of the secondary flow confirmed the phenomenology identified by Poth´erat et al. (2000, 2005) whereby for N (cid:46)
1, the global recirculation associated to the main azimuthalflow induces a flux of angular momentum towards the wall side layer. The ensuing thinning ofthat layer is responsible for an increased dissipation. The intensity of the recirculation matchesclosely the PSM prediction. At statistical equilibrium, the energy associated to average axial flowcomponent was found to scale as E mz ∼ Ha − . both in PSM and the DNS, confirming that it isdriven by this main recirculation. The energy associated to the fluctuating part of this componentwas on the other hand underestimated by PSM compared to the DNS for N (cid:46)
1. The discrepancyoriginates in the small scale turbulence produced in the wall side layer. This effect was foundto incur a significant increase in wall shear stress and in turn a reduction in the global angularmomentum, a phenomenon that previous theories could not capture. The extra dissipation is alsoa possible cause for the shortened transient time observed when the flow is initiated at rest.A major benefit of the 3D DNS was to afford a detailed scrutiny of the three-dimensional effects.The main source of three-dimensionality was found in the small-scale turbulence fed by the meanshear in the free shear layer and the wall side layer. Outside of these regions, turbulent spectrastill exhibit a high frequency range of three-dimensional structures and a low frequency range ofquasi-two dimensional ones, as predicted by Sommeria & Moreau (1982) and found in turbulencedriven by a crystal of vortices (Klein & Poth´erat, 2010). As in this case, the frequency separatingthese two ranges scales with the true interaction parameter N t , albeit with a different exponent . Chen et al. DNS of MATUR as f ∼ N . instead of f ∼ N / . The difference is due to the presence of a strong backgroundflow and when converted into wavenumbers, both experiments exhibit the same cutoff scaling of k c ∼ N / t (Baker et al. , 2018), as predicted by Sommeria & Moreau (1982). This results showsthat the existence of a cutoff wavelength separating three and quasi-two dimensional turbulentfluctuations extends to sheared turbulence. As such, this result and the associated scaling can beexpected to hold in a wider class of flows including duct flows.Simultaneous access to all three velocity components further enabled us to establish a linkbetween the local dimensionality of the turbulence (measured by the cutoff frequency f c ) and itscomponentality, measured by the energy in the axial component of the velocity fluctuations. Whilethe latter decreases monotonically with the former, and both are controlled by the true interactionparameter through simple power laws (cid:104) u (cid:48) z (cid:105) t,z,θ ∼ . N − . t and (cid:104) f c (cid:105) θ (cid:39) . N . t . Unlike forthe dimensionality scaling ( f c ), no prediction exists for the dimensionality scaling ( (cid:104) u (cid:48) z (cid:105) t,z,θ ), sothis raises the question of its applicability to other types of quasi-static MHD turbulence, beyondshear flow turbulence or even beyond this particular experiment. The question is all the morerelevant as in the regime explored here, three-component turbulence was found almost exclusivelywithin a very thin region of the wall side layer. As such, a further step in understanding the linkbetween componentality and dimensionality could target flows where three-dimensional turbulenceis more broadly distributed.The authors acknowledge the support from NSFC under Grants (cid:93) (cid:93) (cid:93) XDB22040201, (cid:93)
QYZDJ-SSW-SLH014. The authors also greatly appreciatethe anonymous reviewers for their comments, which significantly helped to improve the manuscript.Declaration of Interests. The authors report no conflict of interest.
References
Alboussi`ere, T, Uspenski, V & Moreau, R
Experimental Thermal and Fluid Science , 19–24. Alexakis, A. & Biferale, L.
Physics Reports , 1–101.
Baker, Nathaniel T, Poth´erat, Alban & Davoust, Laurent
Journal of Fluid Mechanics , 325–350.
Baker, Nathaniel T, Poth´erat, Alban, Davoust, Laurent & Debray, Franc¸ois
Physical review letters (22), 224502.
Davidson, PA
Journal of Fluid Mechanics , 123–150.
Davidson, PA & Poth´erat, A
European Journalof Mechanics-B/Fluids (5), 545–559. Eckert, S, Gerbeth, G, Witke, W & Langenbrunner, H
International journal ofheat and fluid flow (3), 358–364. Klein, R & Poth´erat, Alban
Physical review letters (3), 034502.
Kljukin, A, A & Kolesnikov, Yu. B
Liquid Metal Magnetohydrodynamics , , vol. 10.Kluwer. . Chen et al. DNS of MATUR
Kobayashi, Hiromichi
Physics of Fluids (4),045107. Kobayashi, Hiromichi
Physics of Fluids (1), 015102. Kolesnikov, Yu B & Tsinober, AB
Fluid Dynamics (4), 621–624. Messadek, Karim & Moreau, Rene
Journal of Fluid Mechanics , 137–159.
Moffatt, H. K.
Journal ofFluid Mechanics (3), 571–592. Moresco, P & Alboussi`ere, T
European Journal of Mechanics-B/Fluids (4), 345–353. Moresco, Pablo & Alboussiere, Thierry
Journal of Fluid Mechanics , 167–181.
Ni, Ming-Jiu, Munipalli, Ramakanth, Huang, Peter, Morley, Neil B & Abdou, Mo-hamed A
Journal of ComputationalPhysics (1), 205–228.
Poth´erat, Alban & Dymkou, Vitali Rm MHDturbulence based on the least dissipative modes.
Journal of Fluid Mechanics , 174–197.
Poth´erat, Alban & Klein, Rico
Journal of Fluid Mechanics , 168–205.
Poth´erat, A. & Kornet, K. Rm . Journal of Fluid Mechanics , 605–636.
Poth´erat, Alban & Schweitzer, Jean-Philippe
Physics of Fluids (5), 055108. Poth´erat, A, Sommeria, J & Moreau, R
Journal of Fluid Mechanics , 75–100.
Poth´erat, Alban, Sommeria, Jo¨el & Moreau, Ren´e
Journal of FluidMechanics , 115–143.
Roberts, Paul Harry
An introduction to magnetohydrodynamics , , vol. 6. Longmans Lon-don.
Schumann, U.
Journal of Fluid Mechanics (1), 31–58. Sommeria, Joel
Journal of Fluid Mechanics , 139–168.
Sommeria, Jo¨eL & Moreau, Ren´e
Journal of Fluid Mechanics , 507–518.
Stelzer, Zacharias, C´ebron, David, Miralles, Sophie, Vantieghem, Stijn, Noir,J´erˆome, Scarfe, Peter & Jackson, Andrew a Experimental and numerical studyof electrically driven magnetohydrodynamic flow in a modified cylindrical annulus. i. base flow.
Physics of Fluids (7), 077101. Stelzer, Zacharias, Miralles, Sophie, C´ebron, David, Noir, J´erˆome, Vantieghem,Stijn & Jackson, Andrew b Experimental and numerical study of electrically drivenmagnetohydrodynamic flow in a modified cylindrical annulus. ii. instabilities.
Physics of Fluids (8), 084108. Tabeling, P & Chabrerie, JP . Chen et al. DNS of MATUR verse pressure gradient.
Physics of Fluids (3), 406–412. Zhao, Yurong & Zikanov, Oleg
Journal of Fluid Mechanics692