Skip to main content

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


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 = ε2Re/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 [58] 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 [1423] and more recently Molokov [2427], 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 ≈ 105) and moderate Hartmann numbers (Ha ≈ 102) 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.

The purpose of this paper is to provide a theoretical basis for practical applications, such as liquid-metal blankets of a fusion reactor, where the usual non-dimensional parameters can reach very high values (e.g. HaO (103) - O(105) and Re ≈ 105), 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 B0. 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 U0 stands for the typical velocity scale (e.g. mean bulk velocity). The magnetic yoke is semi-infinite in the x-direction with a length 2a m in the y-direction and a gap 2b 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 = U0a/ν, the Hartmann number H a = B 0 b σ / ρ ν , and the duct aspect ratio ε = b/a. Here, σ, ρ and ν are the electrical conductivity, density and kinematic viscosity of the fluid, respectively.

Figure 1
figure 1

Formulation of the problem. Geometry and main notations.

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 = Ha2/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

As in most former theories, we assume that the applied magnetic field 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 x = z x B + O ( ε m 3 B 0 ) , B y = z y B + O ( ε m 3 B 0 ) , B z = B ( x , y ) z 2 2 ( x 2 B + y 2 B ) + O ( ε m 4 B 0 ) . }

Here B0 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 x = σ ( x φ + v B ) , j y = σ ( y φ u B ) , j z = σ [ z φ + z ( u y B v x B ) ] . }

Note that the electric charge conservation principle, · j = 0, takes its classic form

Δ φ = B ( × u ) ,

only if the z-dependent component of the electromotive field u × B in equation (2c) 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 ρ ν / σ / 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 Ha1/2 in the case of a two-dimensional short magnet, and b m Ha1/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

Here we derive the expression of a two-dimensional fringing field when the magnetic poles are semi-infinite in the x-direction and infinite in the y-direction. Consider the conformal mapping between complex planes r and t defined by the relation

d r d t = K ( t 1 ) 1 / 2 t ,

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.

As a next step, let us construct the complex potential of the desired two-dimensional magnetic field. Starting from the complex potential of a two-dimensional source flow located at the origin of the t-plane,

w ( t ) = i B o b m π Log ( t ) ,

and, applying transformation (4) with the particular value K = - ib m /π (i2 = - 1), we easily get the x-dependence of the B z = B(x) component, which varies from B0 for x 0 to zero when x →∞. Its expression may be written under the parametric form

B ( x ) B 0 = 1 s , x b m = 2 π [ s + 1 2 Log ( s 1 s + 1 ) ] + C ,

where C is an arbitrary constant, which can be defined from the condition B(x = 0) = B0 /2. After changing the sign of x in order to get x increasing in the flow direction as B(x) from zero to B0, and denoting

X = x a , β = B ( x ) B 0 = 1 s ,

the basic expressions for the magnetic field and its derivative take the following form:

X = 2 ε π [ 2 1 β 1 2 Log ( 3 1 β 1 + β ) ] , β = d β d X = π 2 ε β 2 ( 1 β 2 ) .

Here, 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 (1a). 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.

Figure 2
figure 2

Magnetic field distribution.

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 (2b), 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.

Let us briefly recall the essence of the strong properties of the characteristic surfaces. Consider a small material cylinder of cross-section 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

d Φ d t = C j d s σ ,

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

As stated in the Introduction and will also be accurately checked later, we suggest that in the flow domain any scalar quantity q(x, y, z) can be expressed as q0(x, y) + (z2/2a2)q1(x, y) + This is particularly important for the electric potential

φ = φ 0 ( x , y ) + z 2 2 a 2 φ 1 ( x , y ) + ,

whose parabolic 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:

u = u 0 ( x , y ) + z 2 2 a 2 u 1 ( x , y ) + , v = v 0 ( x , y ) + z 2 2 a 2 v 1 ( x , y ) + , w = z 3 6 a 3 w 1 ( x , y ) +

It is noticeable that the highest order term of the w-component is a small O(ε3) term, so that the incompressibility condition writes:

u 0 = 0 , u 1 + w 1 a = 0.

Relations (11) also yield the following expressions for the vorticity components in the z-direction:

ω 0 = x v 0 y u 0 , ω 1 = x v 1 y u 1 .

In turn, to the leading orders in terms of the small parameter ε, from equations (2), (10) and (11), the components of the current density within the core flow are:

j 0 x = σ ( x φ 0 z 2 2 a 2 2 x φ 1 + v 0 B ) , j 0 y = σ ( y φ 0 z 2 2 a 2 y φ 1 u 0 B ) , j 0 z = σ ( z a φ 1 + z S 0 ) . }

Here, the scalar quantity

S 0 ( x , y ) = ( u 0 × B ) e z = u 0 y B v 0 x B ,

introduces the influence of the non-uniform magnetic field. Notice that the term involving S0 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 S0 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:

Δ φ 0 + φ 1 a 2 = B ω 0 .

Note the presence in (16) of z j z = φ 1 / a 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

Let us now examine the specific properties of the Hartmann layers. Because of symmetry, it is enough to treat the layer along the wall z = -b, whose thickness is δ ( x , y ) = ( 1 / B ) ρ ν / σ . 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:

φ = φ 0 + b 2 2 a 2 φ 1 + φ H ( x , y ) e ξ .

The exponential forms of u, v, j x , j y and j z are then:

u = u 0 ( 1 e ξ ) , v = v 0 ( 1 e ξ ) , j x = σ [ x φ 0 b 2 2 a 2 x φ 1 x φ H e ξ φ H x δ δ ξ e ξ + B v 0 ( 1 e ξ ) ] , j y = σ [ y φ 0 b 2 2 a 2 y φ 1 y φ H e ξ φ H y δ δ ξ e ξ B u 0 ( 1 e ξ ) ] , j z = σ [ b a 2 φ 1 + φ H δ e ξ b S 0 ( 1 e ξ ) ] . }

Notice that the values of the electric potential and current density at the edge of the layer coincide with equations (14), whereas, at the wall (z = -b, ξ = 0), they become:

