Jules Delacroix

Laurent Davoust
email: laurent.davoust@simap.grenoble-inp.fr.fig

On the role of surface rheology in a magnetohydrodynamic swirling Ćow

The original coupling between the surface rheology of a liquid metal surface and a supporting annular Couette magnetohydrodynamic (MHD) flow is theoretically and numerically investigated in this paper, in the general layout of the classical annular viscometer, as developed by Mannheimer and Schechter [ªAn improved apparatus and analysis for surface rheological measurements,º J. Colloid Interface Sci. 32, 195 ±211 (1970)]. The purely hydrodynamic interplay between the main azimuthal flow (induced by a rotating floor) and the secondary overturning flow (generated by centrifugation) is found to be strongly affected by both surface viscous shear and sur-face viscous dilatation. When the flow is subjected to an outer vertical magnetic field, the impact of varying interface boundary conditions (through the surface shear and dilatational viscosities) at the gas/liquid interface profoundly alters the MHD flow topology. Particularly, when centrifugation competes with electromagnetic effects, advection of the main flow by the secondary flow is proved to affect significantly the core MHD flow, leading to a variety of atypical MHD flow patterns.

I. INTRODUCTION

Modeling of magnetohydrodynamic (MHD) two-phase flows has recently become an issue of major interest, given the numerous industrial applications potentially affected:

• metallurgy, with the stirring of molten metal by bubble plumes in crucibles (Haverkort and Peeters [START_REF] Haverkort | ªMagnetohydrodynamic effects on insulating bubbles and inclusions in the continuous casting of steel,º[END_REF] ) or with the entrainment of floating unwetted particles inside a liquid metal bath under a rotating magnetic field, see, e.g., Bojarevics et al.; 2 • microelectronics, where oxidation may affect electrically driven liquid metal processes by conferring non-wetting properties on the fluidÐLiu et al., [START_REF] Liu | ªCharacterization of nontoxic liquid±metal alloy galinstan for applications in microdevices[END_REF] Liu et al., [START_REF] Liu | ªNon-toxic liquid metal microstrip resonators[END_REF] or Sen and Kim; 5 • crystal growth processes and capillarity such as the Czochralski and floating zone methodsÐ see Duffar; 6 • nuclear fusion reactors, with the breeder blanket based cooling loop issues (Helium bubbles in the liquid metal, see, for instance, Morley et al. [START_REF] Morley | ªLiquid magnetohydrodynamics: Recent progress and future directions for fusion[END_REF] or Li et al. [START_REF] Li | ªMHD effect on flow structures and heat transfer characteristics of liquid metal-gas annular flow in a vertical pipe[END_REF] ) or MHD corrosion studies for dual-coolant lead-lithium blanketsÐSmolentsev et al. [START_REF] Smolentsev | ªNumerical study of corrosion of ferritic/martensitic steel in the flowing PbLi with and without a magnetic field[END_REF] Focusing on liquid/gas flows subjected to an outer magnetic field, a typical layout of a stratified two-phase MHD flow is displayed in Fig. 1. In this figure, we can first point out that due to internal motion ⃗ v inside the electro-conductive fluid, in interaction with the outer magnetic field ⃗ B ext , there is a magnetohydrodynamic flow. Moreover, if a gaseous phase is added to form a heterogeneous liquid/gas medium, complex interactions between the two phases are expected to influence the overall MHD flow. More particularly, the role that may be played by the gas/liquid interface is a crucial issue. If the focus is placed on the liquid surface, particles stemming from physico-chemical processes (corrosion, abrasion, oxidation, etc.) can adsorb at this interface. Depending on its level of saturation, the rheological behavior of the interface can be greatly modified, e.g., evolving from a sliding to a non-sliding boundary condition (BC). Consequently, a major impact is expected on the topology of the underlying MHD flow. This particular coupling between MHD and surface rheology needs to be investigated.

To our knowledge, little is actually known about the surface rheology of liquid metals when they are progressively oxidized. The viscoelastic properties of liquid metals have been experimentally investigated by researchers, including Kolevzon et al., [START_REF] Kolevzon | ªLight-scattering spectroscopy of a liquid gallium surface[END_REF][START_REF] Kolevzon | ªLight-scattering study of the mercury liquid-vapor interface[END_REF] Dickey et al., [START_REF] Dickey | ªEutectic gallium-indium (EGaIn): A liquid metal alloy for the formation of stable structures in microchannels at room temperature[END_REF] Larsen et al., [START_REF] Larsen | ªViscoelastic properties of oxide-coated liquid metals[END_REF] Liu et al., [START_REF] Liu | ªCharacterization of nontoxic liquid±metal alloy galinstan for applications in microdevices[END_REF] or Doudrick et al. [START_REF] Doudrick | ªDifferent shades of oxide: From nanoscale wetting mechanisms to contact printing of gallium-based liquid metals[END_REF] These authors have highlighted the preponderant role of the interface level of oxidation on different mechanical properties, such as the response to imposed elastic stresses, surface tension, and contact angles. However, none of these results address the issue of rheological transport, and they are not coupled with MHD.

On the other hand, the MHD of single-phase laminar flows has been extensively studied for many years. The flow can be either confined (see Shercliff [START_REF] Shercliff | ªSteady motion of conducting fluids in pipes under transverse magnetic field[END_REF] ) or may have a free surfaceÐfor a Couette flow with a free surface, see, e.g., Lehnert. [START_REF] Lehnert | ªAn instability of laminar flow of mercury caused by an external magnetic field[END_REF] However, the fundamental issue of varying boundary conditions would seem to have only been considered from an electrical point of view. Thus, in the case of duct flows, the walls can have infinite electrical conductivity (see Shercliff [START_REF] Shercliff | ªSteady motion of conducting fluids in pipes under transverse magnetic field[END_REF] ), no conductivity (Moreau 17 ), mixed infinite and vanishing conductivities (Hunt and Stewartson [START_REF] Hunt | ªMagnetohydrodynamic flows in rectangular ducts. II,º[END_REF] ), or arbitrary conductivity (Tabeling and Chabrerie [START_REF] Tabeling | ªMagnetohydrodynamic secondary flows at high hartmann numbers[END_REF] ). The electrical influence of the walls, governed by the ratio of bulk and wall electrical conductivities σ/σ w , completely modifies the electric circuit and results in a major impact on the topology of the MHD flow.

In this study, the same kind of general approach for the boundary conditions is taken, but this time from a mechanical point of view. In this respect, for flows including liquid/gas interfaces, the competition between bulk MHD and surface rheology (driven by the contamination rate of the liquid surface) may lead to a mechanically varying boundary condition, suspected to affect considerably the overall MHD flow (see, e.g., the oscillatory flow patterns caused by surface oxidation described by Zhang et al., [START_REF] Zhang | ªIntermittent behavior caused by surface oxidation in a liquid metal flow driven by a rotating magnetic field[END_REF] in the case of a liquid metal subjected to a rotating magnetic field). To deal with the practical conditions of molten metals in metallurgy industry, we need to investigate surface mechanics separately from bulk flow, especially when the free surface is gradually oxidizing. The mechanical coupling can therefore be considered between a liquid surface and the underlying bulk, through the emergence of dilatational and shear surface viscosities. These surface viscosities are expected to represent the effects of changing surface boundary conditions on the steady flow of a supporting subphase.

II. OUTLINES

A. The annular MHD viscometer: Description and motivation

To enhance the mechanical coupling between MHD and surface rheology for a liquid metal, we choose to focus on a particular physical layout, i.e., the annular MHD viscometer, shown in Fig. 2.

The annular viscometer, developed by Mannheimer and Schechter, [START_REF] Mannheimer | ªAn improved apparatus and analysis for surface rheological measurements[END_REF] has become a reference solution for determining the rheological properties of a given fluidÐsee, e.g., Hirsa and Lopez, [START_REF] Hirsa | ªDetermination of surface shear viscosity via deep-channel flow with inertia[END_REF][START_REF] Lopez | ªNon-Newtonian behavior of an insoluble monolayer: Effects of inertia[END_REF]

B. Goals of the present work

The first goal of this numerical study is to extend significantly the scope of a companion paper, [START_REF] Delacroix | ªElectrical activity of the hartmann layers relative to surface viscous shearing in an annular magnetohydrodynamic flow[END_REF] focused on the 2-D analytical study of the annular MHD viscometer. This benchmarking analytical study is based on the assumptions Ha ≫ 1 and Ha ≫ Re, where Ha and Re are the Hartmann and Reynolds numbers, respectively. This allowed the authors to ignore the inertial effects and to highlight the competitive effects between surface viscous shearing and a strong transverse uniform magnetic field. Typically, it has been shown that surface rheology actually monitors the electrical activation of Hartmann layers.

For the present study, no particular assumptions have been made concerning the values of Re and Ha (except that Re ≤ 10 4 in order to avoid turbulence issues, see Sec. III A). Consequently, a swirling flow occurs when Ha < Re, leading to inertial centrifugal effects. As a consequence, a new rheological parameter is solicited at the interface, linked to surface viscous dilatation, which might affect core flow. Furthermore, interacting with the outer magnetic field, the centrifugal effects bring an original MHD flow into play, which in turn competes with surface rheology. Their coupling is investigated in this paper.

