# MHD flow in an insulating rectangular duct under a non-uniform magnetic field

- René Moreau
^{1}Email author, - Sergey Smolentsev
^{2}Email author and - Sergio Cuevas
^{3}

**Received: **23 June 2010

**Accepted: **1 October 2010

**Published: **1 October 2010

## Abstract

Followed by a review of previous studies of magnetohydrodynamic (MHD) duct flows in a non-uniform magnetic field at the entry into a magnet (fringing magnetic field), the associated MHD problem is revisited for a particular case of a nonconducting rectangular duct of a small aspect ratio *ε* = *b*/*a* (here, *b* is the duct half-width in the magnetic field direction, and *a* is the half-height). The suggested model includes a realistic three-component *div*- and *curl*-free fringing magnetic field as well as inertia terms and takes into account the mechanism of electric current exchange between the core of the flow and the Hartmann layers. The original three-dimensional flow equations are reduced to a quasi-two-dimensional (Q2D) form for three basic scalar quantities: the vorticity, the streamfunction and the electric potential. This Q2 D formulation implies that the velocity field in the core region between the two Hartmann layers does not change in the magnetic field direction and thus is two-dimensional, while the induced electric current forms both cross-sectional and axial circuits and is essentially three-dimensional. A new parameter *R* = *ε*^{2}*Re*/*Ha* has been identified to characterize the role of inertia in duct flows with insulating walls (*Re* and *Ha* stand for the Reynolds and Hartmann numbers). Computations and analytical studies are performed for inertialess (*R* ≪ 1) and inertial (*R* ≫ 1) flows at *ε* = 0.2 for Re up to 300,000 resulting in new scaling laws for typical lengths, velocities, electric current densities and pressure drops, which provide a new theoretical basis for potential applications.

**PACS Codes:** 47.65.-d, 47, 47.11.-j

## 1 Introduction

Magnetohydrodynamic (MHD) flows in a non-uniform magnetic field, or in ducts of a varying cross-sectional area, have been extensively studied since the seminal paper of Hartmann [1], where the formation of the axial current loops and the origin of the pressure losses associated with the fringing magnetic field either at the entry or exit of a magnet are qualitatively described. Later, Shercliff, in his two books [2, 3] and educational movie [4], showed that the presence of the axial currents promotes the formation of the so-called "M-shaped" velocity profile and associated vorticity. A few important experiments were made in Riga [5–8] and in Saint Petersburg [9], which revealed the spectacular behavior of such non-uniform flows. Incidentally, new phenomena of turbulence suppression in a fringing magnetic field was recently discovered experimentally in Ilmenau, Germany [10, 11]. Attempts to build a theory of those complex flows started with the papers by Kulikovskii [12, 13], where the concept of characteristic surfaces was first introduced. Although this concept is limited to the asymptotic case of a perfectly conducting fluid, it allows for qualitative predictions of the main features of many MHD flows, including those in a fringing field. Subsequently, different research groups, around Walker [14–23] and more recently Molokov [24–27], investigated such flows, using simplifying assumptions, e.g. limiting the analysis to inertialess flows. These authors obtained velocity and electric current distributions, and suggested first explanations of how these flows are organized. In particular, numerical results based on the inertialess flow model [20] have demonstrated a fair agreement with the experimental data for both pipe and rectangular duct flows in a strong magnetic field. However, other experiments for flows in a slotted duct at high Reynolds (*Re* ≈ 10^{5}) and moderate Hartmann numbers (*Ha* ≈ 10^{2}) have demonstrated strong inertial effects [9], which cannot be explained using inertialess models. As a matter of fact, except for Holroyd [28] and recent attempts [29, 30] to model the Ilmenau experiment [10, 11] most of the theoretical efforts are lacking direct comparisons with the experimental data, which often do not comply with the theoretical predictions. Quite recently, Alboussière reexamined the theory, paying a special attention to the full role of the Hartmann layers [31]. His considerations are, however, still based on a few striking assumptions that limit the model applicability.

*e.g*.

*Ha*≈

*O*(10

^{3}) -

*O*(10

^{5}) and

*Re*≈ 10

^{5}), which are sufficiently higher than those in typical laboratory experiments (see [32]). Unlike many previous theoretical studies, we extend our analysis to inertial flows at high Reynolds numbers. In the blanket conditions, the liquid-metal flows in the access ducts at the entry or exit of the blanket module experience strong MHD effects associated with spatial variations of the external plasma-confining magnetic field. This type of flows fits well with the reference problem for MHD flows in a fringing field. In our considerations of this problem, we limit ourselves to the particular case of a non-conducting rectangular duct of a uniform cross-section, located between the poles of a semi-infinite magnet as sketched in figure 1. Such flows are well understood if the magnetic field is uniform. They exhibit typical MHD boundary layers: two Hartmann layers with the thickness

*O*(

*Ha*

^{-1}) at the walls perpendicular to the magnetic field, and two side layers with the thickness

*O*(

*Ha*

^{-1/2}) at the parallel walls, which are also called Shercliff or parallel layers. The internal flow region enclosed by these boundary layers is usually referred to as a core. This terminology is also applied to the flows in a fringing field and is used throughout the paper. We select the rectangular geometry because it allows a clear demonstration of the important role of the Hartmann layers and also elucidates the influence of the side layers, which are almost degenerated in a circular duct. To have a better match with experimental conditions [9] and to simplify derivations we also limit our analysis to ducts with a small aspect ratio

*ε*=

*b*/

*a*(slotted ducts), where the duct half-width

*b*taken in the direction of the magnetic field (Hartmann length) is much smaller than the duct half-height

*a*. In our frame of reference the predominant component of the fringing magnetic field is in the

*z*-direction and varies with

*x*from zero to the constant value

*B*

_{0}. We simplify the analysis by assuming a two-dimensional fringing field produced by magnetic poles that are semi-infinite in the

*x*-direction and infinite in the

*y*-direction. The fluid flows in the

*x*-direction and

*U*

_{0}stands for the typical velocity scale (e.g. mean bulk velocity). The magnetic yoke is semi-infinite in the

*x*-direction with a length 2

*a*

_{ m }in the

*y*-direction and a gap 2

*b*

_{ m }in the

*z*-direction (see figure 1). In addition, it is assumed that the magnet aspect ratio

*ε*=

*b*

_{ m }/

*a*

_{ m }is also small. Under the conditions of the inductionless approximation (when the induced magnetic field is much smaller than the applied one), which is a common approximation in laboratory conditions and even in the blanket flows, the three basic non-dimensional parameters that describe the flow in a fringing magnetic field, are the Reynolds number

*Re*=

*U*

_{0}

*a*/

*ν*, the Hartmann number $Ha={B}_{0}b\sqrt{\sigma /\rho \nu}$, and the duct aspect ratio

*ε*=

*b/a*. Here,

*σ*,

*ρ*and

*ν*are the electrical conductivity, density and kinematic viscosity of the fluid, respectively.

Most of previous theoretical investigations are based on limiting assumptions, which are revisited here. First, since our magnetic field is *div*- and *curl*-free, we do not restrict our considerations to the so-called straight magnetic field approximation, in which only the predominant field component *B*_{
z
}(*x*) is retained (see, for instance, [16, 24, 31]). As a result, the magnetic field derived in the present paper (see Section 2) decays slowly at a large distance from the magnet as *x*^{-1}. As a consequence of this slow decay, a well-known flow structure with the Hartmann and side layers and internal core similar to that in fully developed MHD flows is established at large upstream distances from the magnet.

The exchange of current between the core and the Hartmann layers is another fundamental mechanism that we intend to readdress. Unlike previous theories, where this mechanism is not introduced or included in an approximate way, we derive the three current components directly based on the properties of the Hartmann layers, the vorticity equation and by satisfying the charge conservation equation. This is different, for example, from the approach in Ref. [31] by Alboussière, where the exchange of electric currents between the core flow and the Hartmann layers is included indirectly by integrating (averaging) the current density distribution. Such an approach seems not to be valid in the flow region where the M-shaped velocity profile is formed and the local current density can differ significantly from its averaged value. Despite this limitation, Alboussière's analysis has the merit to renew the basic scaling law in this problem, introducing a new important length scale of the order of *Ha*^{-1/4} in the case of a pipe flow.

A third strong common assumption is the *inertialess approximation* [16] in which inertia in the core flow is neglected in comparison with the Lorentz force and the pressure gradient. This implies that, except the side layers where viscosity is significant, the Lorentz force is purely irrotational and thus cannot generate any vorticity. The rationale for this approximation is the assumption that the ratio of the Lorentz to inertia forces scales as the interaction parameter *N* = *Ha*^{2}/*Re*, which in fusion blanket conditions is much larger than unity. However, in electrically insulating ducts, this scaling law appears not to be correct any longer providing the exchange of current between the core and the Hartmann layers is taken into account. As a consequence, the Lorentz force is in fact *Ha* times smaller, and the relevant parameter is *Ha*/*Re*. As applied to the fusion blanket conditions, this suggests that in the blanket flows inertial effects are not necessarily negligible.

Notice that among various designs of a liquid-metal blanket, the EU Helium-Cooled Lead-Lithium (HCLL) concept requires low velocities (a few mm/s in the blanket ducts and a few cm/s in the feeding duct), whereas in the US Dual-Coolant Lead-Lithium (DCLL) concept much larger velocities are required for heat and tritium removal and proper cooling (of the order of 10 cm/s in the poloidal ducts and about 50 cm/s in the access ducts). Even higher velocities, of the order of 1 m/s or higher, may be needed in a self-cooled blanket. Therefore, to cover the whole velocity range in the blanket applications, our goal is to consider the flows in a fringing magnetic field starting from regimes with zero or weak inertia to highly inertial regimes.