j w x = σ [ x φ 0 b 2 2 a 2 x φ 1 x φ H ] , j w y = σ [ y φ 0 b 2 2 a 2 y φ 1 y φ H ] , j w z = σ [ b a 2 φ 1 + φ H δ ] , φ w = φ 0 + b 2 2 a 2 φ 1 + φ H . }

3.4 Within the Hartmann wall

Since the electric potential is continuous across the fluid-solid interface, the 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

t w ( x j w x + y j w y ) + j w z = 0 ,

where 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 :

σ w t w σ b Δ ( φ 0 + b 2 2 a 2 φ 1 + φ H ) = φ 1 a 2 + φ H b δ .

This relation, which introduces the usual conductance ratio C = σ w t w b as a fourth non-dimensional parameter, can still be significantly simplified when Ha 1, since C ΔφH φ H /. Then (21) becomes:

C Δ ( φ 0 + ε 2 2 φ 1 ) = φ 1 a 2 + φ H b δ .

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

The charge conservation equation has to be written also within the thin Hartmann layer, which, in the case of electrically conducting ducts, exchanges current with the Hartmann wall. First, let us express the components of the excess-current density:

j x j 0 x = σ [ x φ H e ξ φ H x δ δ ξ e ξ B v 0 e ξ ) ] , j y j 0 y = σ [ y φ H e ξ φ H y δ δ ξ e ξ + B u 0 e ξ ) ] . }

The integrals of the quantities on the left-hand-side of equations (23a) and (23b), denoted by g x and g y , must satisfy the global equation

x g x + y g y + j 0 z ( b ) j w z = 0 ,

which can be written as

x 2 ( φ H δ ) x ( φ H δ ) x δ δ + y 2 ( φ H δ ) y ( φ H δ ) y δ δ + φ H δ + ρ ν σ ω 0 + ( b δ ) S 0 = 0.

Clearly, in the limit of high Ha numbers, equation (25) reduces to

φ H δ = ρ ν σ ω 0 b S 0 .

Therefore, φ H can be substituted into (22) to get

C Δ ( φ 0 + ε 2 2 φ 1 ) = φ 1 a 2 ω 0 b ρ ν σ S 0 .

Equation (27) can be understood as a necessary relation between the two-dimensional and three-dimensional parts of the electric potential to feed the Hartmann layers and Hartmann walls, in addition to the condition of charge conservation (16). Eliminating φ1 now yields the equation relating φ0(x, y) to the velocity field:

( 1 + C ) Δ φ 0 C ε 2 a 2 2 [ Δ 2 φ 0 Δ ( B ω 0 ) ] = B ω 0 ρ ν σ ω 0 b S 0 .

Although equations (27) and (28) are complex because of their high differential order, they become rather simple when the duct is electrically insulated (C = 0):

Δ φ 0 = B ω 0 ρ ν σ ω 0 b S 0 ,
φ 1 a 2 = ρ ν σ ω 0 b + S 0 .

In this simple case, we just recover the same expression of j0zas when the magnetic field is uniform:

j 0 z = σ ρ ν ω 0 z b .

3.6 Equation of vorticity

Once equations for the two-dimensional and three-dimensional parts of the electric potential have been obtained, our next step is to derive the two-dimensional form of the equation of motion in the core flow domain. In order to eliminate the pressure and the irrotational part of the electromagnetic force, we limit ourselves to the vorticity budget, where the key ingredient is the z-component of × (j × B), which reads (B · )j0z- (j0 · ) B since both B and j0 are divergence-free. Using equations (1) for the magnetic field, the vorticity budget reads:

ρ ( u ) ω = ( B ) j 0 z + B z j 0 z ( j 0 ) B + ρ ν Δ ω .

For the reference case 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 j0, equation (31) for j0z, 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:

σ ρ ν B ω 0 b σ ρ ν z 2 b ( B ω 0 ) .

The same is true for the third term, which reads:

σ [ φ 0 B + B S 0 + z 2 2 a 2 ( φ 1 B ) ] .

This definitely proves that, if ε 1, the core flow can be modeled with a purely two-dimensional velocity field, by keeping only z-independent terms in equation (32):

ρ ( u ) ω 0 = σ ( x φ 0 B v 0 ) x B + σ ( y φ 0 + B u 0 ) y B σ ρ ν B ω 0 b + ρ ν Δ ω 0 .

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 j0z≠ 0.

Equation (32) could also be used to derive a higher order equation for ω1(x, y) by collecting terms proportional to z2. 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 ( σ ρ ν B ω 0 / b ) . 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 z4 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 j0x(x, y, z = 0) = -σ(∂xφ0 - Bv0), 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 j0xand j0zare derived from equations (35) and (31), while j0ycan be easily obtained from the charge conservation equation.

3.7 Basic set of non-dimensional equations

Let us now introduce dimensionless variables using characteristic scales: the uniform magnetic field B0, the mean bulk velocity 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:

X , Y = x , y a , Z = z b , U , V = u 0 , v 0 U 0 , Ω = ω 0 a U 0 , Ψ = ψ 0 U 0 a , Φ = φ 0 B 0 U 0 a , β = B ( x , y ) B 0 , J X , J Y , J Z = j 0 x , j 0 y , j 0 z σ B 0 U 0 . }

In addition, S denotes the dimensionless expression of S0. 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:

Δ Φ = β Ω Ω H a S ,
ε 2 R e ( U · ) Ω = H a 2 [ β · Φ + β S ] H a β Ω + ε 2 Δ Ω ,
Δ Ψ = Ω .

It is worth noticing that the term Ω/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,

J X = ( X Φ + β X Ψ ) , J Y = ( Y Φ + β Y Ψ ) , J Z = ε H a Ω Z , S = β · Ψ . }

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)).

Here we also consider the case of insulating wall ducts (C = 0). The basic equations become:

Δ Φ = β Ω Ω H a β X Ψ ,
ε 2 R e ( U ) Ω = H a 2 β ( X Φ + β X Ψ ) H a β Ω + ε 2 Δ Ω .

In the most general case, considering inertia and viscous friction, the vorticity budget (42) imposes the following expression for the axial component of the current density:

J X = R Γ H a β β Ω H a β + ε 2 H a 2 β Δ Ω .

where R = ε2Re/Ha and Γ = (U · ) Ω. As discussed below, R results to be the suitable scale for inertia. Once equations (43) and (40c)) are substituted in the electric charge conservation equation, we get