Beyond the goal of highlighting the physical mechanisms of this original coupling, a long-term goal involves measurement. First, thanks to the following numerical calculations, combined with an experimental set-up that is currently being developed, measurement of surface velocity would allow access to the bulk MHD. Second, we would be able to obtain information about oxidized surface rheology, since it would be possible to determine the rheological properties of liquid metal surfaces, for which there is only little information available; see, e.g., Liu et al. [START_REF] Liu | ªCharacterization of nontoxic liquid±metal alloy galinstan for applications in microdevices[END_REF] Finally, depending on the externally applied magnetic field, the imposed angular speed, and the O 2 -rate of the surrounding atmosphere, the model could be used to determine distinct experimental working conditions.

C. Overall method

First, the governing equations as well as the BCs are highlighted, along with the supporting physical phenomena. The coupling between hydrodynamics, electromagnetism, and surface rheology is particularly enhanced, and the scaling parameters of the problem are defined.

Classically, the finite-element method (FEM) is used to discretize Maxwell equations in electromagnetism. However, this is not the case in fluid mechanics, as, contrary to the finite-volume method (FVM), the discrete approximation is a priori not conservative. Consequently, if the FEM method is used to tackle a MHD problem, particular care is required to ensure a conservative solution.

Therefore, numerical modeling is systematically benchmarked with asymptotic cases. For the classical annular viscometer layout (i.e., with no applied magnetic field), the first benchmark takes inertia and surface viscous shearing into account (Hirsa et al. [START_REF] Hirsa | ªDetermination of surface shear viscosity via deep-channel flow with inertia[END_REF] ). The purely hydrodynamic results are then extrapolated to the case where surface viscous dilatation significantly affects the bulk flow. Then, the outer magnetic field is added, and the numerical results are benchmarked with the supporting 2-D analytical study, [START_REF] Delacroix | ªElectrical activity of the hartmann layers relative to surface viscous shearing in an annular magnetohydrodynamic flow[END_REF] which highlights the interactions between surface viscous shearing, electromagnetism, and creeping flow. Once this benchmarking case has been secured, a scaling law is defined to emphasize the interactions between the MHD bulk flow (with inertial effects) and surface rheology (including both surface shearing and dilatation). Their most salient features are then displayed and discussed.

III. PHYSICAL MODELING

A. Notations and assumptions

With respect to purely hydrodynamic assumptions, the annular shear flow considered is assumed to be a permanent, axisymmetric (∂/∂θ = 0), incompressible, and viscous Newtonian flow with no temperature dependence, so that the bulk physical properties of Galinstan are considered to be constant. The flow is laminar with Re = ρΩr 2 o /η, where η is the dynamic viscosity and ρ is the density of Galinstan (indicative values: η = 2.14 × 10 -3 Pa s and ρ = 6.36 × 10 3 kg m -3 at 20 • C, see Plevachuk et al. [START_REF] Plevachuk | ªThermophysical properties of the liquid Ga-In-Sn eutectic alloy[END_REF] ). In this paper, and unlike in the supporting analytical study [START_REF] Delacroix | ªElectrical activity of the hartmann layers relative to surface viscous shearing in an annular magnetohydrodynamic flow[END_REF] where Re is set to be sufficiently small, Re is increased up to the value Re ≤ 10 4 in order to study inertial effects. This relatively high value for the Reynolds number might question the assumption of an axisymmetric laminar flow. However, the Reynolds number defined in the present study (following, e.g., Hirsa et al., [START_REF] Hirsa | ªDetermination of surface shear viscosity via deep-channel flow with inertia[END_REF] Lopez et al., [START_REF] Lopez | ªNon-Newtonian behavior of an insoluble monolayer: Effects of inertia[END_REF] Davoust et al. [START_REF] Davoust | ªFlow-induced melting of condensed domains within a dispersed Langmuir film[END_REF] ) is not really physically consistent, as its definition, only based on the outer radius of the annular channel, does not involve the aspect ratio, h/r o . Our definition clearly overestimates the local Reynolds number which remains to be estimated a posteriori. Therefore, even a value as large as Re = 10 4 remains consistent here with the assumption of a laminar flow (see, e.g., Haddadi and Poncet, [START_REF] Haddadi | ªTurbulence modeling of torsional Couette flows[END_REF] Serre et al., [START_REF] Serre | ªAnnular and spiral patterns in flows between rotating and stationary discs[END_REF] and Launder et al. [START_REF] Launder | ªLaminar, transitional, and turbulent flows in rotor-stator cavities,º[END_REF] for considerations on the onset of turbulence).

Moreover, the quasi-static approximation is made, which consists in ignoring the displacement current when compared with the conduction current. Also, the magnetic Reynolds number Rm = µσr o Ωh is considered to be negligibly small, where µ and σ are the magnetic permeability and the electrical conductivity of the liquid metal, respectively (indicative values: µ = 4π × 10 -7 H m -1 , σ = 3.29 × 10 6 S m -1 , see Plevachuk et al. [START_REF] Plevachuk | ªThermophysical properties of the liquid Ga-In-Sn eutectic alloy[END_REF] ). This allows us to form the classical weakly coupled MHD model, which results in the fact that the electromotive current is approximated as

⃗ v × ⃗ B ∼ ⃗ v × ⃗ B 0
, where ⃗ v is the flow velocity, and that the Lorentz force is approximated as ⃗ j × ⃗ B ∼ ⃗ j × ⃗ B 0 , where ⃗ j is the electric current density. When normalizing the MHD equations, the Hartmann number, Ha = B 0 h  σ/η, is highlighted. Finally, several assumptions are made to describe the behavior of the liquid/gas interface. First, the interface is considered to be flat, with a capillary length l c =  γ/ρg ≪ r o -r i , where g is the gravity and γ is the surface tension of Galinstan (indicative value: γ = 0.534 N m -1 , see Liu et al. [START_REF] Liu | ªCharacterization of nontoxic liquid±metal alloy galinstan for applications in microdevices[END_REF] ). Furthermore, the flow of the liquid surface is modeled through a surface momentum conservation (jump momentum balance (JMB)), in which use is made of two rheological parameters: surface shear viscosity, η S , and surface dilatational viscosity, κ S Ðsee Eqs. ( 14) and (15). It is assumed that these parameters only depend on the O 2 -rate of the surrounding atmosphere (that directly rules interface oxidation), which implicitly requires thermodynamic equilibrium. In particular, and even if the swirling flow were able to induce a radial segregation along the interface, the hypothesis of a uniform oxidation is put forward. This yields that no radial dependence of η S or κ S is taken into account in this paper. The focus is placed here on the coupling with MHD.

Due to centrifugal effects, a swirling flow is generated. A distinction is made between the (main) azimuthal flow v θ and the (secondary) meridian flow (v r , v z ) (also referred to as overturning flow in this paper). When the electroconductive liquid is subjected to an outer magnetic field, electric current densities are induced. Once again, the (main) components ( j r , j z ) (due to the main azimuthal flow) are distinguished from the (secondary) component j θ (due to the secondary meridian flow). Note that even if the velocity has three non-zero components, the physical problem cannot be (strictly speaking) qualified as three-dimensional, since the independent variables in the entire description of the problem are only two, r and z. Consequently, the present study is hereinafter referred as 2-D 1/2. The investigation of the interaction between the main and secondary MHD flows is worthwhile, because, depending on the relative value of Ha and Re, the flow switches from 2-D to 2-D 1/2 topology. This also greatly affects the interface dynamics, because if Ha ≫ Re, the meridian flow vanishes, and only the azimuthal component of JMB is involved. Consequently, the only relevant rheological parameter is surface shear viscosity η S , the value of which may dramatically affect the main MHD flow. [START_REF] Delacroix | ªElectrical activity of the hartmann layers relative to surface viscous shearing in an annular magnetohydrodynamic flow[END_REF] On the contrary, if the inertial effects are significant, the radial component of JMB is brought into play. Thus, the interface tends to stretch, and the surface dilatational viscosity κ S can also affect the surface dynamics. Surface rheology may then modify in turn both the main and the secondary MHD flows.

B. Geometry, governing equations, boundary conditions

Geometry

For symmetry reasons, we can only focus on a cross section of the annular channel. Note that the out-of plane components of MHD quantities must be solved, due to the swirling flow. The geometry of the numerical problem, along with the conditions implemented at the boundaries of the calculation domains, is given in Fig. 3.

Bulk Ćow

a. Electromagnetism. A potential formulation is used to describe the electromagnetic part of the problem. The basic equations derive from Maxwell's equations, and in particular Ampère's law and the continuity equation, completed by Ohm's law generalized to moving matter. Using now the definition of the electric potential φ and the magnetic vector potential ⃗ A, with respect to the electric field ⃗ E and magnetic induction ⃗ B, at steady state,

⃗ E = - ---→ grad (φ) , ⃗ B = --→ curl ( ⃗ A ) . (1)
Ampère's law can be rewritten as

--→ curl ( --→ curl ( ⃗ A ) ) -µσ⃗ v × ( --→ curl ( ⃗ A ) ) + µσ ---→ grad (φ) = 0, (2)
and the continuity equation as div

( µσ⃗ v × ( --→ curl ( ⃗ A ) ) -µσ ---→ grad (φ) ) = 0, (3)
completed by the Coulomb gauge to define uniquely ⃗ A, i.e., div 2) and (3) are the solved electromagnetism equations.

