Three-dimensional simulation of film boiling on a horizontal surface with magnetic field

Abstract

This study conducts a numerical investigation into the three-dimensional film boiling of liquid under the influence of external magnetic fields. The numerical method incorporates a sharp phase-change model based on the volume-of-fluid approach to track the liquid-vapor interface. Additionally, a consistent and conservative scheme is employed to calculate the induced current densities and electromagnetic forces. We investigate the magnetohydrodynamic effects on film boiling, particularly examining the pattern transition of the vapor bubble and the evolution of heat transfer characteristics, exposed to either a vertical or horizontal magnetic field. In single-mode scenarios, film boiling under a vertical magnetic field displays an isotropic flow structure, forming a columnar vapor jet at higher magnetic field intensities. In contrast, horizontal magnetic fields result in anisotropic flow, creating a two-dimensional vapor sheet as the magnetic strength increases. In multi-mode scenarios, the patterns observed in single-mode film boiling persist, with the interaction of vapor bubbles introducing additional complexity to the magnetohydrodynamic flow. More importantly, our comprehensive analysis reveals how and why distinct boiling effects are generated by various orientations of magnetic fields, which induce directional electromagnetic forces to suppress flow vortices within the cross-sectional plane.

1 Introduction

Boiling is widely employed in industrial heat-related processes due to its exceptional heat transfer efficiency (Dhir, 1998; Tomar et al., 2008; Guion et al., 2018; Lavino et al., 2021). Depending on the overheat level, boiling can be categorized into three stages: nucleate boiling, transition boiling, and film boiling. The phenomenon of film boiling is especially common under conditions of significant overheat and is a crucial process in industrial applications. In the context of magnetic confinement fusion, a method used for controllable nuclear fusion, liquid metal lithium can serve as a heat remover within a liquid blanket, facilitating the energy conversion process. Film boiling of this liquid metal occurs under conditions of extremely high heat flux and overheat generated by nuclear fusion. The presence of a magnetic field (MF) induces magnetohydrodynamics (MHD) effects that influence the entire process. Considering this perspective, investigating the film boiling of liquid metal under the influence of a magnetic field emerges as a noteworthy research topic.

The initial studies of film boiling on a horizontal surface were predominantly based on theoretical and experimental investigations. Pioneering work by Chang (1957, 1959) involved analyzing the mechanism of film boiling on a horizontal surface using the theory proposed by Taylor (1950) for the most unstable wavelength. Zuber (1958) extended this theory to derive a model to predict the minimal heat flux during film boiling and demonstrated that Taylor instability predominantly governed the interface with the wavelength satisfying $\lambda_{c}\leq\lambda\leq\lambda_{d}$, where $\lambda_{c}(=2\pi\sqrt{\sigma/(\rho_{l}-\rho_{v})g})$ represents the most critical wavelength, and $\lambda_{d}(=\sqrt{3}\lambda_{c})$ denotes the most dangerous wavelength. Subsequently, Klimenko (1981) and Klimenko & Shelepen (1982) introduced a new correlation with an accuracy of approximately $\pm 25$%. This theoretical model has been widely acknowledged as a foundational basis for subsequent research on film boiling.

The emergence of powerful supercomputers, advancements in numerical simulation techniques, and the development of precise interface tracking methods have steered contemporary research on film boiling over a horizontal surface towards direct numerical simulations. However, many numerical simulations of film boiling have still focused on two-dimensional (2D) simulations and single-mode models, with the computational domain width close to the most dangerous wavelength $\lambda_{d}$. These 2D single-mode film boiling numerical simulations were mainly used to investigate the influence of different physical quantities (e.g., Jakob number, Morton number, density ratio) on morphology and heat transfer in film boiling. For high Jakob numbers, the rapid formation of vapor makes bubble detachment more difficult, resulting in the formation of tall and stable vapor columns (Son & Dhir, 1997; Juric & Tryggvason, 1998; Guo et al., 2011). For low Morton numbers or high-density ratios, it becomes more challenging for bubbles to pinch off (Juric & Tryggvason, 1998; Akhtar & Kleis, 2013; Mortazavi, 2022). Some studies (Agarwal et al., 2004; Tomar et al., 2008; Guo et al., 2011) suggested the periodic nature of bubble growth at nodes and anti-nodes in accordance with wavelength estimated based on Taylor’s theory. Some other researchers have performed multi-mode film boiling simulations, considering factors like the interaction between bubbles and the initial interface. A few studies (Esmaeeli & Tryggvason, 2004b; Hens et al., 2014; Ningegowda & Premachandran, 2019) have explored the impact of the Jakob number on multi-mode film boiling, revealing that, at high Jakob numbers, boiling occurs in the form of vapor columns. The proximity between vapor columns induces mutual influence and potential merging due to the instability of these columns, ultimately resulting in a chaotic process. Several researchers have highlighted the shortcomings in 2D film boiling simulations (Gibou et al., 2007; Akhtar & Kleis, 2013; Zhang & Ni, 2018). These studies observed that bubbles did not detach from the vapor film in their 2D simulations. This issue was attributed to the curvature of the stem connecting the thin film and the bubble approaching zero, thus neglecting the surface tension effect.

Earlier studies on film boiling with large 2D domains have consistently confirmed that the distribution of bubbles corresponds to the most dangerous wavelength in film boiling on a horizontal surface (Esmaeeli & Tryggvason, 2004b; Tomar et al., 2008; Hens et al., 2014; Akhtar, 2018; Ningegowda & Premachandran, 2019). The three-dimensional (3D) numerical simulations of film boiling employing simplified liquid/vapor parameter models primarily focus on single-mode scenarios (Esmaeeli & Tryggvason, 2004a; Tsui & Lin, 2016; Zhang & Ni, 2018; Kumar & Premachandran, 2022). Moreover, limited research has been devoted to exploring the underlying physical mechanisms, with much of the existing work replicating earlier 2D studies. These investigations often revisit aspects like the impact of overheat, Grashof number, and density ratio (Khorram & Mortazavi, 2022). Furthermore, the scope of single-mode film boiling simulations is constrained by computational domain limitations, preventing the complete visualization of bubble generation throughout the boiling process. Notably, three-dimensional (3D) multi-modal film boiling simulations in existing literature have typically employed computational domains of sizes not exceeding 2$\lambda_{d}$ (Esmaeeli & Tryggvason, 2004b; Khorram & Mortazavi, 2022; Shin & Juric, 2002), thus inadequately capturing the actual distribution of bubbles.

From the above discussion, it is evident that research on conventional film boiling has indeed reached a profound level of depth. However, there are still areas where further exploration is warranted, particularly in understanding the underlying physical mechanisms, expanding beyond single-mode simulations, and overcoming computational limitations to capture the full dynamics of multi-modal film boiling. Film boiling coupled with multiphysics, particularly electrohydrodynamics (EHD), has recently emerged as the forefront of film boiling investigation. Extensive experimental studies (Ogata & Yabe, 1993; Wang et al., 2009; Pearson & Seyed-Yagoobi, 2013; Patel et al., 2016; Patel & Seyed-Yagoobi, 2017) and numerical simulations (Tomar et al., 2007, 2009; Pandey et al., 2015, 2016, 2017; Rouzbahani & Mortazavi, 2023; Wang et al., 2023) have shown that electric fields aid in accelerating bubble detachment during the boiling process, thereby enhancing heat transfer efficiency. Furthermore, numerical simulations focusing on 2D multi-model film boiling under electric fields have found that the electric field promotes the occurrence of Rayleigh-Taylor instability at the vapor-liquid interface and increases the number of bubble nucleation sites. However, as another crucial physical phenomenon, magnetic fields offer considerable research potential, particularly in confined nuclear fusion. The effect of MF on fluids can be divided into magnetohydrodynamics (MHD) and ferrohydrodynamics (FHD), where MHD effects result from the action of electromagnetic forces or Lorentz forces, and FHD effects result from Kelvin forces. Regarding the study of MF influence on boiling heat transfer, both experimental and numerical research have focused more extensively on FHD. This is mainly because ferrofluids, compared to magnetic fluids such as liquid lithium, have decent transparency, which brings convenience to experimental investigations. A plethora of experimental (Lee et al., 2012; Shojaeian et al., 2016; Abdollahi et al., 2017; Zonouzi et al., 2019) and numerical (Aminfar et al., 2012, 2014; Malvandi, 2016; Guo et al., 2019) studies have demonstrated that MFs enhance the heat transfer characteristics of ferrofluid boiling, while also increasing the number of bubble nucleation sites during the boiling process.

In contrast, research on magnetohydrodynamic (MHD) effects on film boiling, even today, primarily relies on early theoretical and experimental studies. Faber Jr & Hsu (1967) investigated the effect of a vertical MF on the subcooled boiling of mercury and found that a strong MF suppresses natural convection but promotes nucleate boiling. Fraas et al. (1974) studied the effect of a vertical MF on boiling in potassium salts and metals and concluded that the MF had an insignificant effect. Wagner & Lykoudis (1981) studied the influence of a horizontal magnetic field on nucleate boiling of liquid metallic mercury, finding that the horizontal magnetic field reduced the heat transfer efficiency of boiling and promoted the transition from nucleate to film boiling. Takahashi et al. (1994) examined the impact of a vertical MF on saturated nucleate boiling of mercury on a horizontal surface. It was shown that there was a decrease in the Nusselt number and heat transfer efficiency with increasing MF. Experiments conducted by de Bertodano et al. (1998) showed that the magnetic field slowed down the frequency of bubble release during nucleate boiling. The theoretical research conducted by Arias (2010c, a) corroborated the findings from the aforementioned experimental investigations. However, theoretical investigation of Arias (2010b) within the Taylor–Helmholtz instability framework indicated that the MF did not influence film boiling. The aforementioned literature review suggests that the impact of magnetic fields on boiling does not yield consistent results, and some findings are even contradictory. Additionally, due to the opacity of liquid metals, experimental studies were unable to capture the dynamics of bubbles during the boiling process accurately. Consequently, numerical simulation becomes the only feasible method to explore the mechanisms underlying the influence of magnetic fields on film boiling with liquid metals, a subject that, to our knowledge, has not been investigated previously.