Y J Y = 1 H a [ X ( R Γ β + β Ω β ε 2 H a β Δ Ω ) + Ω ] .

If we also assume that the vorticity is slowly X-dependent, so that it can be approximated as Ω ( X , Y ) = Y 2 Ψ = Y U and Δ Ω = Y 2 Ω , equation (44) can be integrated once. Using the boundary conditions U = J Y = 0 at the insulating side wall Y = -1, yields:

J Y = 1 H a { X [ R β 1 Y Γ d Y β U β + ε 2 H a β ( Y 2 U ( Y 2 U ) w ) ] U } .

The three components (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:

J a X = Y h a = 1 H a Y [ R β 1 Y Γ d Y + β β U ε 2 H a β ( Y 2 U ( Y 2 U ) w ) ] , J a Y = X h a = 1 H a X [ R β 1 Y Γ d Y + β β U ε 2 H a β ( Y 2 U ( Y 2 U ) w ) ] . }

The other family, which is 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:

J c Y = 1 ε Z h c = U H a , J c Z = Y h a = ε Z H a Y U . }

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.

Let us now eliminate the electric potential between the charge conservation equation (41) and the vorticity budget (42). The X-derivative of the charge conservation equation is:

Δ J X = β Y 2 Ψ + X Ω H a .

It is straightforward to substitute equation (43) for J x in (48) to get

Δ ( R Γ + β Ω ε 2 H a Δ Ω β ) + H a β Y 2 Ψ + X Ω = 0 ,

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, Δ(/β'), involves the dimensionless parameter R = ε2Re/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 = ε2Re/Ha2.

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 / H a 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 ≈ 102 [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 × B0, 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 full set of the equations governing the problem includes (48) along with the vorticity equation (42) modified by the inclusion of the time derivative in order to allow the purely numerical time evolution from the initial solution to the final steady or unsteady solution, the condition of incompressibility (39), and the electric charge conservation equation (41), presented in the form

X J X + Y J Y Ω H a = 0.

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 (R1)

As stated previously, 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].

Figure 3
figure 3

Example of numerical results in a regime with negligible inertia. (Ha = 1000, Re = 500, ε = 0.2, and R = 0.02). a) Velocity vectors. b) Vorticity. c) Streamlines. d) Current stream tracers.

For those two Hartmann numbers, all distributions are qualitatively similar. When the liquid proceeds through the fringing field region, the velocity at the axis U0(X) decreases, the near-wall jets form, and the maximum velocity U1(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 U0(X) does not reach a zero value (see figure 4) and the maximum value of U1(X) increases as U1max≈ 0.3Ha1/2 instead of 2.38Ha1/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 U1/δ SL , where δ SL stands for the thickness of the side layer, and yields an estimate of U1. The other peak, of the opposite sign and much smaller value, can be estimated as U1/δ 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.

Figure 4
figure 4

Axial velocity U 0 ( X ). Typical distributions of the axial velocity U0(X) for different combinations of the parameters Ha and Re, illustrating the influence of inertia (in all cases ε = 0.2).

Figure 5
figure 5

Axial velocity U 1 ( X ). Typical distributions of the axial velocity U1(X) for different combinations of the parameters Ha and Re, illustrating the influence of inertia (in all cases ε = 0.2).

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

For these regimes, that include flows with significant or strong inertia, the computational domain is extended to 60 and the number of grid points in the axial direction in all computations was increased to 1024 for better spatial resolution. Also, the time step is decreased to assure that no numerical instability occurs in the course of computations. These changes, along with a slower computational convergence, result in a computational time many times longer compared to that in the inertialess regime (typically 5-7 days of computations using a PC). To reach high values of 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.

Figure 6
figure 6

Example of numerical results in a regime with significant inertia. (Ha = 200, Re = 100, 000, ε = 0.2, and R = 20). a) Velocity vectors. b) Vorticity. c) Streamlines. d) Current stream tracers.

Figure 7
figure 7

Example of numerical results in a regime with strong inertia. R > Ha (Ha = 100, Re = 300, 000, ε = 0.2, and R = 120). a) Velocity vectors. b) Vorticity. c) Streamlines. d) Current stream tracers.

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 X 2 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, X 2 becomes predominant upon Y 2 (see figure 3).

According to the approximation made in Section 3.6, valid for slotted ducts, we neglect the terms proportional to ε2 Z2 so that the components of the Lorentz force read

F X = H a β X h a β U , F Y = H a β Y h a , F Z = H a β Z X h a + β Z U . }

Introducing the non-dimensional pressure defined as

P ( X , Y , Z ) = p ρ U 0 2 R ,

the X- and Y -projections of the motion equation read:

X P = R ( U ) U + ε 2 H a Y 2 U H a β X h a β U , Y P = H a β Y h a . }

Similarly, equation (48) for J X can also be simplified as

Y 2 ( J X β Ψ ) = 1 H a X Ω .

In quasi-parallel flows, the right-hand-side of (54) is negligible when Ha 1, yielding the simple relations

J X = Y h a = β Ψ .

Remarkably, neither the pressure, nor the current density stream function 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 P0(X), and its value within the side layer by P1(X), equations (53) allow deriving the pressure from the following equations

P = P 1 ( X ) H a β h a , P 1 = R ( U ) U β U + ε 2 H a Y 2 U + H a β h a . }

At the side wall, where U = h a = 0, equation (56b) yields the X-dependence of the pressure at this location:

P 1 = ε 2 H a ( Y 2 U ) w .

Subtracting equation (56b) from its form at the edge of the layer, where the velocity reaches its maximum value U1, the relevant equation for the side layer is

R [ X ( U 1 2 U 2 ) Y ( U V ) ] + β ( U 1 U ) + ε 2 H a Y 2 U = 0.

For the inertialess regime (R = 0), in the vicinity of the side wall Y = -1, an elementary solution for the side layer exists, namely,

U = U 1 ( 1 e η ) , η = Y + 1 δ S L , δ S L = ε ( H a β ) 1 / 2 .

However, at inertial regimes, the accelerated side layer does not obey the self-similar solution (59). Nevertheless, let us assume that equation (59) still represents a convenient approximation, except for the layer thickness δ 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