( ⃗ A ) = 0. Equations (
The electromagnetic boundary conditions first consist of an externally applied constant axial magnetic field, imposed through an ªinfiniteº box surrounding the fluid area: ⃗ B ∞ = B 0 ⃗ e z . The assumption Rm ≪ 1 yields: B z ⃗ e z ∼ B 0 ⃗ e z , throughout the entire computational domain. The result is, in terms of the magnetic vector potential, 1 r

∂r A θ ∂r = B 0 . (4)
The second electromagnetic boundary condition is the electrical insulation at the liquid metal boundaries (side-walls, rotating floor, and liquid/gas interface): ⃗ j • ⃗ n fluid walls = 0. Using Ohm's law and the low Rm assumption, this condition is written with respect to electric potential,

σ ( - ---→ grad(φ) + ⃗ v × ⃗ B 0 ) fluid walls • ⃗ n = 0, ( 5
)
where ⃗ n is the unit normal vector at the considered boundaryÐsee Fig. 3.

b. Fluid mechanics. A primitive pressure-velocity formulation is used, based on the complete set of Navier-Stokes equations,

div (⃗ v) = 0, (6)
ρ

( ⃗ v • ---→ grad ) ⃗ v = div ( -pI + η ( grad (⃗ v) + grad ⊤ (⃗ v) )) + ⃗ F. (7)
The coupling term ⃗ F is the Lorentz force, defined as ⃗ F = ⃗ j × ⃗ B, which can be explicitly written as follows, using Ohm's law and the approximation Rm ≪ 1:

⃗ j = σ ( - ---→ grad(φ) + ⃗ v × ⃗ B 0 ) = σ            - ∂φ ∂r + v θ B 0 -v r B 0 - ∂φ ∂z            , (8)
⃗ F = ⃗ j × ⃗ B 0 = σ           -v r B 2 0 B 0 ( ∂φ ∂r -v θ B 0 ) 0           . ( 9
)
Let us now examine the hydrodynamic boundary conditions, starting with the no-slip boundary condition at the motionless side-walls,

⃗ v(r = r i , z) = ⃗ v(r = r o , z) = ⃗ 0. ( 10
)
At the inner and outer side-walls, a no-slip boundary condition is normally imposed for velocity. However, in order to circumvent the boundary discontinuity between the inner or outer side-walls and the rotating floor, two matching functions are introduced for the azimuthal component of velocity (see Fig. 3). These functions apply along segments of typical length δ ≪ r i , r o (indicative value: δ = 2.5 × 10 -4 m) at the inner and outer walls of the annular channel cross section. With these matching functions, a Couette-like profile for velocity is classically assumed at the matching segments, [START_REF] Pintar | ªThe measurement of surface viscosity[END_REF] which leads to

v θ (r = r i , z) =        r i Ω ( 1 - z δ ) if z ∈ [0, δ] , 0 if z ∈ ]δ, h] , (11)
for the inner wall and to

v θ (r = r o , z) =        r o Ω ( 1 - z δ ) if z ∈ [0, δ] , 0 if z ∈ ]δ, h] , (12)
for the outer wall. Thus, velocity is imposed such that lubrication conditions in the gaps are taken into account, the continuity of v θ being warranted at the matching points z = δ at both walls. These lubrication conditions are not simply a numerical technique. They stem from a true experimental layout, where this gap concretely exists between a rotating dish (rotating floor) and an inert cover (side walls), as described Fig. 4. The last remaining velocity boundary conditions at the liquid gas/interface are given by the surface rheology equations, governing surface velocity ⃗ v S dynamicsÐ see Eqs. ( 14) and (15), displayed for a purely hydrodynamic flow, meaning that N = 0, regardless of the value of the Reynolds number. Moreover, the Couette-like layout of the problem, with an imposed velocity at the rotating floor, makes us tend towards the (Ha, Re) description, since the boundary condition at the bottom of the channel is directly linked to Re through the angular speed Ω. Consequently, by monitoring the Reynolds number, different dynamic configurations are described.

