Simulation of magnetohydrodynamic flows of liquid metals with heat transfer or magnetic stirring

We discuss the effects of nonhomogeneous magnetic fields in liquid metal flows in two different configurations. In the first configuration, we briefly report the impact of fringing magnetic fields in a turbulent Rayleigh–Bénard convection setup, where it was shown that the global heat transport decreases with an increase of fringe‐width. The convective motion in regions of strong magnetic fields is confined near the sidewalls. In the second configuration, we numerically study the effects of an oscillating magnetic obstacle with different frequencies of oscillation on liquid metal flow in a duct. The Reynolds number is low such that the wake of the stationary magnetic obstacle is steady. The transverse oscillation of the magnet creates a sinusoidal time‐dependent wake reminiscent of the vortex shedding behind solid obstacles. We examine the behavior of the streamwise and spanwise components of the Lorentz forces as well as the work done by the magnets on the fluid. The frequency of the oscillation of the streamwise component of Lorentz force is twice that of the spanwise component as in the case of lift and drag on solid cylindrical obstacles. The total drag force and the energy transferred from the magnets to the fluid show a nonmonotonic dependence on the frequency of oscillation of the magnetic obstacle indicative of a resonant excitation of the sinusoidal vortex shedding.


Introduction
Magnetohydrodynamic (MHD) flows, i.e., flows of electrically conducting fluids under the influence of magnetic fields, are frequently encountered in engineering and astrophysical applications.In such flows, the fluid is acted upon by the Lorentz force in addition to the force driving the flow.Industrial and technological applications of such flows include heating, pumping, stirring, and levitation of liquid metals, cooling blankets in fusion reactors, and liquid-metal batteries.In the context of astrophysics, magnetic fields strongly influence the flows in the sun and the stars and are responsible for the formation of sunspots and solar flares.
Magnetoconvection has been studied extensively in the past, but most of the studies focused on flows under the influence of a homogeneous magnetic field, which is an idealized approximation.However, in most engineering and astrophysical applications (such liquid metal batteries, cooling blankets in fusion reactors, electromagnetic stirring, solar spots, etc.) the magnetic fields are localized and thus vary in space [1,2].Further, strong homogeneous fields in large regions of space can only be generated by magnets of large size which are difficult to design and very costly to build and operate [3].Thus, it is important to understand the impact of spatially varying magnetic fields on magnetohydrodynamic flows.Recently, Bhattacharya et al. [4] studied the effects of spatially varying magnetic fields on MHD flows driven by buoyancy (magnetoconvection); these effects will be briefly summarized in this paper.There have been several studies on MHD duct flows with different configurations of spatially varying fields (see, for example, Sterl [5] and Prinz et al. [6]); however, in this paper, we focus specifically on duct flows with a localized zone of applied magnetic field (henceforth referred to as magnetic obstacle).Flows past stationary magnetic obstacles have been studied before [7][8][9][10][11].Similarities and differences between stationary magnetic and solid obstacles has been discussed by Votyakov and Kassinos [12].Unsteady wakes were only found for fairly large Reynolds numbers where the flow develops small-scale turbulent eddies [10].In order to realize an unsteady flow past a magnetic obstacle at a rather low Reynolds number it seems necessary to add an additional periodic motion of the magnet.We therefore consider the effects of oscillating magnetic obstacles on MHD duct flow in the present paper, which can be interesting in the context of magnetic stirring.We also remark that oscillating solid obstacles have been studied previously but it appears that such studies are lacking for magnetic obstacles so far.
The outline of the paper is as follows.In Sec. 2, we discuss the mathematical model.Section 3 describes the numerical method used in our simulations.We discuss the results in Sec. 4 and conclude in Sec. 5.