The structure of the paper is the following. Following this introduction, a *div*- and *curl*-free magnetic field is derived in Section 2. This magnetic field is used in Section 3 in the derivations of the basic quasi-two-dimensional (Q2D) equations governing the reference flow in a fringing magnetic field. These equations are obtained from the original three-dimensional equations in terms of three scalar variables, vorticity, stream function and electric potential assuming that the flow variations in the core in the magnetic field direction are negligible, whereas the induced electric current is essentially three-dimensional. Then, in Section 4, typical numerical results are shown. The derivation of explicit asymptotic solutions is the purpose of Section 5, with a special emphasis on regimes with strong inertia and new scaling laws for typical lengths, velocities, electric current densities and pressure drops as compared to previous theories. Finally, a discussion and concluding remarks are given in Section 6.

## 2 A div- and curl-free magnetic field

### 2.1 Beyond the straight magnetic field approximation

**B**is essentially pointing in the

*z*-direction, so that its components can be written in terms of a series expansion of increasing powers of

*z/a*

_{ m }, which is itself of the order of the small parameter

*ε*

_{ m }=

*b*

_{ m }/

*a*

_{ m }where

*b*

_{ m }≪

*a*

_{ m }(see figure 1). Therefore:

*B*

_{0}stands for the uniform magnetic field, which is supposed to exist in the limit

*x*→∞. At the leading order in terms of

*ε*

_{ m }≪ 1, equations (1) satisfy both the conditions ∇ ·

**B**= 0 and ∇ ×

**B**= 0, precisely because the apparently small components

*B*

_{ x }and

*B*

_{ y }are not neglected. Indeed, whereas these components are

*O*(

*ε*

_{ m }) and smaller than the

*B*

_{ z }-component, their

*z*-derivative is of order unity. The same is true for the components of the current density, which from the Ohm's law can be written as:

**j**= 0, takes its classic form

only if the *z*-dependent component of the electromotive field **u** × **B** in equation (2*c*) is taken into account. The term ∂_{
z
} *j*_{
z
} is of order unity even though *j*_{
z
} is much smaller than *j*_{
x
} and *j*_{
y
}. Its contribution appears to be essential as it expresses the leakage of current from the axial loop toward the Hartmann layers.

The aim of this section is to get a correct expression for the function *B*(*x*, *y*), which is in fact one of the main parameters of this problem. Note first that *B*(*x*, *y*) should vary as *x*^{-3} or *y*^{-3} at large distances from the gap, if the magnetic field were a three-dimensional dipole. If this dipole were two-dimensional, as in the case of infinite magnet poles in the *y*-direction, the far field should vary as *x*^{-2}. In the reference case, when the poles are semi-infinite in the *x*-direction and infinite in the *y*-direction, the far potential is logarithmic and the far field must vary as *x*^{-1}. If so, the Hartmann layers still exist at a long distance from the magnetic poles. Indeed, considering that their thickness is proportional to $\sqrt{\rho \nu /\sigma}/{B}_{z}$, the overall length over which the flow exhibits Hartmann layers is *b*_{
m
}*Ha* in the case of a two-dimensional very long magnet. It is *b*_{
m
}*Ha*^{1/2} in the case of a two-dimensional short magnet, and *b*_{
m
}*Ha*^{1/3} in the case of a three-dimensional short magnet.

### 2.2 A correct expression for the two-dimensional fringing field at the entry of a magnet

*x*-direction and infinite in the

*y*-direction. Consider the conformal mapping between complex planes

*r*and

*t*defined by the relation

where *K* is a complex constant which has to be determined. This mapping transforms the half-plane *Im*(*t*) > 0 into the region between the symmetry axis and the yoke of the semi-infinite magnet in the *r*-plane. Due to the exponent 1/2, rotation of the vector representing the complex number *t* - 1 by *π* between the values *t* = 1 + *η* and *t* = 1 - *η* (*η* ≪ 1 is a real number), corresponds to rotation of the vectors (*t* - 1)^{1/2} and *dr* by *π*/2. Therefore, the suggested transformation can reproduce the sharp corner of the yoke, around which the potential lines of the magnetic field turn in the fringing field region.

*t*-plane,

*K*= -

*ib*

_{ m }/π (

*i*

^{2}= - 1), we easily get the

*x*-dependence of the

*B*

_{ z }=

*B*(

*x*) component, which varies from

*B*

_{0}for

*x*≪ 0 to zero when

*x*→∞. Its expression may be written under the parametric form

*C*is an arbitrary constant, which can be defined from the condition

*B*(

*x*= 0) =

*B*

_{0}/2. After changing the sign of

*x*in order to get

*x*increasing in the flow direction as

*B*(

*x*) from zero to

*B*

_{0}, and denoting

*a*and

*b*denote the half-lengths of the duct cross-section, whose thickness

*t*

_{ w }is neglected so that

*b*

_{ m }=

*b*and

*ε*

_{ m }=

*ε*. The slow decay of this magnetic field to zero can be seen in figure 2; note that the asymptotic behavior

*β*= -2

*ε*/π

*X*without the logarithm is justified as soon as

*x/b*

_{ m }> 5. Also, it is easy to show that, in the vicinity of the symmetry plane (

*z*= 0), the

*B*

_{ x }component is proportional to

*z*and to the derivative

*d*

_{ x }

*B*, in agreement with equation (1

*a*). The two-dimensional

*div*- and

*curl*-free magnetic field defined from equations (1) and (8) is used as the reference field for numerical and analytical solutions in Section 4 and Section 5.

## 3 The basic set of equations for the flow in a fringing magnetic field

### 3.1 Comments on the characteristic surfaces

Consider the flow of an incompressible and electrically conducting fluid in the rectangular duct shown in figure 1. Since the duct width is constant (2*b*), and since *ε*_{
m
} ≪ 1, the characteristic surfaces coincide with the cross-sections *x* = const. As shown by Kulikovskii [12, 13] and Holroyd & Walker [16], neither the streamlines, nor the electric current lines can cross them within the core flow.

*dA*and length

*2b*, located between the Hartmann walls and moving along the trajectory of the fluid across the fringing field. Because of mass conservation, its volume

*2bdA*is invariant. Now, if we admit, as Kulikovskii, that the magnetic flux

*BdA*across the area

*dA*is also invariant, taking the ratio of both invariants, we get the property that

*2b/B*and therefore

*B*must also be invariant along a fluid trajectory. In other words, this asymptotic limit only allows the streamlines and the current lines to lie on the planes where

*B*is constant, and does not allow them to cross these planes. However, in fluids of finite electrical conductivity, the generalized form of the Alfvén theorem

