# Modelling of coupled convective and magneto-thermoelectric phenomena for a sphere immersed in the liquid metal using finite element method and adaptive meshes

Barnabeu N., O Budenkova, A Gagnoud, S. Rukolaine, Y. Duterrail-Couvat

<figure><img src="image_1.png" /></figure>

## To cite this version:

Barnabeu N., O Budenkova, A Gagnoud, S. Rukolaine, Y. Duterrail-Couvat. Modelling of coupled convective and magneto-thermoelectric phenomena for a sphere immersed in the liquid metal using finite element method and adaptive meshes. Magnetohydrodynamics c/c of Magnitnaia Gidrodinamika, 2017, 53 (1), pp.149-158. 10.22364/mhd.53.1.16. hal-01897982

HAL Id: hal-01897982

https://hal.science/hal-01897982v1

Submitted on 18 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
The article as submitted to magnetohydrodynamics after the second revision

Modeling of coupled convective and magneto-thermoelectric
phenomena for a sphere immersed in the liquid metal using
finite element method and adaptive meshes

N. Bernabeu¹, O. Budenkova¹,³, A. Gagnoud¹, S. Rukolaine²,
Y. Du Terrail Couvat¹, Y. Fautrelle¹

¹ Univ. Grenoble Alpes, CNRS, Grenoble INP, SIMaP, F-38000 Grenoble, France

² Ioffe Institute, St Petersburg 194021, Russian Federation

³ corresponding author

**Abstract:** Thermoelectric phenomena which occur if an electrically conducting sphere is immersed in a liquid metal and subjected to a thermal gradient and an external static magnetic field are studied numerically and analytically for small Hartmann numbers. Such configuration can happen in solidification of metallic alloys under a static magnetic field, then the sphere corresponds to an equiaxed dendrite grain. Numerical model includes the Ohms law accounting for the thermoelectric and eddy current, stationary Navier-Stokes equations with the Lorentz force and diffusion-advection equation for the temperature. A finite element method with adaptive mesh is used to solve the problem. It is found that for the Hartmann number $Ha < 1$ analytical solutions for a decoupled problem and numerical solutions of a fully coupled one are similar whereas for $Ha = O(1)$ the thermoelectric current in the system is attenuated by the eddy current in the liquid. Consequently, with increase of the magnetic field the intensity of the thermoelectric convection remains bounded not only due to the braking effect of the magnetic field but also because of decrease of the thermoelectric force in the liquid.

**Introduction** Solidification of a liquid metal is usually performed with high thermal gradients that may lead to the appearance of the thermoelectric phenomena if the absolute thermoelectric power of the material is different for the liquid and the solid or varies with concentration [1]. If a static magnetic field **B** is imposed over such system, thermo-electromagnetic forces (TEM forces) may appear both in the solid and in the liquid, generate flow and alter the solidification structure as it was demonstrated for monocrystals [2], for columnar dendrites [3], and for equiaxed grains [4], see references within for other experimental works. Along with numerous experiments some analytical and theoretical studies [1],[5],[6] are far less consistent and may contain approximations which mask essential features of the problem. A better understanding of the fully coupled phenomena presented for simple geometries would be beneficial for the interpretation of the experimental results. The present study continues the work done in [7] and is focused at numerical simulations of the flow generated around a sphere immersed in a liquid metal and driven by a TEM force, accounting for the eddy current [1], [8] and diffusion-advection equation for the temperature. The use of adaptive meshes allowed us to consider a large calculation domain that almost eliminated effect of the walls on the flow field near the sphere.
**1. Statement of the problem and governing equations** A geometry illustrating the problems is shown in Fig. 1a. Calculation domain presents a cube $\Omega$ consisting of two domains, $\Omega_f$ and $\Omega_s$, occupied by a fluid and a solid sphere, respectively, in such a way that $\Omega = \Omega_s \cup \Omega_f$. Both phases are electroconductive. A sphere of a radius $R$ is placed at the cube centre and is significantly less than the cube. The $\partial\Omega$ denotes the boundary of the domain $\Omega$, the top, bottom and sides of the cube are denoted $\partial\Omega_{top}$, $\partial\Omega_{bottom}$, and $\partial\Omega_{side}$, respectively, and $\Gamma$ is used for the interface between $\Omega_s$ and $\Omega_f$. A temperature gradient $G = Ge_z$ along the $z$-axis is imposed over the $\Omega$. Similar to [7], we consider two cases for the directions of a uniform static magnetic field, in the first one it is parallel to the applied temperature gradient $B = Be_z$ and in the second it is orthogonal to $G: B = Be_x$. A spherical coordinate system with the origin at the sphere center used for presentation of results is shown in Fig. 1b, vectors $e_l$ used hereafter are unit vectors, either in spherical or in cartesian system related to the coordinate $l$.