Mathematical model
In this section, we describe the setup and the mathematical formulation of our problems.The study will be conducted under the quasi-static approximation, in which the induced magnetic field is neglected as it is very small compared to the applied magnetic field.This approximation is fairly accurate for MHD flows of liquid metals [2].The governing equations of MHD flows are given by where u and p are the velocity and pressure fields respectively, ν is the kinematic viscosity, and f is the total body force acting on the fluid.For MHD duct flow, f is the specific Lorentz force (i.e.force per unit mass, henceforth denoted as f L ) and is given by where j is the current density, B 0 is the imposed magnetic field strength, σ and ρ are the electric conductivity and mean density of the fluid, respectively, and φ is the electric potential.
In magnetoconvection, f = f L + f b , where f b = αgT ẑ is the buoyancy force, α is the thermal expansion coefficient, g is the gravitational acceleration, and T is the temperature field.Magnetoconvection is additionally governed by the following thermal energy equation which describes the evolution of the temperature field T : where κ is the thermal diffusivity of the fluid.MHD liquid-metal duct flows are governed by two nondimensional parameters: the Reynolds number Re, which is the ratio of the inertial force to the viscous force; and the Hartmann number Ha, which is the ratio of the Lorentz force to the viscous force.Liquid-metal magnetoconvection is governed by three nondimensional parameters: the Rayleigh number Ra, the ratio of the buoyancy force to the dissipative forces; the Prandtl number Pr -the ratio of kinematic viscosity to the thermal diffusivity; and the Hartmann number Ha.These quantities are given by where U , L, and ∆ are the characteristic velocity, length, and temperature scales respectively.For magnetoconvection, we consider the Rayleigh-Bénard setup consisting of fluid enclosed between a cooler top plate and a warmer bottom plate (the temperature difference between the plates being ∆), with the plates separated by a distance H.In this case, H and ∆ respectively are the characteristic length and temperature scales for Ha and Ra.
As discussed in Sec. 1, we describe the effects of spatially varying magnetic fields in liquid metal flow for two configurations -(i) thermal convection in a box, and (ii) pressure-driven duct flow.In the first configuration, we consider a horizontally extended convection box of size l x × l y × H = 16 × 32 × 1 which is influenced by magnetic fields generated by two semiinfinite permanent magnets.The north pole of one magnet faces the bottom of the convection cell and the south pole of the second magnet faces the top of the convection cell.These magnets extend from −∞ to ∞ in the x-direction, l y /2 to ∞ in the y-direction, from near the top wall to ∞ in the positive z-direction, and from near the bottom wall to −∞ in the negative z direction.For a detailed description of the setup, the readers are referred to Bhattacharya et.al [4].In this configuration, the lateral component of the magnetic field (B x ) vanishes, and the longitudinal and vertical components respectively are logarithmic and inverse-tangent functions of the spatial coordinates y and z and the gap δ between the magnetic poles and the horizontal walls.The magnetic field distribution is such that its strength is negligible for 0 < y l y /2, increases steeply at y ∼ l y /2, and saturates close to its maximum value at y l y /2.When δ is increased keeping other parameters same, the total magnetic flux through the convection cell remains the same, but the gradient of the magnetic field at y ∼ l y /2 decreases, thereby increasing the fringe-width of the magnetic field.The aim of the study was to determine the effects of fringe-width on the heat and momentum transport.
The second configuration, which is the main focus of this paper, consists of liquid metal flow in a duct with two oscillating permanent magnetic poles near the top and bottom walls.The magnetic poles are semi-infinite in the z-direction and measure M x = 3 units along the streamwise direction and M y = 4 units along the spanwise direction in agreement with Votyakov and Kassinos [12].The spanwise dimension of the duct is L y = 50 units and the height is L z = 1 unit.The vertical gap between the magnetic poles and the liquid domain (one quarter of the layer height) also corresponds to Ref. [12].A schematic diagram of the setup is shown in Fig. 1.
The magnetic field B = (B x , B y , B z ) generated by the magnetic poles is given by the formula derived by Votyakov et al. [13].The oscillation takes place along the spanwise direction and the y-coordinate of the center of the magnet at time t is given by Copyright line will be provided by the publisher where A and f 0 are the amplitude and frequency of oscillation respectively.The magnets therefore have a velocity u m with respect to the flow domain.Since the induction of currents depends on the relative velocity between conductor and magnet, the difference u − u m must be used in Ohm's law (4b) and in the charge conservation condition (4c).
In our work, the oscillation amplitude is set to A = 1, i.e. the ratio A/M y = 0.25.For the frequency we choose a reference value based on the Strouhal number St 0 = 0.25 in Ref. [12].The nondimensional reference frequency in our work is therefore where U = 1 is the nondimensional mean streamwise velocity.The frequency ratio is defined as F = f /f s , i.e., the ratio of the frequency of oscillation of the magnetic poles to that of vortex shedding for the stationary magnetic obstacle of the same dimensions at a Reynolds number Re = 900 in Ref. [12].