Where Φ is the magnetic flux, implies that the electric current crossing the characteristic surfaces is determined by the Ohmic losses (locally *j*/*σ*). As a consequence, the whole Hartmann layers (not only their portions overlapping the side layers), where these Ohmic losses have the highest values, must be taken into account on both sides of the core. This differs from the Holroyd & Walker [16] ideas, according to which the only possibility for the fluid to cross the characteristic surfaces is within the side layers.

#### 3.2 Core flow and side layers

*q*(

*x*,

*y*,

*z*) can be expressed as

*q*

_{0}(

*x*,

*y*) + (

*z*

^{2}/2

*a*

^{2})

*q*

_{1}(

*x*,

*y*) +⋯ This is particularly important for the electric potential

*z*-dependent part is responsible for the

*j*

_{ z }component of the current density, which feeds the Hartmann layers. As for the velocity field, the difference between its purely two-dimensional part and its complementary part is such that:

*w*-component is a small

*O*(

*ε*

^{3}) term, so that the incompressibility condition writes:

*z*-direction:

*ε*, from equations (2), (10) and (11), the components of the current density within the core flow are:

*S*

_{0}and the derivatives of the function

*B*(

*x*,

*y*), would not be present in the expression of

*j*

_{ z }if the approximation of a

*straight magnetic field*were used. Although we intend to focus on the particular example where

*B*is only

*x*-dependent, it is worth keeping the full expression of

*S*

_{0}in the derivation of our general model equations and to limit ourselves to

*B*(

*x*) later on. The expression of the charge conservation principle, which comes out from (14), is:

Note the presence in (16) of ${\partial}_{z}{j}_{z}={\phi}_{\text{1}}/{a}^{\text{2}}$, which represents the exchange of current between the core and the Hartmann layers, and which, like the other terms, is only a function of *x* and *y*.

### 3.3 Hartmann layers

*z*= -

*b*, whose thickness is $\delta (x,y)=(1/B)\sqrt{\rho \nu /\sigma}$. In addition, let us denote the dimensionless distance from the wall as

*ξ*= (

*z*+

*b*)/

*δ*. It can be easily checked that, at the leading order in terms of

*Ha*

^{-1}, the usual exponential distributions for the

*x*- and

*y*-components of the current density and the velocity are still valid. However, in general, it is necessary to introduce the possibility of having a significant variation of the electric potential across the Hartmann layer,

*φ*

_{ H }. This is essential when the wall has a non-zero electrical conductivity and forces some jump of the electric potential between its value at the edge of the layer and that at the wall:

*u*,

*v*,

*j*

_{ x },

*j*

_{ y }and

*j*

_{ z }are then:

*z*= -

*b*,

*ξ*= 0), they become:

### 3.4 Within the Hartmann wall

*x*- and

*y*-components of the current density within the wall are similar to those given by (19), except that the conductivity is now changed to

*σ*

_{ w }. Using the usual assumption of a thin wall of thickness

*t*

_{ w }and the charge conservation equation

*j*

_{ wz }stands for the current density exiting the wall and entering into the fluid, we get a new important relation between the three potentials

*φ*

_{0},

*φ*

_{1}and

*φ*

_{ H }:

*C*=

*σ*

_{ w }

*t*

_{ w }/σ

*b*as a fourth non-dimensional parameter, can still be significantly simplified when

*Ha*≫ 1, since

*C*Δ

_{⊥φH}≪

*φ*

_{ H }/

*bδ*. Then (21) becomes:

### 3.5 Influence of the Hartmann layer and Hartmann wall on the core electric current

*a*) and (23

*b*), denoted by

*g*

_{ x }and

*g*

_{ y }, must satisfy the global equation

*Ha*numbers, equation (25) reduces to

*φ*

_{ H }can be substituted into (22) to get

*φ*

_{1}now yields the equation relating

*φ*

_{0}(

*x*,

*y*) to the velocity field:

*C*= 0):

*j*

_{0z}as when the magnetic field is uniform:

### 3.6 Equation of vorticity

*z*-component of ∇ × (

**j**×

**B**), which reads (

**B**· ∇)

*j*

_{0z}- (j

_{0}· ∇)

*B*since both

**B**and

**j**

_{0⊥}are divergence-free. Using equations (1) for the magnetic field, the vorticity budget reads:

*C*= 0 (insulated duct), we can further modify equation (32) by using expressions (14) for the

*x*- and

*y*-components of the core current density

**j**

_{0⊥}, equation (31) for j

_{0z}, and equations (29) and (30) for

*φ*

_{0}and

*φ*

_{1}, respectively. The first and second terms on the right-hand-side of this new equation contain a purely two-dimensional term and a parabolic

*z*-dependent term:

*ε*≪1, the core flow can be modeled with a purely two-dimensional velocity field, by keeping only

*z*-independent terms in equation (32):

At the same time, providing that *ω*_{0},*ψ*_{0}, *φ*_{0} and *φ*_{1} are known, the current density given by (14), remains three-dimensional since it contains the *z*-dependent terms and *j*_{0z}≠ 0.

Equation (32) could also be used to derive a higher order equation for *ω*_{1}(*x*, *y*) by collecting terms proportional to *z*^{2}. In this investigation, however, we limit ourselves to a quasi-two-dimensional approximation based on equations (35) for *ω*_{0}(*x*, *y*), (29) for *φ*_{0}(*x*, *y*), (30) for *φ*_{1}(*x*, *y*), and the condition of incompressibility Δ_{⊥}*ψ*_{0} = -*ω*_{0} for *ψ*_{0}(*x*, *y*). In this closed set of equations, the vorticity production due to the non-zero derivatives ∂_{
x
}*B* and ∂_{
y
}*B* is clearly exhibited by the first and second terms on the right-hand-side of equation (35). The Hartmann damping of vorticity is represented by the third term $\left(-\sqrt{\sigma \rho \nu}B{\omega}_{0}/b\right)$. The additional viscous dissipation is represented by the last term, which is significant only in the side layers. This basic equation may be seen as a generalization of the Sommeria & Moreau [33] model equation to the case of non-uniform magnetic fields. It is also worth mentioning that, after an integration of equation (32) between the two Hartmann walls, a more complete model might be obtained, which would additionally take into account the contribution of the parabolic *z*-dependent terms. The derivation of the complete form of the vorticity transport equation would, however, require the introduction of terms proportional to *z*^{4} in the basic expression (10) for the electric potential. Such a model, would be a generalization of the approach presented in [34], allowing for high *ε* values.

Notice that in the particular case *B* = *B*(*x*), equation (35) can also be interpreted as an expression for the axial part of the current density component *j*_{0x}(*x*, *y*, *z* = 0) = -σ(∂_{x}φ_{0} - *Bv*_{0}), parallel to ∇_{⊥}*B*, which is expressed through other unknowns (velocity and vorticity). In other words, this equation explicitly demonstrates that, within the core, electric current lines must cross the characteristic surfaces due, essentially, to two main causes: inertia and Hartmann current sheet. As a consequence, a full three-dimensional electric current distribution can accurately be analyzed (see Section 5), in which *j*_{0x}and *j*_{0z}are derived from equations (35) and (31), while *j*_{0y}can be easily obtained from the charge conservation equation.

### 3.7 Basic set of non-dimensional equations

*B*

_{0}, the mean bulk velocity ${\mathcal{U}}_{0}$, the duct half-width

*a*as a length scale for

*x*and

*y*, and the half-height

*b*between the Hartmann walls for

*z*. The new variables are denoted as follows:

*S*denotes the dimensionless expression of

*S*

_{0}. Since all key functions are only dependent on

*x*and

*y*, in what follows we drop the subscript ⊥. Then, the main set of equations, coming from their dimensional counterparts (29), (35) and (13), is:

*Ha*in (37) represents the exchange of current between the axial loop and the Hartmann layers. It is negligible when

*β*≈ 1, but it may become significant far upstream when

*β*≈

*Ha*

^{-1}.To understand the coming analysis, it is useful to summarize here other dimensionless relations, namely,

The above expressions for *J*_{
X
} and *J*_{
Y
} are, of course, valid to the leading order in terms of *ε*.

### 3.8 The particular case *B* = *B*(*x*)

We now restrict the analysis to the case of a fringing field of the type studied in Section 2.2 where the magnetic field is *y*-independent and given by the function *β*(*X*) (see equation (8)).

*C*= 0). The basic equations become:

*R*=

*ε*

^{2}

*Re*/

*Ha*and Γ = (

**U**· ∇) Ω. As discussed below,

*R*results to be the suitable scale for inertia. Once equations (43) and (40

*c*)) are substituted in the electric charge conservation equation, we get

*X*-dependent, so that it can be approximated as $\Omega (X,Y)=-{\partial}_{Y}^{2}\Psi =-{\partial}_{Y}U$ and $\Delta \Omega ={\partial}_{Y}^{2}\Omega $, equation (44) can be integrated once. Using the boundary conditions

*U*=

*J*

_{ Y }= 0 at the insulating side wall

*Y*= -1, yields:

*J*

_{ X },

*J*

_{ Y },

*J*

_{ Z }) exhibit the noticeable feature that the three-dimensional electric current distribution is a combination of two families of well defined loops. One, the

*z*-independent axial current loops, located in the symmetry (

*X*,

*Y*) plane, can be expressed in terms of the streamfunction

*h*

_{ a }(

*X*,

*Y*) in the form:

*X*-dependent, is located in the cross-sectional (

*Y*,

*Z*) planes. In terms of the stream function

*h*

_{ c }(

*X*;

*Y*;

*Z*), the loops, quoted as cross-sectional loops, are expressed as:

One should notice that equations (46) and (47) show explicitly the exchange of current between the core and the Hartmann and side layers. It is also clear that in the limit of negligible inertia (*R* = 0), the three current density components in the core have the same scaling law, all being proportional to *Ha*^{-1}. Equations (46) have also the merit of exhibiting the contributions of inertia and viscosity on the axial current loop.

*X*-derivative of the charge conservation equation is:

*J*

_{ x }in (48) to get