<figure><img src="image_3.png" /><figcaption>**Figure 1:** (a) Domain of calculation $\Omega$ of size $L \times L \times L$ with a sphere of the radius $R$, ($L \gg R$), subjected to the temperature gradient $G = e_z(T_2 - T_1)/L$ and a uniform magnetic field, either parallel ($B = Be_z$) or perpendicular ($B = Be_x$) to $G$; (b) spherical coordinate system with the origin at the sphere center used in description of results.</figcaption></figure>

A coupled thermoelectric magnetohydrodynamic problem in a steady state is formulated as follows:

$$j = -\sigma_i \operatorname{grad} V - \sigma_i S_i \operatorname{grad} T + \sigma_i (\mathbf{u} \times \mathbf{B}) \text{ in } \Omega, \quad (1)$$

$$\operatorname{div} j = 0 \text{ in } \Omega, \quad (2)$$

$$\rho(\mathbf{u}, \operatorname{div} \mathbf{u}) - \operatorname{div}(2\mu D(\mathbf{u})) + \operatorname{grad} p = \rho \mathbf{g} + j \times \mathbf{B} \text{ in } \Omega_f, \quad (3)$$

$$\operatorname{div} \mathbf{u} = 0 \text{ in } \Omega_f, \text{ and } \mathbf{u} \equiv \mathbf{0} \text{ in } \Omega_s \quad (4)$$

$$\rho_i C_{p,i} \mathbf{u} \operatorname{div} T - \operatorname{div}(\lambda_i \operatorname{grad} T) = 0 \text{ in } \Omega \quad (5)$$

where $j$, $V$ and $T$ in Eqs. 1-2 are the density of the electric current, the electric potential and the temperature, respectively. In Eqs. 3-5, $\mathbf{u}$ is a velocity field in the liquid (zero in the solid), $p$ and $g$ are the pressure and the gravity vector and $D(\mathbf{u})$ denotes the symmetric part of the gradient of the velocity field. Notations and values for thermo-physical properties are given in the Table 1. In Eqs. 1-5 the index $i$ near the physical properties are introduced to stress that their values are different in the fluid ($i = f$) and in the solid ($i = s$).

At $\partial\Omega$ no-slip condition is set for velocity and no electric current through the $\partial\Omega$ is allowed. The $\partial\Omega_{side}$ is adiabatic and two temperatures $T_1 < T_2$ are fixed
at $\partial\Omega_{\text{bottom}}$ and $\partial\Omega_{\text{top}}$. A nominal temperature gradient $G$ and the size of the domain $L$ defined the required temperature difference $T_2 - T_1$. At the solid-liquid interface $\Gamma$, the no-slip condition for velocity and the continuity of heat fluxes and the density of the electric current should be satisfied that gives the following conditions:

$$ u|_{\Gamma} = 0 \qquad (6) $$

$$ (j_s \cdot n)_\Gamma = (j_f \cdot n)_\Gamma \qquad (7) $$

$$ (k_s \text{ grad } T_s \cdot n)_\Gamma = (k_f \text{ grad } T_f \cdot n)_\Gamma \qquad (8) $$

where $n$ is the outward unit normal to $\Gamma$ and $T_i$ and $j_i$ are the temperature and the density current in each medium, $i \in \{s, f\}$.