3 R β 2 d β ( U 1 2 δ S L ) + β U 1 δ S L ε 2 U 1 H a δ S L = 0.

When R is large enough, the Hartmann friction becomes negligible and equation (60) implies the invariance of the missing momentum: U 1 2 δ S L = 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 U1(β):

δ S L = ε H a 1 / 2 U 1 2 .

From these properties of the side layer, it follows that equation (57) for the pressure variation along the duct axis reads:

X P 1 = β U 1 , if R 1 , X P 1 = U 1 5 , if R 1. }

Therefore, the contribution of the fringing magnetic field to the head loss, which incidentally is one of the most important quantities for applications, can be derived from

( Δ P ) 1 = 0 1 ( β U 1 1 ) d β β , if R 1 , ( Δ P ) 1 = 0 1 ( U 1 5 1 ) d β β , if R 1. }

In turn, the stream function of the current density obeys (55) and can be expressed as

h a = β Y 1 Ψ d Y ,

in order to be zero at the insulating side walls. Then, the pressure difference between the axis or the middle of the Hartmann walls, and the middle of the side walls, which is of interest because it has been measured by Tananaev [9], reads:

( Δ P ) 2 = P 0 ( X ) P 1 ( X ) = H a β h a ( 0 ) = H a β β 0 1 Ψ d Y .

It is noticeable that (Δ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:

( Δ P ) 2 = π H a 4 ε β 3 ( 1 β 2 ) .

Its maximum (ΔP)2 = (πHa/10ε)(3/5)3/2 ≈ 0.146(Ha/ε) is reached for β = 3 / 5 0.7746 .

5.2 Inertial regimes

Let us now pay attention to regimes with a strong inertia, where the core velocity remains close enough to unity to justify the use the Oseen's approximation (U · )Ω = X Ω = β'∂ β Ω. Hence (49), which is the basic equation for the vorticity, reads

Y 2 ( R β β Ω + β Ω ε 2 H a Y 2 Ω ) H a β 2 Ω + β 2 β Ω = 0.

The terms of this linear equation express all major effects: from left to right, we have advection of vorticity ( R β Y 2 β Ω ) , Hartmann damping ( β Y 2 Ω ) , damping in the side layers ( ε 2 Y 4 Ω / H a ) , torque (Haβ'2Ω) and exchange of current (β'2βΩ).

Let us express the velocity field in the quasi-parallel core flow in terms of the eigenfunctions sinh α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:

Ψ = Y + F α 3 ( sinh α Y Y sinh α ) , U = 1 + F α 3 ( α cosh α Y sinh α ) , Ω = F α sinh α Y . }

The expression for the vorticity can be simplified as Ω = -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:

R β ( α 2 F ) + β α 2 F β 2 ( H a F F ) = 0.

Note that an additional term of the order of U1δ 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.

The second equation necessary to determine α(β) and F(β) can be derived from a global form of equation (67) integrated between the side wall Y = -1 and the axis Y = 0:

R β β ( Y 2 U ) 0 + β ( Y 2 U ) 0 ε 2 H a ( Y 4 U ) 0 H a β 2 U 0 + β 2 U 0 ' R β β ( Y 2 U ) w β ( Y 2 U ) w + ε 2 H a ( Y 4 U ) w = 0.

The third term of equation (70) is negligible since ε2/Ha 1 and ( Y 4 U ) 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 U0 and maximum velocity in the jet U1, can be written as

U 0 = 1 f F , U 1 = 1 + g F ,

where we denote

f = sinh α α α 3 , g = α cosh α sinh α α 3 .

Equation (70) finally reads:

R β F + β F H a β 2 ( 1 F f ) β 2 ( F f ) = 0.

As soon 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

U U 0 U 1 U 0 = cosh α Y 1 cosh α 1 ,

indicating that values of α 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

U 1 U 0 = F ( g f ) .

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).

To proceed further with a solution of equations (69) and (73), let us substitute the following approximations valid when α ≤ 1, to the general expressions of f(α) and g(α) in equations (72),

f = 1 6 + α 2 120 , g = 1 3 + α 2 30 ,

which keep linear the equations for F (α) and G(α) = α2F:

R β G + β G = H a β 2 F , R β F + β F = H a β 2 ( 1 F 6 G 120 ) . }

The form of the left-hand-side of equations (77a) and (77b) and the expression (8) for β' = (π/2) β2(1 2), suggest the transformation

F = C ( β ) ( 1 β 2 β 2 ) ε / π R , G = D ( β ) ( 1 β 2 β 2 ) ε / π R .

Hence, the differential equations (77a) and (77b) become

D = H a β R C , C + H a β R ( C 6 + D 120 ) = H a β R ( 1 β 2 β 2 ) ε / π R . }

Eliminating C and introducing the rapidly varying variable

t = π H a 2 ε R ( β 3 3 β 5 5 ) ,

such that t' = Haβ'/R, yields the second order equation

d t 2 D + d t D 6 + D 120 = ( 1 β 2 β 2 ) ε / π R .

Since the previous approximations have already put the emphasis on regimes with significant inertia, such that R 1, the right-hand-side of (81) can be assumed to be unity everywhere, except in the inlet vicinity where β → 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:

C = 12 5 sin ( t 12 5 ) e t / 12 , D = 120 { 1 e t / 12 [ cos ( t 12 5 ) + 5 sin ( t 12 5 ) ] } . }

The functions F (β) and G(β) are therefore given by:

F = 12 5 e t / 12 sin ( t 12 5 ) ( 1 β 2 β 2 ) ε / π R , α 2 F = 120 { 1 e t / 12 [ cos ( t 12 5 ) + 5 sin ( t 12 5 ) ] } ( 1 β 2 β 2 ) ε / π R . }

It is now elementary to derive expressions for most characteristic quantities, such as the velocities on the axis, U0 = 1 - f F , and at the edge of the side layer, U1 = 1 + gF :

U 0 = 1 ( 1 β 2 β 2 ) ε / π R { 1 e t / 12 [ cos ( t 12 5 ) 5 sin ( x 12 5 ) ] } , U 1 = 1 + 4 ( 1 β 2 β 2 ) ε / π R { 1 e t / 12 cos ( t 12 5 ) } . }