which is the basic equation for the vorticity in the fringing field. This equation expresses all the relevant contributions to the transport of vorticity, which are, respectively, the advection of vorticity, Hartmann damping, damping in the side layers, torque and exchange of current. Note that the advective term, Δ(*RΓ*/*β'*), involves the dimensionless parameter *R* = *ε*^{2}*Re/Ha*, which results as the suitable measure of the importance of inertia, in contrast with previous theories where the relevant scaling for inertia was *N*^{-1} = *ε*^{2}*Re/Ha*^{2}.

## 4 Numerical results

Although the equations presented in previous sections are quasi-two-dimensional, their accurate numerical solution requires serious efforts, in particular if the Hartmann and Reynolds numbers are high, since computations become more time and memory demanding. Unlike a fully three-dimensional approach, the present one has the advantage of not computing the flow in the Hartmann layers. However, resolution of high-gradient flow regions associated with the side layers of thickness $b/\sqrt{Ha}$ at the walls *y* = ± a is still needed to reproduce correctly the axial and cross-sectional current loops and the current exchange between the bulk flow and the boundary layers. Another important feature of the flow, which should be properly treated in the computations, is the formation of a thin transverse layer at the entry to the magnet (see Section 5). It should also be taken into consideration that the reference flow exhibits high velocity jets once the liquid enters the magnet zone. Such velocity profiles present inflection points, it is thus important that the computed results can indicate whether or not the flow loses its stability and develops a time-dependent behavior. Then, computing time-dependent flows needs small time increments and special approximations of the convective terms to minimize the schematic viscosity and to keep stable computations at the same time.

In the past, many computations of MHD flows demonstrated severe restrictions on the magnitude of the Hartmann number, even for relatively simple and well-understood MHD flows, *e.g*. fully developed flows in a rectangular duct under a uniform magnetic field. Typically, the Hartmann number in such studies was limited to *Ha* ≈ 10^{2} [35, 36]. The limitations were to some extent related to the computational time, which increases roughly exponentially with *Ha*. However, the most severe limitations are caused by a catastrophic loss of accuracy at high *Ha* numbers, which may occur in the course of computations due to accumulation of round-off errors. A potential source of computational errors can be related to the use of the electric potential as the electromagnetic variable. As discussed by Smolentsev *et al.* [37], the loss of accuracy occurs when an induced electric current of small magnitude (*O*(*Ha*^{-1})), is calculated from the Ohm's law as a difference between two large numbers, associated with -∇Φ and **U** × **B**_{0}, which are both *O*(1). In the present formulation, the nature of this problem can be seen directly from the equation for the potential (37). Two terms on the right-hand-side of this equation, *β*Ω and Ω/*Ha*, are different by a factor Ha, once the magnetic field is almost uniform. However, ignoring the term Ω/*Ha*, the whole physical mechanism responsible for the current exchange between the bulk flow and the Hartmann layers would be missed. First numerical tests showed that using the electric potential in the computations for the reference problem is reasonable only for moderate values of *Ha*, limited to a few tens, but the above mentioned accuracy problem becomes extremely critical at higher *Ha*.

Usually the accuracy problem associated with the electric potential does not appear if the alternative formulation based on either the induced magnetic field or electric current is used. For example, when the induced magnetic field is used to calculate the electric current from Ampere's law, the continuity equation for the electric current is automatically satisfied based on the vector identity *div(rot)* ≡ 0. Therefore, the computed electric current always satisfies the charge conservation equation. However, when applying the approach based on the induced magnetic field, the boundary conditions should be formulated far enough from the flow domain, where the induced magnetic field vanishes. It is therefore more convenient to use the electric current, namely its *J*_{
X
} -component, which obeys equation (48), as an electromagnetic variable, since the boundary conditions can be formulated at the interface between the liquid and the wall (insulating duct) or at the external surface of the wall (conducting duct).

The boundary condition at the walls *Y* = ± 1, which are assumed to be non-conducting, is ∂_{
Y
} *J*_{
X
} = 0, as follows directly from the expression for the *Z*-component of ∇ × **J** and taking into account that *J*_{
Y
} = 0 at the walls. The *J*_{
X
} current density component is then used to calculate the Lorentz force term on the right-hand-side of the vorticity equation. In the post-processing phase, the other current component, *J*_{
Y
} , can be derived from (50). After computing the current density components *J*_{
X
} and *J*_{
Y
} , the whole three-dimensional current distribution, including *J*_{
Z
} , can be reconstructed. This formulation was intensively tested recently for both steady and unsteady flows [37]. The tests demonstrated very good match with analogous computations based on the traditional induced magnetic field formulation.

In the code, a finite-difference technique using a non-uniform grid in the *Y* -direction that clusters grid points in the side layers is applied. The whole problem is solved by advancing in time, until a steady-state solution is achieved. The vorticity equation at each time step is solved using the TDMA (Tridiagonal Matrix Algorithm) method in both coordinate directions. The convective terms in the vorticity equation are approximated with the central-difference scheme, which is known to be dissipation free. Using convective term approximations free of numerical dissipation is particularly important when the flow is essentially time-dependent. However, in the chosen range of parameters, including high values of *Re* (up to 300,000), the numerical solution always converges to a steady-state regime without showing any tendency to instability. The elliptic equations for the streamfunction (39) and for the electric current (48) are both solved with the cos fast Fourier transform (FFT) in *X* and the tridiagonal algorithm in the other direction. The use of the FFT technique limits the computational grid in *X* to a uniform spacing, but is more effective compared with much slower relaxation techniques since solution of these elliptic equations may take over 90% of the total computational time. Using the cos FFT in *X* also requires the inlet/outlet boundaries to be positioned far enough from the fringing field region, so that Neuman boundary conditions *∂*_{
X
} = 0 are allowed. Typical calculation domain is 40-60 units. Satisfactory resolution is achieved by using 512 or 1024 grid points in *X*-direction and 201-501 points in *Y* -direction. The numerical tests included grid sensitivity tests, and comparisons with the analytical solution for a fully developed flow and unsteady vortical flows [37, 38]. The comparisons demonstrated coincidence in local values between the analytical and the numerical results in four digits.

### 4.1 Numerical results for flows with weak inertia (*R*≪1)

*R*is the relevant criterion for inertia, so that, when

*R*≪ 1, the flow in the fringing magnetic field can be considered inertialess, while at

*R*≫ 1 inertial effects become important. Figure 3 illustrates computations for high Hartmann number (

*Ha*= 1000). Additional computations in the same range of

*Ha*values without nonlinear terms have demonstrated that for R ≪ 1, the Reynolds number does not affect the results, so that the Hartmann number and the aspect ratio are the only relevant parameters. A similar example, computed for

*Ha*= 2000,

*Re*= 500,

*ε*= 0.2,

*R*= 0.01 can be found in a short paper published elsewhere [39].

*U*

_{0}(

*X*) decreases, the near-wall jets form, and the maximum velocity

*U*

_{1}(

*X*) in these jets increases (typical data are shown in figures 4 and 5). In agreement with the previous results by Lavrent'ev

*et al.*[24], these features are spectacular, although

*U*

_{0}(

*X*) does not reach a zero value (see figure 4) and the maximum value of

*U*

_{1}(

*X*) increases as

*U*

_{1max}≈ 0.3

*Ha*

^{1/2}instead of 2.38

*Ha*

^{1/2}, our theoretical prediction for the jet region (see Section 5). Each jet exhibits two vorticity peaks of opposite signs. The value of the higher peak, which has the same sign as the side layer vorticity, can be estimated as

*U*

_{1}/

*δ*

_{ SL }, where

*δ*

_{ SL }stands for the thickness of the side layer, and yields an estimate of

*U*

_{1}. The other peak, of the opposite sign and much smaller value, can be estimated as

*U*

_{1}/

*δ*

_{ j }, where

*δ*

_{ j }yields an estimate of the jet thickness. The large ratio between the two peaks also indicates the large ratio between the jet thickness

*δ*

_{ j }and that of the side layer

*δ*

_{ SL }, as shown in Section 5.

In the asymptotic case of large *Ha* and small *R*, these computed results suggest that the whole flow domain can be subdivided into five characteristic regions: i) the *far upstream developing Shercliff-like flow*, where the velocity profile is uniform with a moderate over-velocity to compensate the missing flow rate in the side layers; ii) the *upstream region where the M-shape starts to form*, whose position is close to abscissa *X* ≈ -15 and varies very slowly with *Ha*; iii) the *jet-region*, where most of the flow rate is transported by two jets located on both sides of the almost motionless core, and whose length seems to be weakly *Ha*-dependent; iv) a short *transverse layer*, whose length seems significantly smaller than the duct half width, suggesting that the *X*-derivatives should be locally larger than the *Y* -derivatives; and finally v) the *fully-established Shercliff flow* with its familiar plateau-type velocity between the side layers. In Section 5, a theoretical attempt is taken to derive the most important characteristics for the three intermediate regions where the flow is strongly affected by the variation of the magnetic field. The fact that the jet formation occurs over a distance of a few units, which agrees with the present asymptotic theory (see Section 5), seems to be one of the newest ideas arising from this investigation.

The electric current distribution, also shown on figure 3, shows a transition from a three-dimensional structure in the jet region to a purely cross-sectional current distribution in the fully-established Shercliff flow. In turn, in the jet region the three-dimensional current density demonstrates two characteristic families of current loops: an internal one where the current lines are almost quasi-two-dimensional and closed in the core, and an external one where the current lines passing through the jets and side layers are strongly three-dimensional and close themselves within the Hartmann layers. Additional numerical results are presented in Section 5, within the framework of the theoretical ideas that allow the interpretation of these results in terms of scaling laws, specific to each sub-domain.

### 4.2 Numerical results for inertial regimes

*R*characteristic of inertial regimes, the Reynolds number is substantially increased (up to 300,000), compared to the calculations in the inertialess regime, while the maximum Hartmann number is limited to 200. As a result, the highest value of

*R*achieved in this set of computations is 120 (for

*Re*= 300,000,

*Ha*= 100, and

*ε*= 0.2), while the smallest one is 8 (for

*Re*= 50; 000,

*Ha*= 100, and

*ε*= 0.2). The computational results in a regime with significant inertia are shown in figure 6, while those in a regime with strong inertia are presented in figure 7. In both inertial regimes, the flows are similar but pronouncedly different from those in a regime with negligible inertia. Namely, once formed, the jets are convected far downstream, displaying a gradual transition to the Shercliff flow over a long distance of the order of tens of the characteristic channel dimensions. The other particular feature is that the jets themselves are much less pronounced than the high velocity jets in the inertialess regime, computed at higher or even the same Hartmann numbers as shown, for instance, in figure 5. In fact, at the maximum

*R*-value computed (

*R*= 120, figure 7), the jets in the fringing field region can hardly be seen, although the vorticity distribution still demonstrates pronounced variations indicating that the effect of the magnetic field on the flow structure is still detectable. Although small, in the inertial flow the variations of the maximum velocity in the jet can be expressed, as a first approximation, as a function of the parameter

*R*(see Section 5). The induced electric current distribution in the fringing field region is also different from that in the inertialess flows. Namely, the induced current is mostly represented by the axial loop, so that the current lines are remarkably two-dimensional. This indicates that the current exchange mechanism, important in inertialess flows, becomes negligible in flows where inertia is significant or strong.

## 5 Guiding ideas and scaling laws for flows with high *Ha* and small *ε*

### 5.1 General properties

Figures 6 and 7 show that, when inertia is significant, the flow is quasi-parallel all over the fringing region. When inertia is weak, the numerical results suggest that this may still be true, but only in the far upstream region (*β*^{2} ≪ 1) and in the jet region, where the streamlines exhibit a very small curvature (see figure 3). Therefore, let us pay attention to conditions where it can be assumed that ${\partial}_{X}^{2}\ll {\partial}_{Y}^{2}$ in the Laplacians of equations (39), (41), (48) and (49). Of course, there are situations, for instance when inertia is weak, where this approximation is not valid and consequently deviate from the coming analysis. Inertiales regime is studied in Section 5.3.3 where an attempt is presented to analyze the short transition region toward the downstream Shercliff flow. In this case, the flow is far from being quasi-parallel and hence, ${\partial}_{X}^{2}$ becomes predominant upon ${\partial}_{Y}^{2}$ (see figure 3).

*ε*

^{2}

*Z*

^{2}so that the components of the Lorentz force read

*X*- and

*Y*-projections of the motion equation read:

*J*

_{ X }can also be simplified as

*Ha*≫ 1, yielding the simple relations

*h*

_{ a }(

*X*,

*Y*) have any side layer. They are completely determined by equations (53) and (54), when the velocity field in the core is known. In particular, denoting the pressure on the axis by

*P*

_{0}(

*X*), and its value within the side layer by

*P*

_{1}(

*X*), equations (53) allow deriving the pressure from the following equations

*U*=

*h*

_{ a }= 0, equation (56

*b*) yields the

*X*-dependence of the pressure at this location:

*b*) from its form at the edge of the layer, where the velocity reaches its maximum value

*U*

_{1}, the relevant equation for the side layer is

*R*= 0), in the vicinity of the side wall

*Y*= -1, an elementary solution for the side layer exists, namely,

*δ*

_{ SL }(

*β*), which can be derived from a global balance of momentum. If we take

*β*as our new abscissa, with

*∂*

_{ X }=

*β'∂*

_{ β }and

*β'*(

*β*) given by equation (8), and integrate (58) from the wall to the edge, we get

*R*is large enough, the Hartmann friction becomes negligible and equation (60) implies the invariance of the missing momentum: ${U}_{1}^{2}{\delta}_{SL}=\text{Const}$. Consequently, evaluating the Const in the far downstream Shercliff flow, we get the following evolution for the side layer thickness, in terms of the still unknown jet velocity

*U*

_{1}(

*β*):

*P*)

_{2}is weakly dependent on

*R*through the

*R*-influence on the Ψ(

*Y*) profile, whereas the head loss (Δ

*P*)

_{1}significantly depends on

*R*. In the asymptotic case where inertia is strong enough to justify the estimate Ψ ≈

*Y*, this pressure difference across the duct is

*R*-independent and reads:

Its maximum (Δ*P*)_{2} = (*πHa*/10*ε*)(3/5)^{3/2} ≈ 0.146(*Ha*/*ε*) is reached for $\beta =\sqrt{3/5}\approx 0.7746$.

### 5.2 Inertial regimes

**U**· ∇)Ω =

*∂*

_{ X }Ω =

*β'∂*

_{ β }Ω. Hence (49), which is the basic equation for the vorticity, reads

The terms of this linear equation express all major effects: from left to right, we have advection of vorticity $\left(R{\beta}^{\prime}{\partial}_{Y}^{2}{\partial}_{\beta}\Omega \right)$, Hartmann damping $\left(\beta {\partial}_{Y}^{2}\Omega \right)$, damping in the side layers $\left({\epsilon}^{2}{\partial}_{Y}^{4}\Omega /Ha\right)$, torque (*Haβ'*^{2}Ω) and exchange of current (*β*'^{2}∂_{β}Ω).

*αY*and cosh

*αY*, such that the total flow rate is conserved. The functions sin

*αY*and cos

*αY*, which would also conserve the flow rate, cannot be selected because they locate on the axis the maximum curvature of the velocity field. In addition, we introduce two parameters, a shape

*α*(

*β*) factor and an amplitude

*F*(

*β*), which become the main unknowns of the problem:

*FY*as soon as

*αY*≪ 1, what is always true in the vicinity of the axis. Substitution of (68) in the local equation (67) yields the first order differential equation:

Note that an additional term of the order of *U*_{1}*δ*_{
SL
}*Y* might be introduced in the expression for Ψ, to take into account the missing flow rate in the side layers. However, as suggested by (61), this term is small enough to be neglected in this asymptotic approach.

*α*(

*β*) and

*F*(

*β*) can be derived from a global form of equation (67) integrated between the side wall

*Y*= -1 and the axis

*Y*= 0:

*ε*

^{2}/

*Ha*≪ 1 and ${\left({\partial}_{Y}^{4}U\right)}_{0}=O(1)$. Besides, the last three terms cancel exactly, since their addition is just a second derivative of the side layer equation (58). Equation (70) represents the non-trivial contribution of the side layer to the global balance of vorticity. Notice that the axial velocity

*U*

_{0}and maximum velocity in the jet

*U*

_{1}, can be written as

*F*(

*β*) and

*α*(

*β*) are given by the solution of equations (69) and (73), relations (68), (71) and (72) determine completely the velocity distribution at any abscissa

*X*or

*β*. In particular, the normalized velocity profile is given by

*α*significantly higher than unity would imply two thin jets along the side walls, whereas values lower than unity would imply an M-shape extending over the whole width of the duct. The amplitude of the M-shape prifile can be characterized by the quantity

Then, the induced electric current is easily derived from the expressions of the streamfunctions *h*_{
a
} and *h*_{
c
} introduced in (46), (47) and (55), and the pressure follows from equations (56) and (62).

*α*≤ 1, to the general expressions of

*f*(

*α*) and

*g*(

*α*) in equations (72),

*F*(

*α*) and

*G*(

*α*) =

*α*

^{2}

*F*:

*a*) and (77

*b*) and the expression (8) for

*β' =*(

*π*/2)