Our study showed that a special attention should be given to the mesh construction to catch the effect of the convection on the thermoelectric current and to avoid oscillations of $j$ at the solid/liquid interface. We found two empirical criteria, the former of which was based on the value of the Hartmann boundary layer $\delta_H$ related to the Hartmann number $Ha$ and the radius $R$ of a sphere [8]:

$$ \delta_H = \sqrt{\frac{\mu}{\sigma_f} \frac{1}{|\mathbf{B}|}} = \frac{R}{Ha} \qquad (9) $$

This criteria defined a characteristic size $l_1 = O(\delta_H/5)$ for mesh elements located within a spherical layer $r \in [R - 2\delta_H, R + 5\delta_H]$. However, for weak magnetic fields the length $\delta_H$ is not representative. Therefore, another characteristic size $l_2 = O(R/10)$ was defined to assure that the mesh follows the geometry of the particle near its surface. Finally, to construct the mesh around the particle, the minimal value $l = \min\{l_1, l_2\}$ for the size of cells was chosen with the gradual growing of mesh elements toward the walls of the calculation domain and the center of the sphere. To diminish the wall effect on solutions, the length of the cube edge was chosen to be $L = 50R$. The values for thermophysical properties were taken from [7], the values for $S_s$ and $S_f$ were taken in order to get their difference equal to the value used in [7]. Similarly, the nominal thermal gradient $G$ was taken equal to $3000K/m$ as a representative one, effect of this value is discussed below. The radius of the particle was taken equal to $50\mu m$. This choice is a compromise between a reasonable size of the particle and a sufficiently high value of the magnetic field for which calculations could be done with existing computing resource as also discussed below.

<table><caption>**Table 1:** Numerical values of physical properties for a Al-10wt%Cu alloy [7].</caption><thead><tr><th>Physical properties</th><th>value for liquid</th><th>value for solid</th><th>unit</th></tr></thead><tbody><tr><td>Dynamical viscosity μ</td><td>1.3 × 10<sup>-3</sup></td><td>-</td><td>Pa.s</td></tr><tr><td>Thermo-electric constant S<sub>i</sub></td><td>-3.8 × 10<sup>-6</sup></td><td>-5.4 × 10<sup>-6</sup></td><td>V.K<sup>-1</sup></td></tr><tr><td>Density ρ</td><td>2600</td><td>2600</td><td>kg.m<sup>-3</sup></td></tr><tr><td>Electrical conductivity σ<sub>i</sub></td><td>4 × 10<sup>6</sup></td><td>1 × 10<sup>7</sup></td><td>Ω<sup>-1</sup>.m<sup>-1</sup></td></tr><tr><td>Thermal conductivity λ<sub>i</sub></td><td>95</td><td>150</td><td>W.m<sup>-1</sup>.K<sup>-1</sup></td></tr><tr><td>Calorific capacity C<sub>p</sub></td><td>900</td><td>900</td><td>J.kg<sup>-1</sup>.K<sup>-1</sup></td></tr></tbody></table>

The model was realized with the open source C++ FEM library Reolef [9] and the mesh construction was performed with GMSH [10] which had generated tetrahedral cells of variable element size.
## 2. Numerical results compared with analytical solutions for a decoupled problem

If the therm $\sigma_f u \times B$ in Eq. 1 is neglected, the thermo-electrical and convective problems can be decoupled and solved analytically for the infinite domain of liquid. In this case, solution for the thermoelectric current in the solid and liquid phase is given as follows [11], [7]:

$$j_s = \frac{6\sigma_s\sigma_f}{\sigma_s + 2\sigma_f} \frac{\lambda_f}{\lambda_s + 2\lambda_f} (S_s - S_f)Ge_z \quad (10)$$

$$j_f = -\frac{9}{2} \frac{\sigma_s \sigma_f}{\sigma_s + 2\sigma_f} \frac{\lambda_f}{\lambda_s + 2\lambda_f} (S_s - S_f) G \frac{R^3}{r^3} (2 \cos \theta e_r + \sin \theta e_\theta) \quad (11)$$