Thus, the present study aims to comprehend the behavior of film boiling under the influence of applied magnetic fields in different directions. The scarcity of research in this area can be attributed to the increased significance of magnetohydrodynamic (MHD) effects in 3D models. Additionally, the flow and thermal boundary layers of liquid metals are extremely thin, necessitating grids with fine resolution and consequently leading to substantial computational demands. Furthermore, numerical simulations involving incompressible multiphase magnetohydrodynamics with phase change impose high requirements on the robustness of the employed numerical solver. The phase change model utilized in this study was developed in our earlier work (Zhao et al., 2022), which employed a sharp interface method to enforce the temperature and concentration conditions at the liquid/vapor interface, integrating a volume-of-fluid method and an embedded boundary approach for sharp scheme construction. Zhao et al. (2022) extensively validated the phase change model against various benchmark problems and rigorous numerical tests, showcasing its exceptional accuracy and robustness. In addition to this phase change model, we integrate a consistent and conservative scheme (Zhang & Ni, 2014b) for discretizing the magnetohydrodynamic (MHD) equations, enabling the simulation of 3D film boiling on a horizontal surface in the presence of an applied magnetic field. Our investigation comprises two primary components. Firstly, we focus on single-mode film boiling under vertical and horizontal magnetic fields in three dimensions. In contrast to previous studies, we significantly extend the simulation duration to observe the impact of the magnetic field after the complete development of the boiling flow. Additionally, we increase the computational domain height to explore a wider range of boiling dynamics. Secondly, we extend the computational domain size to $4\lambda_{d}\times 4\lambda_{d}\times 4\lambda_{d}$ and analyze the behavior of multi-mode film boiling under both vertical and horizontal magnetic fields.

The rest of the manuscript is structured as follows. In §2, we present the physical model and governing equations and introduce the phase change model developed by Zhao et al. (2022). The results are discussed in §3, where we analyze the influence of horizontal and vertical magnetic fields on both single-mode and multi-mode film boiling phenomena in §3.1 and §3.2, respectively. Finally, concluding remarks are provided in §4.

2 Problem formulation and numerical method

2.1 Physical model

We examine the film boiling phenomenon of liquid under horizontally and vertically applied MFs by performing 3D numerical simulations. Firstly, based on previous studies on film boiling without an imposed MF, a brief description of the entire physical process is provided: during film boiling, there exists a layer of vapor film between the overheated wall and the upper liquid layer, preventing direct contact between the upper fluid and the overheated wall. As the evaporation proceeds, the vapor film gradually thickens, leading to instability at the interface. Eventually, the influence of gravity triggers the Rayleigh-Taylor instability, giving rise to the formation of bubbles. According to this physical process, the computational model set up for this study is illustrated in Fig. 1. The size of the computational domain is $H_{x}\times W_{y}\times W_{z}$, where $H_{x},W_{y}$, and $W_{z}$ are multiples of $\lambda_{d}$. This choice aligns with the observation by Tsui & Lin (2016) that setting the computational domain width as $\lambda_{d}$ in 3D simulations adequately captures the physical phenomenon. A Cartesian coordinate system $(x,y,z)$ is employed, with $x$ representing the vertically upward direction and $y$ and $z$ indicating the horizontal and spanwise directions, respectively. The center of the coordinate system is aligned with the center of the bottom plane $(yz)$. The gravity $\boldsymbol{g}$ acts in the negative $x$ direction. An external uniform magnetic field is applied either in the vertical ($x$) or horizontal ($y$) direction.

We employ the following boundary conditions: no-slip and no-through-flow conditions at the bottom (overheated surface), periodic boundary conditions at the four vertical sides (front/back/left/right), and outflow conditions at the top. The phase interface and the liquid are maintained at the saturation temperature ($T_{\text{sat}}$), which is a valid assumption except when exposed to extremely high heat flux (Juric & Tryggvason, 1998; Akhtar & Kleis, 2013). Additionally, the temperature of the overheated surface is kept constant at $T_{\text{wall}}$ ($T_{\text{wall}}=T_{\text{sat}}+\Delta T$). Initially (at time $t=0$), the interface is defined using sinusoidal functions, and the temperature in the vapor phase linearly decreases from the overheated surface to the phase interface.

\begin{overpic}[scale={0.45},tics=5]{figure/physical-model.jpg} \put(15.0,82.0){$\boldsymbol{g}$} \put(35.0,44.0){\begin{turn}{50.0}$W_{y}$\end{turn}} \put(54.0,48.0){\begin{turn}{-15.0}$W_{z}$\end{turn}} \put(25.0,45.0){$H_{x}$} \put(45.0,82.0){Outflow} \put(38.0,3.0){Overheated surface} \put(73.0,27.0){Phase interface} \put(71.0,12.0){Horizontal MF} \put(73.0,68.0){Vertical MF} \put(15.0,32.0){\begin{turn}{90.0}Periodic boundaries\end{turn}} \put(85.0,32.0){\begin{turn}{90.0}Periodic boundaries\end{turn}} \end{overpic}

Figure 1: Schematic diagram of the film boiling configuration under the applied MFs. The boundary conditions employed along the boundaries of the computational domain are also specified.

2.2 Governing equations

The viscous incompressible MHD two-phase flow, incorporating liquid-vapor phase change, is governed by the mass, momentum, and energy conservation equations. They are expressed as follows:

$$\nabla\cdot\boldsymbol{u}=\ddot{m}\left(\frac{1}{\rho_{v}}-\frac{1}{\rho_{l}}\right),$$ (1)

$$\rho\left(\frac{\partial\boldsymbol{u}}{\partial t}+\boldsymbol{u}\cdot\nabla\boldsymbol{u}\right)=-\nabla p+\nabla\cdot\left(\mu\left(\nabla\boldsymbol{u}+\nabla\boldsymbol{u}^{T}\right)\right)+\rho\boldsymbol{g}+\sigma\kappa\delta_{s}\boldsymbol{n}+\boldsymbol{F_{l}},$$ (2)

$$\frac{\partial T}{\partial t}+\boldsymbol{u}\cdot\nabla T=\nabla\cdot\left(\frac{k}{\rho c_{p}}\nabla T\right)+\Phi.$$ (3)

Here, $\boldsymbol{u}$ signifies the fluid velocity, $\ddot{m}$ represents the volume mass transfer rate, $\rho$ denotes the density, and subscripts $l$ and $v$ designate the physical properties for the liquid and vapor phases, respectively. In Eq. (2), $p$ stands for pressure, $\mu$ is the dynamic viscosity, and $\sigma\kappa\delta_{s}\boldsymbol{n}$ represents the surface tension force, wherein, $\sigma$ is the surface tension coefficient, assumed to be a temperature-independent constant. The interface curvature is denoted by $\kappa$, $\delta_{s}$ is the Dirac distribution function (whose value is equal to 1 at the interface and zero elsewhere), and $\boldsymbol{n}$ denotes the unit normal vector at the local interface.

Besides, $\boldsymbol{F_{l}}$ is the Lorentz force arising from the electromagnetic induction, which is given by

$$\boldsymbol{F_{l}}=\boldsymbol{J}\times\boldsymbol{B},$$ (4)

where $\boldsymbol{J}$ is the induced current density, $\boldsymbol{B}$ represents the applied MF. According to the Ohm’s law, the current density is given by

$$\boldsymbol{J}=\sigma_{e}(-\nabla\varphi+\boldsymbol{u}\times\boldsymbol{B}),$$ (5)

where $\sigma_{e}$ denotes the electric conductivity and $\varphi$ represents the induced electric potential. In Eq. (5), the induced magnetic intensity is neglected because it is significantly smaller than the external field, acknowledging that the magnetic Reynolds number ($Re_{m}=\eta\sigma_{e}Lu^{\prime}$, with $\eta$ the magnetic permeability, $L$ the characteristic length and $u^{\prime}$ the characteristic velocity) is much smaller than 1 for the MHD flow considered in the present study. In order to satisfy the charge conservation, we employ

$$\nabla\cdot\boldsymbol{J}=0.$$ (6)

By combining Eqs. (5) and (6), the electric potential Poisson equation is derived as:

$$\nabla\cdot\left(\sigma_{e}\nabla\varphi\right)=\nabla\cdot\left(\sigma_{e}\boldsymbol{u}\times\boldsymbol{B}\right).$$ (7)

It is to be noted that solving this Poisson equation allows for the determination of the electric potential $\varphi$. Subsequently, the current density $\boldsymbol{J}$ can be calculated using Eq. (5). Finally, the Lorentz force $\boldsymbol{F_{l}}$ can be determined. In Eq. (3) governing the temperature transport, $c_{p}$ represents the volumetric heat capacity, $k$ is the thermal conductivity, $T$ denotes the temperature, and $\Phi$ stands for the internal dissipation term, encompassing viscous dissipation and Joule dissipation. However, in this study, both dissipation effects are considered negligible compared to the temperature terms, a simplification consistent with other studies (Burr & Müller, 2002; Vogt et al., 2018; Chen et al., 2024).

The presence of phase change necessitates that the above equations must satisfy the following interfacial jump conditions:

$$\dot{m}=\rho_{l}\left(\boldsymbol{u}_{\Gamma}-\boldsymbol{u}_{l}\right)\cdot\boldsymbol{n}=\rho_{v}\left(\boldsymbol{u}_{\Gamma}-\boldsymbol{u}_{v}\right)\cdot\boldsymbol{n},$$ (8)

$$\dot{m}=\frac{\left(\boldsymbol{q}_{l}-\boldsymbol{q}_{v}\right)\cdot\boldsymbol{n}}{h_{lg}}=\frac{\left((k\nabla T)_{l}-(k\nabla T)_{\text{v}}\right)\cdot\boldsymbol{n}}{h_{lg}}.$$ (9)

Eq. (8) describes the mass transfer rate from liquid to vapor in accordance with the mass conservation law. Note that $\boldsymbol{u}_{\Gamma}$ is the propagation velocity of the interface, which does not equal to the local fluid velocity $\boldsymbol{u}_{(l,v)}$ in the presence of phase change. Eq. (9) elucidates the source of the mass transfer rate generated by the thermal flux across the interface. Here, $\dot{m}$ denotes the area mass transfer rate, $\boldsymbol{q}$ is the thermal flux at the respective side of liquid and vapor, and $h_{lg}$ is the latent heat. Then the velocity jump condition (Eq. 8) leads to the following divergence-constraint condition on the velocity field

$$\nabla\cdot\boldsymbol{u}=(\frac{1}{\rho_{v}}-\frac{1}{\rho_{l}})\ddot{m}.$$ (10)

The relationship between $\dot{m}$ and $\ddot{m}$ is given by

$$V\ddot{m}=\dot{m}S_{\Gamma},$$ (11)

where $V$ denotes the volume of the discretized cell containing the interface and $S_{\Gamma}$ is the area of that interface. It is to be noted that both $\ddot{m}$ and $\dot{m}$ are exclusively estimated at the interface. In this work, we assume the interface temperature remains constant at the saturation temperature:

$$T_{\Gamma}=T_{sat},$$ (12)

which effectively serves as a Dirichlet-like boundary condition for the solution of the temperature equation in the separated liquid and vapor phases.

2.3 Numerical method

In the current study, our MHD-phase change solver is implemented within the framework of the open-source platform Basilisk (Popinet, 2009, 2015), well-regarded for its capability in simulating complex multiphase flows. A staggered-in-time approximate projection method is employed to discretize the incompressible Navier–Stokes equations, while the volume-of-fluid (VOF) method (Weymouth & Yue, 2010; Scardovelli & Zaleski, 1999) accurately captures and advances the phase interface. The Bell-Colella-Glaz (BCG) second-order scheme (Bell et al., 1989) discretizes the advection term, and a fully implicit scheme is applied to discretize the diffusion term. Additionally, a quad/octree-based adaptive mesh refinement (AMR) technique is utilized for spatial discretization, enabling grid refinement primarily near the interface to optimize computational resources.