*β*

^{ 2 }(1

*-β*

^{ 2 }), suggest the transformation

*a*) and (77

*b*) become

*C*and introducing the rapidly varying variable

*t'*=

*Haβ'/R*, yields the second order equation

*β*→ 0. Nevertheless, it can be checked in the following that the inlet vicinity is indeed irrelevant, because the functions

*C*and

*D*, as well as

*F*and

*G*, rapidly vanish when

*t*→ 0. Thus, (81) reduces to a classical equation, whose solution satisfying the upstream conditions

*C*(0) =

*D*(0) = 0 reads:

*F*(

*β*) and

*G*(

*β*) are therefore given by:

*U*

_{0}= 1 -

*f F*, and at the edge of the side layer,

*U*

_{1}= 1 +

*gF*:

*β →*0:

In particular, in most experimental conditions, these velocities deviate from their asymptotic limit as *Nβ*^{3}/*ε*, where *N* = *Ha*/*R* is the interaction parameter. This shows that the typical length of the upstream region where the M-shape profile is formed, scales as *L*_{
up
} ≈ (*ε*^{2}*Ha/R*)^{1/3} and becomes shorter the higher the *R* values. This can be interpreted as the squeezing of the upstream region against the fringing field by the advection of vorticity, clearly observed in figures 6 and 7.

*β*→ 1, we get 1 -

*β*= (2/3)e

^{-πX/ε}, so that expressions (84

*a*) and (84

*b*) become:

*N*=

*Ha*/

*R*is strong enough, or inertia remains moderate to allow the exponential

*e*

^{-t/12 }to reach its zero limit, then (86

*a*) and (86

*b*) reduce to

*L*

_{ down }≈

*R*and is

*Ha*-independent. Now, in the opposite condition where inertia is strong, such that π

*Ha*/180

*εR*≪ 1, the term

*e*

^{-t/12}can be approximated by 1 -

*t*/12, so that equations (86

*a*) and (86

*b*) reduce to

The correction for the *U*_{0} value is small since *t*/6 ≪ 1, but the correction for *U*_{1} is important and implies a reduction of its maximum value by a factor of the order *Ha/εR*. Nevertheless, the typical length necessary to reach the Shercliff profile remains controlled by e^{-X/R}and scales as *L*_{
down
} ≈ *R*. This downstream transport of the M-shape over a long distance is also clear on figures 6 and 7.

*b*) for flows at large

*R*allows deriving a formula for the fringing field contribution to the head loss. However, the above expressions for

*U*

_{1}(

*X*) are very complex and justify an approximate treatment of the integral in order to get an analytic prediction of this important quantity for potential applications. This approximate treatment is detailed in Appendix 1. The idea consists in replacing the term ${U}_{1}^{5}-1$ by the approximation 5(

*U*

_{1}- 1) since the maximum velocity

*U*

_{1}(

*X*) remains close to unity as shown in figures 6 and 7. Then we use two asymptotic branches for

*U*

_{1}- 1, one for the inlet vicinity (

*β*≪ 1) and one for the downstream transition toward the Shercliff flow (1 -

*β*≪1), instead of the exact expression. Clearly, the maximum of

*U*

_{1}given by the intersection of these two branches is overestimated, but it is reasonable to expect that the full integral yields the correct dependence on the key parameters (

*ε*≪ 1,

*Ha*≫ 1,

*R*≥ 1). This yields the following asymptotic expressions:

