Liquid metal free surface dynamics in rotating permanent magnet stirrer

We study liquid metal stirring by rotating permanent magnets in a laboratory-scale rectangular glass container. The main goal is numerical model validation using experimental free surface shape data. We find reasonable agreement between the experiments and coupled liquid metal magnetohydrodynamics simulations. Since the surface tension forces are not dominant here, free surface is deformed mainly by the dynamic pressure of the bulk flow. Therefore, we can conclude that not only the free surface profile is similar to experiments, but the bulk flow must be also very similar.


Introduction
Liquid metal free surface dynamics is an important aspect in many applications, such as metal casting, stirring etc. Numerical modelling can provide very useful information regarding surface stability under different operating conditions of the system. However, numerical results usually require experimental validation.
Rotating permanent magnet-induced liquid metal flows have been studied in the past, however, acceptable agreement between experiments and simulations regarding flow velocity has not been obtained.
A popular method of liquid metal flow velocity measurements is the ultrasound Doppler velocimetry (UDV). Rotating permanent magnet-induced flow has been studied in [1], where reasonable quantitative agreement between UDV and simulations has been achieved, citing vortical instabilities as a possible reason for some qualitative differences. UDV measurements are known to be sensitive to presence of oxides on the container walls and probe [2], as well as temperature variations [3], so the results may not be consistent.
Another arguably very powerful tool of investigating liquid metal flows is the neutron radiography reported in [4], where the authors have used small tracer particles for neutron beam contrast. Assuming that the particles do not affect the flow directly, flow velocity distribution can be obtained from image sequences using particle image velocimetry, for example. There have been efforts of processing the obtained images and comparing the velocities to simulations, such as in [5], but the results, while qualitatively good, are quantitatively inconclusive, with the high concentration of particles in some areas as the main problem.
The goal of this study is to validate numerical results of rotating permanent magnet-induced liquid metal flow. To achieve this, we take a step back from multi-magnet systems with complicated flow patterns and create a very simple configuration -a single rotating magnet located next to a rectangular liquid metal container. In such a case, the flow pattern is basically a single vortex (apart from possible small secondary vortices which are usually unavoidable in turbulent flows). Rotation direction of the magnet is such that an upward liquid metal jet is induced at the wall closest to it which consequently deforms the free surface.

Experiments
Experimental setup is shown in Fig.1. All experiments are conducted at room temperature of approximately 20°C. The magnet is a nickel-plated NdFeB cylinder magnetized perpendicular to its axis (N52, remanence B r = 1.42 − 1.47T ) enclosed in a non-magnetic shell fixed to a shaft connected to a motor.
Magnet rotation frequency is varied using a frequency converter of the motor up to 50Hz in the present configuration. The frequency is set via a computer program.
Liquid metal is Galinstan with density 6440 kg/m 3 , viscosity 2.4 mPa·s, electrical conductivity 3.46 MS/m, surface tension coefficient 0.535 N/m and melting point 11°C. Initial experiments were performed with pure gallium (melting point 29°C), but it was difficult to keep it from solidifying during experiments. To prevent free surface oxidation and avoid complicated wetting phenomena, we pour HCl solution on the surface.
The main experimental measurement tool is video recording. Free surface shape is extracted from the images using a self-written python script. To achieve better contrast, we use a bright background behind the container, which then makes the liquid metal dark. Examples of raw as well as postprocessed images will be shown in the results section.