From the asymptotic behavior of these typical velocities near the inlet, it can be assured that no singularity appears when β → 0:

U 0 = 1 π H a 36 ε R β 3 2 ε π R , U 1 = 1 + π H a 18 ε R β 3 2 ε π R .

In particular, in most experimental conditions, these velocities deviate from their asymptotic limit as 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 ≈ (ε2Ha/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.

The analogous derivation in the outlet vicinity, which determines how the M-shape profile evolves toward the downstream Shercliff profile, deserves some care. It must be noticed that, according to (8), near the limit β → 1, we get 1 - β = (2/3)eX/ε, so that expressions (84a) and (84b) become:

U 0 = 1 ( 4 3 e π X / ε ) ε / π R { 1 e t / 12 [ cos ( t 12 5 ) 5 sin ( t 12 5 ) ] } , U 1 = 1 + 4 ( 4 3 e π X / ε ) ε / π R [ 1 e t / 12 cos ( t 12 5 ) ] . }

If N = Ha/R is strong enough, or inertia remains moderate to allow the exponential e-t/12 to reach its zero limit, then (86a) and (86b) reduce to

U 0 = 1 ( 4 3 ) ε / π R e X / R , U 1 = 1 + 4 ( 4 3 ) ε / π R e X / R .

These equations show that the downstream length scale between the abscissa where the M-shape profile is maximum and the Shercliff flow, scales as 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/12can be approximated by 1 - t /12, so that equations (86a) and (86b) reduce to

U 0 = 1 ( 4 3 ) ε / π R e X / R ( 1 t 6 ) , U 1 = 1 + ( 4 3 ) ε / π R e X / R t 3 .

The correction for the U0 value is small since t/6 1, but the correction for U1 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/Rand scales as L down R. This downstream transport of the M-shape over a long distance is also clear on figures 6 and 7.

In principle, equation (63b) for flows at large R allows deriving a formula for the fringing field contribution to the head loss. However, the above expressions for U1(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(U1 - 1) since the maximum velocity U1(X) remains close to unity as shown in figures 6 and 7. Then we use two asymptotic branches for U1 - 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 U1 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:

( Δ P ) 1 = 5 ( ε 2 N ) 1 / 3 , if ( ε 2 N ) 1 / 3 R , ( Δ P ) 1 = π H a / 2 ε , if N 1 . }

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

Let us now consider the asymptotic condition R 1 (see figure 3). As mentioned in Section 4.1, appart from the far upstream region where a Shercliff-like flow develops and the far downstream region where a fully established Shercliff flow exists, we distinguish three intermediate regions strongly affected by the fringing field: an upstream region where the M-shape velocity profile starts to form, a jet region where the flow rate is concentrated in two jets on both sides of an almost motionless core, and a transition region toward the downstream fully established flow. The main novelty in the vorticity equation (49), or (67), in comparison with previous investigations, (see for instance [16, 24, 31]), is the presence of the last term, which represents the exchange of current. To appreciate its potential influence, let us use the transformation

Ω = ω ( X , Y ) e H a β ( X ) ,

which allows combining the last two terms of (49) into only one of order unity, so that the inertialess form now reads

ε 2 H a Y 4 ω β Y 2 ω = β X ω .

The implications of that equation in the core flow appear particularly strong, since (91) reduces to the simple form

X ω = β β Y 2 ω .

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 eHaβ . 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

In this region where β 1, providing that the assumption X 2 Y 2 is still justified, the local equation (69) reads

F ( R β α 2 + β 2 ) + F ( R β ( α 2 ) + β α 2 H a β 2 ) = 0.

Since F and F' are O(1), the coefficient of F must be much smaller than Ha, so that:

π R 2 ε β ( α 2 ) + α 2 = π 2 H a β 3 4 ε 2 .

The solution of equation (94) such that α (0) = 0 is

α 2 = π 2 H a β 3 4 ε 2 ( 1 + 3 π R 2 ε ) ,

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: Lup ≈ (ε2Ha/R)1/3It 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.

Besides, using the asymptotic value f = 1/6 and neglecting the exchange of current since β « 12ε2R/π in this region, the global equation (73) can be reduced to the simple form

π R 2 ε β F + ( 1 + π 2 H a β 3 24 ε 2 ) F = π 2 H a β 3 4 ε 2 .

In the inertialess limit the solution is straightforward:

F = π 2 H a β 3 4 ε 2 1 + π 2 H a β 3 24 ε 2 , U 0 = 1 1 + π 2 H a β 3 24 ε 2 , U 1 = 1 + π 2 H a β 3 8 ε 2 1 + π 2 H a β 3 24 ε 2 . }

For non-zero but weak inertia, it is also simple to derive the small influence of the vorticity advection in the vicinity of the limit β →0:

F = π 2 H a β 3 4 ε 2 1 + 3 π R 2 ε , U 0 = 1 π 2 H a β 3 24 ε 2 1 + 3 π R 2 ε , U 1 = 1 + π 2 H a β 3 12 ε 2 1 + 3 π R 2 ε . }

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 (56a) and (56b), 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

Let us now depart from the initial upstream region and consider abscissas where, as suggested by figure 3, quasi-parallel jets exist on both sides of the almost motionless core, being separated from the walls by the side layers. The relevant form of equation (67) in the jet near the wall is:

β Y 2 Ω β 2 ( H a Ω β Ω ) = 0.

If the two first terms of equation (99) are predominant upon the exchange of current, a self-similar solution can be found such that U = U0 + (U1-U0) f(η), with η = (y +1)/δj(β)where δ j stands for the jet thickness. It implies that:

δ j 2 = 4 ε 2 π 2 H a 1 β 3 ( 1 β 2 ) , f = e η .