In our previous study (Zhao et al., 2022), we developed a numerical method for the phase change model capable of simulating boiling and evaporating flows. This method imposes jump conditions at the interface using a sharp scheme while exactly conserving both mass and energy. The relative conservation error is found to be consistently smaller than $10^{-3}$. In the present study, we adopt this numerical method to address the phase-change model for film boiling. Additionally, in the presence of an external magnetic field (MF), we utilize a consistent and conservative scheme, which was also developed earlier (Zhang & Ni, 2014b, a), for computing electric currents and Lorentz forces. Detailed information about the numerical technique and validations regarding the MHD-phase change solver can be found in previous works (Zhao et al., 2022; Zhang & Ni, 2014b, a). In view of this, we will only present some crucial information featuring the progression of the numerical solution from time level $n$ to level $n+1$:

1. A split advection method proposed by Weymouth & Yue (2010) is employed to discretize the VOF advection equation, and the volume fraction $C^{n+\frac{1}{2}}$ is updated accordingly. However, the presence of phase change alters the form of the VOF advection equation to

   $$\frac{\partial C}{\partial t}+\nabla\cdot\left(\boldsymbol{u}C\right)=-\frac{\ddot{m}}{\rho_{l}},$$ (13)

   which is consistent with the formulation of the mass conservation equation. Due to the discontinuity in velocity at the interface, it is not possible to solve the VOF advection equation directly, as is the case with ordinary two-phase flows. To address this, an extended velocity field $\boldsymbol{u}_{E}$ is introduced, $\boldsymbol{u}_{E}$=$\boldsymbol{u}$-$\boldsymbol{u}_{J}$. $\boldsymbol{u}_{J}$ is the velocity jump, which could be obtained by solving for the velocity potential $\psi$ using the following equations:

   $$\left\{\begin{array}[]{l}\nabla^{2}\psi=\left(\frac{1}{\rho_{v}}-\frac{1}{\rho_{l}}\right)\ddot{m},\\ \boldsymbol{u}_{J}=\nabla\psi.\end{array}\right.$$ (14)

   To avoid the negative effect of the source term, the interfacial velocity ($\boldsymbol{u}_{\Gamma}$) is chosen to advect the VOF advection equation:

   $$\frac{\partial C}{\partial t}+\boldsymbol{u}_{\Gamma}\cdot\nabla C=0,$$ (15)

   where $\boldsymbol{u}_{\Gamma}$=$\boldsymbol{u}_{E}$- $\frac{\ddot{m}}{\rho_{l}}\boldsymbol{n}$, wherein $\boldsymbol{n}$ directs towards the vapor phase in normal direction at the interface.

2. The physical properties $\phi^{n+\frac{1}{2}}$ are updated using

   $$\phi^{n+\frac{1}{2}}=C^{n+\frac{1}{2}}\phi_{l}+\left(1-C^{n+\frac{1}{2}}\right)\phi_{v},\quad{\rm where}~{}~{}\phi=\rho,~{}\mu,~{}\sigma_{e}.$$ (16)

3. The temperature $T^{n+\frac{1}{2}}$ is then updated. Note that the temperature fields of the liquid and vapor phases are solved separately to ensure sharp temperature gradients at the interface for accurate heat flux calculations. Specifically, the temperature equation is solved only in the vapor phase for saturated film boiling, maintaining the liquid temperature at $T_{\Gamma}$. The interfacial temperature ($T_{\Gamma}$) is imposed as the interfacial boundary condition, which is realized via an embedded boundary method. Besides, the advection term of the temperature is discretized by using a geometrical scheme, which is consistent with that applied in Eq. (15).

   $$\frac{T^{n+\frac{1}{2}}-T^{n-\frac{1}{2}}}{\Delta t}+\nabla\cdot(T\boldsymbol{u})^{n+\frac{1}{2}}=\nabla\cdot\left(\left(\frac{\lambda}{\rho c_{p}}\right)\nabla_{f}T\right)^{n+\frac{1}{2}}+\left(\frac{\ddot{m}}{\rho}T_{\Gamma}\right)^{n+\frac{1}{2}}.$$ (17)

4. The Lorentz force $\boldsymbol{F_{l}}^{n+\frac{1}{2}}$ is updated by solving the electric potential Poisson equation and Ohm’s law using

   $$\nabla\cdot\left(\sigma_{e}\nabla\varphi\right)^{n+\frac{1}{2}}=\nabla\cdot\left(\sigma_{e}(\boldsymbol{u}\times\boldsymbol{B})\right)^{n+\frac{1}{2}},$$ (18)

   $$\boldsymbol{J}^{n+\frac{1}{2}}=\sigma_{e}^{n+\frac{1}{2}}\left(-\nabla\varphi+\boldsymbol{u}\times\boldsymbol{B}\right)^{n+\frac{1}{2}},$$ (19)

   $$\boldsymbol{F_{l}}^{n+\frac{1}{2}}=\boldsymbol{J}^{n+\frac{1}{2}}\times\boldsymbol{B}^{n+\frac{1}{2}}.$$ (20)

5. The mass transfer rate $\dot{m}^{n+\frac{1}{2}}$ is updated according to Eq. (9), while the embedded boundary method is used to calculate the temperature gradient, ensuring that the interpolation template only uses those cell values located in the vapor phase. This method avoids errors caused by artificially thickening the interface and ensuring accurate mass transfer rate calculations.

6. Finally, the pressure field ${p}^{n+1}$ and the velocity field $\boldsymbol{u}^{n+1}$ are updated using the approximate projection method, which has been incorporated in Basilisk (Popinet, 2009, 2015).

2.4 Dimensionless parameters

In this study, the length, velocity and time scales are defined as $\lambda^{\prime}$=$\sqrt{\sigma/(\rho_{l}-\rho_{v})g}$, $t^{\prime}$=$\sqrt{\lambda^{\prime}/g}$ and $u^{\prime}$=$\sqrt{\lambda^{\prime}g}$. The critical wavelength is given by $\lambda_{c}$=$2\pi$$\lambda^{\prime}$. The most dangerous wavelength is given by $\lambda_{d}$=$\sqrt{3}$$\lambda_{c}$. The relevant dimensionless parameters in this problem are listed below. The Grashof number $(Gr=\rho_{v}(\rho_{l}-\rho_{v})g\lambda^{\prime 3}/\mu_{v}^{2})$ describes the ratio of the buoyancy to viscous force. The Prandlt number $(Pr=c_{p,v}\mu_{v}/k_{v})$ represents the ratio of the momentum diffusion to the thermal diffusion in vapor. The Jakob number $(Ja=c_{p,v}\Delta T/h_{lg})$ represents the ratio of overheat to latent heat.

To investigate the heat transfer characteristics in film boiling studies, the Nusselt number $(Nu)$ is the most relevant dimensionless number, representing the ratio between heat convection and heat conduction. It measures the strength of boiling heat transfer during film boiling. In the current study, the space-averaged Nusselt number is defined as

$$Nu=\left.\frac{1}{W_{y}}\frac{1}{W_{z}}\int_{0}^{W_{y}}\int_{0}^{W_{z}}\frac{\lambda^{\prime}}{T_{\text{wall }}-T_{\text{sat }}}\frac{\partial T}{\partial n}\right|_{\text{wall }}\mathrm{~{}d}y\mathrm{~{}d}z.$$ (21)

Further, the space- and time-averaged Nusselt number can be calculated as

$$\overline{Nu}=\frac{1}{t_{1}-t_{0}}\int_{t_{0}}^{t_{1}}Nu\mathrm{~{}d}t,$$ (22)

where $t_{0}$ is the initial time and $t_{1}$ is the end time.

For MHD flow, the interaction number ($N$) describes the ratio between the Lorentz and gravitational forces. This is represented as

$$N=\sigma_{e(l)}B^{2}\lambda^{\prime}/\rho_{l}\sqrt{g\lambda^{\prime}}.$$ (23)

Additionally, we note that another dimensionless parameter, the Hartmann number ($Ha$), is used to describe the magnitude of the Lorentz force. It is expressed as $Ha=\sqrt{ReN}$, where $Re$ is the Reynolds number, defined as $Re=\rho_{l}u^{\prime}\lambda^{\prime}/\mu_{l}$. In this study, we primarily use the interaction number ($N$) to characterize the strength of the external magnetic field (MF), while also providing the corresponding $Ha$ values.

In the present study, we explore the film-boiling phenomenon based on the liquid metal lithium. However, in the case of liquid metal lithium, both the thermal and flow boundary layers are very thin, posing a significant challenge to achieving the required high spatial resolution for resolving the extremely thin vapor film which has a thickness of $\delta\propto\left[\frac{\mu_{v}k_{v}\Delta T}{h_{lg}\rho_{v}(\rho_{l}-\rho_{v})g}\left(\frac{\sigma}{(\rho_{l}-\rho_{v})g}\right)^{1/2}\right]^{1/4}$ (Berenson, 1961; Kim et al., 2015). Typically, we found that when the computational domain is confined to $H_{x}\times W_{y}\times W_{z}=3\lambda_{d}\times\lambda_{d}\times\lambda_{d}$, the minimum vapor film thickness cannot be resolved even when grid resolution is 512 grids within a most dangerous wavelength, which represents the diameter of the vapor bubble in the single-mode film boiling. It implies more than fifty million grids are required even if we use adaptive mesh refinement technique to reduce the computational cost. More seriously, the time step is highly restrained by the size of the smallest mesh. To solve this problem, we manually adjust some of the physical properties of liquid lithium in our numerical simulations, especially the viscosity and conductivity of the vapor phase so that the vapor layer is now almost eight times thicker than that of the original liquid lithium. The physical properties of the fluid under investigation are listed in Table 1, and we will show that the grid resolution of $\lambda_{d}/128$ can accurately resolve the thin vapor layer. Furthermore, in Appendix A, we demonstrate that even when using more realistic physical properties, the MHD effects on the rising behavior of the vapor bubble remain very similar. This supports the notion that adjustments to the fluid properties do not have a qualitative impact on the conclusions drawn in the present study.

According to Table 1, the values of characteristic scales and dimensionless parameters are presented in Table 2. The physical quantities described in the rest of the manuscript have been non-dimensionalized using the characteristic scales. And the temperature field is normalized based on the saturation temperature $T_{\mathrm{sat}}$, given that $(T-T_{\mathrm{sat}})/(T_{\mathrm{wall}}-T_{\mathrm{sat}})$ standardizing the relative magnitude of the vapor phase temperature field to values between 0 and 1. As a result, the temperature of the liquid phase is set to 0, while the temperature of the overheated wall is fixed at 1.

	$\rho$(kg/m3)	$\mu$(Ps$\cdot$s)	$k$(W/(m$\cdot$K))	$c_{p}$(J/(kg$\cdot$K))	$h_{lg}$(J/kg)	$\sigma$(N/m)	$\Delta$$T$(K)
Vapor	0.546	0.049	6.528	9259	1.933$\times$ $10^{6}$	0.23	20
Liquid	401.3	0.159	68.5	4250

Table 1: Physical properties of fluids used in the present study.

$\lambda^{\prime}$	$t^{\prime}$	$u^{\prime}$	$\lambda_{c}$	$\lambda_{d}$	$Gr$	$Pr$	$Ja$
0.0076 m	0.0279 s	0.273 m/s	0.048 m	0.083 m	0.392	69.499	0.096

Table 2: Characteristic scales and values of the dimensionless numbers used in the present study.

3 Results and discussion

3.1 Single-mode film boiling

In the case of single-mode film boiling with and without MFs, a computational domain with dimensions $3\lambda_{d}\times\lambda_{d}\times\lambda_{d}$ is used. An initial perturbation in the following form is applied to the phase interface:

$$x=\frac{\lambda_{d}}{128}\left(8.0+\mathrm{cos}\left(\frac{2\pi y}{\lambda_{d}}\right)+\mathrm{sin}\left(\frac{2\pi z}{\lambda_{d}}\right)\right).$$ (24)

A grid convergence test is conducted for film boiling without MF, considering three different grids: $192\times 64\times 64$ (referred to as grid-level 6), $384\times 128\times 128$ (grid-level 7), and $768\times 256\times 256$ (grid-level 8). The vapor bubble shapes on the vertical cross-section ($XOY$ plane) at $t=14.8$ and $t=16.7$ are illustrated in Fig. 2(a) and Fig. 2(b), respectively. The temporal evolution of the space-averaged Nusselt number ($Nu$) is also depicted in Fig. 2(c). It is evident that the grid resolution at grid-level 7 is adequate to achieve the desired computational accuracy. Consequently, this grid size is adopted for the remainder of our study.

The computational results presented in this study have been obtained using workstations equipped with two AMD^® EPYC 7742 CPUs, each with 64 cores operating at 2.20 GHz and 256 GB of RAM. For single-mode boiling flows without an external magnetic field, as discussed in §3.1, a typical run with a grid refined to Level 7 and covering a dimensionless time span of $0<t<75$ required more than two months using 24 processors. When a vertical magnetic field was applied at $N=6.3$, the computational time extended to over three months due to the additional Poisson equation that needed to be solved. For multi-mode film boiling simulations discussed in §3.2, the computations were carried out using 128 AMD processors. In this case, with a grid refined to Level 7 and a dimensionless computational time span of $0<t<100$ for $N=11.1$, each simulation took nearly six months to complete. Moreover, test cases indicated that refining the grid from Level 7 to Level 8 would increase the computational time by approximately threefold, primarily due to the increased time constraints imposed by CFL condition. For example, the dimensionless time step was reduced from 2.7 $\times 10^{-4}$ to 1.0 $\times 10^{-4}$. Based on a comprehensive evaluation of computational efficiency and accuracy, we concluded that a spatial resolution of Level 7 is the most practical choice for our numerical simulations.

Figure 2: The grid independence test performed at different grid resolutions. The bubble shapes on the vertical cross-section ($XOY$ plane) at (a) $t=14.8$ and (b) $t=16.7$ for different grid resolutions. (c) The temporal evolution of the space-averaged Nusselt number ($Nu$) obtained using different grid resolutions.

3.1.1 Film boiling without a MF

Figure 3: Temporal evolution of the space-averaged $Nu$ for the single-mode film boiling without a MF. The dashed line represents the theoretical results of Klimenko (1981).

We begin the presentation of our results by analyzing the film boiling phenomenon without MF and compare with existing boiling mechanisms and empirical models. Fig. 3 shows the temporal evolution of the space-averaged $Nu$. It can be seen that $Nu$ oscillates periodically around the correlation proposed by Klimenko (1981), given that $Nu_{k}=0.19(GrPr)^{\frac{1}{3}}\times(0.89Ja^{-\frac{1}{3}})$, this agreement indicates the credibility of the numerical simulation results. This periodic oscillation corresponds to the shedding of bubbles. The evolution of the liquid/vapor interface for the fifth bubble cycle (from $t=51.9$ to $t=63.0$), corresponding to the blue region in Fig. 3, is depicted in Fig. 4, while Fig. 5 presents the corresponding streamline and temperature field distributions within the vertical cross-section ($XOY$). It is noteworthy that during the fifth bubble cycle, the bubble generated in the fourth bubble cycle has just detached, resulting in the presence of a rising bubble in Fig. 4(a-d). The film boiling process within one bubble cycle can be divided into four stages: Stage one: From $t=51.9$ to $t=53.7$, the phase change causes the expansion of the vapor phase volume, and the bulge grows rapidly. Stage two: From $t=53.7$ to $t=55.6$, the Rayleigh-Taylor instability occurs, leading to the formation of a significant bulge and vortices. Stage three: From $t=55.6$ to $t=59.3$, the neck connecting the thin film and bubble starts to contract, leading to a mushroom-shaped bubble and the vortex structure. Stage four: From $t=59.3$ to $t=63.0$, the bubble pinches off, and a new bubble cycle initiates. Furthermore, within one bubble cycle, the vapor film thickness initially increases due to the expansion of the vapor phase volume and the growth of the bulge. Later, it decreases when the neck phenomenon occurs. Conversely, the Nusselt number decreases first and then increases within one bubble cycle, as explained by Eq. (21). The temperature profile reveals a linear distribution from the overheated wall to the liquid/vapor interface throughout the film boiling process. A smaller distance between the liquid/vapor interface and the overheated wall corresponds to a more significant temperature gradient, resulting in a higher $Nu$. In contrast, a larger distance leads to a smaller $Nu$. These findings align with previous studies (Juric & Tryggvason, 1998; Zhang & Ni, 2018). Next, we investigate the influence of a vertical MF on the film boiling dynamics.

Figure 4: Vapor bubble generation during the film boiling flow without a MF at different dimensionless times, which correspond to the blue region in Fig. 3. (a) $t=51.9$, (b) $t=53.7$, (c) $t=55.6$, (d) $t=57.4$, (e) $t=59.3$ and (f) $t=63.0$.

Figure 5: Temperature profile (contours) and velocity streamlines (black arrow lines) in the vertical cross-section ($XOY$) of the computational domain. The time instants from panels (a-f) correspond to those displayed in Fig. 4(a-f).

3.1.2 Effect of the vertical MF

\begin{overpic}[width=298.75394pt,tics=1]{figure/single-ver-shape.eps} \put(-3.5,98.1){\large(a)} \put(-3.5,73.3){\large(b)} \put(-3.5,48.6){\large(c)} \put(-3.5,23.8){\large(d)} \end{overpic}

Figure 6: Temporal evolution of the bubble shapes during the film boiling flow without and with the vertical magnetic fields. (a) $N=0$ (B = 0 T), (b) $N=6.3$ (B = 0.15 T), (c) $N=11.1$ (B = 0.2 T), and (d) $N=44.5$ (B = 0.4 T).

We examine the influence of the vertical MF applied along the $x$ direction on film boiling flows. The horizontal directions perpendicular to the applied MF are represented by $y$ and $z$. Thus, it is essential to highlight that under the vertical MF, the MHD effect is expected to be isotropic in any horizontal plane parallel to $YOZ$. Consequently, the results in two vertical cross-sections, namely the $XOY$ and $XOZ$ planes, are the same. Therefore, we only illustrate findings for the $XOY$ cross-section. We vary the MF intensities with B = 0.05 T, B = 0.1 T, B = 0.15 T, B = 0.2 T and B = 0.4 T. According to Eq. (23), the corresponding interaction numbers (Hartmann numbers) are $N=0.7$ ($Ha=1.9$), $N=2.8$ ($Ha=3.8$), $N=6.3$ ($Ha=5.7$), $N=11.1$ ($Ha=7.6$), and $N=44.5$ ($Ha=15.3$), respectively. The MF intensity choice considers a spectrum ranging from inertia being predominant to the Lorentz force being dominant.

Figure 7: (a-c) Comparison of the temperature field and streamline patterns in the vertical cross-section ($XOY$) under different vertical magnetic fields. The right-hand side of panels (b) and (c) shows the comparison of the Lorentz force (black arrows) and horizontal velocity distribution ($u_{y}$, the positive direction: from left to right) (contours). (d-e) The distribution of the current line. (a) $N=0$, $t=57.4$. (b) and (d) $N=11.1$, $t=53.7$. (c) and (e) $N=44.5$, $t=40.7$.

To illustrate the MHD effect, we compare the liquid/vapor interface patterns within four bubble generation cycles, and the results are presented in Fig. 6. From the first and second columns, it can be seen that the vertical MF inhibits the necking phenomenon and delays the detachment of the bubble. For instance, at $t=16.7$, the bubble is about to fall off when $N=0$, while the necking phenomenon does not appear when $N=44.5$. Besides, as the MF strength increases, the periodic detachment of the bubble is gradually suppressed, and the vapor jet emerges. At $N=0$, the periodic detachment of the bubble caused by Rayleigh-Taylor instability can be observed, and the bubble detachment locations for each bubble cycle are almost the same. However, when $N=6.3$ and $N=11.1$, this phenomenon can only be observed in the previous short period. After the boiling flow completely develops, the vapor jet arises, and a small bubble is released from the top of the jet, corresponding to the onset of vapor jet instability. Also, the bubble detachment location increases, implying that the MF suppresses the onset of Rayleigh-Taylor instability. When $N$ increases to $44.5$, the detachment of the bubble is barely observable, and boiling takes place in the form of a steady vapor column after the flow completely develops. In this stage, small amplitude surface waves can still be observed at the interface, indicating that the high-intensity MF further inhibits the occurrence of jet instability.

To investigate the impact of a vertical MF on film boiling flow and its consequent influence on the boiling pattern, we initially compare the streamlines for cases $N=0$, $N=11.1$ and $N=44.5$, as depicted in Fig. 7(a-c). It can be seen that with an increase in MF intensity, the suppression of the vortex structure, owing to the MHD effect, becomes more pronounced. In Fig. 7(a), where $N=0$, vortices are observed both at locations where the bubble has detached and where the bubble is still connected to the thin film. However, in the left-half panel of Fig. 7(b), where $N=11.1$, vortices are only present at the location where the bubble has detached. Vortices no longer appear when $N$ increases to 44.5, as depicted in the left-half panel of Fig. 7(c). Additionally, close inspection of the bubble shapes reveals that the interface curvature at the bottom of the bulge decreases, and the necking phenomenon is thus weakened.

Now, we aim to elucidate the mechanisms underlying the aforementioned phenomenon. First, the distribution of the current line surrounding the liquid/vapor interface is plotted in Fig. 7(d) and 7(e). It can be seen that the current lines form closed loops in two cases, indicating that the computation of the induced current is conserved. The right-hand panels of Fig. 7(b) and 7(c) present the distribution of the Lorentz force (black arrows) and the horizontal component of velocity (contours, $u_{y}$ in the $XOY$ cross-section). It is essential to clarify that, due to a significantly large ratio of electric conductivities ($\sigma_{e,l}/\sigma_{e,v}$), the Lorentz forces in the vapor phase can be neglected so that they only manifest in the liquid phase. It is evident that the Lorentz force induced by the vertical MF acts in the opposite direction to $u_{y}$. Notably, under moderate MF intensity (Fig. 7(b)), the Lorentz force near the upper part of the bubble points towards the inside of the bubble, while at the lower part, it points towards the outside. This results in a torque opposite to the direction of the local vortex, as indicated by the purple circle with arrows in the diagram. Consequently, the horizontal velocity is weakened, consistent with the color scale results in Fig. 7(a-c). This explains the suppression of vortex formation, and such decay in vortices counteracts the pressure difference across the vapor-liquid interface. Thus, the surface tension, which balances this pressure difference in the form of $\sigma\kappa=p_{v}-p_{l}$, should be reduced correspondingly. The magnitude of surface tension is positively correlated with the curvature of the interface, and thus, the presence of the Lorentz force reduces interface curvature, making the bubble more challenging to pinch off. Simultaneously, due to the rapid evaporation of the liquid, the continuous generation of vapor converges towards the bulge, ultimately leading to the formation of a vapor jet under the influence of buoyancy. For vapor jets, bubble detachment is triggered by the amplification of instabilities and surface waves (Lin & Reitz, 1998). Then, in the scenario with a large MF intensity ($N=44.5$, Fig. 7(c)), the Lorentz forces entirely suppress the vortices and lead to uniformly distributed $u_{y}$, with the flow directing toward the center of the vapor column and upward flow. In this context, the Lorentz force also plays a role in suppressing the growth of surface waves while concurrently sustaining boiling in the form of a vapor column. In §3.1.4, an additional case based on the large MF intensity case $(N=44.5)$ will be further introduced to validate this finding.

Figure 8: (a) Temporal evolution of the space-averaged $Nu$ in the single-mode film boiling under different vertical MF intensities. The dashed line represents the theoretical results of Klimenko (1981). (b) Variation of the space- and time-averaged $Nu$ number ($\overline{Nu}$) with vertical MF intensity in the case of the single-mode film boiling. The dimensional parameter $B$ is used on the horizontal axis to enhance clarity in presenting the value distribution.

The temporal evolution of the space-averaged Nusselt number ($Nu$) under different vertical MFs is shown in Fig. 8(a). It can be observed that as the vertical MF intensity increases, the $Nu$ transitions from approximately constant amplitude oscillation ($N\leq 2.8$) to damped oscillation ($N\geq 6.3$). This observation is consistent with the variations in the bubble detachment pattern under the influence of a vertical MF, suggesting that the vertical MF serves as a damping force, restraining bubble pinch-off behavior. Specifically, for cases with $N=6.3$, $N=11.1$, and $N=44.5$, the initially oscillated $Nu$ eventually stabilizes at a constant value. This stabilization occurs because, at these MF intensities, film boiling tends to maintain a stable jet form after the full development of boiling flow, as depicted in Fig. 6. Bubble detachment occurs at locations far from the overheated wall or even without bubble detachment, keeping the interface near the overheated wall relatively stable. Given that, during saturated film boiling, the temperature field in the vapor phase exhibits a nearly linear distribution, the temperature distribution remains unchanged, as observed in Figs. 6 and 7.

Fig. 8(b) illustrates the space- and time-averaged Nusselt number ($\overline{Nu}$) calculated using Eq. (22) for different MF intensities. In cases of constant amplitude oscillation ($N=0$, $N=0.7$, and $N=2.8$), the $Nu$ is averaged over time from the second trough to the seventh trough in Fig. 8(a), in order to eliminate the influence of the first bubble cycle. The final constant value is utilized in the case of damped oscillation, where $N=6.3$, $N=11.1$, and $N=44.5$. The dimensional parameter $B$ is used on the horizontal axis to enhance clarity in presenting the distribution of $\overline{Nu}$ under different MF intensities. The chart is demarcated with an orange dashed line based on the pattern of bubble detachment and $Nu$ variation. The left side of the dashed line represents the periodic bubble detachment area, where the $Nu$ exhibits constant amplitude oscillation. In this region, the vertical MF results in a slightly higher $\overline{Nu}$ than the case without a MF. On the right side of the dashed line is the vapor column boiling area, where the $Nu$ exhibits damped oscillation. The $\overline{Nu}$ sharply declines in the transition zone. After that, the $\overline{Nu}$ gradually increases with the increase in MF intensity.

Figure 9: (a) Comparison of the bubble shapes at $t=70.4$ under three MF intensities ($N=6.3$, $N=11.1$, and $N=44.5$). (b) Temporal evolution of the space-averaged Nusselt number ($Nu$) and the minimum distance between the liquid/vapor interface and the overheated wall over time after boiling reaches a steady state under these three MF intensities.

Figure 10: (a) Variation of the vertical component of velocity ($u_{x}$, the positive direction: from bottom to top) with dimensionless height ($x/\lambda_{d}$) on the axis line near the overheated wall (see black dashed line in panel b) under three different magnetic field intensities ($N=6.3$, $N=11.1$, and $N=44.5$) at $t=70.4$. Panels (b)-(d) show the distribution of $u_{x}$ for $N=6.3$, $N=11.1$, and $N=44.5$ at $t=70.4$, respectively. (e) The distribution of $u_{x}$ for $N=0$ at $t=57.4$.

To understand the reason behind the increase in the $\overline{Nu}$ with the rise in MF intensity in this region ($N\geq 6.3$), we compare the bubble shapes at $t=70.4$ under three different MF intensities, as shown in Fig. 9(a). Note that in the same picture, the embedded panel illustrates the vapor interface in the vapor film region ($y>0.2$), and Fig. 9(b) presents the temporal evolution of the minimum distance between the interface and the overheated wall, and the development of the corresponding $Nu$ is also depicted. It can be observed that the vapor thickness decreases with the increase in MF intensity, leading to an increase in the $Nu$ and $\overline{Nu}$. Now, the question arises: why does the vapor thickness decrease by generating a vapor jet as the MF intensity rises? Fig. 10(a) explains the reasons by plotting the variation of the vertical component of velocity $(u_{x})$ with height on the axisymmetric line near the overheated wall (as shown by the black dashed line in Fig. 10(b-d)) under three different MF intensities ($N=6.3$, $N=11.1$, and $N=44.5$) at $t=70.4$. This figure shows that the greater the MF intensity, the higher the vertical velocity component near the overheated wall on the axisymmetric line. This allows vapor generated through evaporation to flow quickly from the near-wall region into the vapor jet, resulting in a smaller minimum distance between the vapor interface and the overheated wall. Analyzing the vertical component of the velocity contour ($u_{x}$) plot (Fig. 10(b-d)), it is observed that the vertical velocity is highest at the necking point of the jet, leading to a blocking effect that there is a sharp decrease in vertical velocity at its lower end (marked in the red box in the pictures). This decrease suppresses the process of upstream vapor entering the jet interior. Previous analyses have indicated that an increase in MF intensity suppresses jet instability, causing the necking point of the jet to move farther from the overheated wall. Therefore, the necking phenomenon has a weaker inhibitory effect on the velocity near the wall on the symmetry axis. In summary, when boiling occurs in the form of a jet, the location of bubble detachment influences the velocity at which vapor enters the center of the jet. The higher the bubble detachment position, the greater the vapor velocity in the jet, making vapor overflow more quickly. Thus, a higher MF intensity decreases the vapor film thickness near the overheated wall, promoting the $Nu$ and $\overline{Nu}$.

The abrupt cliff-like decrease in the $\overline{Nu}$ during the transition region in Fig. 8(b) can be explained through the same analytical method. Fig. 10(e) presents the distribution of $u_{x}$ for $N=0$ at $t=57.4$, when the necking phenomenon caused by Rayleigh-Taylor instability is distinctly visible. When film boiling occurs in the form of a periodic bubble detachment, the location of bubble necking is near the vapor film and close to the overheated wall. Thus, the necking phenomenon brings about a sharp increase in vertical velocity, directly translating into an abrupt rise in vapor overflow velocity, as shown by the color scale in Fig. 10(e). This facilitates the rapid influx of vapor generated by evaporation into the bulge, significantly reducing the vapor film thickness during the bubble necking stage. As a result, the $Nu$ increases, enhancing the $\overline{Nu}$.

3.1.3 Effect of the horizontal MF

\begin{overpic}[width=298.75394pt,tics=1]{figure/single-hor-shape.eps} \put(-3.5,98.1){\large(a)} \put(-3.5,73.3){\large(b)} \put(-3.5,48.6){\large(c)} \put(-3.5,23.8){\large(d)} \end{overpic}

Figure 11: Temporal evolution of the vapor bubble shapes during the film boiling flow without and with the horizontal MFs. (a) $N=2.8$ (B = 0.1 T), (b) $N=6.3$ (B = 0.15 T), (c) $N=11.1$ (B = 0.2 T) and (d) $N=44.5$ (B = 0.4 T).

Next, we investigate the effect of a horizontal MF on the single-mode boiling configuration. Here, the MF is applied in the $y$ direction. Thus, the vertical $XOY$ plane is parallel to the MF, while the vertical $XOZ$ plane is perpendicular to the MF. The choice of MF strengths is the same as those used in the case of vertical MF configuration, namely $N=0.7$ (B = 0.05 T), $N=2.8$ (B = 0.1 T), $N=6.3$ (B = 0.15 T), $N=11.1$ (B = 0.2 T) and $N=44.5$ (B = 0.4 T).

Fig. 11 depicts the evolution of the liquid/vapor interface patterns under different horizontal MF intensities. Anisotropy is observed in the boiling flow when a horizontal MF is applied. For instance, in the latter part of the process ($t>30.6$) for $N=2.8$, columnar bulges out at the thin film grow along the MF direction, while the detaching bubble is also non-axisymmetric. With further increases in the MF intensity ($N=6.3$, $N=11.1$, and $N=44.5$), this characteristic manifests much earlier that a two-dimensional boiling flow forms and the vapor features as a thin sheet along the direction of the MF after complete development. Furthermore, after the vapor sheet grows to a certain height (see those scenarios at $t=74.1$), the amplification of minor disturbances can lead to the instability of the vapor sheet. In the study of sheet instability (Hagerty & Shea, 1955), two modes of instability for sheets have been identified: the sinuous mode and the varicose mode. From Fig. 11, it is evident that the flow corresponds to the sinuous mode.

Figure 12: (a-b) Comparison of temperature field and streamline patterns in two vertical cross-sections ($XOY$: parallel to the magnetic field, on the left-hand side; $XOZ$: perpendicular to the magnetic field, on the right-hand side) under different horizontal magnetic field intensities. (c-d) Comparison of the Lorentz forces (black arrows) and the distribution of the vertical component of velocity ($u_{x}$, contours) corresponding to panels (a-b). (e) Distribution of the current lines. (a) and (c) $N=6.3$, $t=27.8$. (b), (d), and (e) $N=11.1$, $t=27.8$.

To understand those scenarios caused by the horizontal MFs, we examine the temperature and streamline distribution in two vertical cross-sections: the $XOY$ plane parallel to the MF and the $XOZ$ plane perpendicular to the MF. These distributions are plotted for two different MF intensities ($N=6.3$ and $N=11.1$) at $t=27.8$, as illustrated in Fig. 12(a) and 12(b), respectively. The corresponding distribution of Lorentz forces and vertical velocity component ($u_{x}$) are presented in Fig. 12(c-d). Fig. 12(e) shows the distribution of current lines under the horizontal MF. It can be observed that the current lines are not closed. This occurs because we set periodic boundary conditions around the computational domain, causing the current lines to pass through the boundaries. By comparing the velocity streamline distribution (see Fig. 12(a) and 12(b)), it becomes evident that the horizontal MF primarily influences the flow vortices in the direction perpendicular to the MF, which suppresses the vortex on the side of the bubble near the wall. When $N$ increases to 11.1, the vortex completely disappears, and no downward flow is observed on the $XOY$ plane. Then, the distribution of Lorentz force, as displayed in Fig. 12(c) and 12(d), explains the disappearance of the vortex perpendicular to the MF. From left-half panels of Fig. 12(c) and 12(d), it is observed that in the $XOY$ plane, the direction of the Lorentz force is upward on the side of the bubble near the wall, while the Lorentz force is directed vertically downward at the top of the bubble. This results in a torque opposite to the direction of the local vortex, as indicated by the orange circle with arrows in the diagrams. So, the horizontal MF tends to suppress vortices perpendicular to the MF. As a result, vortices gradually disappear, and the necking phenomenon weakens, causing the interface to move away from the overheated wall under the influence of continuously generated vapor. Eventually, a columnar bulge along the MF direction is formed. In contrast, for the $XOZ$ plane perpendicular to the MF (see right-half panels of Fig. 12(c) and 12(d)), the direction of the Lorentz force is uniformly vertical downward, and the flow field structure does not exhibit significant changes with variations in MF intensity. Thus, the necking phenomenon still occurs in this plane, forming a vapor sheet. These conclusions align with the theoretical prediction of MHD effects(Sommeria & Moreau, 1982; Davidson, 1995), demonstrating that the Lorentz force tends to decrease the velocity gradient parallel to the MF to reduce Joule dissipation. Furthermore, the temperature distribution depicted in Fig. 12(a) and 12(b) exhibits a pattern similar to that of the vertical MF cases. It can be seen that as the MF intensity increases, the temperature inside the bulge becomes lower. This trend is in line with the distribution of the vertical component of velocity $(u_{x})$ shown in Fig. 12(c) and 12(d), and the underlying reason is analogous to the vertical MF situation. Specifically, with the increase in MF intensity, the weakening of the necking phenomenon leads to a reduction in the velocity of vapor entering the bulge, making heat transfer more challenging. Consequently, the internal temperature becomes lower.

Fig. 13(a) presents the temporal evolution of the $Nu$ under varying horizontal MF intensities. At lower MF intensities $(N<6.3)$, the $Nu$ variation demonstrates approximate periodic oscillations, aligning with the periodic detachment of bubbles. Notably, when $N=2.8$, a considerable reduction in amplitude occurs, attributed to the initiation of the horizontal MF effect. This effect induces the formation of a columnar bulge along the MF direction at the vapor film’s bottom, leading to an increased vapor layer thickness and subsequently causing a decline in $Nu$. As the MF intensity increases $(N\geq 6.3)$, $Nu$ undergoes a brief period of periodic oscillations in the initial stages due to bubble detachment. However, once a stable vapor sheet is established, a transition to a more stable state occurs, and higher MF intensities accelerate this transition.

Figure 13: (a) Temporal evolution of the space-averaged $Nu$ under different horizontal MF intensities in single-mode film boiling. The dashed line represents the theoretical results of Klimenko (1981). (b) Variation of the space and time averaged $Nu$ number ($\overline{Nu}$) with horizontal MF intensity in the case of the single-mode film boiling. The dimensional parameter B is used on the horizontal axis to enhance clarity.

Figure 14: Comparison of the bubble shapes in the vertical cross-section ($XOZ$ plane) perpendicular to the magnetic field under three different MF intensities ($N=6.3$, $N=11.1$, and $N=44.5$). In panel (a), the red, green, and blue lines depict the results for $N=6.3$ at $t=42.6$, $N=11.1$ at $t=33.3$, and $N=44.5$ at $t=22.2$, respectively. In panel (b), the red, green, and blue lines depict the results for $N=6.3$ at $t=57.4$, $N=11.1$ at $t=50.0$, and $N=44.5$ at $t=40.7$, respectively.

Figure 15: Comparison of the vertical component of the Lorentz force ($F_{l,x}$) on the liquid/vapor interface under different MF intensities corresponding to Fig. 14. In panel (a), the red, green, and blue lines depict the results for $N=6.3$ at $t=42.6$, $N=11.1$ at $t=33.3$, and $N=44.5$ at $t=22.2$, respectively. In panel (b), the red, green, and blue lines depict the results for $N=6.3$ at $t=57.4$, $N=11.1$ at $t=50.0$, and $N=44.5$ at $t=40.7$, respectively.

The $\overline{Nu}$ for the aforementioned cases is plotted in Fig. 13(b). Similar to the vertical MF, for small MF intensity cases $(N<6.3)$, the average of the $Nu$ from the second trough to the seventh trough in Fig. 13(a) is selected, while for large MF intensity cases, the final value is employed. Based on the final character of the film boiling, the chart is divided into two regions. The left side represents the periodic detachment region of the bubble. In this region, the $\overline{Nu}$ experiences a sharp decrease after the MF reaches a certain intensity. This is consistent with the previously mentioned reduction in the amplitude of the $Nu$, both attributed to the appearance of the elongated bulge in the vapor film. The right side is the vapor sheet boiling region, where the $\overline{Nu}$ increases as the MF intensity increases. To investigate the reason, we compare the shapes of the interface in the $XOZ$ plane at the initial formation of the vapor sheet (Fig. 14(a)) and when the vapor sheet stably exists (Fig. 14(b)) under these three different MF intensities. It is clear that there is a significant difference in the thickness of the thin film near the overheated wall, particularly away from the vapor sheet side, as indicated by the red square in the diagram. The higher the MF strength, the smaller the thickness, resulting in a larger $Nu$ and $\overline{Nu}$. From Fig. 12(c) and 12(d), we know that in the $XOZ$ plane, the Lorentz force is consistently directed downward. Therefore, we hypothesize that the Lorentz force might suppress the expansion of the vapor film near the overheated wall. To verify this hypothesis, we output the vertical component of the Lorentz force ($F_{l,x}$) on the near wall interface corresponding to Fig. 14, as depicted in Fig. 15. It can be observed that as the MF intensity increases, the downward effect of the Lorentz force becomes more pronounced. In other words, the suppression of the expansion of the vapor film near the overheated wall is stronger, resulting in a smaller thickness of the vapor film and a larger $\overline{Nu}$.

After examining the MHD effect on the growth of vapor bubbles, it is useful to compare our findings with those of previous studies. The most relevant work is by Chandrasekhar (2013), who used linear stability analysis to investigate the effects of vertical and horizontal magnetic fields on Rayleigh-Taylor instability in inviscid interfacial flows. Specifically, the exponential growth rate $n$ of disturbances for a flat interface between inviscid fluids in ideal MHD, with a vertical MF applied to the interface, is given by: $n^{3}+n^{2}kB_{0}\left(\sqrt{\rho_{l}/\rho_{0}}+\sqrt{\rho_{v}/\rho_{0}}\right)+n\left(k^{2}B_{0}^{2}-(\rho_{l}-\rho_{v})/(\rho_{l}+\rho_{v})gk\right)=2(\rho_{l}-\rho_{v})/(\rho_{l}+\rho_{v})gk^{2}B_{0}/(\sqrt{\rho_{l}/\rho_{0}}+\sqrt{\rho_{v}/\rho_{0}})$, where $g$ is the gravitational acceleration, $k$ is the wavenumber, $B_{0}=B/\sqrt{\rho_{0}\mu_{0}}$, $\rho_{0}=(\rho_{l}+\rho_{v})/2$, and $\mu_{0}$ is the vacuum magnetic permeability. Notably, there is no critical wavelength, so all modes are unstable, but a larger $B_{0}$ enhances the stability of the interface. This aligns with our results, which show that the vapor bubble grows more slowly under a stronger vertical MF. However, the RT instability in boiling flow cannot be completely eliminated, even with a very strong MF. For the horizontal MF parallel to the interface, linear stability analysis also reveals anisotropic effects induced by the external MF. Horizontal disturbances transverse to the $\mathbf{B}$-lines are unaffected by the horizontal MF, while disturbances parallel to the $\mathbf{B}$-lines have an exponential growth rate $n$ given by: $n^{2}=gk(\rho_{l}-\rho_{v})/(\rho_{l}+\rho_{v})-k^{2}{B_{0}}^{2}$, with $k$ denoting the wavenumber of the disturbances in that direction. This shows that a sufficiently strong MF can suppress Rayleigh-Taylor instability in the direction parallel to the $\mathbf{B}$-lines but has little impact in the direction transverse to the $\mathbf{B}$-lines, consistent with our numerical results. We also acknowledge the importance of future studies addressing these challenges by incorporating a phase change model into the framework of linear stability analysis for a more comprehensive understanding of boiling coupled RT instability under MHD conditions.

3.1.4 What happens after MF removal?

To further understand the MHD effect on film boiling, in this section, we investigate what happens when the MF is suddenly removed after film boiling reaches a steady state under MF (vapor column boiling for vertical MF scenario and vapor sheet boiling for horizontal MF scenario).

Figure 16: Temporal evolution of the temperature profile (color contours) and velocity streamline pattern (black arrow lines) in the vertical cross-section ($XOY$) after removing the magnetic field at $t=44.4$ for $N=44.5$. (a) $t=44.4$, (b) $t=50.0$, (c) $t=51.9$, (d) $t=54.6$, and (e) $t=56.5$ and $t=59.3$. It can be seen that the removal of the magnetic field leads to vapor jet instability, resulting in the reappearance of bubble detachment.

Figure 17: Comparison of the temporal evolution of (a) the space-averaged Nusselt number ($Nu$) and (b) the evaporation rate ($\dot{m}S_{\Gamma}$) with the vertical MF and after the removal of the MF at $t=44.4$ for $N=44.5$. Note that the oscillation of the curves, corresponding to the magnetic field removal, is due to the reappearance of vapor bubble detachment.

Initially, we examine the scenario based on the vertical MF with large intensity $(N=44.5)$. The MF is suddenly removed when the film boiling reaches a relatively stable state $(t=44.4)$. The results are depicted in Fig. 16. It is observed that, after the removal of the MF, the amplitude of surface waves progressively increases, and vortices appear, ultimately leading to the detachment of small bubbles. This observation supports the findings discussed in §3.1.2 that the Lorentz force plays a vital role in suppressing the growth of surface waves and the subsequent detachment of the vapor bubble. The comparison of the temporal evolution of the $Nu$ with the MF and upon its removal is illustrated in Fig. 17(a). After removing the vertical MF, a noticeable decrease in the $Nu$ can be observed. This is attributed to the fact that, upon the removal of the vertical MF, the stable columnar boiling pattern can no longer be sustained, leading to the detachment of bubbles at the top of the vapor column (see Fig. 16). As analyzed in §3.1.2, the bubble detachment behavior inhibits the vapor influx into the columnar bulge, causing the near wall vapor film to thicken and reducing $Nu$. Fig. 17(b) compares the evaporation rate (here, we measure the evaporation rate using the mass flow at the interface ($\dot{m}S_{\Gamma}$)) between the cases, revealing a significant decrease in the evaporation rate after the removal of the MF, which is also caused by the increase in thickness of the near wall vapor film. Note that the oscillation of $\dot{m}S_{\Gamma}$ after removing the MF is caused by the reappeared detachment of the vapor bubble. The change in evaporation rate provides another perspective on the reason for bubble detachment. After removing the MF, the decrease in evaporation rate prevents film boiling from supplying sufficient vapor for stable columnar bulge.

Figure 18: (a-c) The temporal evolution of the liquid/vapor interface after removing the horizontal MF at $t=44.4$ for $N=44.5$. (d-f) Temperature profile (color contours) and velocity streamline pattern (black arrow lines) in the vertical cross-section ($XOZ$) corresponding to panels (a-c). Panels (a,d), (b,e), and (c,f) correspond to $t=59.3$, $t=63.0$, and $t=66.7$, respectively.

Figure 19: Comparison of the temporal evolution of (a) the space-averaged Nusselt number ($Nu$) and (b) the evaporation rate ($\dot{m}S_{\Gamma}$) with the horizontal MF and after the removal of the magnetic field at $t=57.4$ for $N=22.2$. We observe that in both plots, the values are almost identical before and after removing the horizontal magnetic field.

Now, we examine the scenario based on a horizontal MF with the intensity of $N=22.2$. When film boiling reaches a relatively stable vapor sheet boiling at $t=57.4$, the MF is suddenly removed. However, intriguingly, we observe a very interesting result: film boiling in the vapor sheet state does not appear to be affected after removing the horizontal MF, as depicted in Fig. 18. Furthermore, by comparing the $Nu$ and evaporation rate (Fig. 19), it is observed that with and without the presence of a horizontal MF, both $Nu$ and evaporation rate are almost identical. This indicates that the vapor sheet formed under horizontal MF is initially stable, in contrast to the vapor column generated under vertical MF. To explain this, we should remember that the vertical MF eliminates the vortices around the vapor bubble in an axisymmetrical manner, and such a controlling effect leads to the generation of an isotropic vapor column. Therefore, the axisymmetrical Lorentz forces, directing radial-wards, maintain the force balance around the vapor column. Consequently, if the Lorentz forces are removed to break such force balance, the downward flow must develop due to any disturbance that emerges on the vapor column, resulting in an expanding vortex region and, hence, a detached bubble. Nevertheless, the horizontal MF suppresses the vortices in the parallel plane ($XOY$) but has little influence on the vortices in the perpendicular plane ($XOZ$). Besides, Fig. 12(c) and 12(d) reveals that corresponding to the stable vapor sheet, the Lorentz forces direct almost downwards in terms of $F_{l,x}\propto u_{x}\cdot B_{y}^{2}$, because any $z-$ component of $\boldsymbol{F}_{l}$ requires the participation of $u_{z}$, which almost disappears for a stable and two-dimensional vapor sheet. It indicates that for the vapor sheet scenario under horizontal MF, the Lorentz forces play little role in the force balance in radial direction, and thus, removing it causes little influence on the vapor pattern.

3.2 Multi-mode film boiling

For multi-mode film boiling scenarios with and without MFs, the computational domain size is increased to $4\lambda_{d}\times 4\lambda_{d}\times 4\lambda_{d}$. To eliminate the influence of artificial disturbances, a randomly initial perturbation is introduced and defined by the following equation:

$$x=\frac{\lambda_{d}}{16}+0.05/n\sum_{i=1}^{n}A(i)\left(\cos\left(\frac{2\pi iy}{3\lambda_{d}}\right)+\sin\left(\frac{2\pi iz}{3\lambda_{d}}\right)+\cos\left(\frac{2\pi iy}{3\lambda_{d}}\right)+\sin\left(\frac{2\pi iz}{3\lambda_{d}}\right)\right),$$ (25)

where $n=10$ denotes the number of random waves, and $0\leq A(i)\leq 1$ represents a random number. We vary the MF intensities to B = 0 T, B = 0.1 T, B = 0.2 T, B = 0.4 T and B = 0.6 T. The corresponding interaction numbers (Hartmann numbers) for both vertical and horizontal magnetic fields are as follows: $N=0$ ($Ha=0$), $N=2.8$ ($Ha=3.8$), $N=11.1$ ($Ha=7.6$), $N=44.5$ ($Ha=15.3$), and $N=100.1$ ($Ha=22.9$), respectively.

3.2.1 Effect of the vertical MF

\begin{overpic}[width=349.96898pt,tics=1]{figure/mut-ver-shape.jpg} \put(-4.0,76.0){\large(a)} \put(-4.0,50.5){\large(b)} \put(-4.0,24.8){\large(c)} \end{overpic}

Figure 20: Temporal evolution of the bubble shapes during multi-mode film boiling without and with the vertical MFs. (a) $N=0$ (B=0.0 T), (b) $N=11.1$ (B=0.2 T), and (c) $N=100.1$ (B=0.6 T).

Fig. 20 depicts the evolution of the liquid/vapor interface under three different vertical MFs. In the initial stage $(t<60.0)$ of film boiling, as observed in the first two columns of figures, the distribution of bubbles aligns with the theory proposed by Zuber (1958) for the most dangerous wavelength ($\lambda_{d}$). The distance between adjacent bubbles correlates with $\lambda_{d}$ even if the initial disturbances are randomly imposed. Additionally, in the early stage, the interaction between bubbles is not pronounced, allowing for a clear observation that under the multi-bubble model, the impact of a vertical MF on bubble shapes follows the conclusions drawn from the single-mode film boiling: as the MF strength increases, the film boiling transitions from the periodic detachment of bubbles from the bottom to the vapor jet releasing small bubbles from the top. This evolution ultimately leads to vapor columnar boiling without any bubble detachment, a pattern confirmed by Fig. 21, which presents the vertical cross-section of a single bubble in the multi-mode, illustrating the distribution of the Lorentz force, the $z-$ component of velocity $(u_{z})$, the temperature field, and the velocity streamlines. A comparison with Fig. 7 reveals that in the multi-mode boiling, the effect of the vertical MF aligns with the observations in the single-mode boiling. The evolution of $Nu$ also mirrors the single-mode cases in the early stages. As the magnetic field strength increases, the Nusselt number shifts from approximate periodic oscillations to damped oscillations. A slight increase in the averaged value of $Nu$ for large $N$ is also observed, as shown in Fig. 22.

Next, let us focus on the later stages of the film boiling process under vertical MFs. In the last two columns ($t>64$) of Fig. 20, it becomes evident that the positions of bubbles during boiling undergo displacement in both scenarios—with and without MF. Particularly, such boiling pattern comprising vapor columns under strong vertical MF becomes unstable too. This displacement is primarily attributed to the non-uniform growth of bubbles generated by the initial random perturbation, leading to amplified interactions between bubbles over time. This phenomenon is depicted in Fig. 23, chosen with $N=11.1$ as an example since this behavior is observed consistently regardless of the presence of the MF. In Fig. 23, two adjacent bubbles (A and B) are selected for observation, and the evolution of their shapes and velocity streamlines in the vertical plane ($XOY$) at different time instants is presented. At $t=51.9$, the detachment behaviors of bubbles A and B differ, resulting in asymmetry in the flow field on either side of bubble B, as indicated by regions R1 and R2 in the figure. Vortices are present in region R1, while region R2 lacks vortices. As is well-known, regions with vortices exhibit lower pressure. Consequently, bubble B tends to approach bubble A. At $t=63.0$, the two bubbles are very close together, and there is almost no horizontal flow in region R1, while horizontal flow persists in region R2. This causes bubble B to approach bubble A further, eventually leading to their fusion and resulting in a chaotic flow state. Additionally, when bubbles merge, vapor overflows into new positions to form bubbles. The entire boiling system is no longer in a state where bubbles simultaneously detach and generate, but rather, the formation, detachment, and merging of bubbles occur simultaneously. Consequently, at this stage, the height fluctuation of the near-wall thin film is no longer significant. Therefore, the $Nu$ number variation curve shown in Fig. 22 exhibits small amplitude oscillations around the correlation proposed by Klimenko (1981) during the latter stage.

Figure 21: Distribution of the temperature field, the streamline patterns, the Lorentz forces, and the horizontal component of velocity ($u_{z}$, the positive direction: from left to right), in the vertical cross-section of a single bubble in multi-mode film boiling under different vertical MFs. (a) $N=11.1$, $t=47.2$ and (b) $N=100.1$, $t=37.0$.

Figure 22: Temporal evolution of the space-averaged $Nu$ under different vertical MF intensities in the case of the multi-mode film boiling. The dashed line represents the theoretical results of Klimenko (1981).

Figure 23: Interaction of two neighboring bubble sites and the evolution of the streamline patterns at different time instants from $t=51.9$ to $t=66.7$ under the vertical MF with $N=11.1$.

3.2.2 Effect of the horizontal MF

\begin{overpic}[width=349.96898pt,tics=1]{figure/mut-hor-shape.jpg} \put(-4.0,76.0){\large(a)} \put(-4.0,50.5){\large(b)} \put(-4.0,24.8){\large(c)} \end{overpic}

Figure 24: Temporal evolution of the bubble shapes during the multi-mode film boiling flow without and with the horizontal MFs. From the top to the bottom row ((a)-(c)): $N=2.8$ (B=0.1 T), $N=11.1$ (B=0.2 T), $N=100.1$ (B=0.6 T).

Fig. 24 presents the evolution of the liquid/vapor interface under three different horizontal MFs. In the initial stage of boiling flow, as shown in the first column of Fig. 24, the relatively weak interaction between bubbles results in the effect of the horizontal MF on each bubble in the multi-mode cases being consistent with the conclusion drawn for the single-mode situation. With an increase in MF strength, the boiling flow tends to become two-dimensional along the MF direction. For small MF strength ($N=2.8$ in Fig. 24), boiling still occurs in the form of bubble detachment. For large MF strength ($N=11.1$ and 100.1 in Fig. 24), boiling takes place in the form of a vapor sheet. Fig. 25, using $N=2.8$ and $N=11.1$ as examples, illustrates the distribution of the flow field, temperature, Lorentz forces, and vertical velocity component $(u_{x})$ in two cross-sections of a single bubble in the multi-mode case during the initial stage of film boiling. The results align with the single-mode situation (Fig. 12). In other words, the mechanism of the horizontal MF cases is consistent with that of single-mode situations. Additionally, the variation of $Nu$, as shown in Fig. 26, is also similar to the single-mode cases before $t=60.0$.

Figure 25: Distribution of the temperature field, the streamline patterns, the Lorentz forces, and the vertical component of velocity $(u_{x})$, in two vertical cross-sections ($XOY$: parallel to the MF and $XOZ$: perpendicular to the MF) of a single bubble in the multi-mode film boiling under the horizontal MFs. (a) $N=2.8$, $t=42.6$. (b) $N=11.1$, $t=31.5$.

Figure 26: Temporal evolution of the space-averaged $Nu$ under different horizontal MF intensities in the case of the multi-mode film boiling. The dashed line represents the theoretical results of Klimenko (1981).

Similar to the case with a vertical MF, as boiling progresses to a certain extent, the flow becomes chaotic due to the influence of bubble interactions, as shown in the last two columns ($t>57.0$) of Fig. 24. Moreover, the flow exhibits more chaotic behavior in horizontal MF cases. This is attributed to the fact that, under the horizontal MF, each single column of bubbles along the direction of the MF forms a vapor sheet. As boiling proceeds, the interaction between vapor sheets intensifies compared to bubbles and vapor columns, ultimately resulting in larger-scale rupture or fusion of vapor sheets, exerting a more significant influence on the entire flow system. Fig. 27 illustrates the evolution of streamline patterns in the vertical cross-section perpendicular to the MF under the condition of $N=100.1$. It can be observed that, as boiling progresses, the asymmetry of vortices on both sides of the four vapor sheets (A, B, C, D) becomes more pronounced. The vapor sheets are then twisted, eventually leading to the fusion of vapor sheets B, C, and D at the base and fragmentation at the top, as shown in Fig. 27 at $t=48.2$. The fragments produce small bubbles, making the flow system more chaotic. After fusion, new vapor sheets E and F are generated, as shown in Fig. 27 at $t=66.7$, corresponding to the sudden increase in $Nu$ number in Fig. 26. For other cases where $N<100.1$, the variation in $Nu$ exhibits smaller fluctuations in the later stages. This is because, in these situations, the merging and generation of bubbles or vapor sheets occur almost simultaneously, resulting in a smaller fluctuation in the distance between the interface and the overheated wall.

Figure 27: Dynamics of vapor sheets and the evolution of the streamline pattern for different time instants range from $t=22.2$ to $t=66.7$, with the horizontal MFs $(N=100.1)$.

Moreover, it is essential to note that bubbly flow in multi-mode film boiling eventually transitions into turbulence or pseudo-turbulence, making the accurate resolution of bubble coalescence and collisions a significant area of ongoing research. The complexity of these interactions is beyond the scope of the present work. For example, Volume of Fluid (VOF) methods often simplify the coalescence process (Scardovelli & Zaleski, 1999). Additionally, recent work by Innocenti et al. (2021) suggests that a spatial resolution of $\Delta=2d/Ar$ (where $Ar=\rho_{l}g^{1/2}d^{3/2}/\mu_{l}$, and $d$ is the typical size of the bubble) is sufficient for convergence in the spectra of liquid velocity fluctuations in bubbly flows. In our study, the multi-mode film boiling flows are characterized by an Archimedes number of $Ar\approx 120$. The spatial resolution at grid level 7, which corresponds to $\Delta=d/96$, meets the convergence criteria suggested by Innocenti et al. While this does not constitute a full grid convergence verification, it provides confidence in the resolution used to capture the essential physics of heat and mass transfer during multi-mode film boiling. Future research will focus on a more detailed study of vapor bubbly flows under external magnetic fields, including a comprehensive statistical analysis and examination of spectra of kinetic energy and velocity fluctuations on various grids, to rigorously address these complex interactions in the presence of external MFs.

4 Conclusions

Direct numerical simulations of the 3D film boiling process under magnetic fields (MFs) have been conducted using a numerical method based on the combined magnetohydrodynamics (MHD) and phase change model developed in our previous studies Zhao et al. (2022); Zhang & Ni (2018). The present study analyzes the effects of vertical and horizontal MFs on the morphology and heat transfer characteristics of liquid film boiling. The results demonstrate that the MF primarily suppresses vortices within the cross-section parallel to the MF direction. Consequently, the boiling pattern and heat transfer are altered depending on the MF direction.

In the presence of a vertical magnetic field, we observe that the flow structure remains isotropic in single-mode scenarios. As the strength of the MF increases, film boiling gradually transitions from periodic bubble detachment to a columnar jet formation. This transition occurs because the induced Lorentz forces counteract the horizontal velocity component, thereby suppressing the generation of vortices near the bubbles and inhibiting their detachment. Regarding the heat transfer characteristics corresponding to the change in boiling form, the time variation of the Nusselt number $(Nu)$ changes from approximately periodic oscillations to damped oscillations. This transition occurs because, during boiling in the form of a columnar jet, bubble detachment happens away from the wall, resulting in a relatively stable bottom vapor film. In terms of overall heat transfer efficiency, the impact of small MF strength in the bubble periodic detachment area on the space- and time-average Nusselt number ($\overline{Nu}$) is not substantial. However, under high MF strength in the vapor column boiling area, $\overline{Nu}$ increases with the rising MF strength. Additionally, a cliff-like sharp decrease in $\overline{Nu}$ occurs during the transition from the region of periodic bubble detachment to the columnar boiling region. This drastic change is attributed to the different boiling forms resulting in varying bubble detachment positions, consequently affecting the thickness of the bottom vapor film. Furthermore, it is observed that after removing the MF in the case of complete columnar jetting boiling, the boiling flow tended to revert to the unstable state observed without a magnetic field.

In the presence of a horizontal magnetic field, film boiling gradually transitions from periodic bubble detachment to a two-dimensional planar vapor film sheet as the MF strength increases. This shift occurs because the induced Lorentz force opposes the vertical velocity component within the cross-section parallel to the MF, suppressing vortex generation on both sides of the bubbles and preventing necking phenomena. Meanwhile, the flow field within the cross-section perpendicular to the MF remains largely unaffected, prompting the boiling flow to adopt a two-dimensional flow along the MF direction. Regarding heat transfer characteristics, the time variation of $Nu$ shifts from approximately periodic oscillations to nearly linear changes. This transition is attributed to the formation of a relatively stable vapor film during boiling in the form of a planar vapor film sheet. Concerning overall heat transfer efficiency, in the region of bubble periodic detachment, the initial effects of the horizontal MF lead to a decrease in $\overline{Nu}$ with increasing MF strength. However, $\overline{Nu}$ increases with rising MF strength in the planar vapor sheet boiling area. This is because, during boiling in the form of a planar vapor sheet, the Lorentz force acting vertically downward on the bottom vapor film becomes more pronounced with increasing MF strength, resulting in a thinner bottom vapor film thickness. Unlike the vertical MF, it was observed that under complete planar vapor film jetting boiling conditions, removing the MF still maintained the original state of the boiling flow.

In multi-mode film boiling, irrespective of whether it is under vertical or horizontal MF conditions, the early stages of boiling exhibit similar patterns observed in single-mode film boiling, with consistent effects of the MF on boiling morphology and heat transfer characteristics. However, once the boiling flow reaches complete development, the interaction between bubbles increases, leading to the complexity of the boiling flow and the emergence of numerous scattered small bubbles. Furthermore, under horizontal MF conditions, the flow structure tends to be more chaotic compared to vertical MF conditions. During this stage, bubble formation, detachment, and merging occur simultaneously, causing the $Nu$ to deviate from its original variation pattern with a significantly reduced amplitude.

5 Appendix A.

Here, we demonstrate that even when using more realistic physical properties for the liquid metal, the MHD effects on the rising characteristics of the vapor bubbles remain qualitatively consistent. As detailed in §2.4, the value of $Gr$ in the present study is set at $0.392$. For the additional test cases, we have simulated scenarios with $Gr$ values of 1.745 and 24.265, which are one and two orders of magnitude larger than the original value, respectively. Correspondingly, the grid resolutions were increased to 1536 $\times$ 512 $\times$ 512 and 3072 $\times$ 1024 $\times$ 1024, representing 8 and 64 times the grid numbers of the $Gr=0.392$ case. Figure 28 shows snapshots of the liquid/vapor interface in the $XOY$ vertical plane: panel (a) represents $Gr=0.392$, panel (b) represents $Gr=1.745$, and panel (c) represents $Gr=24.265$. In each case, the two sub-figures depict the first two bubble detachment cycles; the blue line represents the non-MHD case, while the red line shows the MHD case at $N=100.08$. The results indicate that even with higher $Gr$ numbers, the influence of the vertical magnetic field on film boiling remains consistent: the vertical magnetic field suppresses the Rayleigh-Taylor instability, causing a delay or even complete suppression of vapor bubble detachment. These additional simulations support the conclusion that the adjustments in physical properties do not significantly alter the fundamental physical processes observed in the original simulations.

Figure 28: Comparison of the vapor-liquid interface without (red line) and with (blue line) a vertical MF for (a) $Gr=0.392$, (b) $Gr=1.745$ and (c) $Gr=24.265$, while 1 and 2 are the time moments just after bubble detachment in the absence of a MF. (a-1) $t=18.9$ and (a-2) $t=37.6$; (b-1) $t=17.5$ and (b-2) $t=29.8$; (c-1) $t=17.9$ and (c-2) $t=35.8$.

Acknowledgement: The authors gratefully acknowledge the supports from National Key R $\&$ D Program of China (2023YFA1011000, 2022YFE03130000) and from the NSFC (12222208, 51927812). K.C.S. thanks IIT Hyderabad for the financial support through grant IITH/CHE/F011/SOCH1.

Declaration of interests: The authors report no conflict of interest.