Followed by a review of previous studies of magnetohydrodynamic (MHD) duct flows in a nonuniform 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 halfwidth in the magnetic field direction, and a is the halfheight). The suggested model includes a realistic threecomponent div and curlfree 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 threedimensional flow equations are reduced to a quasitwodimensional (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 twodimensional, while the induced electric current forms both crosssectional and axial circuits and is essentially threedimensional. A new parameter R = ε^{2}Re/Ha has been identified to characterize the role of inertia in duct flows with insulating walls (Re and Ha stand for the Reynolds and Hartmann numbers). Computations and analytical studies are performed for inertialess (R ≪ 1) and inertial (R ≫ 1) flows at ε = 0.2 for Re up to 300,000 resulting in new scaling laws for typical lengths, velocities, electric current densities and pressure drops, which provide a new theoretical basis for potential applications.
PACS Codes: 47.65.d, 47, 47.11.j
1 Introduction
Magnetohydrodynamic (MHD) flows in a nonuniform magnetic field, or in ducts of a varying crosssectional 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 socalled "Mshaped" velocity profile and associated vorticity. A few important experiments were made in Riga [5–8] and in Saint Petersburg [9], which revealed the spectacular behavior of such nonuniform flows. Incidentally, new phenomena of turbulence suppression in a fringing magnetic field was recently discovered experimentally in Ilmenau, Germany [10, 11]. Attempts to build a theory of those complex flows started with the papers by Kulikovskii [12, 13], where the concept of characteristic surfaces was first introduced. Although this concept is limited to the asymptotic case of a perfectly conducting fluid, it allows for qualitative predictions of the main features of many MHD flows, including those in a fringing field. Subsequently, different research groups, around Walker [14–23] and more recently Molokov [24–27], investigated such flows, using simplifying assumptions, e.g. limiting the analysis to inertialess flows. These authors obtained velocity and electric current distributions, and suggested first explanations of how these flows are organized. In particular, numerical results based on the inertialess flow model [20] have demonstrated a fair agreement with the experimental data for both pipe and rectangular duct flows in a strong magnetic field. However, other experiments for flows in a slotted duct at high Reynolds (Re ≈ 10^{5}) and moderate Hartmann numbers (Ha ≈ 10^{2}) have demonstrated strong inertial effects [9], which cannot be explained using inertialess models. As a matter of fact, except for Holroyd [28] and recent attempts [29, 30] to model the Ilmenau experiment [10, 11] most of the theoretical efforts are lacking direct comparisons with the experimental data, which often do not comply with the theoretical predictions. Quite recently, Alboussière reexamined the theory, paying a special attention to the full role of the Hartmann layers [31]. His considerations are, however, still based on a few striking assumptions that limit the model applicability.
The purpose of this paper is to provide a theoretical basis for practical applications, such as liquidmetal blankets of a fusion reactor, where the usual nondimensional parameters can reach very high values (e.g. Ha ≈ O (10^{3})  O(10^{5}) and Re ≈ 10^{5}), which are sufficiently higher than those in typical laboratory experiments (see [32]). Unlike many previous theoretical studies, we extend our analysis to inertial flows at high Reynolds numbers. In the blanket conditions, the liquidmetal 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 plasmaconfining 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 nonconducting rectangular duct of a uniform crosssection, located between the poles of a semiinfinite 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 halfwidth b taken in the direction of the magnetic field (Hartmann length) is much smaller than the duct halfheight a. In our frame of reference the predominant component of the fringing magnetic field is in the zdirection and varies with x from zero to the constant value B_{0}. We simplify the analysis by assuming a twodimensional fringing field produced by magnetic poles that are semiinfinite in the xdirection and infinite in the ydirection. The fluid flows in the xdirection and U_{0} stands for the typical velocity scale (e.g. mean bulk velocity). The magnetic yoke is semiinfinite in the xdirection with a length 2a_{
m
} in the ydirection and a gap 2b_{
m
} in the zdirection (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 nondimensional parameters that describe the flow in a fringing magnetic field, are the Reynolds number Re = U_{0}a/ν, the Hartmann number , and the duct aspect ratio ε = b/a. Here, σ, ρ and ν are the electrical conductivity, density and kinematic viscosity of the fluid, respectively.
Most of previous theoretical investigations are based on limiting assumptions, which are revisited here. First, since our magnetic field is div and curlfree, we do not restrict our considerations to the socalled 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 wellknown 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 Mshaped velocity profile is formed and the local current density can differ significantly from its averaged value. Despite this limitation, Alboussière's analysis has the merit to renew the basic scaling law in this problem, introducing a new important length scale of the order of Ha^{1/4} in the case of a pipe flow.
A third strong common assumption is the inertialess approximation [16] in which inertia in the core flow is neglected in comparison with the Lorentz force and the pressure gradient. This implies that, except the side layers where viscosity is significant, the Lorentz force is purely irrotational and thus cannot generate any vorticity. The rationale for this approximation is the assumption that the ratio of the Lorentz to inertia forces scales as the interaction parameter N = Ha^{2}/Re, which in fusion blanket conditions is much larger than unity. However, in electrically insulating ducts, this scaling law appears not to be correct any longer providing the exchange of current between the core and the Hartmann layers is taken into account. As a consequence, the Lorentz force is in fact Ha times smaller, and the relevant parameter is Ha/Re. As applied to the fusion blanket conditions, this suggests that in the blanket flows inertial effects are not necessarily negligible.
Notice that among various designs of a liquidmetal blanket, the EU HeliumCooled LeadLithium (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 DualCoolant LeadLithium (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 selfcooled 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 curlfree magnetic field is derived in Section 2. This magnetic field is used in Section 3 in the derivations of the basic quasitwodimensional (Q2D) equations governing the reference flow in a fringing magnetic field. These equations are obtained from the original threedimensional 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 threedimensional. 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 curlfree 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 zdirection, 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:
Here B_{0} stands for the uniform magnetic field, which is supposed to exist in the limit x →∞. At the leading order in terms of ε_{
m
} ≪ 1, equations (1) satisfy both the conditions ∇ · B = 0 and ∇ × B = 0, precisely because the apparently small components B_{
x
} and B_{
y
} are not neglected. Indeed, whereas these components are O(ε_{
m
}) and smaller than the B_{
z
}component, their zderivative 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:
Note that the electric charge conservation principle, ∇ · j = 0, takes its classic form
only if the zdependent 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 threedimensional dipole. If this dipole were twodimensional, as in the case of infinite magnet poles in the ydirection, the far field should vary as x^{2}. In the reference case, when the poles are semiinfinite in the xdirection and infinite in the ydirection, 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 , the overall length over which the flow exhibits Hartmann layers is b_{
m
}Ha in the case of a twodimensional very long magnet. It is b_{
m
}Ha^{1/2} in the case of a twodimensional short magnet, and b_{
m
}Ha^{1/3} in the case of a threedimensional short magnet.
2.2 A correct expression for the twodimensional fringing field at the entry of a magnet
Here we derive the expression of a twodimensional fringing field when the magnetic poles are semiinfinite in the xdirection and infinite in the ydirection. Consider the conformal mapping between complex planes r and t defined by the relation
where K is a complex constant which has to be determined. This mapping transforms the halfplane Im(t) > 0 into the region between the symmetry axis and the yoke of the semiinfinite magnet in the rplane. 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 twodimensional magnetic field. Starting from the complex potential of a twodimensional source flow located at the origin of the tplane,
and, applying transformation (4) with the particular value K =  ib_{
m
}/π (i^{2} =  1), we easily get the xdependence of the B_{
z
} = B(x) component, which varies from B_{0} for x ≪ 0 to zero when x →∞. Its expression may be written under the parametric form
where C is an arbitrary constant, which can be defined from the condition B(x = 0) = B_{0} /2. After changing the sign of x in order to get x increasing in the flow direction as B(x) from zero to B_{0}, and denoting
the basic expressions for the magnetic field and its derivative take the following form:
Here, a and b denote the halflengths of the duct crosssection, 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 twodimensional div and curlfree magnetic field defined from equations (1) and (8) is used as the reference field for numerical and analytical solutions in Section 4 and Section 5.
3 The basic set of equations for the flow in a fringing magnetic field
3.1 Comments on the characteristic surfaces
Consider the flow of an incompressible and electrically conducting fluid in the rectangular duct shown in figure 1. Since the duct width is constant (2b), and since ε_{
m
} ≪ 1, the characteristic surfaces coincide with the crosssections 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 crosssection dA and length 2b, located between the Hartmann walls and moving along the trajectory of the fluid across the fringing field. Because of mass conservation, its volume 2bdA is invariant. Now, if we admit, as Kulikovskii, that the magnetic flux BdA across the area dA is also invariant, taking the ratio of both invariants, we get the property that 2b/B and therefore B must also be invariant along a fluid trajectory. In other words, this asymptotic limit only allows the streamlines and the current lines to lie on the planes where B is constant, and does not allow them to cross these planes. However, in fluids of finite electrical conductivity, the generalized form of the Alfvén theorem
Where Φ is the magnetic flux, implies that the electric current crossing the characteristic surfaces is determined by the Ohmic losses (locally j/σ). As a consequence, the whole Hartmann layers (not only their portions overlapping the side layers), where these Ohmic losses have the highest values, must be taken into account on both sides of the core. This differs from the Holroyd & Walker [16] ideas, according to which the only possibility for the fluid to cross the characteristic surfaces is within the side layers.
3.2 Core flow and side layers
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 q_{0}(x, y) + (z^{2}/2a^{2})q_{1}(x, y) +⋯ This is particularly important for the electric potential
whose parabolic zdependent 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 twodimensional part and its complementary part is such that:
It is noticeable that the highest order term of the wcomponent is a small O(ε^{3}) term, so that the incompressibility condition writes:
Relations (11) also yield the following expressions for the vorticity components in the zdirection:
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:
Here, the scalar quantity
introduces the influence of the nonuniform magnetic field. Notice that the term involving S_{0} and the derivatives of the function B(x, y), would not be present in the expression of j_{
z
} if the approximation of a straight magnetic field were used. Although we intend to focus on the particular example where B is only xdependent, it is worth keeping the full expression of S_{0} in the derivation of our general model equations and to limit ourselves to B(x) later on. The expression of the charge conservation principle, which comes out from (14), is:
Note the presence in (16) of , 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 . 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 ycomponents 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 nonzero electrical conductivity and forces some jump of the electric potential between its value at the edge of the layer and that at the wall:
The exponential forms of u, v, j_{
x
}, j_{
y
} and j_{
z
} are then:
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:
3.4 Within the Hartmann wall
Since the electric potential is continuous across the fluidsolid interface, the x and ycomponents 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
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
} :
This relation, which introduces the usual conductance ratio C = σ_{
w
}t_{
w
}/σb as a fourth nondimensional parameter, can still be significantly simplified when Ha ≫ 1, since C Δ_{⊥φH
}≪ φ_{
H
}/bδ. Then (21) becomes:
3.5 Influence of the Hartmann layer and Hartmann wall on the core electric current
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 excesscurrent density:
The integrals of the quantities on the lefthandside of equations (23a) and (23b), denoted by g_{
x
} and g_{
y
}, must satisfy the global equation
which can be written as
Clearly, in the limit of high Ha numbers, equation (25) reduces to
Therefore, φ_{
H
} can be substituted into (22) to get
Equation (27) can be understood as a necessary relation between the twodimensional and threedimensional 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:
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):
In this simple case, we just recover the same expression of j_{0z
}as when the magnetic field is uniform:
3.6 Equation of vorticity
Once equations for the twodimensional and threedimensional parts of the electric potential have been obtained, our next step is to derive the twodimensional 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 zcomponent of ∇ × (j × B), which reads (B · ∇)j_{0z
} (j_{0} · ∇) B since both B and j_{0⊥} are divergencefree. Using equations (1) for the magnetic field, the vorticity budget reads:
For the reference case C = 0 (insulated duct), we can further modify equation (32) by using expressions (14) for the x and ycomponents of the core current density j_{0⊥}, equation (31) for j_{0z
}, and equations (29) and (30) for φ_{0} and φ_{1}, respectively. The first and second terms on the righthandside of this new equation contain a purely twodimensional term and a parabolic zdependent term:
The same is true for the third term, which reads:
This definitely proves that, if ε ≪1, the core flow can be modeled with a purely twodimensional velocity field, by keeping only zindependent terms in equation (32):
At the same time, providing that ω_{0},ψ_{0}, φ_{0} and φ_{1} are known, the current density given by (14), remains threedimensional since it contains the zdependent terms and j_{0z
}≠ 0.
Equation (32) could also be used to derive a higher order equation for ω_{1}(x, y) by collecting terms proportional to z^{2}. In this investigation, however, we limit ourselves to a quasitwodimensional 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 nonzero derivatives ∂_{
x
}B and ∂_{
y
}B is clearly exhibited by the first and second terms on the righthandside of equation (35). The Hartmann damping of vorticity is represented by the third term . 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 nonuniform 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 zdependent terms. The derivation of the complete form of the vorticity transport equation would, however, require the introduction of terms proportional to z^{4} in the basic expression (10) for the electric potential. Such a model, would be a generalization of the approach presented in [34], allowing for high ε values.
Notice that in the particular case B = B(x), equation (35) can also be interpreted as an expression for the axial part of the current density component j_{0x
}(x, y, z = 0) = σ(∂_{x}φ_{0}  Bv_{0}), parallel to ∇_{⊥}B, which is expressed through other unknowns (velocity and vorticity). In other words, this equation explicitly demonstrates that, within the core, electric current lines must cross the characteristic surfaces due, essentially, to two main causes: inertia and Hartmann current sheet. As a consequence, a full threedimensional electric current distribution can accurately be analyzed (see Section 5), in which j_{0x
}and j_{0z
}are derived from equations (35) and (31), while j_{0y
}can be easily obtained from the charge conservation equation.
3.7 Basic set of nondimensional equations
Let us now introduce dimensionless variables using characteristic scales: the uniform magnetic field B_{0}, the mean bulk velocity , the duct halfwidth a as a length scale for x and y, and the halfheight b between the Hartmann walls for z. The new variables are denoted as follows:
In addition, S denotes the dimensionless expression of S_{0}. Since all key functions are only dependent on x and y, in what follows we drop the subscript ⊥. Then, the main set of equations, coming from their dimensional counterparts (29),(35) and (13), is:
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,
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 yindependent 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:
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:
where R = ε^{2}Re/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
If we also assume that the vorticity is slowly Xdependent, so that it can be approximated as and , equation (44) can be integrated once. Using the boundary conditions U = J_{
Y
} = 0 at the insulating side wall Y = 1, yields:
The three components (J_{
X
} , J_{
Y
} , J_{
Z
} ) exhibit the noticeable feature that the threedimensional electric current distribution is a combination of two families of well defined loops. One, the zindependent 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:
The other family, which is Xdependent, is located in the crosssectional (Y, Z) planes. In terms of the stream function h_{
c
}(X; Y; Z), the loops, quoted as crosssectional loops, are expressed as:
One should notice that equations (46) and (47) show explicitly the exchange of current between the core and the Hartmann and side layers. It is also clear that in the limit of negligible inertia (R = 0), the three current density components in the core have the same scaling law, all being proportional to Ha^{1}. Equations (46) have also the merit of exhibiting the contributions of inertia and viscosity on the axial current loop.
Let us now eliminate the electric potential between the charge conservation equation (41) and the vorticity budget (42). The Xderivative of the charge conservation equation is:
It is straightforward to substitute equation (43) for J_{
x
} in (48) to get
which is the basic equation for the vorticity in the fringing field. This equation expresses all the relevant contributions to the transport of vorticity, which are, respectively, the advection of vorticity, Hartmann damping, damping in the side layers, torque and exchange of current. Note that the advective term, Δ(RΓ/β'), involves the dimensionless parameter R = ε^{2}Re/Ha, which results as the suitable measure of the importance of inertia, in contrast with previous theories where the relevant scaling for inertia was N^{1} = ε^{2}Re/Ha^{2}.
4 Numerical results
Although the equations presented in previous sections are quasitwodimensional, 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 threedimensional approach, the present one has the advantage of not computing the flow in the Hartmann layers. However, resolution of highgradient flow regions associated with the side layers of thickness at the walls y = ± a is still needed to reproduce correctly the axial and crosssectional 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 timedependent behavior. Then, computing timedependent 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 wellunderstood MHD flows, e.g. fully developed flows in a rectangular duct under a uniform magnetic field. Typically, the Hartmann number in such studies was limited to Ha ≈ 10^{2} [35, 36]. The limitations were to some extent related to the computational time, which increases roughly exponentially with Ha. However, the most severe limitations are caused by a catastrophic loss of accuracy at high Ha numbers, which may occur in the course of computations due to accumulation of roundoff errors. A potential source of computational errors can be related to the use of the electric potential as the electromagnetic variable. As discussed by Smolentsev et al. [37], the loss of accuracy occurs when an induced electric current of small magnitude (O(Ha^{1})), is calculated from the Ohm's law as a difference between two large numbers, associated with ∇Φ and U × B_{0}, which are both O(1). In the present formulation, the nature of this problem can be seen directly from the equation for the potential (37). Two terms on the righthandside 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
The boundary condition at the walls Y = ± 1, which are assumed to be nonconducting, is ∂_{
Y
}J_{
X
} = 0, as follows directly from the expression for the Zcomponent 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 righthandside of the vorticity equation. In the postprocessing 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 threedimensional 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 finitedifference technique using a nonuniform 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 steadystate 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 centraldifference scheme, which is known to be dissipation free. Using convective term approximations free of numerical dissipation is particularly important when the flow is essentially timedependent. However, in the chosen range of parameters, including high values of Re (up to 300,000), the numerical solution always converges to a steadystate 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 4060 units. Satisfactory resolution is achieved by using 512 or 1024 grid points in Xdirection and 201501 points in Y direction. The numerical tests included grid sensitivity tests, and comparisons with the analytical solution for a fully developed flow and unsteady vortical flows [37, 38]. The comparisons demonstrated coincidence in local values between the analytical and the numerical results in four digits.
4.1 Numerical results for flows with weak inertia (R ≪1)
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].
For those two Hartmann numbers, all distributions are qualitatively similar. When the liquid proceeds through the fringing field region, the velocity at the axis U_{0}(X) decreases, the nearwall jets form, and the maximum velocity U_{1}(X) in these jets increases (typical data are shown in figures 4 and 5). In agreement with the previous results by Lavrent'ev et al. [24], these features are spectacular, although U_{0}(X) does not reach a zero value (see figure 4) and the maximum value of U_{1}(X) increases as U_{1max
}≈ 0.3Ha^{1/2} instead of 2.38Ha^{1/2}, our theoretical prediction for the jet region (see Section 5). Each jet exhibits two vorticity peaks of opposite signs. The value of the higher peak, which has the same sign as the side layer vorticity, can be estimated as U_{1}/δ_{
SL
}, where δ_{
SL
} stands for the thickness of the side layer, and yields an estimate of U_{1}. The other peak, of the opposite sign and much smaller value, can be estimated as U_{1}/δ_{
j
}, where δ_{
j
} yields an estimate of the jet thickness. The large ratio between the two peaks also indicates the large ratio between the jet thickness δ_{
j
} and that of the side layer δ_{
SL
}, as shown in Section 5.
In the asymptotic case of large Ha and small R, these computed results suggest that the whole flow domain can be subdivided into five characteristic regions: i) the far upstream developing Sherclifflike flow, where the velocity profile is uniform with a moderate overvelocity to compensate the missing flow rate in the side layers; ii) the upstream region where the Mshape starts to form, whose position is close to abscissa X ≈ 15 and varies very slowly with Ha; iii) the jetregion, 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 Hadependent; iv) a short transverse layer, whose length seems significantly smaller than the duct half width, suggesting that the Xderivatives should be locally larger than the Y derivatives; and finally v) the fullyestablished Shercliff flow with its familiar plateautype 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 threedimensional structure in the jet region to a purely crosssectional current distribution in the fullyestablished Shercliff flow. In turn, in the jet region the threedimensional current density demonstrates two characteristic families of current loops: an internal one where the current lines are almost quasitwodimensional and closed in the core, and an external one where the current lines passing through the jets and side layers are strongly threedimensional 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 subdomain.
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 57 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 Rvalue 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 twodimensional. This indicates that the current exchange mechanism, important in inertialess flows, becomes negligible in flows where inertia is significant or strong.
5 Guiding ideas and scaling laws for flows with high Ha and small ε
5.1 General properties
Figures 6 and 7 show that, when inertia is significant, the flow is quasiparallel 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 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 quasiparallel and hence, becomes predominant upon (see figure 3).
According to the approximation made in Section 3.6, valid for slotted ducts, we neglect the terms proportional to ε^{2}Z^{2} so that the components of the Lorentz force read
Introducing the nondimensional pressure defined as
the X and Y projections of the motion equation read:
Similarly, equation (48) for J_{
X
} can also be simplified as
In quasiparallel flows, the righthandside of (54) is negligible when Ha ≫ 1, yielding the simple relations
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 P_{0}(X), and its value within the side layer by P_{1}(X), equations (53) allow deriving the pressure from the following equations
At the side wall, where U = h_{
a
} = 0, equation (56b) yields the Xdependence of the pressure at this location:
Subtracting equation (56b) from its form at the edge of the layer, where the velocity reaches its maximum value U_{1}, the relevant equation for the side layer is
For the inertialess regime (R = 0), in the vicinity of the side wall Y = 1, an elementary solution for the side layer exists, namely,
However, at inertial regimes, the accelerated side layer does not obey the selfsimilar 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
When R is large enough, the Hartmann friction becomes negligible and equation (60) implies the invariance of the missing momentum: . Consequently, evaluating the Const in the far downstream Shercliff flow, we get the following evolution for the side layer thickness, in terms of the still unknown jet velocity U_{1}(β):
From these properties of the side layer, it follows that equation (57) for the pressure variation along the duct axis reads:
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
In turn, the stream function of the current density obeys (55) and can be expressed as
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:
It is noticeable that (ΔP)_{2} is weakly dependent on R through the Rinfluence 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 Rindependent and reads:
Its maximum (ΔP)_{2} = (πHa/10ε)(3/5)^{3/2} ≈ 0.146(Ha/ε) is reached for .
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
The terms of this linear equation express all major effects: from left to right, we have advection of vorticity , Hartmann damping , damping in the side layers , torque (Haβ'^{2}Ω) and exchange of current (β'^{2}∂_{β}Ω).
Let us express the velocity field in the quasiparallel 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:
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:
Note that an additional term of the order of U_{1}δ_{
SL
}Y might be introduced in the expression for Ψ, to take into account the missing flow rate in the side layers. However, as suggested by (61), this term is small enough to be neglected in this asymptotic approach.
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:
The third term of equation (70) is negligible since ε^{2}/Ha ≪ 1 and . 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 nontrivial contribution of the side layer to the global balance of vorticity. Notice that the axial velocity U_{0} and maximum velocity in the jet U_{1}, can be written as
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
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 Mshape extending over the whole width of the duct. The amplitude of the Mshape prifile can be characterized by the quantity
Then, the induced electric current is easily derived from the expressions of the streamfunctions h_{
a
} and h_{
c
} introduced in (46), (47) and (55), and the pressure follows from equations (56) and (62).
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),
which keep linear the equations for F (α) and G(α) = α^{2}F:
The form of the lefthandside of equations (77a) and (77b) and the expression (8) for β' = (π/2) β^{
2
}(1 β^{
2
}), suggest the transformation
Hence, the differential equations (77a) and (77b) become
Eliminating C and introducing the rapidly varying variable
such that t' = Haβ'/R, yields the second order equation
Since the previous approximations have already put the emphasis on regimes with significant inertia, such that R ≫ 1, the righthandside 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:
The functions F (β) and G(β) are therefore given by:
It is now elementary to derive expressions for most characteristic quantities, such as the velocities on the axis, U_{0} = 1  f F , and at the edge of the side layer, U_{1} = 1 + gF :
From the asymptotic behavior of these typical velocities near the inlet, it can be assured that no singularity appears when β → 0:
In particular, in most experimental conditions, these velocities deviate from their asymptotic limit as Nβ^{3}/ε, where N = Ha/R is the interaction parameter. This shows that the typical length of the upstream region where the Mshape profile is formed, scales as L_{
up
} ≈ (ε^{2}Ha/R)^{1/3} and becomes shorter the higher the R values. This can be interpreted as the squeezing of the upstream region against the fringing field by the advection of vorticity, clearly observed in figures 6 and 7.
The analogous derivation in the outlet vicinity, which determines how the Mshape 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)e^{πX/ε
}, so that expressions (84a) and (84b) become:
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
These equations show that the downstream length scale between the abscissa where the Mshape profile is maximum and the Shercliff flow, scales as L_{
down
} ≈ R and is Haindependent. Now, in the opposite condition where inertia is strong, such that πHa/180εR ≪ 1, the term e^{t/12}can be approximated by 1  t /12, so that equations (86a) and (86b) reduce to
The correction for the U_{0} value is small since t/6 ≪ 1, but the correction for U_{1} is important and implies a reduction of its maximum value by a factor of the order Ha/εR. Nevertheless, the typical length necessary to reach the Shercliff profile remains controlled by e^{X/R
}and scales as L_{
down
} ≈ R. This downstream transport of the Mshape 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 U_{1}(X) are very complex and justify an approximate treatment of the integral in order to get an analytic prediction of this important quantity for potential applications. This approximate treatment is detailed in Appendix 1. The idea consists in replacing the term by the approximation 5(U_{1}  1) since the maximum velocity U_{1}(X) remains close to unity as shown in figures 6 and 7. Then we use two asymptotic branches for U_{1}  1, one for the inlet vicinity (β ≪ 1) and one for the downstream transition toward the Shercliff flow (1  β ≪1), instead of the exact expression. Clearly, the maximum of U_{1} given by the intersection of these two branches is overestimated, but it is reasonable to expect that the full integral yields the correct dependence on the key parameters (ε≪ 1, Ha ≫ 1, R ≥ 1). This yields the following asymptotic expressions:
5.3 The particular case of flow regimes with zero or weak inertia
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 Sherclifflike 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 Mshape 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
which allows combining the last two terms of (49) into only one of order unity, so that the inertialess form now reads
The implications of that equation in the core flow appear particularly strong, since (91) reduces to the simple form
This is a heatlike equation with the nonuniform 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 Mshape profile is formed, nor the jet thickness. The Hartmann number influences only the local value of Ω(β) within this region through the amplitude coefficient e^{
Haβ
} . It is worth noticing that the two terms involving ∂_{
β
}Ω in (67) illustrate the competition between the downstream advection of vorticity, predominant in the case analyzed in Section 5.2 and the upstream diffusion due to the exchange of current, predominant in the present analyzed case.
5.3.1 The upstream Mshape profile formation region
In this region where β ≪ 1, providing that the assumption is still justified, the local equation (69) reads
Since F and F' are O(1), the coefficient of F must be much smaller than Ha, so that:
The solution of equation (94) such that α (0) = 0 is
and shows that α deviates from its zero value in the far upstream region when β^{3} ≈ (ε^{2}/Ha)(1+3πR/2ε). Therefore, the length over which the Mshape profile develops is L_{
Msh
} ≈ (εHa)^{1/3} in the inertialess limit. In the case of flows with weak inertia such that R » 1, we just recover the previous estimate obtained in Section 5.2: L_{up} ≈ (ε^{2}Ha/R)^{1/3}It is worth noticing that this estimate combines two effects: the upstream diffusion in a geometry independent of H a, and the vorticity distribution, which is itself H adependent.
Besides, using the asymptotic value f = 1/6 and neglecting the exchange of current since β « 12ε^{2}R/π in this region, the global equation (73) can be reduced to the simple form
In the inertialess limit the solution is straightforward:
For nonzero but weak inertia, it is also simple to derive the small influence of the vorticity advection in the vicinity of the limit β →0:
The components of the current density in this region immediately follow from equation (46) for the axial loop and (47) for the crosssectional 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, quasiparallel 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:
If the two first terms of equation (99) are predominant upon the exchange of current, a selfsimilar solution can be found such that U = U_{0} + (U_{1}U_{0}) f(η), with η = (y +1)/δ_{j}(β)where δ_{
j
} stands for the jet thickness. It implies that:
The conservation of the flow rate requires that the stream function Ψ =U_{0}Y(U_{1}U_{0})δ_{
j
}e^{η} must take the value at the edge of the side layer (η = 0). This yields one equation to determine the two unknowns U_{0} and U_{1}. The complementary equation comes from the global equation
which requires U_{0} = 0. Note that the numerical results suggest values significantly smaller than unity, but nevertheless positive, even for Ha = 10^{3} or larger. This small discrepancy is certainly a consequence of the quasitwodimensional approximation that neglects the exchange of current and derivatives 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 U_{1} is completely determined by the conservation of the flow rate:
This expression shows that U_{1}(β) has a maximum at the abscissa where δ_{
SL
}(β) is minimum, namely, β_{
min
} = (3/7)^{1/2} ≈ 0.655, or X_{
min
} ≈ 0.366 for ε = 0.2. In turn, U_{1max
}must vary as Ha^{1/2}/ε. This behavior can be checked on figure 8, although the numerical value of U_{1max
}, close to , is definitely smaller than the predicted value, close to . 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 Mshape formation region, the typical length of the jet region must be L_{
jet
} ≈(εHa)^{1/3}. In the computed cases (ε= 0.2, Ha ≈ 10^{3}), this yields L_{
jet
} ≈ 6.
It is straightforward to derive the expressions of the current density and pressure from the selfsimilar solution:
In turn, the Y component of the current density in the motionless core becomes:
In fact, the current density component has also a maximum at an abscissa slightly different from that where δ_{
j
} is minimum and U_{1} maximum. This behavior is in fair agreement with the numerical data shown on figure 9(a). The maximum for the positive quasiuniform J_{
Y
}component in the almost motionless core, as well as a strong concentration of the current lines in the transverse layer, where their sign is changed, are exhibited.
5.3.3 The transverse layer
Let us now consider the vicinity of β = 1, and denote β = 1  ν with ν ≪1. In this region, the assumption and the heatlike equation (92) are not valid since they would yield a singularity at β = 1. This means that the Xderivatives 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 lowvelocity core and the downstream fully developed Shercliff flow. The side layer with the thickness still exists along the side walls. Based on assumption , a relevant equation for the transverse layer can be derived from equation (49). Keeping only the predominant terms, this equation reads
If we take β'^{2}/β = (πν/ε)^{2} and introduce the transformed abscissa ζ = 2(εv/π)^{1/2}to avoid the Xdependent coefficient, equation (105) can be transformed in the elementary equation
This is a wellknown equation in MHD duct flows in the presence of a uniform magnetic field although, usually, a coefficient Ha^{2} appears instead of Ha. This equation characterizes all kinds of free shear layers parallel to the applied magnetic field, often named inertialess Ludford layers, located between cores having different properties. The main novelty introduced by (106) is the thickness of the transverse layer, of the order Ha^{1/4} in ξunits, instead of Ha^{1/2} in the case of Ludford layers. This layer has a clear similarity with the free shear layer discovered by Alboussière [31] along a singular characteristic surface "which splits" in ducts having a nonuniform 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 crosssections of the rectangular duct). Indeed, it appears that this layer exists as soon as the slowly Xdependent solution exhibits a singularity. In our case, the origin of the singularity lies in the solution for the Mshape profile region, when β →1 and β' →0. Hence, the transverse layer provides the smooth transition between the Mshape 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. 138150; Müller & Bühler [42], pp. 8589). 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 Ydirection and diffusing in the ± Xdirections, and the other issuing from Y = +1 with a zero thickness, propagating in the negative Ydirection and diffusing also in the ± Xdirections. 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 wellknown technique can be applied for its solution, which can be epxressed as a combination of two selfsimilar solutions:
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 midheight of the upper half duct computed for a very high value of the Hartmann number.
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 , or δ_{
ζ
} = Ha^{1/4}, its thickness in β and Xunits must be
This prediction of a Haindependent thickness in Xunits 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 Rvalues to maintain inertia negligible. This actually demonstrates that the properties of this flow in a fringing field are not generic enough to be analyzed without making the choice of a specific β(X) distribution. For instance, would we use the distribution β' = (πν/ε)^{
n
}, the above derivation would yield δ_{
X
} ≈ (Ha^{
n1}ε^{2n
})^{1/(42n)}. 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
together with the condition of charge conservation, which yields
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:
It can be verified that the global current is Yindependent and takes the remarkably simple value . 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 , in the symmetry plane (Z = 0) of the transverse layer are
From the specific properties of the transverse layer, we get:
The pressure profile finally reads
We observe that pressure differences seem to scale as Ha^{1/2}/ε and are, consequently, much higher than pressure differences in both the Mshape 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:
which yields the expected linear dependence with the mean velocity . 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:
This a priori surprising result for the inertial regime is a consequence of the required balance between inertia (proportional to ) and Lorentz force (proportional to ).
6 Discussion and conclusion
In this paper, the problem of the flow in a nonuniform 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 curlfree magnetic field which models the inhomogeneities of the actual field at the entrance of a large yoke. A second important point is the explicit consideration of the Hartmann layers and the electric current exchange mechanism between them and the core flow. One of the key consequences of this mechanism is that the relevant scaling for the ratio of inertial and Lorentz forces in the case of insulating walls is R =ε^{2}Re/Ha instead of N^{1} = ε^{2}Re/Ha^{2}. Another important consequence is that the threedimensional electric current distribution is formed by a pair of twodimensional current loops, namely, one in the axial (x, y), and another in the crosssectional (y, z) planes. In the observed threedimensional 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 nonzero 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 quasiparallel flow assumption () and Oseen's approximation to flows with strong inertia we were able to construct scaling laws for the length over which the Mshape profile develops, which is found to be (ε^{2}Ha/Re)^{1/3}, and for the velocity maximum in the form U_{1}  U_{0} ≈ Ha/εR. These scaling laws are in qualitative agreement with the numerical predictions. For comparison, in the inertialess regime this length scale was found as (εHa)^{1/3} and the maximum velocity of the Mshape profile as U_{1}  U_{0} ≈ 1. Besides, the predicted downstream advection of the Mshape 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 threedimensional, with more currents closing in the axial planes, as inertia increases. Hence, for the inertial regime the exchange of current between the core and the Hartmann layers seems to be much less important than for inertialess flows.
Another remarkable property of the regimes with strong inertia is that the computational results remain quite stable, even at large Reynolds numbers (Re ≈ 3 ×10^{5}). The stability of the computations is in fact in agreement with recent observations of turbulence suppression by Andreev et al. [10, 11] in the experiments with a short magnet. This striking feature can be considered as one of the most noticeable properties of the flows in fringing fields. While in a strong uniform field quasitwodimensional flow disturbances are promoted and damped slowly for a typical Hartmann time (see [33, 43]), in a region of nonuniform 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 xdirection. In the nonuniform 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 nonlinear interactions between vortices and the mean flow, and the instability mechanisms based on inflectional velocity profiles. Secondly, the quasitwodimensional vortices typical of MHD turbulence, which travel across the characteristic surfaces, are more rapidly damped than vortices elongated in the ydirection, 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 crosssection 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 Mshape velocity profile develops, provided that the side friction is neglected, the basic vorticity equation reduces to a heatlike 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.
Let us emphasize that the choice of a particular div and curlfree twodimensional magnetic field is particularly important, instead of assumptions on the length scale over which this field varies. Indeed, many crucial results, such as the thickness of the transverse layer and the final head losses, would be different with another magnetic field distribution. There remains some questions where the asymptotic analysis does not provide any clear answer, such as the dependence of the maximum velocity U_{1max
}reached in the jets on the parameters Ha and ε. Additional important points, such as the flow properties when ε ≥ 1, or when C ≠ 0 would deserve to be reexamined. 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 U_{1} ≈ 1 (see figures 5 and 7), the expression can be approximated by 5(U_{1} 1) in the integral of the equation (63b). Instead of the full expression (84) for U_{1}, let us use the two asymptotic branches, easily derived from (85), (87), (88) when ε ≪R:
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:
Then the integral of (63) can be approximated as:
Substituting the values of β_{*} and X_{*} yields: . In the exponent of the exponential, the ratio Ha/R^{4} = N/R^{3} can be assumed quite small, so that the exponential can be taken as unity. Similarly, the factor (8/9)^{1/3} can be taken as unity. This yields the following estimate for the head loss:
In most cases where N ≫ 1, the number R = Ha/ε^{2}Re is small enough to justify neglecting the second term at the righthandside of (120). Then the estimate for the head loss becomes:
Consider now the opposite case where inertia is strong: N ≪1. In the expression of U_{1}  1 for β ≫ 1, N ≪1, substitute to the varying variable t its limit when β = 1: t_{
max
} = πHa/15εR.
The asymptotic branches intersect at . 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:
The last term at the righthandside is predominant, yielding the following estimate for the head loss:
It is noticeable that, if N ≫ 1 the total head loss is controlled by the initial asymptotic branch where U_{1} increases toward its maximum, whereas it is controlled by the final asymptotic branch where U_{1} decreases from its maximum if N ≪1.
Appendix 2. Solving the differential equation (106)
Let us look for a selfsimilar solution of the form (107). In ξunits, this forms combines two symmetric selfsimilar layers with thicknesses:
Substituting (124) into (106) yields the following 4th order ordinary differential equation for either F_{+}(η_{+}) or F_{}(η_{

}):
We first disregard the obvious solution F = constant and denote F' = η f(η). Equation (125) becomes:
Using the transformation
reduces once more the differential order of equation (126) and provides the modified Bessel equation of order 2:
The general solution reads y(t) = AI_{1/2}(t) + BK_{1/2}(t). It leads to the general expression for the basic function F (η):
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 This requires
It namely eliminates the Bessel function I_{1/2}(η^{2}/4). Reaching a uniform velocity in the downstream Shercliff flow requires F_{+}(0) = F_{} (0) = 1/2, which leads to
Finally, we get the rather simple explicit solution:
Declarations
Acknowledgements
We are indebted to Mohamed Abdou, Thierry Alboussière, Neil Morley, Mingjiu Ni, Ramakanth Munipali and the group at Hypercomp Inc., for fruitful discussions. RM is indebted to the UCLA fusion group for supporting his visits during the accomplishment of this work. SS acknowledges support from the U.S. Department of Energy via grant DEFG0286ER52123A040. SC thankfully acknowledges support from CONACYT under project 59977.
Authors’ Affiliations
(1)
Laboratoire SIMAP/EMP, GrenobleINP
(2)
Mechanical and Aerospace Eng. Dept, University of California
(3)
Centro de Investigación en Energía, Universidad Nacional Autónoma de México, A. P. 34
References
Hartmann J, Lazarus F: K Den Vidensk Selsk MatFys Meddl. 1937, 15:1–27.
Shercliff JA: The theory of electromagnetic flowmeasurement. Cambridge, UK: Cambridge University Press; 1962.
Shercliff JA: A textbook of Magnetohydrodynamics. Oxford, UK: Pergamon Press; 1965.
Shercliff JA: Magnetohydrodynamics (a 30 min. educational film). Educational Services Inc. for the National Committee on Fluid Mechanic Films, USA; 1965.
Slyusarev N, Shilova E, Shcherbinin EV: Magnetohydrodynamics. 1970, 6:497–502.
Lielausis O: Atomic Energy Review. 1975, 13:527–581.
Gelfgat YM, Lielausis O, Shcherbinin EV: Liquid metal under the action of electromagnetic forces. Riga, USSR: Zinatne Press; 1976.
Tananaev AV: Flows in MHD ducts. Moscow, USSR: Atomizdat; 1979.
Andreev O, Kolesnikov Y, Thess A: MHDchannel 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.
Andreev O, Kolesnikov Y, Thess A: Phys Fluids. 2007, 19:039902.View ArticleADS
Kulikovskii A: Isv Akad Nauk SSSR Mekh Zhidk Gaza. 1968, 3:3–10.
Kulikovskii A: Isv Akad Nauk SSSR Mekh Zhidk Gaza. 1973, 3:144–150.
Walker JS: Threedimensional laminar MHD flows in rectangular ducts with thin conducting walls and strong transverse nonuniform magnetic fields. In Liquid Metal Flows and Magnetohydrodynamics, Proc. Third Int. BerSheva Seminar in the MHD Flows and Turbulence Series. Edited by: Branover H, Lykoudis P, Yakhot A. AIAA; 1983:3–19.
Walker JS, Petrykowski JC: Approximate side layer solutions for a liquid metal flow in a rectangular duct with a strong nonuniform magnetic field. In Proc Fourth Int BerSheva Seminar on MHD Flows and Turbulence. Edited by: Branover H, Yakhot A. AIAA; 1985:17–31.
Molokov S, Shishko A: Fully developed magnetohydrodynamic flows in rectangular ducts with insulating walls.Tech Rep 5247, Kernforschungszentrum Karlsruhe 1993.
Molokov S, Reed CB: Liquid metal flow in an insulated circular duct in a strong nonuniform magnetic field. Part 2. Inclination of the field gradient to the duct axis.Tech rep., ANL/TD/TM01–19, Argonne National Laboratory 2001.
Votyakov EV, Zienicke E, Thess A: MHDchannel 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.
Votyakov EV, Zienicke E: Fluid Dyn Mat Process. 2007, 3:97–113.
Moreau R, Smolentsev S, Cuevas S: Why and how turbulence can be suppressed in a nonuniform magnetivc field.Proc 6th Int. Conf. on Electromagnetic Processing of Materials 2009, 37–40.
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.