where $e_r$ and $e_\theta$ are unit vectors in spherical coordinate system (Fig. 1b). Further, analytical solutions for the flow field arising due to the Lorentz force $\sigma_f j_f \times B$ can be obtained for weak flows where inertia terms can be neglected. In the next sections numerical solutions of a fully coupled problem given by Eqs.1-5 are compared with the analytical solutions for flow fields in a case of fully decoupled problem.

### 2.1. Magnetic field imposed along z-axis

For $B = Be_z$ in case of a decoupled problem only the azimuthal force acts in the liquid [7] that leads to formation of two axially symmetric azimuthal vortices in the upper ($z > 0$) and lower ($z < 0$) hemisphere, which rotate in opposite directions (Fig. 2a). Analytical solution in case of weak flow is presented in Appendix A and gives

$$U_{\varphi}(r, \theta) = 4\sqrt{12} U_{\varphi, \max} \left( \frac{R}{r} - \frac{R^3}{r^3} \right) \sin \theta \cos \theta \quad (12)$$

where $U_{\varphi,max}$ is an absolute maximal value of the velocity located at $r = \sqrt{3}R \approx 1.732R$ and angles $\theta = 45^\circ$ and $\theta = 135^\circ$ and equal to

$$U_{\varphi, \max} = \frac{\sqrt{12}}{8\mu} \frac{\sigma_s \sigma_f}{\sigma_s + 2\sigma_f} \frac{\lambda_f}{\lambda_s + 2\lambda_f} (S_s - S_f) GBR^2 \quad (13)$$

In Fig. 2b the module of the velocity $|U_\varphi(r, \theta)|$ given by Eq. 12 and the module of the velocity $|u(r)|$ obtained numerically by resolving the Eqs.1-5 for a particle with $R = 50 \ \mu m$ are plotted along the line $\theta = 45^\circ$ passing through the sphere center. A good agreement between the analytical and numerical solutions can be observed for magnetic fields $B_z = 0.01 \ T$ ($Ha \approx 0.03$) and $B_z = 0.1 \ T$ ($Ha \approx 0.3$) in a region near the particle ($r/R < 5$). For a larger distance from the particle the difference between the two solutions increases because of the finite size of the calculation domain. However, for $B_z = 0.5 \ T$ ($Ha \approx 1.4$) the maximal velocity obtained numerically is smaller than the one estimated analytically and location of this maximum is closer to the surface of the sphere than $\sqrt{3}R$.

The contours of the calculated value of velocity $|u|$ visualised in a vertical meridian plane for $B_z = 0.1 \ T$ and $B_z = 0.5 \ T$ are presented in Fig.3a,b. According to Fig.3, the magnetic field tends to "shrink" the azimuthal vortices radially and to stretch them along the direction of the imposed magnetic field, that was observed earlier in two-dimensional axisymmetric simulations [13]. The same tendency can be found in [1] where a similar coupled problem was considered analytically but without a sharp solid-liquid interface.
<figure><img src="image_4.png" /><figcaption>**Figure 2:** (a) Scheme for a flow around a sphere for $B = Be_z$; (b) Velocity profiles along the line $\theta = 45^\circ$ passing through the sphere center for $B = 0.01\,\text{T}$ ($Ha \approx 0.03$), $B = 0.1\,\text{T}$ ($Ha \approx 0.3$) and $B = 0.5\,\text{T}$ ($Ha \approx 1.4$): solid lines are for $|U_\varphi(r)|$ (Eq.12) and dashed, dotted and dashed-dotted lines are for numerical solution of Eqs.1-5 $|u|$ for the indicated values of $B$, respectively. Vertical dashed line at $r/R = \sqrt{3}$ shows the location of maximal velocities in analytical solution.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 3:** Contours for $|u|$ in a vertical meridional plane (a plane $(xz)$ was chosen) containing the sphere center for $B = 0.1\,\text{T}$ (a) and $B = 0.5\,\text{T}$ (b), velocity scale is the same for both figures.</figcaption></figure>