The conservation of the flow rate requires that the stream function Ψ =U0Y-(U1-U0 j e must take the value ( 1 + ε / H a β ) at the edge of the side layer (η = 0). This yields one equation to determine the two unknowns U0 and U1. The complementary equation comes from the global equation

β ( Y Ω ) 0 + β 2 ( H a U 0 d β U 0 ) = 0 ,

which requires U0 = 0. Note that the numerical results suggest values significantly smaller than unity, but nevertheless positive, even for Ha = 103 or larger. This small discrepancy is certainly a consequence of the quasi-two-dimensional approximation that neglects the exchange of current and derivatives 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.

Hence, the maximum velocity U1 is completely determined by the conservation of the flow rate:

U 1 = 1 + δ S L δ j = β ( H a β ) 1 / 2 ( 1 + δ S L ) π 2 ε [ ( H a β 3 ) 1 / 2 ( 1 β 2 ) ] .

This expression shows that U1(β) 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, U1maxmust vary as Ha1/2/ε. This behavior can be checked on figure 8, although the numerical value of U1max, close to 0.3 H a , is definitely smaller than the predicted value, close to 2.38 H a . 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 ≈ 103), this yields L jet ≈ 6.

Figure 8
figure 8

Maximum velocity in the jet. Numerical results for the maximum velocity in the jet as a function of Ha in a regime with negligible inertia (Re = 500, ε = 0.2) and their approximation.

It is straightforward to derive the expressions of the current density and pressure from the self-similar solution:

h a = β U 1 H a β ' ( 1 e η ) , P 0 = P m a x β 2 β ' 1 + δ S L δ j e η .

In turn, the Y -component of the current density in the motionless core becomes:

J Y c o r e = X h a c o r e = π 4 ε β 3 / 2 ( 1 β 2 ) .

In fact, the current density component has also a maximum at an abscissa slightly different from that where δ j is minimum and U1 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.

Figure 9
figure 9

Induced electric current. Induced electric current (J x , J y ) stream tracers in the symmetry plane (Z = 0). a) In a regime with negligible inertia (Ha = 1000, Re = 500, ε = 0.2, and R = 0.02). b) In a regime with significant inertia (Ha = 100, Re = 20, 000, ε = 0.2, and R = 8).

5.3.3 The transverse layer

Let us now consider the vicinity of β = 1, and denote β = 1 - ν with ν 1. In this region, the assumption X 2 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 ε / H a still exists along the side walls. Based on assumption Y 2 X 2 , a relevant equation for the transverse layer can be derived from equation (49). Keeping only the predominant terms, this equation reads

X 4 Ψ = H a β 2 β Y 2 Ψ .

If we take β'2/β = (πν/ε)2 and introduce the transformed abscissa ζ = 2(εv/π)1/2to avoid the X-dependent coefficient, equation (105) can be transformed in the elementary equation

ζ 4 Ψ = H a Y 2 Ψ .

This is a well-known equation in MHD duct flows in the presence of a uniform magnetic field although, usually, a coefficient Ha2 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.

Due to the analogy between equation (106) and the corresponding equation for a flow in a uniform magnetic field, a well-known technique can be applied for its solution, which can be epxressed as a combination of two self-similar solutions:

Ψ ( x , Y ) = ( Y + 1 ) F + ( η + ) + ( Y 1 ) F ( η ) , η + = H a 1 / 4 ζ Y + 1 , η = H a 1 / 4 ζ Y 1 F ( η ) = 1 2 ( 1 + η 2 2 ) ( 1 erf η 2 ) η e η 2 / 4 2 π . }

Specific details on the algebra are given in Appendix 2. The transverse layer, whose organization is sketched on figure 10(a), is the superposition of two symmetric layers. Their thicknesses are parabolic functions of 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.

Figure 10
figure 10

The transverse layer. a) Schematic view of the double layer with some streamlines in the bottom half and a typical profile of the V -component of the velocity in the upper half; b) derived V -profile computed for Ha = 105 and ε = 0.2 at Y = 0.5 and Z = 0.

To show the existence of this transverse region and determine its scaling laws are, certainly, the main results of this attempt to clarify the properties of the inertialess flow in the fringing field of a magnet. The way the thickness depends on the key parameters Ha and ε still deserves some comments. Since equation (106) requires ζ 4 H a , or δ ζ = Ha -1/4, its thickness in β- and X-units must be

δ β = ( δ ζ ) 2 / ε = H a 1 / 2 ε 1 , δ X = ε .

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 ≈ (Han-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.

The current density in this transverse layer can now be derived from the following relations

J X T L = β H a β X 2 Ψ , J Z T L = ε Z H a X 2 Ψ ,

together with the condition of charge conservation, which yields

Y J Y = β H a β ' X 3 Ψ .

As expected above, equation (110) determines the way the electric circuit is closing via the transverse layer, where the most important component of the current density is given by:

J Y T L = π 2 ε H a 1 / 2 [ ( l + Y ) η + G + + ( 1 Y ) η G ] , G ( η ) = n 2 [ 1 erf ( η 2 ) ] 1 π e η 2 / 4 . }

It can be verified that the global current I Y T L = + J Y T L d X is Y-independent and takes the remarkably simple value I Y T L = 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.

The relevant equations to derive the dimensionless pressure distribution P = R p / ρ U 0 2 , in the symmetry plane (Z = 0) of the transverse layer are

X P = H a β J Y T L , Y P = H a β J X T L .

From the specific properties of the transverse layer, we get:

Y P = ζ 2 Ψ = H a 1 / 2 ( F + F ) = 1 2 [ erf ( η + 2 ) erf ( η 2 ) ] .

The pressure profile finally reads

P = π H a 1 / 2 ε ( 1 + Y 2 K + + 1 Y 2 K ) , K ( η ) = erf ( η 2 ) η 2 2 [ 1 erf ( η 2 ) ] + η π e η 2 / 4 . }

We observe that pressure differences seem to scale as Ha1/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:

Figure 11
figure 11

Numerical results for the pressure distribution. a) In a regime with negligible inertia (Ha = 200, Re = 2000, ε = 0.2, and R = 0.4). b) In a regime with significant inertia (Ha = 200, Re = 10, 000, ε = 0.2, and R = 2).

p o u t p i n = π H a 1 / 2 ε R ρ U 0 2 = π H a 3 / 2 ρ ν U 0 a ( a b ) 3 ,

which yields the expected linear dependence with the mean velocity U 0 . For the sake of comparison, let us remind that according to equation (89b), in regimes with strong inertia (N 1) the pressure drop also depends linearly on the mean velocity:

p o u t p i n = π H a 2 ε R ρ U 0 2 = π 2 σ B 0 2 U 0 ( a b ) 2 .