⃗ v(r, z = h) = ⃗ v S (r) = v r S (r)⃗ e r + v θS (r)⃗ e θ . ( 13
The bulk dimensionless quantities of interest, superscripted ⋆ , are then defined as follows:

• the radial coordinate r ⋆ = r/r o ;

• the axial coordinate z ⋆ = z/h;

• the velocity ⃗ v ⋆ = ⃗ v/V , where V = r o Ω. The azimuthal component v ⋆ θ and the meridian compo- nents (v ⋆ r , v ⋆ z
) allow for analysis of the main and secondary MHD bulk flows, respectively;

• the magnetic induction: b ⋆ r = b r /B, b ⋆ θ = b θ /B, where B = µV √ ση and B ⋆ z = B z /B 0 ;
• and finally, for the electric current densities: ⃗ j ⋆ = ⃗ j/J, where J = σB 0 V . The meridian components ( j ⋆ r , j ⋆ z ) and the azimuthal component j ⋆ θ allow for analysis of the main and secondary MHD bulk flows, respectively.

Surface Ćow

a. Surface rheology. Boundary condition (13) brings a new unknown into play, namely, in-plane surface velocity ⃗ v S . This stands as the first coupling term of the two-way coupling between the MHD bulk flow and surface equations. The latter conditions can be derived from the momentum balance written on an elementary heterogeneous volume straddling a liquid surface of zero thickness, in accordance with the Gibbs approach (see Slattery et al. [START_REF] Slattery | Interfacial Transport Phenomena[END_REF] ). This allows us to circumvent the issue of the strong and often unknown variability of physical variables across a diffuse interface. In return, the surface ªexcessº quantities, e.g., surface viscosities, must be introduced to balance the jump in bulk shear stress at the singular interface (Edwards et al. [START_REF] Edwards | Interfacial Transport Processes and Rheology[END_REF] ). To model the relationship between surface stress and surface strain, the Boussinesq-Scriven constitutive law is used (Aris 39 ) to write the two following components of the JMB:

(η S + κ S ) ( d 2 v r S dr 2 + 1 r dv r S dr - v r S r 2 ) = η ∂v r ∂z z=h , (14)
η S ( d 2 v θS dr 2 + 1 r dv θS dr - v θS r 2 ) = η ∂v θ ∂z z=h . (15)
The left-hand term of Eq. ( 14) represents a combination of surface viscous shear and surface viscous dilatation along the liquid surface, while only the surface viscous shear is present in Eq. (15), where η S is the surface (in-excess) shear viscosity, and κ S is the surface (in-excess) dilatation viscosity. They stand as the second coupling term of the aforementioned two-way coupling. The right-hand terms of both equations are the liquid shears vertically imposed from the supporting subphase flow, with η, the Newtonian bulk shear viscosity. Equations ( 14) and ( 15) are discretized along the interface by calculating their respective weak forms. Finally, to solve for the JMB, the following Dirichlet end-point boundary conditions for surface flow are required:

v r S (r = r i ) = v r S (r = r o ) = v θS (r = r i ) = v θS (r = r o ) = 0. ( 16
)
b. Dimensionless quantities. Normalizing Eqs. ( 14) and ( 15) leads to the definition of two rheological scaling parameters:

Bo η S = η S ηh , Bo κ S = κ S ηh , (17)
where Bo η S and Bo κ S are the surface shear and surface dilatational Boussinesq numbers, respectively. The Bo ηS number describes the balance between bulk and surface viscous shears, while Bo κ S expresses the ratio between the dilatational stress along the interface and the bulk viscous shear. It should be noted that, compared to the supporting article, [START_REF] Delacroix | ªElectrical activity of the hartmann layers relative to surface viscous shearing in an annular magnetohydrodynamic flow[END_REF] Bo κ S is a new feature emerging due to the inertial effects. The dimensionless quantities of interest for surface dynamics consist only of surface velocity:

⃗ v ⋆ S = ⃗ v S /V ; v ⋆ r S
and v ⋆ θS allow for analysis of the main and secondary surface flows, respectively.

Strong coupling and numerical implementation

Overall coupling between subphase and surface flows is made evident in Fig. 5. The electromagnetism and hydrodynamics modules interact by means of the coupling terms, i.e., the electromotive current and the Lorentz force. The surface rheology module provides the velocity boundary condition for solving bulk flow through surface velocity, while the hydrodynamics module provides in turn for the vertical bulk shears at the interface, to solve the surface flow.

With respect to numerical methods, due to the basic layout of the computational domain, a fully coupled approach is implemented, based on the full Jacobian matrix as one entity. This approach is based on the Newton-Raphson method, which linearizes the non-linear problem based on the current solution, at each iteration. A linear stationary direct solver is implemented to solve the linearized problem, i.e., the MUltifrontal Massively Parallel sparse direct Solver (MUMPS), based on Lower-Upper (LU) factorization (see MUMPS support [START_REF]ªMultifrontal massively parallel solver (mumps 4.10.0) user's guide,º[END_REF] for further details).

Finally, the implemented mesh is displayed in Fig. 6(a). It consists of 27 524 elements, mainly triangular, with mesh refinement on the Galinstan domain. As shown in Fig. 6(b), a specific rectangular boundary layer mesh is set up at the boundaries of the fluid domain. Typically, the relative thickness of the first layer is set so as to be much lower than the reciprocal of the Hartmann and Reynolds numbers, which both monitor the thickness of the physical boundary layers.

IV. RESULTS AND INTERPRETATION

A. The annular hydrodynamic viscometer (Ha = 0)

In this section, the interaction between the purely hydrodynamic subphase flow and the surface rheology is analyzed. Let us first validate the 2-D 1/2 numerical modeling with a benchmarking asymptotic study.

The deep-channel viscometer benchmark: Inertia and surface viscous shearing

The first selected benchmark is the numerical study of the classical deep-channel annular viscometer (with no outer magnetic field, i.e., Ha = 0), performed by Hirsa et al., [START_REF] Hirsa | ªDetermination of surface shear viscosity via deep-channel flow with inertia[END_REF] from a numerical code based on a second order centered differences scheme. In this case, only the hydrodynamic issue is brought into focus, taking inertia into account, but considering that the interfacial condition along the radial direction (Eq. ( 14) of the present study) reduces to no-slip, which is similar to the condition Bo κ S → ∞. As a result, the role of the surface viscous shearing alone is enhanced, through the parameter Bo η S .

In this part, relevant quantities must be temporarily redefined to fit with Hirsa et al. notations and typical values and are subscripted HL (corresponding to Hirsa layout). Thus, the height and the inner and outer radii of the channel are noted h HL , r iHL , and r oHL , respectively. Note that h HL ∼ r oHL -r iHL , hence the ªdeep-channelº denomination. The surface shear Boussinesq number used in this part is Bo η SHL = η S /ηr oHL . The definitions of the other relevant quantities are unchanged (particularly the Reynolds number Re).

Some results obtained for the same geometrical layout as Hirsa et al. are then shown in Fig. 7. The quantities displayed are

• the contours of the azimuthal vorticity ω ⋆ = ∂v ⋆ r /∂z ⋆ -∂v ⋆ z /∂r ⋆ , which rules the meridian flow, in Fig. 7 The agreement between both studies is quite satisfying. Both models predict a secondary overturning flow due to inertial effects, intensifying with Re. At low Reynolds numbers, the angular momentum diffuses from the rotating floor into the liquid metal (see Fig. 7(b), top). For higher Re, an Ekman layer appears near the rotating floor, concentrating the vortex lines. Due to the outer wall, this layer is turned vertically into a jet-like flow, and then turned inwards at the interface, advecting fluid with large angular momentum radially inwards. This results in a spin-up of the fluid interface increasingly closer to the inner wall [START_REF] Hirsa | ªDetermination of surface shear viscosity via deep-channel flow with inertia[END_REF] (see Fig. 7(b), bottom, and Fig. 7(c)).

With respect to the impact of surface viscous shearing, it can be seen in Fig. 7(b) that if Bo η SHL is low (top), the vortex lines are normal to the liquid surface, and if it is high (bottom), then the lines are tangential to the liquid/gas interface, evolving consistently from a free surface (Neumann) to no-slip (Dirichlet) boundary conditions. Consequently, v ⋆ θS decreases as Bo η SHL increases, as seen when comparing Figs. 7(c) and 7(d). When Bo η SHL is sufficiently high, the interface is then only governed by surface viscous shearing, and the bulk shear taken at the vicinity of the interface does not affect it. This leads to a vanishing Poiseuille-like velocity profile along the interface, as shown in Fig. 7(d).

The shallow annular channel: Inertia, surface viscous shearing, and dilatation

Now that the relevance of hydrodynamic numerical modeling is proved, let us return to the geometrical layout of the shallow annular channel, as shown in Fig. 2. In this section, the interaction between inertia and surface rheology is particularly emphasized for bulk flow, while the previous benchmark mainly focuses on surface phenomena. The results are displayed in Fig. 8.

a. Bo η S impact. The way in which surface viscous shearing affects the hydrodynamic flow through Bo η S is investigated in this section, by comparing Figs. 8(a) and 8(b) (Bo κ S is negligibly FIG. 7. Numerical benchmark for the purely hydrodynamic case. [START_REF] Hirsa | ªDetermination of surface shear viscosity via deep-channel flow with inertia[END_REF] small: Bo κ S = 10 -4 ). For the azimuthal flow v ⋆ θ , velocity profile evolution with respect to Bo η S is similar to Sec. IV A 1, with slight changes due to the difference in geometrical layouts between the deep and the shallow channel viscometer.

With respect to the overturning flow (v ⋆ r , v ⋆ z ), an increase in Bo η S leads to flow homogenization, as shown in Figs. 8(a) and 8(b). When Bo η S = 10 -4 , i.e., when liquid surface acts similarly to a free surface, the main vortex governing the overturning flow is mainly located in the outer part of the channel. When Bo η S = 10 4 , this vortex expands radially inwards throughout the whole cross section of the channel. This difference can be accounted for by the nature of the interface boundary condition. When the latter is similar to a free surface, the momentum injected from the rotating floor is dissipated in the bulk and at the sliding interface, with a significant surface radial velocity (see Fig. 8(a)). When the liquid surface is rigid (vanishing surface velocity, see Fig. 8(b)), it no longer participates in viscous damping of the injected momentum. Bulk viscous damping is therefore enhanced, leading to expansion of the main vortex inside the subphase flow.

Note that when Re varies at fixed Bo η S , overturning flow magnitude increases with Re, and the main vortex governing the secondary flow is enlarged. Momentum is then increasingly confined within the Ekman layers, similarly to what is observed in Sec. IV A 1 for the main flow. However, the qualitative impact of Bo η S on the meridian flow remains identical irrespective of Re. Bo η S = 10 -4 ). For the azimuthal flow, an increase in Bo κ S does not lead to a significant change in the velocity profile v ⋆ θ . It seems only that the main flow is less efficiently advected by the secondary flow. This observation is confirmed by turning to analysis of the secondary flow (v ⋆ r , v ⋆ z ). Contrary to what is observed for Bo η S , in this case, the main vortex governing the secondary flow decreases when Bo κ S increases in Fig. 8(c). This phenomenon is quite unexpected, because for the secondary flow, the boundary condition at the liquid surface along the radial direction is the same both in the case of an infinite value for Bo η S and for Bo κ S , i.e., a vanishing v ⋆ r S (where the interface is similar to a non-sliding wall, see Eq. ( 14)). However, as Bo η S also governs the azimuthal boundary condition for the main flow, when the surface viscous shear is very high, the main azimuthal flow must also match with a vanishing v ⋆ θS (see Eq. ( 15)). Consequently, there is a complete deflection of momentum at the interface, favoring expansion of the main overturning vortex as mentioned above (Fig. 8(b)). For surface dilatation, even when Bo κ S greatly increases, it has no direct consequence on the boundary condition v ⋆ θS for the main flow. Therefore, for Bo κ S = 10 4 , the overturning flow must match with a vanishing value of v ⋆ r S , whereas the main flow remains unchanged at the vicinity of the interface (with a fixed Bo η S = 10 -4 ). This favors damping of the main vortex governing the secondary flow near the interface observed in Fig. 8(c).

Note that, as for Bo η S , changing the value of Re does not affect the qualitative impact of Bo κ S on the centrifugal flow.

B. The annular MHD viscometer (Ha 0)

The electroconductive flow is now subjected to an outer magnetic field B 0 ⃗ e z . This time, the interaction between the MHD subphase flow and the surface rheology is investigated. We shall begin with a benchmarking asymptotic study to validate 2-D 1/2 MHD numerical modeling.

MHD, no inertia, and surface viscous shearing

The second benchmark is the analytical study of the annular MHD viscometer performed in a companion paper, [START_REF] Delacroix | ªElectrical activity of the hartmann layers relative to surface viscous shearing in an annular magnetohydrodynamic flow[END_REF] within the limit of ignoring inertial effects when Ha ≫ Re and Ha ≫ 1 (see Eq. ( 19)). If the meridian flow is not taken into account, the governing equations are considerably simplified, thus allowing implementation of a (⃗ v, ⃗ B) formulation for the MHD bulk flow. Then, without the overturning flow, the interface dynamics is only governed by Eq. ( 15), meaning that only surface viscous shearing is brought into play through the Boussinesq number Bo η S . A two-way coupling is then set up, allowing solving for the surface dynamics, followed by determination of the whole MHD flow by means of an asymptotic matched expansion based on the small parameter 1/Ha ≪ 1. Some results are shown in Fig. 9. Note that cutting lines (A) and (B) are defined in Fig. 3 and allow for analysis of the MHD core flow and the Shercliff or the Hartmann layers, respectively. In addition to this analytical benchmark, a second 2-D numerical study within the asymptotic limit N ≫ 1 is added for the sake of multiple benchmarking. This 2-D numerical study is based on the same (⃗ v, ⃗ B) formulation as the one developed in the analytical study.

The general agreement between all studies is once again quite satisfying. They predict radically different topologies for the MHD flows, depending on the relative values of Ha and Bo η S . If Ha ≫ Bo η S , the v ⋆ θ contours demonstrate a 2-D tendency with exclusively radial velocity gradients (except for near the side-walls), as seen in Fig. 9(a). This rigid-body motion, expressing the electromagnetic blocking of the flow, first observed by Lehnert, [START_REF] Lehnert | ªAn instability of laminar flow of mercury caused by an external magnetic field[END_REF] is caused by the well-known two-dimensionality FIG. 9. 2-D analytical and numerical benchmarks for the annular MHD viscometer, with Ha ≫ Re (no inertia). [START_REF] Delacroix | ªElectrical activity of the hartmann layers relative to surface viscous shearing in an annular magnetohydrodynamic flow[END_REF] The electric current densities are normalized with respect to the maximum electric current J max = 6.8 × 10 2 A m -2 reached in all cases, i.e., for Ha = 50. tendency of magnetic induction. Consequently, the interface is perfectly aligned with the bulk, and the bulk viscous shear at the interface is no longer significant. Therefore, the electric current density is essentially confined to the Shercliff layers, with two electric loops closing up near the side-walls (see Figs. 9(c) and 9(d), case (Ha, Bo η S ) = (50, 0.01)). Now, when Ha ≪ Bo η S , the three modelings lead to a quite homogeneous ªmotionlessº configuration, where the momentum is mainly concentrated near the right-hand corner at the bottom (Fig. 9(b)). This singular phenomenon is partially explained by the fact that, in this case, the surface dynamics is governed by surface viscous shear and behaves as a non-sliding membrane. Thus, v ⋆ θ must match with the vanishing component v ⋆ θS along the surface. However, this cannot solely account for the motionless layout across the whole cross section. Other reasons are found by focusing on the electric current densities. Due to strong velocity gradients near the liquid/gas surface, and to current continuity, electric current densities are now found to flow within the top and bottom Hartmann layers, which are therefore electrically activeÐsee Fig. 9(d). The presence of a strong radial component of electric current density, combined with the imposed magnetic field B 0 ⃗ e z , leads to the emergence of a Lorentz forcej r B 0 along the azimuthal direction. As deduced from Fig. 9(d), this Lorentz force is negative at the bottom and positive at the top of the channel cross section. Consequently, this leads to an electromagnetic damping of the momentum injected from the rotating floor at the bottom, while it enhances momentum in the top part of the channel. Both contributions lead to a homogenization of the flow, which accounts for the overall flow patterns.

To conclude, we confirm that for the asymptotic case Ha ≫ Re, the surface shear Boussinesq number drives electrical activation of the Hartmann layers. Thus, surface rheology can really monitor the magnitude of both velocity and electric currents in the MHD core and the boundary layers, highlighting the competitive effects between MHD tendency towards two-dimensionality, and surface viscous shearing.

Scaling laws for inertial MHD and surface rheology interaction

2-D 1/2 numerical modeling, the reliability of which is demonstrated from the previous asymptotic case, can now be confidently extrapolated to conditions including MHD of an annular swirling flow, coupled with surface rheology. For this purpose, a scaling law is required to determine the relevant values of Ha and Re, under which the interactions between the (main and secondary) subphase MHD and surface flows are the most insightful. We choose the following criterion, based on the ratio of surface velocities for each (Re, Ha) pair:

v ⋆ r S v ⋆ θS max ≥ 0.01. ( 18
)
The choice of this particular limit is explained in Fig. 10, which represents the magnitude of this ratio compared with the values of (Re, Ha) pairs. In this figure, three main areas can be distinguished: first, the hatched one represents the area where inertia can be ignored. In this part, which includes both 2-D analytical and numerical studies, the impact of surface viscous shear only can be investigated. The validity limit of this area stems from the criterion taken from Tabeling and Chabrerie, [START_REF] Tabeling | ªMagnetohydrodynamic secondary flows at high hartmann numbers[END_REF] stating that inertial effects can be ignored as soon as

h 3 Re 2 Rr 2 o Ha 2 √ Ha < 1, (19)
with R = r o -r i (this criterion is better adapted to the annular channel layout than the generic criterion N = Ha 2 /Re ≫ 1). Note that even if the matched asymptotic expansion performed in the analytical study is theoretically only valid for Ha → ∞, the results are in perfect agreement with the 2-D numerical study as soon as Ha ≥ 10.

The second area is referred to as ªweakly swirling flow.º In this transitional zone, inertia can no longer be ignored. However, the overturning flow is not strong enough to stretch significantly the interface, and the magnitude of the ratio v ⋆ r S /v ⋆ θS max continues to be very small. Therefore, surface dilatation impact on bulk MHD is likely to be negligible, and this ªweakly swirling flowº area is left apart from the present analysis. The color-mapped area represents the working conditions for which v ⋆ r S /v ⋆ θS max ≥ 0.01. Note that the purely hydrodynamic study by Hirsa et al. [START_REF] Hirsa | ªDetermination of surface shear viscosity via deep-channel flow with inertia[END_REF] is also indicated as the abscissa axis Ha = 0. In this part, surface radial velocity tends to become the same order of magnitude as the surface azimuthal velocity. As a result, surface viscous dilatation can become significant and may strongly interact with bulk flow through the Bo κ S number, in addition to surface viscous shear.

The goal of the following study is to highlight the secondary MHD flow due to inertia and emphasize the coupling between bulk MHD and surface viscous shear and dilatation. To maximise the most salient interactions, we choose the working conditions matching with the third area. Fig. 10 gives us the relevant scaling laws, i.e., for Re ≤ 10 4 , in the following results, the value of the Hartmann number does not exceed the threshold value Ha = 10.

MHD, inertia, and surface viscous shear

In this section, the influence of surface viscous shearing only (through Bo η S ) on main and secondary MHD bulk flows is investigated, in order to complete the analytical analysis. [START_REF] Delacroix | ªElectrical activity of the hartmann layers relative to surface viscous shearing in an annular magnetohydrodynamic flow[END_REF] For this purpose, the surface dilatation Boussinesq number is set to be negligibly small: Bo κ S = 10 -4 . a. Main MHD bulk flow. In this part, there is no need to produce any graphs either for v ⋆ θ or for ( j ⋆ r , j ⋆ z ) (note that an overall view for ⃗ j is shown in Fig. 12). Depending on the relative values of Ha and Re, the qualitative effects of surface viscous shearing on the velocity and electric current density profiles are essentially the same as in the two previously detailed benchmarking asymptotic cases. For instance, when Re = 10 3 and Ha = 5, if Bo η S ≪ Ha, then the fluid flow tends towards a rigid body motion (though not really pronounced as Ha is not high enough) where the electric current densities are confined within the (thick) Shercliff layers. If Bo η S ≫ Ha, an electromagnetic damping of the flow is also observed (again not so marked), and the Hartmann layers become electrically active, as described in Sec. IV B 1. When the Reynolds number is increased (for instance, Re = 10 4 ), the velocity profile evolution with respect to Bo η S is similar to Sec. IV A. The influence of Ha is less significant, as the flow is almost purely hydrodynamic.

b. Secondary MHD bulk flow. Let us now analyse the overturning MHD flow, beginning with the meridian velocities (v ⋆ r , v ⋆ z ), for which no figures are displayed. The outer magnetic field does not yield a significant qualitative difference with respect to the purely hydrodynamic case (Sec. IV A 2), regarding the Bo η S impact. Indeed, for given (Re, Bo η S ) values, an increase in Ha affects both the meridian velocities magnitude and topology. The main vortex governing the meridian flow is increasingly diminished, due to MHD tendency towards two-dimensionality. However, for given (Re, Ha) values, increasing Bo η S still leads to the stretching of the overturning flow, as in Figs. 8(a However, this is very different when the azimuthal electric current j ⋆ θ is considered. The typical order of magnitude for electric current densities is in fact J = σB 0 V (V is the order of magnitude of the main flow velocity v ⋆ θ , directly imposed through the rotating floor). Consequently, an increase in Ha or Re leads to higher orders of magnitude for the meridian electric current densities ( j ⋆ r , j ⋆ z ) (linked to the main MHD flow). However, as j θ = -σv r B 0 (see Eq. ( 8)), there is competition between the electromagnetic and inertial effects concerning the secondary flow: when Ha (and therefore B 0 ) is increased, the radial velocity v r is damped due to the MHD tendency towards two-dimensionality. Figs. 11(a) and 11(b) show the evolution of j ⋆ θ with respect to different (Re, Ha, Bo η S ) triplets, along cutting lines (A) and (B), defined in Fig. 3.

For given (Ha, Bo η S ) values, an increase in Re causes a significant increase of j ⋆ θ along (A) and (B), directly linked to the expansion and strengthening of the centrifugal flow highlighted in Sec. IV A 2 (see for instance curves (blue solid line) and (blue dashed-dotted line), or (green square) and (green plus)).

For given (Re, Bo η S ) values, variation of Ha results in many evolutions for j ⋆ θ , from which no general principle can be deduced. Thus, it appears clear from Fig. 11(b) that along (B), an increase in Ha leads to a drop of j ⋆ θ magnitude for low Bo η S values (see, for instance, blue curves (blue solid line) and (blue dashed line)), whereas for high Bo η S values (green symbols, for instance, (green diamond) and (green square)), j ⋆ θ increases with Ha. This seems to be globally the contrary in Fig. 11(a), along cutting line (A). This complex interaction expresses the aforementioned competition between electromagnetic and inertial effects, which have opposite influence on the evolution of j ⋆ θ . Finally, for given (Re, Ha) values, an increase in Bo η S helps to expand the main vortex governing the secondary flow, similarly to what is shown in Figs. 8(a) and 8(b). Consequently, along (A), the curves are flattened with high Bo η S values, because j ⋆ θ is proportional to -v ⋆ r (Fig. 11(a)). Concerning (B), for Re = 4 × 10 3 , Ha = 10 and Bo η S = 10 -4 (Fig. 11(b)), the main vortex closes up without including the vertical cutting line, leading to a vanishing j ⋆ θ due to vanishing v ⋆ r (curve (blue dashed line)). At higher Bo η S = 10 4 , expansion of the main vortex leads to the emergence of j ⋆ θ along (B) (curve (green square)). The overturning shape of the centrifugal flow is consistently found, for the latter case as well as for other (Re, Ha) values. This centrifugal profile is increasingly symmetric with respect to the z ⋆ = 0.5 axis with increasing Bo η S , with the asymptotic value Bo η S = 10 4 imposing a vanishing j ⋆ θ at the liquid surface (see for instance curves (blue dashed-dotted line) and (green plus)).

We shall now make one final comment about this part, concerning the role played by the Shercliff or Hartmann boundary layers. In the benchmarking case of Sec. IV B 1 linked to the main MHD flow, the influence of these layers is significantly affected by the surface viscous shear. In turn, they dramatically impact bulk MHD. For the secondary MHD flow analyzed in this section, this crucial role of these layers has not been found. This can be accounted for by the nature of the secondary FIG. 12. ⃗ j ⋆ for Ha = 5. ⃗ j ⋆ is log-scaled by the magnitude exp

( ( ln ( ∥ ⃗ j ⋆ ∥/∥ ⃗ j ⋆ ∥ max ) ) /(1 + p) )
; p = 1.5 for (a)±(c), and 2 for (d). Black (greyscale: darkest) arrows are essentially meridian, while beige (greyscale: lightest) and blue (greyscale: intermediate) arrows correspond to significantly (i.e., j ⋆ θ /∥ ⃗ j ⋆ ∥ ≥ 0.01) positive and negative azimuthal components, respectively.

MHD flow: j ⋆ θ is indeed purely electromotive, being only linked to radial velocity. Consequently, it naturally closes within the MHD core all around the channel with no need for Hartmann or Shercliff boundary layers. The fact that there is no requirement for any electric fields leads to a interplay between j ⋆ and the meridian flow (v ⋆ r , v ⋆ z ). c. Overall view of electric current densities. To conclude this section about the impact of surface viscous shearing on bulk MHD, let us display some figures allowing an overview of the overall MHD bulk flow. Note that only graphs for electric current densities ⃗ j ⋆ are shown in Fig. 12, for a given Ha = 5. Displaying figures for the velocity field ⃗ v ⋆ is not worthwhile, as the ratio v ⋆ r /v ⋆ θ is relatively low throughout the cross section of the annular channel (except in the outer side layer, where it may be locally significant), ⃗ v remains essentially azimuthal. Besides, the most interesting features are similar to those highlighted in Fig. 8 for the purely hydrodynamic case.

Concerning electric current densities, very interesting features appear depending on the values of (Re, Bo η S ). At low Re, the main electric loops are located near the side walls in the Shercliff layers for low Bo η S values (see Fig. 12(a)). An increase in Bo η S magnifies the role of the Hartmann layers (see Fig. 12(b)). The impact of the overturning flow can be already seen, with quite significant azimuthal perturbation throughout the cross section (increasing as Bo η S increases). This effect is even more dramatic for very high Re numbers. When Bo η S = 10 -4 , the azimuthal component j ⋆ θ is very strong, although the helical pattern of electric current density is not well-structured, because the main vortex governing the secondary flow is still not expanded (Fig. 12(c)). For Bo η S = 10 4 , this expansion takes place throughout the entire cross section, and the current density closes up into a kind of well-developed helical pattern in the core flow (Fig. 12(c)). Since the secondary flow is radially inwards along the surface, the electric current loops are counter-clockwise when the channel is seen from above. Along the rotating floor, the electric current loops are found to be clockwise. This is the reason why the electric current is so twisted at the mid-height of the cavity.

MHD, inertia, and surface viscous dilatation

Following analysis of the impact of surface viscous shear on the overall MHD bulk flow, we shall now investigate the interactions between surface viscous dilatation and the subphase flow. As we focus only on the parameter Bo κ S , the surface shear Boussinesq number is set to be negligibly small in this section: Bo η S = 10 -4 . a. Main MHD bulk flow. When an outer magnetic field is applied, some interesting new features appear depending on the Ha value. Concerning azimuthal velocity v ⋆ θ , the results are shown in Fig. 13, for Re = 10 4 , Ha = 10, Bo κ S = 10 4 (contours). These results can be compared to the purely hydrodynamic case (Ha = 0) of Fig. 8(c). Contrary to the latter, the cross section can be divided into two sub areas in this case. For the inner part of the cross section, where the velocity magnitude is lower, the inertial effects are weak compared to the electromagnetic effects. Consequently, the electromagnetic blocking observed previously in the benchmarking case Sec. IV B 1 predominates. For the outer part, inertia is preponderant, and the advection v ⋆ θ by the overturning flow is strong, as Sec. IV A 2. In between, a transitional area is found, governed by the balance between inertial and electromagnetic effects.

With respect to meridian electric current densities, there is no need to display some figures (note that an overall view for ⃗ j is shown in Fig. 15). For given Re and Ha values, an increase in Bo κ S does not significantly affect the electric circuit: the current densities remain confined within the Shercliff layers. Particularly, it is important to highlight that, contrary to Bo η S , Bo κ S does not have a noticeable influence over the boundary layers. Thus, an increase in Bo κ S does not lead to electrical activation of the Hartmann layers, which considerably restricts the impact on the main MHD bulk flow.

b. Secondary MHD bulk flow. The interaction between surface dilatation and the meridian flow is now investigated, by focusing on the quantities (v ⋆ r , v ⋆ z ) and j ⋆ θ . We shall begin with the meridian velocities, for Re = 10 4 , Ha = 10, Bo κ S = 10 4 (arrows), in Fig. 13. The first obvious impact of the outer magnetic field, when compared to the purely hydrodynamic case of Fig. 8(c), is the decay of the main vortex governing the secondary flow. This vortex increasingly confined to the outer part of the cross section either with increasing Ha, due to MHD tendency towards two-dimensionality, or with increasing Bo κ S , as observed in Fig. 8(c) for the purely hydrodynamic case.

Moreover, depending on the (Re, Ha, Bo κ S ) values, the emergence of new flow patterns is observed. In the transitional area mentioned above, a counter-rotating vortex appears. Such a complex overturning flow pattern greatly depends on the values of the three parameters (Re, Ha, Bo κ S ): it is not observed with Re = 10 3 , or Ha = 5, or Bo κ S = 10 -4 . Therefore, contrary to the Bo η S investigation, the meridian velocities (v ⋆ r , v ⋆ z ) exhibit different flow patterns depending on the value of the Ha number. The emergence of this new vortex seems to stem from a purely hydrodynamic phenomenon. The radial Lorentz force in this area is found to be always oriented in the opposite direction to the velocity field, due to Lenz's law. The purely electro-motive azimuthal component of the electric current density j ⋆ θ only results in electromagnetic braking of this counter-rotating flow, and as such, cannot generate this new vortex. Rather, it is the coincidence of the values of the three parameters (Re, Ha, Bo κ S ) that leads to this original flow pattern: sufficiently high Ha and Bo κ S values enhance electromagnetic blocking of the flow in the inner part of the channel (by damping the main vortex governing the secondary flow, thus promoting the rigid-body motion), whereas a high Re number allows inertial effects to predominate in the outer part of the cross section. The geometrical configuration of the flow is therefore greatly modified, with the centrifugal flow concentrated against the outer wall, leading to the hydrodynamic development of this new counter-rotating vortex. Note that this new flow pattern is not found in the corresponding part concerning the Bo η S effect: for the same Re = 10 4 , Ha = 10 values, an increase in Bo η S leads to damping of the main azimuthal flow, through electrical activation of the Hartmann layers, and to stretching of the secondary flow. Flow topology is significantly different from that observed here for the Bo κ S effect, and as such, does not lead to the emergence of the counter-rotating flow.

We shall now examine the evolution of j ⋆ θ , along the cutting lines (A) and (B), in Figs. 14(a) and 14(b). As explained in the part concerning the influence of Bo η S , it is worthwhile investigating the dependence of j ⋆ θ not only with respect to Bo κ S but also to Re and Ha. For given (Ha, Bo κ S ), the impact of varying Re is the same as for Bo η S , i.e., j ⋆ θ increases with Re due to the expansion of the overturning flow (e.g., curves (blue dashed line) and (blue dotted line), or (green diamond) and (green circle)).

A change in the Ha value leads to complex evolutions along (A) and (B), similar to the Bo η S effect, relating to the competition between inertial and electromagnetic effects. However, there is a major difference between surface dilatation and shear. For Bo κ S , along a given cutting line, the impact of increasing Ha remains qualitatively the same at low or at high Bo κ S values, whereas it depends on the value of Bo η S for the previous case. The effect is nonetheless the opposite between the two cutting lines. For instance, along (A), the magnitude of j ⋆ θ globally increases with Ha both at low (see curves (blue dashed line) and (blue dotted line)) and at high (see curves (green plus) and (green circle)) Bo κ S values. Along (B), the magnitude of j ⋆ θ is lower when there is an increase in Ha, whatever the value of Bo κ S (same curves).

Finally, for given (Re, Ha) values, along (A) (Fig. 14(a)), an increase in Bo κ S does not lead to flattening of the profiles previously observed with Bo η S (Fig. 11(a)). Rather, at high Bo κ S values (green symbols), j ⋆ θ oscillates around the j ⋆ θ = 0 axis, where the amplitude and the number of peaks are maximum for (Re, Ha, Bo κ S ) = (10 4 , 10, 10 4 ) (curve (green plus)). These oscillations are linked to the emergence of the counter-rotating vortex for high Bo κ S values ( j ⋆ θ is proportional to -v ⋆ r ). Note that near the inner wall, the curves for low (blue lines) and high (green symbols) Bo κ S values are superimposed. This is because an increase in Bo κ S does not lead to a significant change for the overturning flow in this part of the channel cross section. Along (B), and contrary to what is noticed for the Bo η S effect (Fig. 11(b)), the overturning profile is not recovered with an increase in Bo κ S (Fig. 14(b)). Rather, the curves are increasingly flattened, due to the progressive damping of the main vortex governing the secondary flow. Consequently, the cutting line (B) is no longer included in the area where v ⋆ r is still significant. As for the Bo η S counterpart, the interaction between the Shercliff or Hartmann layers and the surface dilatation (through Bo κ S ) does not provide any insights into the swirling flow topology. c. Overall view of the electric current densities. To complete this section on the impact of surface viscous dilatation on bulk MHD, we shall now display some results allowing an overview of overall MHD bulk flow in Fig. 15. Note that only graphs for the electric current densities ⃗ j ⋆ are shown, for the same reason as for the Bo η S effect: the most interesting features for ⃗ v were previously discussed in relation to Fig. 13. Contrary to the Bo η S impact, it was previously shown that the flow pattern qualitatively depends on the value of the Ha number as well as on Re and Bo κ S (with the emergence of the counter-rotating vortex). Therefore, the Bo κ S impact is analyzed for Ha = 5 (left part) and for Ha = 10 (right part).

In this section, we shall begin with the Bo κ S influence over the whole MHD bulk flow. For this purpose, the left parts of To highlight some interesting new features linked to a high Bo κ S value, the Ha number is increased to 10 (Fig. 15, right part). When Re = 10 3 , the electric current densities are essentially confined inside the Shercliff layers (Fig. 15(b)), similarly to what is observed in the benchmarking case of Sec. IV B 1 for a low Bo η S = 10 -4 . For a higher Re = 10 4 , the new patterns previously observed for the secondary flow are recovered. The helical pattern for ⃗ j, previously observed in the core for the Bo η S effect (Fig. 12(d)), is recovered in Fig. 15(d). However, the electric current is mostly located near the outer wall, with the inner half of the flow being current-free, due to the confinement of the main vortex governing the secondary flow for a high Bo κ S . Besides, ⃗ j is significantly perturbed by the emergence of the counter-rotating vortex. This vortex induces radially outwards (at the top) and radially inwards (near the bottom) flows, whereas along the rotating floor, the flow must remain radially outwards (see inset in Fig. 13). Consequently, the electric current densities are twisted at two locations, essentially because their azimuthal orientation is governed by the radial component of the meridian flow: a first twist is observed in the middle of the secondary vortex, while a second twist is observed near the rotating floor. Such an original pattern is a typical example of how surface dilatation may affect the MHD bulk flow.

V. CONCLUSIONS AND OUTLOOK

To complete the mathematical approach presented in a companion paper, [START_REF] Delacroix | ªElectrical activity of the hartmann layers relative to surface viscous shearing in an annular magnetohydrodynamic flow[END_REF] which focused on the interaction between surface viscous shear and an annular MHD flow with no inertia, the overall coupling mechanism between the full viscous rheology of the liquid surface (including surface dilatation) and a supporting annular MHD flow (taking inertial effects into account) is successfully investigated in this paper.

This paper shows how a change in the mechanical properties of a fluid interface can greatly influence a MHD core flow. Moreover, it has proved that shear and dilatational viscosities of the surface do not generate the same changes. On the one hand, viscous shearing of the interface actively modifies the main annular flow by means of the Hartmann layers, which become electrically active. A damping is clearly demonstrated, as is also a 2-D MHD tendency. On the other hand, the dilatational viscosity of the interface is only responsible for damping the meridian flow driven by centrifugation, with a new pattern if inertia and electromagnetic blocking are both significant.

In addition to this analysis, the relevant working conditions highlighting the competitive effects of inertia, electromagnetism, and surface rheology have been identified. The corresponding values of the dimensionless numbers allow us to establish some scaling laws for an upcoming experimental facility.

One concluding remark relates to the activity of the electrical boundary layers. Given the relatively low Ha values required for development of the swirling flow, the influence of the Hartmann layers is in any case reduced. However, regarding the possible influence of surface viscous dilatation on the main and secondary flows, the Shercliff layers and, in particular, the Hartmann layers do not especially become electrically active throughout the many configurations tested in this paper. Consequently, their impact on the core flow is not as dramatic as when a strong surface viscous shear is applied, the latter case promoting the damping of the main azimuthal flow via the electrical activation of the Hartmann layers. [START_REF] Delacroix | ªElectrical activity of the hartmann layers relative to surface viscous shearing in an annular magnetohydrodynamic flow[END_REF]

FIG. 1 .

 1 FIG.1. Typical layout of a MHD liquid/gas flow.

FIG. 3 .

 3 FIG. 3. Geometry and boundary conditions of the channel cross section used for numerical computation. Note the presence of cutting lines ( A) : z = h 0 /2 and (B) : r = (r i + r o )/2 used hereafter for interpretation of the results.

  FIG.4. Lubrication gaps between the rotating floor and the motionless side-walls with respect to the experimental layout (see, e.g., Mannheimer and Schechter[START_REF] Mannheimer | ªAn improved apparatus and analysis for surface rheological measurements[END_REF] ).

FIG. 5 .

 5 FIG. 5. Overall coupling.

FIG. 6 .

 6 FIG. 6. Mesh used for the numerical computation (27 524 elements): (a) global view, (b) zoom on the boundary layer mesh.

  (a), for different Bo η SHL and Re values; • the contours of the axial angular momentum α ⋆ = r ⋆ v ⋆ θ , in Fig. 7(b), for different Bo η SHL and Re values; • the surface azimuthal velocity v ⋆ θS in Figs. 7(c) and 7(d), for Bo η SHL = 0.01 or 10, respectively, and various Re values.

b.

  FIG. 8. Bo η S and Bo κ S impact on the velocity field for the classical annular viscometer layout (Ha = 0). (v ⋆ r , v ⋆ z ) is log-scaled by the magnitude exp ln ∥(v ⋆ r , v ⋆ z )∥/∥(v ⋆ r , v ⋆ z )∥ max /(1 + p) ; p = 0.5 for (a)±(c).

FIG. 10 .

 10 FIG.10. Scaling law for determining the most insightful working conditions, in view of overall coupling between MHD with inertia and surface rheology with both surface shear and dilatation.

  ) and 8(b). Consequently, the most salient phenomena qualitatively remain the same irrespective of the Re and Ha values.

FIG. 11 .

 11 FIG. 11. Azimuthal electric current density j ⋆ θ , along cutting lines ( A) for (a) and B for (b), for various Re, Ha, and Bo η S values. j ⋆ θ is normalized with respect to the maximum electric current J max = 3.96 × 10 3 A m -2 reached in all cases, i.e., for the highest (Re, Ha) = (10 4 , 10) values.

FIG. 13 .

 13 FIG. 13. Bo κ S impact on the velocity field for the annular MHD viscometer: main (contours: v ⋆ θ ) and secondary (arrows: (v ⋆ r , v ⋆ z )) flows, Re = 10 4 , Ha = 10, Bo κ S = 10 4 . (v ⋆ r , v ⋆ z ) is log-scaled by the magnitude exp ((ln(∥(v ⋆ r , v ⋆ z )∥/∥(v ⋆ r , v ⋆ z )∥ max ))/ (1 + p)), where p = 1. Inset: streamlines of the counter-rotating vortex.

FIG. 14 .

 14 FIG. 14. Azimuthal electric current density j ⋆ θ , along cutting lines ( A) for (a) and B for (b), for various Re, Ha, and Bo κ S values. j ⋆ θ is normalized with respect to the maximum electric current J max = 3.96 × 10 3 A m -2 reached in all cases, i.e., for the highest (Re, Ha) = (10 4 , 10) values.

FIG. 15 .

 15 FIG. 15. ⃗ j ⋆ for Re = 10 4 . ⃗ j ⋆ is log-scaled by the magnitude exp ( ( ln ( ∥ ⃗ j ⋆ ∥/∥ ⃗ j ⋆ ∥ max ) ) /(1 + p) ) ; p = 2 for (a) and (b) and 1.5 for (c) and (d). Black (greyscale: darkest) arrows are essentially meridian, while beige (greyscale: lightest) and blue (greyscale: intermediate) arrows correspond to significantly (i.e., j ⋆ θ /∥ ⃗ j ⋆ ∥ ≥ 0.01) positive and negative azimuthal components, respectively.

Fig. 12 (

 12 corresponding to the case Bo κ S = 10 -4 for given (Re, Ha, Bo η S ) values) and of Fig. 15 (Bo κ S = 10 4 for the same (Re, Ha, Bo η S ) values) are compared. As expected from the previous results, for a given Ha = 5, the increase in Bo κ S does not result in an obvious impact (Figs. 12(a) and 15(a), Figs. 12(c) and 15(c)), irrespective of Re value (Re = 10 3 for (a), Re = 10 4 for (c)). The main difference occurs at the liquid surface, where ⃗ j becomes purely radial for Bo κ S = 10 4 .

ACKNOWLEDGMENTS

This work was supported by a national Ph.D. fellow from the French Ministry for Higher Education and Research. The laboratory SIMAP is part of the LabEx Tec 21 (Investissements d'AvenirÐ Grant Agreement No. ANR-11-LABX-0030).


Jesus Argote-Aguilar

Florin-Doru Hutu

Guillaume Villemaud

Matthieu Gautier

Olivier Berder

Alimentation par énergie RF d'une Wake-Up Radio ultra faible consommation

Abstract-Nowadays, the efficiency of the Radio-Frequency (RF) energy harvesting circuits is continuously increasing and, at the same time, the energy consumption of connected devices is drastically decreasing. Despite that, collecting, storing and delivering such kind of harvested energy to the device in an appropriate manner is still a challenge. This paper focuses on a strategy of harvesting energy in an efficient way from both low and high levels of the RF field. The objective is to power an ultralow power consumption wake-up radio requiring a regulated voltage. The starting point is the characterization of a commercial power management integrated circuit which operates with ultralow power levels. Then, some design guidelines of the RF power rectifiers, specifically designed for our envisaged application are given.

Index Terms-rectifier, RF/DC, wireless power transfer, power management integrated circuit, ultra-low power wake-up radio.

I. INTRODUCTION

In a context of an increasing deployment of the Internet of Things (IoT), replacing the battery of billions of connected devices is an expensive and exhausting task. Powering low power connected devices with energy available in the environment avoids the use of a battery or increases its lifetime. For example, [START_REF] Aoudia | Leveraging energy harvesting and wake-up receivers for long-term wireless sensor networks[END_REF] allows a Wireless Sensor Node (WSN) to efficiently use the energy harvested from ambient sources (e.g. solar, wind, vibration). Moreover, the WSN power consumption is reduced by using ultra-low power Wake-up radios (WURx), avoiding the unnecessary activation of the main transceiver. The present paper has the main goal to power a WURx consuming hundreds of nW and requiring a Direct Current (DC) regulated input voltage by RF energy harvested from different power levels of the RF field. An architecture combining two RF rectifiers, each one optimized for a specific power level, with a Power Management Integrated Circuit (PMIC) is considered. The RF/DC conversion is performed by the former and the supply of a regulated voltage is assured by the latter. The paper is organized as follows: Section II introduces the proposed architecture, Section III presents the features of a PMIC and presents measurement results of the commercial PMIC BQ25570 from Texas Instruments and Section IV describes the rectifier design considerations for our application.

II. PROPOSED ARCHITECTURE

Aiming to harvest RF power from Wireless Power Transfer (WPT) and from Self Energy Recycling (S-ER) during the transmission state of the main transceiver [START_REF] Demir | Optimum qos-aware beamformer design for full-duplex relay with self-energy recycling[END_REF], knowing that each of the mentioned sources of RF energy have different power levels (lower than -20 dBm for WPT and around -10 dBm for S-ER) the design of a rectifier for each power level is considered. The received power levels may considerably vary (depending on the propagation channel). Since the output voltage of a rectifier depends on its input power, the use of a voltage regulator is mandatory and the use of a PMIC appears as an appropriate solution. A PMIC provides an output regulated voltage independently of its input power and power consumed by its load. The schematic of the proposed architecture is depicted in Fig. 1. The proposed architecture consists in the association of two rectifiers, each one optimized for a specific input power level, in series with a PMIC. The latter, allows to appropriately load a rectifier as explained in Section III. Moreover, as explained in [START_REF] Piñuela | Ambient rf energy harvesting in urban and semi-urban environments[END_REF], in this architecture, the PMIC efficiently loads the rectifier with the highest output voltage at a given moment.

III. POWER MANAGEMENT INTEGRATED CIRCUIT MEASUREMENTS

In the context of the envisaged application, three functionalities offered by a PMIC are identified. The Maximum Power Point Tracking (MPPT) which enables a PMIC to efficiently load a rectifier for maximizing its RF/DC conversion efficiency, the management of the collected energy (among the WURx and a storage element) and the functionality that ensures to supply a regulated output voltage to the load. These functions are enabled when the PMIC is in its Normal Mode (NM) state after exiting its Cold Start (CS) state which harvests energy in a less efficient way, although necessary to collect the minimal energy required by the former state. In this paper the measured features of the BQ25570 in NM are presented. In [START_REF] Adami | Self-powered ultra-low power dc-dc converter for rf energy harvesting[END_REF], the authors claim that a rectifier may be modeled by a voltage source in series with its internal resistor R int and conclude that a maximum transfer of power occurs when the rectifier load has the same value as R int . Aiming to study the influence of R int in the minimal voltage V pmic-in and power P pmic-in required by a PMIC in NM while powering a load P pmic-load emulating the consumption of an ultra-low power WURx, the BQ25570 is connected to a voltage source in series with two different values of resistors. These elements are modeling the RF rectifier. To load the rectifier model with the same value of R int , MPPT is set to 50 % implying the output voltage of the rectifier model is fixed at half its open circuit voltage. Furthermore, a ceramic capacitor C store of 100 µF is used as storage element and its voltage V store is set to 2.5 V (to limit leakage currents). The output voltage of the PMIC is set to a regulated voltage of 1.8 V. Results are summarized in Table III. Using (1) the equivalent load R rect-load at which the PMIC loads the rectifier model is calculated as

R rect-load = V 2 pmic-in P pmic-in . (1)
R rect-load of 2.4 kΩ and 52.6 kΩ are obtained for a R int of 2.2 kΩ and 47 kΩ respectively. Furthermore, it is observed that higher R int require higher voltages, however decreasing the power absorbed by the PMIC. Actually, a minimal P pmic-in of 4.12 µW and 2.7 µW and a minimal V pmic-in of 99.1 mV and 376.8 mV were obtained for a R int of 2.2 kΩ and 47 kΩ respectively. These results give design considerations for the RF rectifier at input power levels of -20 dBm, since at low powers rectifier efficiencies are reduced and thus obtaining high voltages is difficult.

R int (kΩ) V pmic-in (

IV. LOW POWER RECTENNAS

In [START_REF] Fumtchum | A survey of rf energy harvesting circuits[END_REF], the authors show that Schottky SMS7630 diodes are well suited for low power rectifiers. Using this technology, voltage doubler, series and shunt topologies are compared in simulation at 868 MHz using Keysight ADS software. Design considerations are mainly given at -20 dBm, since a rectifier design is more critical at low powers. The strategy adopted is to design a rectifier having a high R int in order to limit P pmic-in and able to supply the minimum V pmic-in and P pmic-in required. Consequently, the three aforementioned topologies were optimized to have a maximum efficiency at 47 kΩ. The matching circuit employed for all the cases is a L-C topology where the reflection coefficient S 11 was optimized to values inferior to -25 dB for a power of -20 dBm, a central frequency of 868 MHz and a bandwith of 6 MHz. Fig. 2 shows the voltages and efficiencies of the rectifiers in function of their input powers P rect-in when loaded with their optimal loads. At -20 dBm, single diode topologies (Series and Shunt) slightly outperform the Voltage Doubler, having an output voltage of 492.6 mV and an efficiency of 49.5 % which implies an output power of 4.9 µW and a R int of 49 kΩ. With regards to the rectifier optimized at -10 dBm, since the minimal V pmic-in and P pmic-in are not a constraint, the optimization was performed favoring efficiency. At -10 dBm, it was found the shunt rectifier has higher performance with an efficiency of 77.5 % and an output voltage of 859.4 mV.

CONCLUSIONS

The use of a PMIC is an appropriate solution to power ultralow power consumption devices requiring regulated voltage supply. Although, in our application, the quiescent power of a PMIC is one order of magnitude higher than the targeted load power consumption, current voltage regulators have even or higher quiescent power consumption when operating at low powers. Moreover, a PMIC allows to maximise the RF/DC power conversion efficiency of a rectifier. After characterizing the PMIC features for the architecture proposed in this paper, it was determined that at low powers, the design considerations are mainly oriented by the rectifier impedance to limit the quiescent power of a PMIC, and at relatively higher powers, efficiency needs to be fostered.

Fig. 1 .

 1 Fig. 1. Proposed architecture for energy harvesting from WPT (-20 dBm) and for harvesting during the transmission of the main transceiver (-10 dBm).

Fig. 2 .

 2 Fig. 2. Rectifier Efficiencies and output voltages when loaded with their optimal loads.