Numerical simulations showed as well that the density of the thermoelectric current $j$ remains uniform inside the particle but its value strongly decreases (Fig. 4) due to the term $\sigma_i u \times B$ in Eq.1. The latter acts in the liquid but affects the current in the solid through the continuity condition expressed by Eq. 7. Apart from already indicated [13] and [1], to our knowledge, no other article reported this phenomena. Therefore, increase of the magnetic field has a twofold effect in the system since the thermoelectric force acting in the liquid increases but the eddy current increases as well and tends to extinguishes the thermoelectric current. A certain balance can be expected between these phenomena. The
simulations showed that for the particle with the radius $R = 50$ µm the intensity of the convective flow seems to stabilise for $B \ge 2$ T ($Ha \ge 6.32$) around the value 0.8 mm/s.

<figure><img src="image_6.png" /><figcaption>**Figure 4:** Dependence of the maximal value of the density of the electric current $|j|$ (circle and dashed line) and maximal value of the velocity in the liquid $|u|$ on the value of the imposed magnetic field. Symbols (square for $|j|$ and circles for $|u|$) correspond to the values of $B$ for which the calculations were made.</figcaption></figure>

Because of a too big number of mesh cells required for obtaining of a non-oscillating solution it was not possible to perform calculations for higher values of the magnetic field for this size of the particle. The size of the particle does not affect the $j$, but, according to Eq.13 affects the maximal value of the velocity: a more intense flow can be expected for a larger particle. Consequently, a lower magnetic filed would lead to the stabilization of the velocity but in this case the number of required mesh cells grows significantly because of increase of the area of sphere surface. A lower magnetic filed at which the liquid velocity is stabilized can be expected as well for a higher thermal gradient. Because of vast variation of parameters used in simulations the comparison of the presented results with other reported elsewhere ([5], [15], [16]) is hardly possible.

## 2.2. Magnetic field imposed along x-axis
For the magnetic field $B = Be_x$ the resulting Lorentz force acting on the particle is not zero, but in the present study the particle is considered to be fixed. Analytical solution for a convective problem decoupled from the thermoelectric one for such configuration can be found in [14] for the problem of electromagnetophoresis. Because of the symmetry of the solution with respect to $x$ and $z$ axis, the eq.(23) from [14] is also a solution for the present problem and is written as

$$ \mathbf{u}_{an} = \frac{3B \sigma_s \sigma_f \lambda_f}{2\mu (\sigma_s + 2\sigma_f) (\lambda_s + 2\lambda_f)} (S_s - S_f) G \frac{R^3}{r^3} \left( \frac{R^2}{r^2} - 1 \right) \times [xye_x + (z^2 - x^2) e_y - yze_z] \quad (14) $$

Velocity flow filed for $B = 0.1$T ($Ha \approx 0.3$) and $B = 0.5$T ($Ha \approx 1.4$) are shown in Fig. 5 and Fig. 6, respectively. One can observe that the velocity in $(yz)$ plane is weaker than in $(xy)$ plane and that with increase of the magnetic field the difference between velocity in indicated planes increases.
<figure><img src="image_7.png" /><figcaption>**Figure 5:** Vector field of **u** in half-planes (xy) and (yz) and on a spherical surface with the radius 1.2R. Calculations are made for R = 50µm for B=0.1T ($Ha \approx 0.3$)</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 6:** Vector field of **u** in half-planes (xy) and (yz) and on a spherical surface with the radius 1.2R. Calculations are made for R = 50µm for B=0.5T ($Ha \approx 1.4$)</figcaption></figure>