### 5.3 The particular case of flow regimes with zero or weak inertia

This is a heat-like equation with the non-uniform diffusivity *β*/*β'*, the minus sign implying that diffusion is directed upstream. Remarkably, the Hartmann number has disappeared from equation (92), what implies that it does not affect the geometry, neither the typical length scale over which the M-shape profile is formed, nor the jet thickness. The Hartmann number influences only the local value of Ω(*β*) within this region through the amplitude coefficient *e*^{
Haβ
} . It is worth noticing that the two terms involving *∂*_{
β
}Ω in (67) illustrate the competition between the downstream advection of vorticity, predominant in the case analyzed in Section 5.2 and the upstream diffusion due to the exchange of current, predominant in the present analyzed case.

#### 5.3.1 The upstream M-shape profile formation region

*β*≪ 1, providing that the assumption ${\partial}_{X}^{2}\ll {\partial}_{Y}^{2}$ is still justified, the local equation (69) reads

*F*and

*F*' are

*O*(1), the coefficient of

*F*must be much smaller than

*Ha*, so that:

*α*(0) = 0 is

and shows that *α* deviates from its zero value in the far upstream region when *β*^{3} ≈ (*ε*^{2}/*Ha*)(1+3*πR*/2*ε*). Therefore, the length over which the M-shape profile develops is *L*_{
M-sh
} ≈ (*εHa*)^{1/3} in the inertialess limit. In the case of flows with weak inertia such that *R* » 1, we just recover the previous estimate obtained in Section 5.2: *L*_{up} ≈ (*ε*^{2}*Ha*/R)^{1/3}It is worth noticing that this estimate combines two effects: the upstream diffusion in a geometry independent of *H a*, and the vorticity distribution, which is itself *H a*-dependent.

*f*= 1/6 and neglecting the exchange of current since

*β*« 12

*ε*

^{2}

*R*/π in this region, the global equation (73) can be reduced to the simple form

*β*→0:

The components of the current density in this region immediately follow from equation (46) for the axial loop and (47) for the cross-sectional loop, while expression (64) yields the stream function *h*_{
a
}. Similarly, the pressure variation follows from equations (56*a*) and (56*b*), whereas equations (63) yield the head loss. We do not give details of these expressions since, as will be shown below, the predominant contribution does not come from the upstream region, but from the short transition toward the downstream Shercliff flow.

#### 5.3.2 The jet region

*U*=

*U*

_{0}+ (

*U*

_{1}-

*U*

_{0})

*f*(

*η*), with

*η =*(

*y*+1)/δ

_{j}(

*β*)where δ

_{ j }stands for the jet thickness. It implies that:

*U*

_{0}

*Y*-(

*U*

_{1}-

*U*

_{0})δ

_{ j }

*e*

^{-η}must take the value $-(1+\epsilon /\sqrt{Ha\beta})$ at the edge of the side layer (

*η*= 0). This yields one equation to determine the two unknowns

*U*

_{0}and

*U*

_{1}. The complementary equation comes from the global equation

which requires *U*_{0} = 0. Note that the numerical results suggest values significantly smaller than unity, but nevertheless positive, even for *Ha* = 10^{3} or larger. This small discrepancy is certainly a consequence of the quasi-two-dimensional approximation that neglects the exchange of current and derivatives ${\partial}_{X}^{2}$ in the Laplacians. Notice that these quantities, neglected in both the upstream and jet regions, are certainly signficant between these regions, as suggested by figure 3.

*U*

_{1}is completely determined by the conservation of the flow rate:

*U*

_{1}(

*β*) has a maximum at the abscissa where

*δ*

_{ SL }(

*β*) is minimum, namely,

*β*

_{ min }= (3/7)

^{1/2}≈ 0.655, or

*X*

_{ min }≈ 0.366 for

*ε*= 0.2. In turn,

*U*

_{1max}must vary as

*Ha*

^{1/2}/

*ε*. This behavior can be checked on figure 8, although the numerical value of

*U*

_{1max}, close to $0.3\sqrt{Ha}$, is definitely smaller than the predicted value, close to $2.38\sqrt{Ha}$. Equations (89) can also be used to get a relevant scaling for the typical length of the jet region, which must be such that

*β*≈

*ε*

^{2/3}

*Ha*

^{-1/3}. In agreement with the typical length obtained for the M-shape formation region, the typical length of the jet region must be

*L*

_{ jet }≈(

*εHa*)

^{1/3}. In the computed cases (

*ε*= 0.2,

*Ha*≈ 10

^{3}), this yields

*L*

_{ jet }≈ 6.

*Y*-component of the current density in the motionless core becomes:

*δ*

_{ j }is minimum and

*U*

_{1}maximum. This behavior is in fair agreement with the numerical data shown on figure 9(a). The maximum for the positive quasi-uniform

*J*

_{ Y }-component in the almost motionless core, as well as a strong concentration of the current lines in the transverse layer, where their sign is changed, are exhibited.

#### 5.3.3 The transverse layer

*β*= 1, and denote

*β*= 1 -

*ν*with

*ν*≪1. In this region, the assumption ${\partial}_{X}^{2}\ll {\partial}_{Y}^{2}$ and the heat-like equation (92) are not valid since they would yield a singularity at

*β*= 1. This means that the

*X*-derivatives must become predominant in order to smooth out the physically unacceptable singularity. Therefore, as noticed on the streamlines plotted on figure 3, the transition toward the downstream Schercliff flow takes place over a transverse layer located between the upstream low-velocity core and the downstream fully developed Shercliff flow. The side layer with the thickness $\epsilon /\sqrt{Ha}$ still exists along the side walls. Based on assumption ${\partial}_{Y}^{2}\ll {\partial}_{X}^{2}$, a relevant equation for the transverse layer can be derived from equation (49). Keeping only the predominant terms, this equation reads

*β'*

^{2}/

*β*= (

*πν*/

*ε*)

^{2}and introduce the transformed abscissa

*ζ*= 2(

*εv/π*)

^{1/2}to avoid the

*X*-dependent coefficient, equation (105) can be transformed in the elementary equation

This is a well-known equation in MHD duct flows in the presence of a uniform magnetic field although, usually, a coefficient *Ha*^{2} appears instead of *Ha*. This equation characterizes all kinds of free shear layers parallel to the applied magnetic field, often named inertialess Ludford layers, located between cores having different properties. The main novelty introduced by (106) is the thickness of the transverse layer, of the order *Ha* ^{-1/4} in ξ-units, instead of *Ha* ^{-1/2} in the case of Ludford layers. This layer has a clear similarity with the free shear layer discovered by Alboussière [31] along a singular characteristic surface "which splits" in ducts having a non-uniform width in the magnetic field direction. We now find that this layer appears along one of the characteristic surfaces whose shape is not at all singular, since all surfaces have the same shape (plane cross-sections of the rectangular duct). Indeed, it appears that this layer exists as soon as the slowly *X*-dependent solution exhibits a singularity. In our case, the origin of the singularity lies in the solution for the M-shape profile region, when *β* →1 and *β'* →0. Hence, the transverse layer provides the smooth transition between the M-shape flow region and the downstream fully developed Shercliff flow.

Several solutions of equations similar to (106) for flows in a uniform magnetic field have been found, for instance by Moffatt [40] (see also Moreau [41], pp. 138-150; Müller & Bühler [42], pp. 85-89). They have a remarkable structure determined by the superposition of two layers, one issuing from *Y* = -1 with a zero initial thickness, propagating in the positive *Y*-direction and diffusing in the ± *X*-directions, and the other issuing from *Y* = +1 with a zero thickness, propagating in the negative *Y*-direction and diffusing also in the ± *X*-directions. Here the term "propagating" reminds that these layers can be interpreted as diffusive Alfvén waves propagating in the ± **B** directions. The analogy between the previously studied free shear layers and the present transverse layer that appears in rectangular ducts at the entry of a strong magnetic field seems to stop at this point, since the analogous propagation is not any more in the direction of the magnetic field. This transverse layer seems to be the only possibility to feed the downstream fully developed flow, complying with the requirement that the characteristic surfaces cannot be crossed by the streamlines.

*Y*± 1 in analogy with the Ludford layers, which are parabolic functions of

*Z*± 1. Typical streamlines, drawn in the bottom half duct, show how this transverse layer feeds the downstream core flow from the jets. Figure 10(b) shows the velocity profile

*V*(

*X*,

*Y*= 1/2) at mid-height of the upper half duct computed for a very high value of the Hartmann number.

*Ha*and

*ε*still deserves some comments. Since equation (106) requires ${\partial}_{\zeta}^{4}\approx Ha$, or

*δ*

_{ ζ }=

*Ha*

^{-1/4}, its thickness in

*β*- and

*X*-units must be

This prediction of a *Ha*-independent thickness in *X*-units is *a priori* surprising. However, it has been confirmed through additional numerical results (not shown in this paper), obtained when the parameters used in figure 3 are modified within the range of small *R*-values to maintain inertia negligible. This actually demonstrates that the properties of this flow in a fringing field are not generic enough to be analyzed without making the choice of a specific *β*(*X*) distribution. For instance, would we use the distribution *β'* = (*πν*/*ε*)^{
n
}, the above derivation would yield *δ*_{
X
} ≈ (*Ha*^{n-1}*ε*^{2n})^{1/(4-2n)}. This clearly illustrates that the choice *n* = 1, corresponding to the magnetic field distribution (8), is the one which makes the transverse layer thickness independent of *Ha*.