Numerical model
The present problem involves coupled electromagnetism and two-phase flow, in other words, we are dealing with liquid metal magnetohydrodynamics (MHD). Numerical solution of such problems is quite time-consuming task. Recently developed open-source solution EOF-Library [6][7] efficiently couples Elmer (finite element method) [8] and OpenFOAM (finite volume method) [9]. In this work, we use these programs and solve electromagnetics in Elmer and two-phase flow in OpenFOAM. Liquid metal shape data is transferred from OpenFOAM to Elmer and electromagnetic forces back from Elmer to OpenFOAM. All simulations are 3D. Scheme of the system (magnet and liquid metal volume) is shown in Fig.2. Distance from the magnet edge to liquid metal is fixed at 7 mm and the bottom of the magnet coincides with the bottom of liquid metal.  Large Re means that the flow is fully turbulent and turbulence models are needed. Small Re m means that the fluid flow velocity v doesn't affect the external magnetic field much, i.e. the magnetic field of the σ v × B current is negligible. This would mean that fluid flow and magnetic field can be decoupled and v can be neglected in the electromagnetics simulation. This does not mean that the v × B term can be dropped from the Lorentz force, though.
Large W e means that surface tension forces are much weaker than fluid inertia. This is very important here. It allows us to propose that the free surface deformation is caused and affected almost exclusively by the fluid flow or, in other words, the dynamic pressure p v = 0.5ρv 2 . This is the main point of simulation validation in this paper -if the simulated free surface profile is the same as in experiments, we can conclude that the flow velocity must also be the same.

Electromagnetics
Using potential formulation, i.e. writing magnetic flux density as B = ∇× A, where A -magnetic vector potential with ∇ · A = 0, relevant electromagnetic field equations in this problem are the following: where µ 0 -vacuum permeability, j -current density, σ -electrical conductivity, v -fluid velocity, M -magnetization. Equation (1) is Ampere's law and is the main equation solved to obtain magnetic field distribution. Equation (2) is the Ohm's law and Eq. (3) is the current continuity.
Electromagnetism is solved in frequency domain by expressing all oscillating fields using complex exponents, e.g. A = A 0 cos(ωt) = Re{ A 0C e iωt }, where subscripts 0 and C denote amplitude and a complex field, respectively, ω is angular frequency, t is time. In the complex form, the time-dependent exponents drop out of equations which are then solved for amplitudes. The magnet rotation is modelled as a rotating magnetization vector and in complex representation that means we must set amplitudes of real and imaginary parts of magnetization. For counter-clockwise rotation time-dependent form is M x = Br µ 0 cos (ωt) , M y = Br µ 0 sin (ωt) and in complex form this is M x,re = Br µ 0 , M y,im = Br µ 0 . Result of the electromagnetic simulation in Elmer is the complex amplitude of A from which current density j and magnetic field B can be calculated, which are then sent to OpenFOAM. The fluid dynamics simulation is transient and we convert the complex amplitudes back to time-dependent sinusoidal form.
Two important simplifications in electromagnetics model must be mentioned here. Firstly, the fluid velocity is neglected in the electromagnetic simulation, but the v × B term is considered in fluid flow simulation. This approximation is justified by the small magnetic Reynolds number. Secondly, electromagnetic simulation in the whole coupled simulation process is done only once at the initial stage with undeformed free surface. This is justified by noticing that the rotating magnet is relatively far from the free surface and so the induced forces for the deformed free surface shape would be almost the same. This reduces simulation time considerably.
Electromagnetics model consists of the whole system -liquid metal, magnet and air around them. Container is non-conductive and non-magnetic, so it is modelled as air. There is only one boundary condition -A = 0 on the external boundary of the domain, which means that the magnetic field is parallel to the boundary.