In Figs. 5-6 the lines marked as 1 and 2 are the lines along the x-axis, i.e. (y = 0, z = 0) and along the z-axis, i.e. (x = 0, y = 0). According to Eq.14, velocity distribution along the lines 1 and 2 should be the same if the eddy current is not taken into account. However, in a coupled thermoelectric problem one can observe a difference in velocity distribution along these lines which increases with the value of the imposed magnetic field (Fig. 7). Similar to the previous case, analytical and numerical solutions are remarkably close for $B_x = 0.1$ T but differ significantly in terms of value of the velocity for a higher magnetic field, this difference is more important for the velocity distribution along z-axis shown with a dashed-dotted line. Note also that location of the maximal value of the velocity becomes closer to the sphere surface with increase of B. Similar to the previous case, magnetic field tends to "shrink" the flow field in a plane perpendicular to the direction of B, i.e. plane (yz) and to stretch it along the direction parallel to the magnetic field (x-axis).

Similar to the previous case, the eddy current which appeared in the liquid due to convection, affects the value of the thermoelectric current **j**. Examination of the dependence of **j** on the value of $B_x$ revealed the same behavior as shown in Fig. 4.
<figure><img src="image_9.png" /><figcaption>**Figure 7:** Profile of $|u_{an}|$ given by Eq. 14 (solid lines) and of calculated velocity $|u|$ along the line 1 (dashed) and 2 (dashed-dotted) shown in Figs. 5-6. Calculations are made for a sphere with $R = 50$ µm and different values of $B = Be_x$. Vertical dashed line at $r/R = \sqrt{3}$ shows the location of maximal velocities in analytical solution.</figcaption></figure>

**3. Conclusion** Numerical study performed in the present article demonstrated that the eddy current arising in the liquid attenuates the thermoelectric current generated in the system. Analytical solutions obtained for a decoupled problem where the term $\sigma_f u \times B$ in the Eq. 1 is neglected provided a good reference for results in case of $Ha < 1$ but overestimated significantly both the density of the thermoelectric current and intensity of convection already for $Ha = O(1)$.

No numerical simulations for magnetic filed higher than $2T$ for a particle with the radius $R = 50$µm were performed because the characteristic scale $\delta_H$ (Eq.9) became very small, that lead to severe conditions for the mesh construction within the Hartmann layer around the sphere and numerical solution required too much memory for computations.

A high sensitivity of the calculation to the mesh used in calculation was revealed and it can be stated that a special precaution in mesh construction should be made to avoid errors in modelling. According to our experience, the existence of such errors can be identified via oscillation or significant change of the values of the thermoelectric current $|j|$ near the solid-liquid interface.