Numerical method
We conduct direct numerical simulations of our setups using a second-order finite difference code developed by Krasnov et al. [14,15].For the magnetoconvection setup, a non-uniform grid of resolution 4800 × 9600 × 300 was used.All the walls were rigid and electrically insulated such that the electric current density j formed closed field lines inside the cell.The top and bottom walls were fixed at T = −0.5 and T = 0.5 respectively, and the sidewalls were adiabatic.The Rayleigh number, Prandtl number, and the Hartmann number based on the maximum value of the vertical magnetic field were fixed at Ra = 10 5 , Pr = 0.021, and Ha z,max = 120.The gap δ between the magnetic poles and the conducting plates was varied from δ = 0.01H to δ = 9H, where H is the cell height.
For the configuration of flow past oscillating magnetic obstacle, we employ a rectangular domain of dimensions L x × L y × L z = 200 × 50 × 1 with a grid resolution of 1024 × 384 × 32.The fluid enters the domain at x = 0 with a nearly fully-developed laminar flow profile that is approximated by the analytical expression The fluid leaves the domain at x = L x where ∂u/∂x = 0.The magnetic poles are located x = 50.The mesh is non-uniform in y and z-directions.The top, bottom, and sidewalls are rigid (no-slip) and electrically insulated.We fix Re = 400 and the Hartmann number based on the maximum vertical magnetic field as Ha z,max = 70, and vary the frequency ratio from F = 0.2 to F = 0.8.It must be noted that the characteristic length and velocity for the above nondimensional quantities are L z /2 (that is, half of the duct height) and the bulk horizontal velocity at the inlet (U ), respectively.For both the simulation setups, the elliptic equations for pressure, electric potential, and the temperature were solved based on applying cosine transforms in along the directions with uniform grid-stretching and using a tridiagonal solver along the direction with non-uniform grid stretching.The diffusive term in the temperature transport equation is treated implicitly.The time discretization of the momentum equation uses the fully explicit Adams-Bashforth/Backward-Differentiation method of second order.

Results
In this section, we first briefly summarize the results on the magnetoconvection simulations and then describe in detail the results on the flow past oscillating magnetic obstacle.

Results on magnetoconvection
A schematic of the magnetoconvection setup is illustrated in Fig. 2(a).The magnetic field generated by the magnets is strong enough to cease the flow in high magnetic flux region of the convection cell.We observe that as the local vertical magnetic field strength increases, the large scale structures become thinner and align themselves perpendicular to the longitudinal sidewalls.
The dependence of the local Reynolds and Nusselt numbers on the local Hartmann number (based on the vertical component of the magnetic field) was determined; this dependence was observed to be independent of the fringe-width.The global heat transport was observed to decrease with increasing fringe-width for strong magnetic fields but decrease with increasing fringe-width for weak magnetic fields.The convective motion became confined to the vicinity of the sidewalls in the regions of strong magnetic fields as shown in Fig. 2(b).The amplitudes of these wall modes were shown to exhibit a non-monotonic dependence on the fringe-width.
For further details on the results, the readers are referred to Bhattacharya et al. [4].In the next section, we discuss the results for the MHD duct flow setup.