Fluid flow
The incompressible two-phase flow is governed by the Navier-Stokes equations with body forces: where ρ -fluid density, v -velocity, µ -viscosity, µ t -turbulent viscosity, p -pressure, f -body forces. To account for turbulence effects and calculate µ t , the k − ω SST [10] model is used.
To model the two immiscible phases (liquid metal and HCl solution (assumed to be pure water)), we use the Volume of Fluid method [11], implemented in the interFoam solver, which allows solving only one set of Navier-Stokes equations for both phases. This method introduces additional equation: with α = 0 for one fluid and α = 1 for the other. The interface between the fluids is usually assumed to be in numerical cells where α = 0.5. Density and viscosity are then where indices 1 and 2 designate fluid 1 and 2, respectively. Using VOF model, surface tension forces are implemented as a continuum surface force f s = γκ∇α (9) where κ = −∇ n -surface curvature, n -surface normal, γ -surface tension coefficient. Note that even though surface tension forces are not dominant in this particular case, we cannot neglect surface tension completely because we need to set contact angle different from 90 degrees and the surface must be rounded near the container walls as in experiments. As mentioned in the previous subsection, the σ v × B term is included in the fluid flow simulation. Since this is a part of electric current density, it must obey the continuity equation ∇ · j = ∇ · (σ v × B) = 0. Since such an equation may not hold as written, additional correction φ is introduced and equation is implemented in OpenFOAM to find φ. In the end, the total current density used for Lorentz force calculation is j = j elmer + σ v × B − ∇φ. Body force in equation (4) is then where j elmer represents the current density received from Elmer that does not include the v × B or ∇φ terms. Fluid flow model consists of only the liquid metal and some water above it, enough for free surface deformation. Boundary conditions are v = 0, ∂p ∂n = 0 on container walls and ∂ v ∂n = 0 (for outflow) or v = 0 (for inflow) and static pressure p s = 0 on the top boundary. Boundary condition for α is a fixed gradient corresponding to 160 degree contact angle (the value is estimated experimentally) on walls and α = 0 (for inflow) or ∂α ∂n (for outflow) on the top boundary.

Results
Let us first look at typical Lorentz force and velocity distributions in this system, shown in Fig.3. As expected from the magnet rotation, the force is directed mainly upwards and so is the flow in that zone, deforming the free surface.  1.7 mm at 20 Hz. Note that the vertical axis scales are different for each plot. The rectangular bump on the right side for experimental data is due to a black tape present in the background (see Fig.4). Similarities and differences between simulations and experiments are immediately clear. There is one pronounced maximum near the magnet, followed by a minimum and an increased height on the other end. A quantitative difference is that the maximum surface height is smaller in simulations than in experiments for all considered frequencies. This can simply be due to turbulence model -such two-parameter models tend to overestimate turbulent viscosity [12] leading to increased flow damping. Another reason could be that the distance between the magnet and melt is not matched absolutely exactly. Since the magnetic field decreases rapidly by increasing distance to the magnet and the force is proportional to the magnetic field squared, even a small distance error can lead to considerable induced force difference.
The main qualitative difference is the presence of a local maximum in the middle of container in simulations. This is due to unphysical water flow above the free surface. Test simulations with water flow damping (we multiply the velocity v with the volume fraction α which is 0 in water and 1 in liquid metal) confirms this -the surface is then flatter, more resembling the experimental one. However, such artificial second phase damping can lead to slight primary phase damping, so another solution to this problem must be sought.
As already mentioned above, the free surface oscillates. To demonstrate that, Fig.6 shows simulated free surface height approximately at the maximum position in time.
Two distinct oscillations can be distinguished here. The smallest one has frequency equal to double the magnet rotation frequency or equal to the Lorentz force frequency. The large oscillations are of hydrodynamic nature -surface waves and oscillations determined by material properties and container geometry.  Figure 6. Simulated free surface height oscillations at 10Hz magnet rotation frequency.

Conclusions
The simulated free surface shape is similar to the experimental one for different magnet rotation frequencies. As mentioned above, differences can arise because of turbulence model, some system parameter differences (distance from magnet to melt), as well as unphysical flow in the fluid above the liquid metal.
Further numerical work will focus on the solution to the unphysical secondary phase flow above the melt. Improved turbulence models, such as Large Eddy Simulation, which have been shown in literature to describe liquid metal free surface dynamics more realistically, will also be employed. Experimental improvements should include more precise control of crucial system parameters. Video camera with higher frame-rate will be used to capture the faster free surface oscillations.