It can be verified that the global current ${I}_{Y}^{TL}={\displaystyle {\int}_{-\infty}^{+\infty}{J}_{Y}^{TL}dX}$ is *Y*-independent and takes the remarkably simple value ${I}_{Y}^{TL}=-H{a}^{-1/2}$. As shown in figure 9, this explains that the axial current issuing from the jets, completes its trajectory through the transverse layer.

*Z*= 0) of the transverse layer are

*Ha*

^{1/2}/

*ε*and are, consequently, much higher than pressure differences in both the M-shape and jet regions. This explains the quite pronounced pressure drop across this layer exhibited on figure 11(a), which is a clear consequence of the scaling of the current density component

*J*

_{ Y }. Thus, it appears that the predominant contribution to the total pressure drop in the inertialess flow regime is originated in this transverse layer. This is important in the context of the potential applications to fusion blankets. At high

*Ha*numbers, the difference between the dimensional pressure outside

*p*

_{ out }and inside

*p*

_{ in }the magnet, obeys the scaling law:

*b*), in regimes with strong inertia (

*N*≪1) the pressure drop also depends linearly on the mean velocity:

This *a priori* surprising result for the inertial regime is a consequence of the required balance between inertia (proportional to ${\mathcal{U}}_{0}^{2}$) and Lorentz force (proportional to ${\mathcal{U}}_{0}$).

## 6 Discussion and conclusion

In this paper, the problem of the flow in a non-uniform magnetic field is revisited by paying special attention to accurate representation of the applied fringing magnetic field, electric current exchange between the core and the Hartmann layers, and the role of inertia. The starting point is the derivation of a consistent *div*- and *curl*-free magnetic field which models the inhomogeneities of the actual field at the entrance of a large yoke. A second important point is the explicit consideration of the Hartmann layers and the electric current exchange mechanism between them and the core flow. One of the key consequences of this mechanism is that the relevant scaling for the ratio of inertial and Lorentz forces in the case of insulating walls is *R* =*ε*^{2}*Re*/*Ha* instead of *N*^{-1} = *ε*^{2}*Re*/*Ha*^{2}. Another important consequence is that the three-dimensional electric current distribution is formed by a pair of two-dimensional current loops, namely, one in the axial (*x*, *y*), and another in the cross-sectional (*y*, *z*) planes. In the observed three-dimensional electric current distribution, in contrast with previous theories, a part of the induced current can cross the Kulikovskii characteristic surfaces within the core flow, in spite of the constraint imposed by the Alfvén theorem. As a result, the core flow can exhibit a non-zero vorticity generated by the rotational Lorentz force, which can be balanced by inertia, local Hartmann damping, and friction in the side layers.

The quantification of inertia in the reference flows is another novelty. After adapting the asymptotic theory based on the quasi-parallel flow assumption (${\partial}_{X}^{2}\ll {\partial}_{Y}^{2}$) and Oseen's approximation to flows with strong inertia we were able to construct scaling laws for the length over which the M-shape profile develops, which is found to be (*ε*^{2}*Ha*/*Re*)^{1/3}, and for the velocity maximum in the form *U*_{1} - *U*_{0} ≈ *Ha*/*εR*. These scaling laws are in qualitative agreement with the numerical predictions. For comparison, in the inertialess regime this length scale was found as (*εHa*)^{1/3} and the maximum velocity of the M-shape profile as *U*_{1} - *U*_{0} ≈ 1. Besides, the predicted downstream advection of the M-shape profile over a distance of order of *R* into the region where the magnetic field is uniform, also seems to be a distinguished feature of inertial flows. We have also observed that the electric current circuits in the fringing field region become less and less three-dimensional, with more currents closing in the axial planes, as inertia increases. Hence, for the inertial regime the exchange of current between the core and the Hartmann layers seems to be much less important than for inertialess flows.

Another remarkable property of the regimes with strong inertia is that the computational results remain quite stable, even at large Reynolds numbers (*Re* ≈ 3 ×10^{5}). The stability of the computations is in fact in agreement with recent observations of turbulence suppression by Andreev *et al.* [10, 11] in the experiments with a short magnet. This striking feature can be considered as one of the most noticeable properties of the flows in fringing fields. While in a strong uniform field quasi-two-dimensional flow disturbances are promoted and damped slowly for a typical Hartmann time ${\tau}_{H}=Ha\rho /\sigma {B}_{0}^{2}$ (see [33, 43]), in a region of non-uniform field, disturbances seem to be much more rapidly damped. As shown elsewhere [44], the main damping mechanism in the fringing field zone is due to electric currents crossing the Kulikovski's characteristic surfaces which create extra pressure losses in the mean flow and extra Joule dissipation affecting turbulent pulsations. These currents are induced by any vortical disturbance with a significant velocity component in the *x*-direction. In the non-uniform field zone where *β'* and *β* are of order unity, the additional damping is *Ha* times larger than the Hartmann damping. This has major consequences. In the first place, the extra damping overcomes the Hartmann braking effect, as well as both the non-linear interactions between vortices and the mean flow, and the instability mechanisms based on inflectional velocity profiles. Secondly, the quasi-two-dimensional vortices typical of MHD turbulence, which travel across the characteristic surfaces, are more rapidly damped than vortices elongated in the *y*-direction, parallel to these surfaces. This explains that, in spite of its simplicity, the linear model presented in [44] provides predictions in fair agreement with measurements reported in [10, 11]. In short, the characteristic surfaces progressively laminarize the turbulent upstream disturbances and stabilize the flow more rapidly and more efficiently than any other mechanism.

*β*→1. The most important findings presented in Sections 4 and 5 are summarized in the form of a flow sketch in figure 12. This sketch shows the most important lengths and thicknesses along with scaling laws derived in the paper for both inertialess and inertial flows.

Let us emphasize that the choice of a particular *div*- and *curl*-free two-dimensional magnetic field is particularly important, instead of assumptions on the length scale over which this field varies. Indeed, many crucial results, such as the thickness of the transverse layer and the final head losses, would be different with another magnetic field distribution. There remains some questions where the asymptotic analysis does not provide any clear answer, such as the dependence of the maximum velocity *U*_{1max}reached in the jets on the parameters *Ha* and *ε*. Additional important points, such as the flow properties when *ε* ≥ 1, or when *C* ≠ 0 would deserve to be re-examined. The general guideline of this investigation might easily be adapted to such additions, by adjusting the general formulation to new assumptions.

Finally, the results presented in this paper for flows with a strong or moderate inertia agree rather well with Tananaev's experiments performed in a similar geometry. The orders of magnitude of typical velocities and pressures agree fairly well. It is however difficult to quantitatively benchmark these data, because the exact values of the parameters (aspect ratio, definitions of the *Re* and *Ha* numbers) are not given with enough care in Tananaev's book [9].

## Appendix 1. Estimate for the head loss in inertial regimes (*R*≫ 1)

*U*

_{1}≈ 1 (see figures 5 and 7), the expression ${U}_{1}^{5}-1=({U}_{1}-1)({U}_{1}^{4}+{U}_{1}^{3}+{U}_{1}^{2}+{U}_{1}+1)$ can be approximated by 5(

*U*

_{1}-1) in the integral of the equation (63

*b*). Instead of the full expression (84) for

*U*

_{1}, let us use the two asymptotic branches, easily derived from (85), (87), (88) when

*ε*≪

*R*:

*N*≫ 1 (moderate inertia). The values of

*β*and

*X*where the two asymptotic branches intersect, denoted as

*β*

_{ * }or

*X**, can be estimated as follows:

*β*

_{*}and

*X*

_{*}yields: $5{(8/9)}^{1/3}{({\epsilon}^{2}N)}^{1/3}+20R{e}^{({\epsilon}^{2}Ha/9{\pi}^{2}{R}^{4}}{)}^{1/3}$. In the exponent of the exponential, the ratio

*Ha*/

*R*

^{4}=

*N*/

*R*

^{3}can be assumed quite small, so that the exponential can be taken as unity. Similarly, the factor (8/9)

^{1/3}can be taken as unity. This yields the following estimate for the head loss:

*N*≫ 1, the number

*R*=

*Ha*/

*ε*

^{2}

*Re*is small enough to justify neglecting the second term at the right-hand-side of (120). Then the estimate for the head loss becomes:

Consider now the opposite case where inertia is strong: *N* ≪1. In the expression of *U*_{1} - 1 for *β* ≫ 1, *N* ≪1, substitute to the varying variable *t* its limit when *β* = 1: *t*_{
max
} = π*Ha*/15*εR*.

*N*≪ 1, the value of

*X*

_{ * }is small enough to take the exponential as unity and to suggest

*β** = (2/5)

^{1/3}≈ 0.74 Then the relevant approximation for integral of (63) is:

It is noticeable that, if *N* ≫ 1 the total head loss is controlled by the initial asymptotic branch where *U*_{1} increases toward its maximum, whereas it is controlled by the final asymptotic branch where *U*_{1} decreases from its maximum if *N* ≪1.

## Appendix 2. Solving the differential equation (106)

*F*

_{+}(

*η*

_{+}) or

*F*

_{-}(

*η*

_{ - }):

*F*= constant and denote

*F'*=

*η f*(

*η*). Equation (125) becomes:

*y*(

*t*) =

*AI*

_{1/2}(

*t*) +

*BK*

_{1/2}(

*t*). It leads to the general expression for the basic function