**Acknowledgement** This work is supported by the Labex Tec21 (Investments d'Avenir -grant agreement ANR-11-LABX-0030).

**A. Analytical solutions for the velocity field** In spherical coordinate system the radial and polar components of the imposed field $B = Be_z$ are given as $B_r = B \cos\theta$ and $B_\theta = -B \sin\theta$, respectively and only the azimuthal force is acting in the liquid:

$$F_f = e_\varphi \frac{9}{2} \frac{\sigma_s \sigma_f}{\sigma_s + 2\sigma_f} \frac{\lambda_f}{\lambda_s + 2\lambda_f} (S_s - S_f) GB \frac{R^3}{r^3} \sin(2\theta) \quad (15)$$

For weak flows only the azimuthal component of the velocity $U_\varphi$ can be taken into account and described by the following equation in spherical coordinate system:

$$\frac{\partial}{\partial r} \left( r^2 \frac{\partial U_\varphi}{\partial r} \right) + \frac{\partial}{\partial \theta} \left( \frac{1}{\sin \theta} \frac{\partial (U_\varphi \sin \theta)}{\partial \theta} \right) + \frac{r^2}{\mu} F_{f,\varphi} = 0 \quad (16)$$
with boundary conditions

$$U_{\varphi(r,\theta)|r=R} = 0, \quad U_{\varphi(r,\theta)|r\to\infty} = 0 \qquad (17)$$

$$U_{\varphi(r,\theta)|\theta=0} = 0, \quad U_{\varphi(r,\theta)|\theta=\pi} = 0 \qquad (18)$$

Solution for Eq. 16 can be looked for in a form

$$U_{\varphi} = \sum_{n=1}^{\infty} u_{r,n}(r) P_{n}^{1}(\cos \theta) \qquad (19)$$

where $P_n^1(\cos \theta)$ are associated Legendre functions of the order 1 and degree $n$, and which satisfy the orthogonality condition with a weight ($\sin \theta$):

$$\int_{0}^{\pi} P_{m}^{1}(\cos \theta) P_{n}^{1}(\cos \theta) \sin \theta d\theta = \frac{2(n+1)}{(2n+1)(n-1)!} \delta_{m,n} \qquad (20)$$

where $\delta_{m,n}$ is the Kronecker delta. Note that the force $F_{f,\phi}$ given by Eq. 15 can be expressed in terms of $P_2^1(\cos \theta) = -3 \sin \theta \cos \theta$:

$$\frac{r^2}{\mu} F_{f,\varphi} = -\frac{3}{\mu} \frac{\sigma_s \sigma_f}{\sigma_s + 2\sigma_f} \frac{\lambda_f}{\lambda_s + 2\lambda_f} (S_s - S_f) GB \frac{R^3}{r} P_2^1(\cos \theta)$$

Eq. 16 is multiplied with a function $P_m^1(\cos \theta)$ weighted with the ($\sin \theta$). Using Eq. 20 and Eq. 18, obtain for $u_{r,n}(r)$:

$$\frac{2(n+1)!}{(2n+1)(n-1)!} \left[ \frac{d}{dr} \left( r^2 \frac{du_{r,n}}{dr} \right) - n(n+1)u_{r,n} \right] = \frac{1}{\mu} \frac{\sigma_s \sigma_f}{\sigma_s + 2\sigma_f} (S_s - S_f) \tilde{G} B \frac{R^3}{r} \frac{2(n+1)!}{(2n+1)(n-1)!} \delta_{2,n} \qquad (21)$$

with boundary conditions

$$u_{r,n}(r)|_{r=R} = 0, \quad u_{r,n}(r)|_{r\to\infty} = 0 \qquad (22)$$

In Eq.21, one can distinguish $n=1, n=2$ and $n \ge 3$ that gives:

$$\frac{d}{dr} \left( r^2 \frac{du_{r,1}}{dr} \right) - 2u_{r,1} = 0 \qquad (23)$$

$$\frac{d}{dr} \left( r^2 \frac{du_{r,2}}{dr} \right) - 6u_{r,2} = \frac{3}{\mu} \frac{\sigma_s \sigma_f}{\sigma_s + 2\sigma_f} \frac{\lambda_f}{\lambda_s + 2\lambda_f} (S_s - S_f) GB \frac{R^3}{r} \qquad (24)$$

$$\frac{d}{dr} \left( r^2 \frac{du_{r,l}}{dr} \right) - l(l+1)u_{r,l} = 0, \quad l \ge 3 \qquad (25)$$

Accounting for the boundary conditions (Eq. 22), one can see that Eq. 23 and Eq. 25 allows only trivial solutions $u_{r,1} \equiv 0$ and $u_{r,l} \equiv 0$ for $l \ge 3$. A general solution for a homogeneous equation corresponding to Eq. 24 is

$$u_{r,2,h} = \frac{C_1}{r^3} + C_2 r^2$$

As a particular solution of a non-homogeneous equation one can accept $u_{r,2,p} = C_3/r$ and putting it into equation, find

$$u_{r,2,p} = -\frac{1}{2\mu} \frac{\sigma_s \sigma_f}{\sigma_s + 2\sigma_f} \frac{\lambda_f}{\lambda_s + 2\lambda_f} (S_s - S_f) GB \frac{R^3}{r}$$
Then the general solution for $u_{r,2}$ is written as a sum $u_{r,2,h} + u_{r,2,p}$ where $C_1$ is found from the condition at the sphere surface while $C_2 \equiv 0$. That gives for $u_{r,2}$:

$$u_{r,2} = \frac{1}{2\mu} \frac{\sigma_s \sigma_f}{\sigma_s + 2\sigma_f} \frac{\lambda_f}{\lambda_s + 2\lambda_f} (S_s - S_f) GBR^2 \left( \frac{R^3}{r^3} - \frac{R}{r} \right)$$

Finally, solution has a form written in the article in Eqs. 12-13.

## REFERENCES

1. J.A. SHERCLIFF. Thermoelectric magnetohydrodynamics. *Journal of fluid mechanics*, vol. 91(2), p. 231-251, 1979

2. A.E. MIKELSON AND Y.K. KARKLIN. Control of crystallization processes by means of magnetic fields. *Journal of Crystal Growth*, vol. 52, p. 524-529, 1981

3. X. LI, Y. FAUTRELLE, Z.-M. REN, A. GAGNOUD, R. MOREAU, Y. ZHANG, C. ES-LING. Effect of a high magnetic field on the morphological instability and irregularity of the interface of a binary alloy during directional solidification. *Acta Materialia*, vol. 57(5), p. 1689-1701, 2009

4. G. SALLOUM-ABOU-JAOUDE, J. WANG, L. ABOU-KHALIL, G. REINHART, Z. REN, N. MANGELINCK-NOEL, X. LI, Y. FAUTRELLE, H. NGUYEN-THI. Motion of equiaxed grains during directional solidification under static magnetic field. *Journal of Crystal Growth*, vol. 417, p. 25-30, 2015

5. A. KAO. Analytic solutions to determine critical magnetic fields for thermoelectric magnetohydrodynamics in alloy solidification. *Metallurgical and Materials Transactions A*, vol. 46(9), p. 4215-4233, 2015

6. O. BUDENKOVA, N. BERNABEU, S. RUKOLAINE, Y DU TERRAIL COUVAT, A. GAGNOUD, R. TARPAGKOU, Y. FAUTRELLE. Thermoelectric problem for an axisymmetric ellipsoid particle in the liquid metal: analytical solution and numerical modeling. *Applied Mathematical Modelling*, vol. 45, p. 590-605, 2017

7. F. BALTARETU, J. WANG, S. LETOUT, Z.-M. REN, X. LI, O. BUDENKOVA AND Y. FAUTRELLE. Thermoelectric effects on electrically conducting particles in liquid metal. *Magnetohydrodynamics*, vol. 51(1), 2015.

8. P. A. DAVIDSON. *An introduction to magnetohydrodynamics* (Cambridge university press, 2001)

9. P. SARAMITO. *Efficient C++ finite element computing with Rheolef* (CNRS-CCSD ed., 2016).

10. C. GEUZAINE, J.-F. REMACLE. Gmsh: a three-dimensional finite element mesh generator with built-in pre- and post-processing facilities. *International Journal for Numerical Methods in Engineering*, vol. 79(11), p. 1309-1331, 2009

11. P. B. NAGY, A. H. HAYFEH. On the thermoelectric magnetic field of spherical and cylindrical inclusions. *Journal of Applied Physics*, vol. 87, p. 7481-7490, 2000.

12. G. A. KORN AND T. M. KORN. *Mathematical Handbook for Scientists and Engineers* (McGraw-Hill, New York, 1968)

13. Y. DU TERRAIL COUVAT, A. GAGNOUD, D. BRASILIANO, O. BUDENKOVA, AND Y. FAUTRELLE. Numerical modelling of thermoelectric magnetic effects in solidification. *Proceeding 8th International Conference on Electromagnetic Processing of Materials EPM2015* Oct 2015, Cannes, France.

14. D. LEENOV AND A. KOLIN. Theory of electromagnetophoresis. I. Magnetohydrodynamic forces experienced by spherical and symmetrically oriented cylindrical particles. *The Journal of Chemical Physics*, vol. 22(4), p. 683-688, 1954

15. I. KALDRE, Y. FAUTRELLE, J. ETAY, A. BOJAREVICS, L. BULIGINS. Investigation of liquid phase motion generated by the thermoelectric current and magnetic field interaction. *Magnetohydrodynamics*, vol. 46(4), p. 371-380, 2010

16. JIANG WANG. *PhD Thesis of University of Grenoble and University of Shanghai*, 2013