Results on flow past oscillating magnetic obstacle
The simulations of the flow past magnetic obstacles are run for 300 convective time units after reaching a fully-developed state.The contour plots of instantaneous streamwise velocity are exhibited in Figs.3(a-c) and those with time-averaging in Figs.3(d-f) for F = 0.2, F = 0.5, and F = 0.8.The figures show regions of reduced and even reversed streamwise velocity in the regions of strong magnetic field and also in the wake of the magnetic obstacle.The regions of reduced instantaneous velocity exhibit a wavy pattern.It can be visually observed from Figs. 3(a-c) that as the magnets oscillate faster, the wavelength of spatial oscillation decreases.There is an increase in the amplitude of this path from F = 0.2 to F = 0.5, but the amplitude decreases with a further rise in F .For F = 0.5, the wake of the magnetic obstacle comprises of small-scale eddies, indicating that the flow becomes turbulent.The time-averaged streamwise velocity contours show that the length of the reversed flow region first decreases as F is increased to 0.5, and then increases with a further increase of F .
We examine the components of the total Lorentz force in the streamwise (f L,x ) and spanwise (f L,y ) directions.Note that f L,x and f L,y are the analog of the drag and lift forces, respectively, in flow past solid cylinders.These quantites are Figures 4(a,b,c) exhibit the plots of the above quantities versus the convective time t for F = 0.5, F = 0.6, and F = 0.8, respectively.The figures show a periodic sinusoidal time dependence.The magnitude of f L,x is higher than f L,y ; however, f L,y oscillates with a higher amplitude than f L,x .The amplitude of oscillation increases with an increase of F .The response frequency of f L,y is equal to the excitation frequency f 0 of the magnets; however, the response frequency of f L,x is twice of f 0 .We further compute f L,x t , the streamwise component of Lorentz force averaged over 300 timeframes, and plot it versus the frequency ratio in Fig. 5(a).It can be seen that f L,x t increases rapidly from F = 0.2 to F = 0.55.On further increase of F , f L,x t decreases sharply till F = 0.7 above which f L,x t saturates close to a constant value.Interestingly, the aforementioned behaviors of f L,x and f L,y closely resemble that of the drag and lift forces, respectively, in flows past an oscillating cylinder [16,17].
For flows past an oscillating cylinder, the non-dimensional mechanical energy transferred from the cylinder to the fluid is expressed as where t P is the motion period, d is the diameter of the cylinder, y is the spanwise position of the cylinder's axis, and f lift is the magnitude of the lift force [16].In our work, the energy transferred from the oscillating magnets to the fluid can be expressed similarly as follows: where t P = 300 convective time units for our case.We compute E for different frequency ratios and plot it versus F in Fig. 5(b).The figure shows that E is always positive, implying that the magnets perform work on the fluid for all frequencies.
The figure further shows that there is a gradual growth of E until F = 0.45 and then it sharply decreases to a minimum value at F = 0.53.The energy transfer increases monotonically on further increase of F .Interestingly, the point of minimum energy transfer almost coincides with the point of maximum average streamwise Lorentz force.This point corresponds to resonance where the velocity field exhibits stronger fluctuations compared to other frequency ratios.We finally examine the trends of the phase angle between the spanwise component of Lorentz force and the spanwise displacement of the magnets.This parameter is used as an indicator of energy transfer from the magnets to the fluids [16] where a phase angle between 0 and 180 degrees indicates positive energy transfer.We compute the phase angle using our data by fitting it with a sinusoidal function using the method of least squares.The computed phase angle is plotted versus the frequency ratio in Fig. 5(b).The figure shows that the phase angle lies between 0 and 180 degrees, consistent with the fact that the energy is transferred from the magnets to the fluid.The maximum phase angle at F = 0.55 reaches about 170 degrees.
Copyright line will be provided by the publisher (a) (b) × 10 −4 Fig. 5: For flow past oscillating magnetic obstacle: (a) Variation of the time-averaged Lorentz force in the streamwise direction with the frequency ratio, and (b) plots of energy transferred from the magnets to the fluid and the phase angle between the spanwise Lorentz force and displacement of the magnets.The above quantities exhibit a nonmonotonic dependence on F .

Conclusions
In this paper, we numerically examined the effects of non-homogeneous magnetic fields in liquid metal flows using a finitedifference fluid solver.We briefly summarized the results of Bhattacharya et al. [4] in which the influence of fringing magnetic fields on turbulent convection was studied.An important finding was that for strong magnetic fields, the global heat transport decreases with an increase of fringe-width, whereas for weak magnetic fields, the heat transport marginally increases with an increase of fringe-width.The convective motion gets confined near the sidewalls in regions of strong magnetic fields.
We numerically examined the effects of oscillating magnetic obstacle with different frequencies on liquid metal flow in a duct.We showed the presence of reduced and reversed streamwise velocity in the regions of strong magnetic field and in the wake of the magnetic obstacle.The regions of reduced velocity exhibit a wavy pattern with the wavelength of spatial oscillation decreasing with the excitation frequency of the magnets.The amplitude of wake oscillation shows a nonmonotonic dependence on the frequency of the magnets and exhibits a maximum at a particular frequency f max , which appears to correspond to the point of maximum Lorentz force in the streamwise direction and the minimum work done by the magnets on the fluid.The total streamwise and spanwise components of the Lorentz force oscillate sinusoidally with time with the same frequency and twice the frequency of the magnet's oscillation.The mean of the spanwise Lorentz force is zero.Its amplitude increases with an increase of the frequency of oscillation of the magnets.The frequency f max ≈ 0.5f s is considerably smaller than the reference value f s taken from Ref. [12].Although the stationary magnet does not produce vortex shedding in our case, it seems plausible that our value f max is indicative of the intrinsic shedding frequency when Re (and possibly Ha) are increased further.Lower values than St 0 = 0.25 of the Strouhal number of stationary magnetic obstacles were also reported by Kenjereš et al. [10].

Fig. 1 :
Fig. 1: Schematic of the setup for the duct flow with a localized oscillating magnetic field.The magnetic poles are semi-infinite in z-direction and oscillate along y-direction.The height Lz of the duct is unity.

Fig. 2 :
Fig. 2: (a) Schematic diagram of the magnetoconvection setup, and (b) isosurface contours of vertical velocity uz = 0.01 (red) and uz = −0.01(blue) for magnetoconvection with δ/H = 3 [4].The magnetic poles are semi-infinite in y and z-directions, and infinite along x-direction.The fluid motion in the region of strong magnetic fields is restricted to narrow zones adjacent to the sidewalls.