*F*(

*η*):

*A*

_{+},

*B*

_{+},

*C*

_{+},

*D*

_{+}and

*A*

_{-},

*B*

_{-},

*C*

_{-},

*D*

_{-}. To get the upstream core flow motionless the general solution (129) must satisfy

*F*

_{+}(∞) =

*F*

_{-}(∞) = 0 as well as ${{F}^{\prime}}_{+}(\infty )={{F}^{\prime}}_{-}(\infty )=0$ This requires

*I*

_{1/2}(

*η*

^{2}/4). Reaching a uniform velocity in the downstream Shercliff flow requires

*F*

_{+}(0) =

*F*

_{-}(0) = 1/2, which leads to

## Declarations

### Acknowledgements

We are indebted to Mohamed Abdou, Thierry Alboussière, Neil Morley, Mingjiu Ni, Ramakanth Munipali and the group at Hypercomp Inc., for fruitful discussions. RM is indebted to the UCLA fusion group for supporting his visits during the accomplishment of this work. SS acknowledges support from the U.S. Department of Energy via grant DE-FG02-86ER52123-A040. SC thankfully acknowledges support from CONACYT under project 59977.

## Authors’ Affiliations

## References

- Hartmann J, Lazarus F: K Den Vidensk Selsk Mat-Fys Meddl. 1937, 15: 1-27.Google Scholar
- Shercliff JA: The theory of electromagnetic flow-measurement. 1962, Cambridge, UK: Cambridge University PressGoogle Scholar
- Shercliff JA: A textbook of Magnetohydrodynamics. 1965, Oxford, UK: Pergamon PressGoogle Scholar
- Shercliff JA: Magnetohydrodynamics (a 30 min. educational film). 1965, Educational Services Inc. for the National Committee on Fluid Mechanic Films, USAGoogle Scholar
- Branover H, Vasil'ev AS, Gelfgat YM: Magnetohydrodynamics. 1967, 3: 61-65.Google Scholar
- Slyusarev N, Shilova E, Shcherbinin EV: Magnetohydrodynamics. 1970, 6: 497-502.Google Scholar
- Lielausis O: Atomic Energy Review. 1975, 13: 527-581.Google Scholar
- Gelfgat YM, Lielausis O, Shcherbinin EV: Liquid metal under the action of electromagnetic forces. 1976, Riga, USSR: Zinatne PressGoogle Scholar
- Tananaev AV: Flows in MHD ducts. 1979, Moscow, USSR: AtomizdatGoogle Scholar
- Andreev O, Kolesnikov Y, Thess A: MHD-channel flow of liquid metal under an inhomogeneous magnetic field. Part 1: experiment. 15th Riga and 6th Pamir Int. Conf. on Fundamental and Applied MHD. 2005, 1: 185-190.Google Scholar
- Andreev O, Kolesnikov Y, Thess A: Phys Fluids. 2007, 19: 039902-10.1063/1.2718400.View ArticleADSGoogle Scholar
- Kulikovskii A: Isv Akad Nauk SSSR Mekh Zhidk Gaza. 1968, 3: 3-10.Google Scholar
- Kulikovskii A: Isv Akad Nauk SSSR Mekh Zhidk Gaza. 1973, 3: 144-150.Google Scholar
- Walker JS, Ludford GSS, Hunt JCR: J Fluid Mech. 1972, 56: 121-141. 10.1017/S0022112072002228.View ArticleADSMATHGoogle Scholar
- Walker JS, Ludford GSS: J Fluid Mech. 1972, 56: 481-496. 10.1017/S0022112072002460.View ArticleADSMATHGoogle Scholar
- Holroyd RJ, Walker JS: J Fluid Mech. 1978, 84: 471-495. 10.1017/S0022112078000282.View ArticleADSGoogle Scholar
- Walker JS: Three-dimensional laminar MHD flows in rectangular ducts with thin conducting walls and strong transverse nonuniform magnetic fields. Liquid Metal Flows and Magnetohydrodynamics, Proc. Third Int. Ber-Sheva Seminar in the MHD Flows and Turbulence Series. Edited by: Branover H, Lykoudis P, Yakhot A. 1983, AIAA, 3-19.Google Scholar
- Petrykowski JC, Walker JS: J Fluid Mech. 1984, 139: 309-324. 10.1017/S0022112084000379.View ArticleADSMATHGoogle Scholar
- Walker JS, Petrykowski JC: Approximate side layer solutions for a liquid metal flow in a rectangular duct with a strong nonuniform magnetic field. Proc Fourth Int Ber-Sheva Seminar on MHD Flows and Turbulence. Edited by: Branover H, Yakhot A. 1985, AIAA, 17-31.Google Scholar
- Reed CB, Picologlou BF, Hua TQ, Walker JS: IEEE. 1987, 2 (87CH2507-2): 1267-1270.Google Scholar
- Hua T, Walker JS: Numerical simulations of three-dimensional MHD flows in strong non-uniform transverse magnetic fields. Liquid Metal Magnetohydrodynamixcs, FMIA. Kluwer Acad. Pub. Edited by: Lielpeteris J, Moreau R. 1988, 13-19.Google Scholar
- Hua T, Walker JS: Int J Engng Sci. 1989, 27: 1079-1091. 10.1016/0020-7225(89)90086-4.View ArticleMATHGoogle Scholar
- Sellers CC, Walker JS: Int J Engng Sci. 1999, 37: 541-552. 10.1016/S0020-7225(98)00088-3.View ArticleMATHGoogle Scholar
- Lavrent'ev IV, Molokov SY, Sidorenkov SI, Shishko AR: Magnetohydrodynamics. 1990, 26: 328-338.MATHGoogle Scholar
- Molokov S, Shishko A: Fully developed magnetohydrodynamic flows in rectangular ducts with insulating walls. Tech Rep 5247, Kernforschungszentrum Karlsruhe. 1993Google Scholar
- Molokov S, Reed CB: Liquid metal flow in an insulated circular duct in a strong non-uniform magnetic field. Part 2. Inclination of the field gradient to the duct axis. Tech rep., ANL/TD/TM01-19, Argonne National Laboratory. 2001Google Scholar
- Molokov SY, Reed CB: Fusion Sci Tech. 2003, 32: 200-216.Google Scholar
- Holroyd RJ: J Fluid Mech. 1979, 93: 609-630. 10.1017/S0022112079001956.View ArticleADSGoogle Scholar
- Votyakov EV, Zienicke E, Thess A: MHD-channel flow of liquid metal under an inhomogeneous magnetic field. Part 2: Direct numerical simulation. 15th Riga and 6th Pamir Int. Conf. on Fundamental and Applied MHD. 2005, 1: 215-218.Google Scholar
- Votyakov EV, Zienicke E: Fluid Dyn Mat Process. 2007, 3: 97-113.Google Scholar
- Alboussière T: J Fluid Mech. 2004, 521: 125-154. 10.1017/S0022112004001740.View ArticleADSMathSciNetMATHGoogle Scholar
- Wong CP, Malang S, Sawan M, Dagher M, Smolentsev S: Fusion Engng Design. 2006, 81: 461-67. 10.1016/j.fusengdes.2005.05.012.View ArticleGoogle Scholar
- Sommeria J, Moreau R: J Fluid Mech. 1982, 118: 507-518. 10.1017/S0022112082001177.View ArticleADSMATHGoogle Scholar
- Pothérat A, Sommeria J, Moreau R: J Fluid Mech. 2000, 424: 75-100. 10.1017/S0022112000001944.View ArticleADSMATHGoogle Scholar
- Araseki H, Kotake S: J Comput Phys. 1994, 110: 301-309. 10.1006/jcph.1994.1027.View ArticleADSMathSciNetMATHGoogle Scholar
- Leboucher L: J Comput Phys. 1999, 150: 181-198. 10.1006/jcph.1998.6170.View ArticleADSMathSciNetMATHGoogle Scholar
- Smolentsev S, Cuevas S, Beltran A: J Comp Phys. 2010, 229: 1558-1572. 10.1016/j.jcp.2009.10.044.View ArticleADSMathSciNetMATHGoogle Scholar
- Cuevas S, Smolentsev S, Abdou MA: J Fluid Mech. 2006, 553: 227-252. 10.1017/S0022112006008810.View ArticleADSMathSciNetMATHGoogle Scholar
- Moreau R, Smolentsev S, Cuevas S: Magnetohydrodynamics. 2009, 45: 181-192.Google Scholar
- Moffatt HK: Electrically driven steady flows in magnetohydrodynamics. Proc 11th Int Cong Appl Mech. Edited by: Görtler H. 1994, Springer-Verlag, 946-953.Google Scholar
- Moreau R: Magnetohydrodynamics. 1990, Dordrecht, The Netherlands: Kluwer Acad PubView ArticleMATHGoogle Scholar
- Müller U, Bühler L: Magnetofluiddynamics in channels and containers. 2001, Berlin, Germany: SpringerView ArticleGoogle Scholar
- Knaepen B, Moreau R: Ann Rev Fluid Mechanics. 2008, 40: 25-45. 10.1146/annurev.fluid.39.050905.110231.View ArticleADSMathSciNetGoogle Scholar
- Moreau R, Smolentsev S, Cuevas S: Why and how turbulence can be suppressed in a non-uniform magnetivc field. Proc 6th Int. Conf. on Electromagnetic Processing of Materials. 2009, 37-40.Google Scholar

## Copyright

This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.