This a priori surprising result for the inertial regime is a consequence of the required balance between inertia (proportional to U 0 2 ) and Lorentz force (proportional to 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 =ε2Re/Ha instead of N-1 = ε2Re/Ha2. 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 ( X 2 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 (ε2Ha/Re)1/3, and for the velocity maximum in the form U1 - U0Ha/ε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 U1 - U0 ≈ 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 ×105). 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 τ H = H a ρ / σ 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.

The important and subtle role of the electric current exchange in the inertialess regime, seems to be analyzed in detail for the first time in this paper. The inclusion of the exchange of current gives the vorticity equation (91) a parabolic character. This means that the flow in each cross-section is affected by the flow in sections where the magnetic field is stronger. Therefore, in the case of a flow entering into a magnet diffusion is directed upstream. In particular, in the region where the M-shape velocity profile develops, provided that the side friction is neglected, the basic vorticity equation reduces to a heat-like equation which, as the main important consequence, predicts the existence of the transverse layer when β →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.

Figure 12
figure 12

Organization of the flow. Organization of the flow in different regions with the most typical length scales (not to scale): a) inertialess regime, b) inertial regime.

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 U1maxreached 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)

Since U1 ≈ 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(U1 -1) in the integral of the equation (63b). Instead of the full expression (84) for U1, let us use the two asymptotic branches, easily derived from (85), (87), (88) when ε R:

* for β 1 : U 1 1 = π H a 18 ε R β 3 , *for β 1 : { U 1 1 = 4 e X / R , if N 1 , U 1 1 = t 3 e X / R , if N 1. }

Let us start with the case N 1 (moderate inertia). The values of β and X where the two asymptotic branches intersect, denoted as β * or X*, can be estimated as follows:

β * 3 = 72 ε π N , X * = ( ε 2 H a 9 π 2 R ) 1 / 3 .

Then the integral of (63) can be approximated as:

0 1 5 ( U 1 1 ) d β β ' = 0 β * 5 H a 9 R β d β + X * 2 0 e X / R d X = 5 N 18 β * 2 + 20 R e X * / R .

Substituting the values of β* and X* yields: 5 ( 8 / 9 ) 1 / 3 ( ε 2 N ) 1 / 3 + 20 R e ( ε 2 H a / 9 π 2 R 4 ) 1 / 3 . In the exponent of the exponential, the ratio Ha/R4 = N/R3 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:

( Δ P ) 1 5 ( ε 2 N ) 1 / 3 + 20 R , if N 1.

In most cases where N 1, the number R = Ha/ε2Re is small enough to justify neglecting the second term at the right-hand-side of (120). Then the estimate for the head loss becomes:

( Δ P ) 1 5 ( ε 2 N ) 1 / 3 , if ( ε 2 N ) 1 / 3 R .

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

The asymptotic branches intersect at β * 3 = ( 2 / 5 ) e X / R . Since 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:

0 1 5 ( U 1 1 ) d β β ' = 0 β * 5 H a 9 R β d β + X * π H a 9 ε e X / R d X = 5 N 18 β * 2 + π H a 2 ε .

The last term at the right-hand-side is predominant, yielding the following estimate for the head loss:

( Δ P ) 1 = π H a 2 ε , if N 1.

It is noticeable that, if N 1 the total head loss is controlled by the initial asymptotic branch where U1 increases toward its maximum, whereas it is controlled by the final asymptotic branch where U1 decreases from its maximum if N 1.

Appendix 2. Solving the differential equation (106)

Let us look for a self-similar solution of the form (107). In ξ-units, this forms combines two symmetric self-similar layers with thicknesses:

δ + = ( Y + 1 ) 1 / 2 H a 1 / 4 , δ = ( Y 1 ) 1 / 2 H a 1 / 4 .

Substituting (124) into (106) yields the following 4th order ordinary differential equation for either F+(η+) or F-(η - ):

4 F I V η 2 F + η F = 0.

We first disregard the obvious solution F = constant and denote F' = η f(η). Equation (125) becomes:

4 ( η f ) η 3 f = 0.

Using the transformation

f ( η ) = y ( t ) η , t = η 2 4 ,

reduces once more the differential order of equation (126) and provides the modified Bessel equation of order 2:

t 2 y + t y ( t 2 + 1 4 ) y = 0.

The general solution reads y(t) = AI1/2(t) + BK1/2(t). It leads to the general expression for the basic function F (η):

F ( η ) = ( η [ A I 1 / 2 ( η 2 4 ) + B K 1 / 2 ( η 2 4 ) ] d η η ) d η + C η 2 2 + D .

There are height constants to determine from the relevant boundary conditions: 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 + ( ) = F ( ) = 0 This requires

A + = A = 0 , B + = B = 1 π 2 , C + = C = 1 2 .

It namely eliminates the Bessel function I1/2(η2/4). Reaching a uniform velocity in the downstream Shercliff flow requires F+(0) = F- (0) = 1/2, which leads to

D + = D = 1 2 .

Finally, we get the rather simple explicit solution:

F ( η ) = 1 2 ( 1 + η 2 2 ) ( 1 erf η 2 ) η e η 2 / 4 2 π , F ' ( η ) = η 2 ( 1 erf η 2 ) e η 2 / 4 π , F ' ' ( η ) = 1 2 ( 1 erf η 2 ) . }


  1. Hartmann J, Lazarus F: K Den Vidensk Selsk Mat-Fys Meddl. 1937, 15: 1-27.

    Google Scholar 

  2. Shercliff JA: The theory of electromagnetic flow-measurement. 1962, Cambridge, UK: Cambridge University Press

    Google Scholar 

  3. Shercliff JA: A textbook of Magnetohydrodynamics. 1965, Oxford, UK: Pergamon Press

    Google Scholar 

  4. Shercliff JA: Magnetohydrodynamics (a 30 min. educational film). 1965, Educational Services Inc. for the National Committee on Fluid Mechanic Films, USA

    Google Scholar 

  5. Branover H, Vasil'ev AS, Gelfgat YM: Magnetohydrodynamics. 1967, 3: 61-65.

    Google Scholar 

  6. Slyusarev N, Shilova E, Shcherbinin EV: Magnetohydrodynamics. 1970, 6: 497-502.

    Google Scholar 

  7. Lielausis O: Atomic Energy Review. 1975, 13: 527-581.

    Google Scholar 

  8. Gelfgat YM, Lielausis O, Shcherbinin EV: Liquid metal under the action of electromagnetic forces. 1976, Riga, USSR: Zinatne Press

    Google Scholar 

  9. Tananaev AV: Flows in MHD ducts. 1979, Moscow, USSR: Atomizdat

    Google Scholar 

  10. 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 

  11. Andreev O, Kolesnikov Y, Thess A: Phys Fluids. 2007, 19: 039902-10.1063/1.2718400.

    Article  ADS  Google Scholar 

  12. Kulikovskii A: Isv Akad Nauk SSSR Mekh Zhidk Gaza. 1968, 3: 3-10.

    Google Scholar 

  13. Kulikovskii A: Isv Akad Nauk SSSR Mekh Zhidk Gaza. 1973, 3: 144-150.

    Google Scholar 

  14. Walker JS, Ludford GSS, Hunt JCR: J Fluid Mech. 1972, 56: 121-141. 10.1017/S0022112072002228.

    Article  ADS  MATH  Google Scholar 

  15. Walker JS, Ludford GSS: J Fluid Mech. 1972, 56: 481-496. 10.1017/S0022112072002460.

    Article  ADS  MATH  Google Scholar 

  16. Holroyd RJ, Walker JS: J Fluid Mech. 1978, 84: 471-495. 10.1017/S0022112078000282.

    Article  ADS  Google Scholar 

  17. 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 

  18. Petrykowski JC, Walker JS: J Fluid Mech. 1984, 139: 309-324. 10.1017/S0022112084000379.

    Article  ADS  MATH  Google Scholar 

  19. 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 

  20. Reed CB, Picologlou BF, Hua TQ, Walker JS: IEEE. 1987, 2 (87CH2507-2): 1267-1270.

    Google Scholar 

  21. 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 

  22. Hua T, Walker JS: Int J Engng Sci. 1989, 27: 1079-1091. 10.1016/0020-7225(89)90086-4.

    Article  MATH  Google Scholar 

  23. Sellers CC, Walker JS: Int J Engng Sci. 1999, 37: 541-552. 10.1016/S0020-7225(98)00088-3.

    Article  MATH  Google Scholar 

  24. Lavrent'ev IV, Molokov SY, Sidorenkov SI, Shishko AR: Magnetohydrodynamics. 1990, 26: 328-338.

    MATH  Google Scholar 

  25. Molokov S, Shishko A: Fully developed magnetohydrodynamic flows in rectangular ducts with insulating walls. Tech Rep 5247, Kernforschungszentrum Karlsruhe. 1993

    Google Scholar 

  26. 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. 2001

    Google Scholar 

  27. Molokov SY, Reed CB: Fusion Sci Tech. 2003, 32: 200-216.

    Google Scholar 

  28. Holroyd RJ: J Fluid Mech. 1979, 93: 609-630. 10.1017/S0022112079001956.

    Article  ADS  Google Scholar 

  29. 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 

  30. Votyakov EV, Zienicke E: Fluid Dyn Mat Process. 2007, 3: 97-113.

    Google Scholar 

  31. Alboussière T: J Fluid Mech. 2004, 521: 125-154. 10.1017/S0022112004001740.

    Article  ADS  MathSciNet  MATH  Google Scholar 

  32. Wong CP, Malang S, Sawan M, Dagher M, Smolentsev S: Fusion Engng Design. 2006, 81: 461-67. 10.1016/j.fusengdes.2005.05.012.

    Article  Google Scholar 

  33. Sommeria J, Moreau R: J Fluid Mech. 1982, 118: 507-518. 10.1017/S0022112082001177.

    Article  ADS  MATH  Google Scholar 

  34. Pothérat A, Sommeria J, Moreau R: J Fluid Mech. 2000, 424: 75-100. 10.1017/S0022112000001944.

    Article  ADS  MATH  Google Scholar 

  35. Araseki H, Kotake S: J Comput Phys. 1994, 110: 301-309. 10.1006/jcph.1994.1027.

    Article  ADS  MathSciNet  MATH  Google Scholar 

  36. Leboucher L: J Comput Phys. 1999, 150: 181-198. 10.1006/jcph.1998.6170.

    Article  ADS  MathSciNet  MATH  Google Scholar 

  37. Smolentsev S, Cuevas S, Beltran A: J Comp Phys. 2010, 229: 1558-1572. 10.1016/

    Article  ADS  MathSciNet  MATH  Google Scholar 

  38. Cuevas S, Smolentsev S, Abdou MA: J Fluid Mech. 2006, 553: 227-252. 10.1017/S0022112006008810.

    Article  ADS  MathSciNet  MATH  Google Scholar 

  39. Moreau R, Smolentsev S, Cuevas S: Magnetohydrodynamics. 2009, 45: 181-192.

    Google Scholar 

  40. 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 

  41. Moreau R: Magnetohydrodynamics. 1990, Dordrecht, The Netherlands: Kluwer Acad Pub

    Book  MATH  Google Scholar 

  42. Müller U, Bühler L: Magnetofluiddynamics in channels and containers. 2001, Berlin, Germany: Springer

    Book  Google Scholar 

  43. Knaepen B, Moreau R: Ann Rev Fluid Mechanics. 2008, 40: 25-45. 10.1146/annurev.fluid.39.050905.110231.

    Article  ADS  MathSciNet  Google Scholar 

  44. 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 

Download references


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.

Author information

Authors and Affiliations


Corresponding authors

Correspondence to René Moreau or Sergey Smolentsev.

Authors’ original submitted files for images

Rights and permissions

Open Access This is an open access article distributed under the terms of the Creative Commons Attribution Noncommercial License (, which permits any noncommercial use, distribution, and reproduction in any medium, provided the original author(s) and source are credited.

Reprints and Permissions

About this article

Cite this article

Moreau, R., Smolentsev, S. & Cuevas, S. MHD flow in an insulating rectangular duct under a non-uniform magnetic field. PMC Phys B 3, 3 (2010).

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI:


  • Vorticity
  • Head Loss
  • Hartmann Number
  • Core Flow
  • Fringe Field