1. Introduction
Boiling is widely employed in industrial heat-related processes due to its exceptional heat transfer efficiency (Dhir Reference Dhir1998; Tomar et al. Reference Tomar, Biswas, Sharma and Welch2008; Guion et al. Reference Guion, Afkhami, Zaleski and Buongiorno2018; Lavino et al. Reference Lavino, Smith, Magnini and Matar2021). 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 MF 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 (Reference Chang1957, Reference Chang1959) involved analysing the mechanism of film boiling on a horizontal surface using the theory proposed by Taylor (Reference Taylor1950) for the most unstable wavelength. Zuber (Reference Zuber1958) 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{\rm \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. Here, 
$\sigma$ denotes the interfacial tension, 
$g$ represents the acceleration due to gravity and 
$\rho_l$ and 
$\rho_v$ are the density of the liquid and vapour phases, respectively. Subsequently, Klimenko (Reference Klimenko1981) and Klimenko & Shelepen (Reference Klimenko and Shelepen1982) 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 (2-D) simulations and single-mode models, with the computational domain width close to the most dangerous wavelength 
$\lambda _d$. These 2-D 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 the morphology and heat transfer in film boiling. For high Jakob numbers, the rapid formation of vapour makes bubble detachment more difficult, resulting in the formation of tall and stable vapour columns (Son & Dhir Reference Son and Dhir1997; Juric & Tryggvason Reference Juric and Tryggvason1998; Guo et al. Reference Guo, Sun, Li and Tao2011). For low Morton numbers or high density ratios, it becomes more challenging for bubbles to pinch off (Juric & Tryggvason Reference Juric and Tryggvason1998; Akhtar & Kleis Reference Akhtar and Kleis2013; Mortazavi Reference Mortazavi2022). Some studies (Agarwal et al. Reference Agarwal, Welch, Biswas and Durst2004; Tomar et al. Reference Tomar, Biswas, Sharma and Welch2008; Guo et al. Reference Guo, Sun, Li and Tao2011) 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 Reference Esmaeeli and Tryggvason2004b; Hens, Biswas & De Reference Hens, Biswas and De2014; Ningegowda & Premachandran Reference Ningegowda and Premachandran2019) 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 vapour columns. The proximity between vapour 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 2-D film boiling simulations (Gibou et al. Reference Gibou, Chen, Nguyen and Banerjee2007; Akhtar & Kleis Reference Akhtar and Kleis2013; Zhang & Ni Reference Zhang and Ni2018). These studies observed that bubbles did not detach from the vapour film in their 2-D 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 2-D domains have consistently confirmed that the distribution of bubbles corresponds to the most dangerous wavelength in film boiling on a horizontal surface (Esmaeeli & Tryggvason Reference Esmaeeli and Tryggvason2004b; Tomar et al. Reference Tomar, Biswas, Sharma and Welch2008; Hens et al. Reference Hens, Biswas and De2014; Akhtar Reference Akhtar2018; Ningegowda & Premachandran Reference Ningegowda and Premachandran2019). The three-dimensional (3-D) numerical simulations of film boiling employing simplified liquid/vapour parameter models primarily focus on single-mode scenarios (Esmaeeli & Tryggvason Reference Esmaeeli and Tryggvason2004a; Tsui & Lin Reference Tsui and Lin2016; Zhang & Ni Reference Zhang and Ni2018; Kumar & Premachandran Reference Kumar and Premachandran2022). Moreover, limited research has been devoted to exploring the underlying physical mechanisms, with much of the existing work replicating earlier 2-D studies. These investigations often revisit aspects like the impact of overheat, Grashof number and density ratio (Khorram & Mortazavi Reference Khorram and Mortazavi2022). 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, 3-D multi-modal film-boiling simulations in existing literature have typically employed computational domains of sizes not exceeding 2
$\lambda _d$ (Shin & Juric Reference Shin and Juric2002; Esmaeeli & Tryggvason Reference Esmaeeli and Tryggvason2004b; Khorram & Mortazavi Reference Khorram and Mortazavi2022), 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, has recently emerged as the forefront of film-boiling investigation. Extensive experimental studies (Ogata & Yabe Reference Ogata and Yabe1993; Wang et al. Reference Wang, Lewin, Swaffield and Chen2009; Pearson & Seyed-Yagoobi Reference Pearson and Seyed-Yagoobi2013; Patel et al. Reference Patel, Seyed-Yagoobi, Sinha-Ray, Sinha-Ray and Yarin2016; Patel & Seyed-Yagoobi Reference Patel and Seyed-Yagoobi2017) and numerical simulations (Tomar et al. Reference Tomar, Gerlach, Biswas, Alleborn, Sharma, Durst, Welch and Delgado2007, Reference Tomar, Biswas, Sharma and Welch2009; Pandey, Dalal & Biswas Reference Pandey, Dalal and Biswas2015; Pandey, Biswas & Dalal Reference Pandey, Biswas and Dalal2016, Reference Pandey, Biswas and Dalal2017; Rouzbahani & Mortazavi Reference Rouzbahani and Mortazavi2023; Wang et al. Reference Wang, Pérez, Huang, Guan and Wu2023) have shown that electric fields aid in accelerating bubble detachment during the boiling process, thereby enhancing heat transfer efficiency. Furthermore, numerical simulations focusing on 2-D multi-model film boiling under electric fields have found that the electric field promotes the occurrence of Rayleigh–Taylor instability at the vapour–liquid interface and increases the number of bubble nucleation sites. However, as another crucial physical phenomenon, MFs offer considerable research potential, particularly in confined nuclear fusion. The effect of a MF on fluids can be divided into 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 has focused more extensively on FHD. This is mainly because ferrofluids, compared with magnetic fluids such as liquid lithium, have decent transparency, which brings convenience to experimental investigations. A plethora of experimental (Lee, Lee & Jeong Reference Lee, Lee and Jeong2012; Shojaeian et al. Reference Shojaeian, Yildizhan, Coşkun, Ozkalay, Tekşen, Gulgun, Acar and Koşar2016; Abdollahi, Salimpour & Etesami Reference Abdollahi, Salimpour and Etesami2017; Zonouzi, Aminfar & Mohammadpourfard Reference Zonouzi, Aminfar and Mohammadpourfard2019) and numerical (Aminfar, Mohammadpourfard & Sahraro Reference Aminfar, Mohammadpourfard and Sahraro2012; Aminfar, Mohammadpourfard & Maroofiazar Reference Aminfar, Mohammadpourfard and Maroofiazar2014; Malvandi Reference Malvandi2016; Guo et al. Reference Guo, Li, Feng, Wang and Zhao2019) 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 MHD effects on film boiling, even today, primarily relies on early theoretical and experimental studies. Faber & Hsu (Reference Faber and Hsu1967) 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, Lloyd & MacPherson (Reference Fraas, Lloyd and MacPherson1974) 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 (Reference Wagner and Lykoudis1981) studied the influence of a horizontal MF on nucleate boiling of liquid metallic mercury, finding that the horizontal MF reduced the heat transfer efficiency of boiling and promoted the transition from nucleate to film boiling. Takahashi, Inoue & Kaneko (Reference Takahashi, Inoue and Kaneko1994) 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, Leonardi & Lykoudis (Reference de Bertodano, Leonardi and Lykoudis1998) showed that the MF slowed down the frequency of bubble release during nucleate boiling. The theoretical research conducted by Arias (Reference Arias2010a,Reference Ariasc) corroborated the findings from the aforementioned experimental investigations. However, theoretical investigation of Arias (Reference Arias2010b) within the Taylor–Helmholtz instability framework indicated that the MF did not influence film boiling. The aforementioned literature review suggests that the impact of MFs 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 MFs 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 behaviour of film boiling under the influence of applied MFs in different directions. The scarcity of research in this area can be attributed to the increased significance of MHD effects in 3-D 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 MHD 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, Zhang & Ni Reference Zhao, Zhang and Ni2022), which employed a sharp interface method to enforce the temperature and concentration conditions at the liquid/vapour interface, integrating a volume-of-fluid method and an embedded boundary approach for sharp scheme construction. Zhao et al. (Reference Zhao, Zhang and Ni2022) 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 Reference Zhang and Ni2014b) for discretizing the MHD equations, enabling the simulation of 3-D film boiling on a horizontal surface in the presence of an applied MF. Our investigation comprises two primary components. Firstly, we focus on single-mode film boiling under vertical and horizontal MFs in three dimensions. In contrast to previous studies, we significantly extend the simulation duration to observe the impact of the MF 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 analyse the behaviour of multi-mode film boiling under both vertical and horizontal MFs.
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. (Reference Zhao, Zhang and Ni2022). The results are discussed in § 3, where we analyse the influence of horizontal and vertical MFs 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 3-D 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 vapour 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 vapour 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 figure 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 (Reference Tsui and Lin2016) that setting the computational domain width as 
$\lambda _d$ in 3-D 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 centre of the coordinate system is aligned with the centre of the bottom plane 
$(yz)$. Gravity 
$\boldsymbol {g}$ acts in the negative 
$x$ direction. An external uniform MF is applied either in the vertical (
$x$) or horizontal (
$y$) direction.
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.
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_{{sat}}$), which is a valid assumption except when exposed to extremely high heat flux (Juric & Tryggvason Reference Juric and Tryggvason1998; Akhtar & Kleis Reference Akhtar and Kleis2013). Additionally, the temperature of the overheated surface is kept constant at 
$T_{{wall}}$ (
$T_{{wall}} = T_{{sat}} + \Delta T$). Initially (at time 
$t=0$), the interface is defined using sinusoidal functions, and the temperature in the vapour phase linearly decreases from the overheated surface to the phase interface.
2.2. Governing equations
The viscous incompressible MHD two-phase flow, incorporating liquid–vapour phase change, is governed by the mass, momentum and energy conservation equations. They are expressed as follows:
\begin{gather} \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u}=\ddot{m} \left(\frac{1}{\rho_{v}}-\frac{1}{\rho_{l}}\right), \end{gather}
\begin{gather} \rho\left(\frac{\partial \boldsymbol{u}}{\partial t}+\boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u}\right)={-}\boldsymbol{\nabla} p+\boldsymbol{\nabla} \boldsymbol{\cdot}\left(\mu\left(\boldsymbol{\nabla} \boldsymbol{u}+\boldsymbol{\nabla} \boldsymbol{u}^{T}\right)\right)+\rho \boldsymbol{g}+\sigma \kappa \delta_{s} \boldsymbol{n}+ \boldsymbol{F_l}, \end{gather}
\begin{gather} \frac{\partial T}{\partial t}+\boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} T=\boldsymbol{\nabla} \boldsymbol{\cdot}\left(\frac{k}{\rho c_{p}} \boldsymbol{\nabla} T\right)+\varPhi. \end{gather}
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 vapour phases, respectively. In (2.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
\begin{equation} \boldsymbol{F_l}=\boldsymbol{J} \times \boldsymbol{B}, \end{equation}
where 
$\boldsymbol {J}$ is the induced current density and 
$\boldsymbol {B}$ represents the applied MF. According to the Ohm's law, the current density is given by
\begin{equation} \boldsymbol{J}=\sigma_{e}(-\boldsymbol{\nabla} \varphi+\boldsymbol{u} \times \boldsymbol{B}), \end{equation}
where 
$\sigma _{e}$ denotes the electric conductivity and 
$\varphi$ represents the induced electric potential. In (2.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 L u'$, with 
$\eta$ the magnetic permeability, 
$L$ the characteristic length and 
$u'$ 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
\begin{equation} \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{J}=0. \end{equation}By combining (2.5) and (2.6), the electric potential Poisson equation is derived as
\begin{equation} \boldsymbol{\nabla} \boldsymbol{\cdot}\left(\sigma_{e} \boldsymbol{\nabla} \varphi\right)=\boldsymbol{\nabla} \boldsymbol{\cdot}\left(\sigma_{e} \boldsymbol{u} \times \boldsymbol{B}\right). \end{equation}
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 (2.5). Finally, the Lorentz force 
$\boldsymbol {F_l}$ can be determined. In (2.3) governing the temperature transport, 
$c_{p}$ represents the volumetric heat capacity, 
$k$ is the thermal conductivity, 
$T$ denotes the temperature and 
$\varPhi$ stands for the internal dissipation term, encompassing viscous dissipation and Joule dissipation. However, in this study, both dissipation effects are considered negligible compared with the temperature terms, a simplification consistent with other studies (Burr & Müller Reference Burr and Müller2002; Vogt et al. Reference Vogt, Ishimi, Yanagisawa, Tasaka, Sakuraba and Eckert2018; Chen, Yang & Ni Reference Chen, Yang and Ni2024).
The presence of phase change necessitates that the above equations must satisfy the following interfacial jump conditions:
\begin{gather} \dot{m}=\rho_{l}\left(\boldsymbol{u}_{\varGamma}-\boldsymbol{u}_{l}\right) \boldsymbol{\cdot} \boldsymbol{n}=\rho_{v }\left(\boldsymbol{u}_{\varGamma}-\boldsymbol{u}_{v}\right) \boldsymbol{\cdot} \boldsymbol{n}, \end{gather}
\begin{gather} \dot{m}=\frac{\left(\boldsymbol{q}_{l}-\boldsymbol{q}_{v}\right) \boldsymbol{\cdot} \boldsymbol{n}}{h_{lg}}=\frac{\left((k \boldsymbol{\nabla} T)_{l}-(k \boldsymbol{\nabla} T)_{{v}}\right) \boldsymbol{\cdot} \boldsymbol{n}}{h_{lg}}. \end{gather}
Equation (2.8) describes the mass transfer rate from liquid to vapour in accordance with the mass conservation law. Note that 
$\boldsymbol {u}_{\varGamma }$ is the propagation velocity of the interface, which does not equal the local fluid velocity 
$\boldsymbol {u}_{(l, v)}$ in the presence of phase change. Equation (2.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 vapour and 
$h_{lg}$ is the latent heat. Then the velocity jump condition (2.8) leads to the following divergence-constraint condition on the velocity field:
\begin{equation} \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u} = \left(\frac{1}{\rho_v} - \frac{1}{\rho_l}\right)\ddot{m}. \end{equation} The relationship between 
$\dot {m}$ and 
$\ddot {m}$ is given by
\begin{equation} V \ddot{m}=\dot{m} S_{\varGamma}, \end{equation}
where 
$V$ denotes the volume of the discretized cell containing the interface and 
$S_{\varGamma }$ 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
\begin{equation} T_{\varGamma}=T_{sat}, \end{equation}which effectively serves as a Dirichlet-like boundary condition for the solution of the temperature equation in the separated liquid and vapour 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 Reference Popinet2009, Reference Popinet2015), 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 (Scardovelli & Zaleski Reference Scardovelli and Zaleski1999; Weymouth & Yue Reference Weymouth and Yue2010) accurately captures and advances the phase interface. The Bell–Colella–Glaz second-order scheme (Bell, Colella & Glaz Reference Bell, Colella and Glaz1989) discretizes the advection term, and a fully implicit scheme is applied to discretize the diffusion term. Additionally, a quad/octree-based adaptive mesh refinement technique is utilized for spatial discretization, enabling grid refinement primarily near the interface to optimize computational resources.
In our previous study (Zhao et al. Reference Zhao, Zhang and Ni2022), 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 MF, we utilize a consistent and conservative scheme, which was also developed earlier (Zhang & Ni Reference Zhang and Ni2014a,Reference Zhang and Nib), 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 (Zhang & Ni Reference Zhang and Ni2014a,Reference Zhang and Nib; Zhao et al. Reference Zhao, Zhang and Ni2022). 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$:
(i) A split advection method proposed by Weymouth & Yue (Reference Weymouth and Yue2010) is employed to discretize the VOF advection equation, and the volume fraction
$C^{n+{1}/{2}}$ is updated accordingly. However, the presence of phase change alters the form of the VOF advection equation to
(2.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\begin{equation} \frac{\partial C}{\partial t}+\boldsymbol{\nabla} \boldsymbol{\cdot}\left(\boldsymbol{u} C \right)={-}\frac{\ddot{m}}{\rho_l}, \end{equation}$\boldsymbol {u}_E$ is introduced, $\boldsymbol {u}_E=\boldsymbol {u}-\boldsymbol {u}_J$. Here, $\boldsymbol {u}_J$ is the velocity jump, which could be obtained by solving for the velocity potential $\psi$ using the following equations:
(2.14)To avoid the negative effect of the source term, the interfacial velocity (
\begin{equation} \left. \begin{gathered} \nabla^{2} \psi=\left(\frac{1}{\rho_{v }}-\frac{1}{\rho_{l}}\right) \ddot{m}, \\ \boldsymbol{u}_{J}=\boldsymbol{\nabla} \psi. \end{gathered} \right\}\end{equation}
$\boldsymbol {u}_{\varGamma }$) is chosen to advect the VOF advection equation
(2.15)where\begin{equation} \frac{\partial C}{\partial t} + \boldsymbol{u}_{\varGamma} \boldsymbol{\cdot} \boldsymbol{\nabla} C =0, \end{equation}$\boldsymbol {u}_{\varGamma }=\boldsymbol {u}_E- ({\ddot {m}}/{\rho _l})\boldsymbol {n}$, wherein $\boldsymbol {n}$ directs towards the vapour phase in the normal direction at the interface.(ii) The physical properties
$\phi ^{n+{1}/{2}}$ are updated using
(2.16)\begin{equation} \phi^{n+{1}/{2}}=C^{n+{1}/{2}} \phi_{l}+\left(1-C^{n+{1}/{2}}\right) \phi_{v}, \quad {\rm where} \ \phi=\rho, \mu, \sigma_{e}. \end{equation}(iii) The temperature
$T^{n+{1}/{2}}$ is then updated. Note that the temperature fields of the liquid and vapour 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 vapour phase for saturated film boiling, maintaining the liquid temperature at $T_\varGamma$. The interfacial temperature ($T_{\varGamma }$) 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 (2.15)
(2.17)
\begin{align}
\frac{T^{n+{1}/{2}}-T^{n-{1}/{2}}}{\Delta
t}+\boldsymbol{\nabla} \boldsymbol{\cdot}(T
\boldsymbol{u})^{n+{1}/{2}}=\boldsymbol{\nabla}
\boldsymbol{\cdot}\left(\left(\frac{\lambda}{\rho
c_{p}}\right) \boldsymbol{\nabla}_{f}
T\right)^{n+{1}/{2}}+\left(\frac{\ddot{m}}{\rho}
T_{\varGamma}\right)^{n+{1}/{2}}.
\end{align}
(iv) The Lorentz force
$\boldsymbol {F_l}^{n+{1}/{2}}$ is updated by solving the electric potential Poisson equation and Ohm's law using
(2.18)\begin{gather} \boldsymbol{\nabla} \boldsymbol{\cdot}\left(\sigma_{e} \boldsymbol{\nabla} \varphi \right)^{n+{1}/{2}}=\boldsymbol{\nabla} \boldsymbol{\cdot}\left(\sigma_{e}(\boldsymbol{u} \times \boldsymbol{B}) \right)^{n+{1}/{2}}, \end{gather}(2.19)\begin{gather} \boldsymbol{J}^{n+{1}/{2}}=\sigma_{e}^{n+{1}/{2}}\left(-\boldsymbol{\nabla} \varphi+\boldsymbol{u} \times \boldsymbol{B}\right)^{n+{1}/{2}}, \end{gather}(2.20)\begin{gather} \boldsymbol{F_l}^{n+{1}/{2}}=\boldsymbol{J}^{n+{1}/{2}} \times \boldsymbol{B}^{n+{1}/{2}}. \end{gather}(v) The mass transfer rate
$\dot {m}^{n+{1}/{2}}$ is updated according to (2.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 vapour phase. This method avoids errors caused by artificially thickening the interface and ensuring accurate mass transfer rate calculations.(vi) 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 Reference Popinet2009, Reference Popinet2015).
2.4. Dimensionless parameters
In this study, the length, velocity and time scales are defined as 
$\lambda '=\sqrt {\sigma /(\rho _l-\rho _v)g}$, 
$t'=\sqrt {\lambda '/g}$ and 
$u'=\sqrt {\lambda 'g}$. The critical wavelength is given by 
$\lambda _c=2{\rm \pi} \lambda '$. 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 '^3/\mu _v^2)$ describes the ratio of the buoyancy to viscous force. The Prandtl number 
$(Pr=c_{p,v}\mu _v/k_v)$ represents the ratio of the momentum diffusion to the thermal diffusion in vapour. 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
\begin{equation} Nu=\left.\frac{1}{W_y} \frac{1}{W_z} \int_{0}^{W_y} \int_{0}^{W_z} \frac{\lambda^{\prime}}{T_{{wall }}-T_{{sat }}} \frac{\partial T}{\partial n}\right|_{{wall }} \,\mathrm{d} y \,\mathrm{d} z. \end{equation}Further, the space- and time-averaged Nusselt number can be calculated as
\begin{equation} \overline{Nu}=\frac{1}{t_{1}-t_{0}} \int_{t_{0}}^{t_{1}} Nu \,\mathrm{d} t, \end{equation}
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
\begin{equation} N=\sigma_{e,l}B^2\lambda'/\rho_l\sqrt{g\lambda'}. \end{equation} 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'\lambda '/\mu _l$. In this study, we primarily use the interaction number (
$N$) to characterize the strength of the external 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 vapour film which has a thickness of
\begin{equation} \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}\end{equation}
(Berenson Reference Berenson1961; Kim, Lee & Kim Reference Kim, Lee and Kim2015). 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 vapour film thickness cannot be resolved even when grid resolution is 512 grids within a most dangerous wavelength, which represents the diameter of the vapour bubble in single-mode film boiling. It implies more than fifty million grids are required even if we use an 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 vapour phase so that the vapour 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 vapour layer. Furthermore, in the Appendix, we demonstrate that, even when using more realistic physical properties, the MHD effects on the rising behaviour of the vapour 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.
Table 1. Physical properties of fluids used 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_{{sat}}$, given that 
$(T-T_{{sat}})/(T_{{wall}} - T_{{sat}})$, standardizing the relative magnitude of the vapour 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.
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:
\begin{equation} x=\frac{\lambda_d}{128} \left( 8.0 + \cos\left( \frac{2 {\rm \pi}y}{\lambda_d} \right)+ \sin\left( \frac{2 {\rm \pi}z}{\lambda_d} \right) \right). \end{equation}
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 vapour bubble shapes on the vertical cross-section (
$XOY$ plane) at 
$t=14.8$ and 
$t=16.7$ are illustrated in figures 2(a) and 2(b), respectively. The temporal evolution of the space-averaged Nusselt number (
$Nu$) is also depicted in figure 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.
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.
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 MF, 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 MF 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 the 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 the Courant–Friedrichs–Lewy (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.
3.1.1. Film boiling without a MF
We begin the presentation of our results by analysing the film-boiling phenomenon without MF and compare with existing boiling mechanisms and empirical models. Figure 3 shows the temporal evolution of the space-averaged 
$Nu$. It can be seen that 
$Nu$ oscillates periodically around the correlation proposed by Klimenko (Reference Klimenko1981), given that 
$Nu_k=0.19(Gr\,Pr)^{{1}/{3}}\times (0.89Ja^{-{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/vapour interface for the fifth bubble cycle (from 
$t=51.9$ to 
$t=63.0$), corresponding to the blue region in figure 3, is depicted in figure 4, while figure 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 figure 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 vapour 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 vapour film thickness initially increases due to the expansion of the vapour 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 (2.21). The temperature profile reveals a linear distribution from the overheated wall to the liquid/vapour interface throughout the film boiling process. A smaller distance between the liquid/vapour 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 Reference Juric and Tryggvason1998; Zhang & Ni Reference Zhang and Ni2018). Next, we investigate the influence of a vertical MF on the film-boiling dynamics.
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 (Reference Klimenko1981).
Figure 4. Vapour bubble generation during the film-boiling flow without a MF at different dimensionless times, which correspond to the blue region in figure 3. Panels show (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 figure 4(a–f).
3.1.2. Effect of the vertical MF
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\,\textrm {T}}$, 
${B = 0.1\,\textrm {T}}$, 
${B = 0.15\,\textrm {T}}$, 
${B = 0.2\,\textrm {T}}$ and 
${B = 0.4\,\textrm {T}}$. According to (2.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 dominant to the Lorentz force being dominant.
To illustrate the MHD effect, we compare the liquid/vapour interface patterns within four bubble generation cycles, and the results are presented in figure 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 vapour 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 vapour jet arises, and a small bubble is released from the top of the jet, corresponding to the onset of vapour 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 vapour 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.
Figure 6. Temporal evolution of the bubble shapes during the film-boiling flow without and with the vertical MFs. Note that in the pictures, the long arrow represents the direction of the MF, while the same holds in the following figures of the present study. Panels show (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).
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 figure 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 figure 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 figure 7(b i), 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 figure 7(c i). 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.
Figure 7. (a–c) Comparison of the temperature field and streamline patterns in the vertical cross-section (
$XOY$) under different vertical MFs. Panels (b ii,c ii) show 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. Panels show (a) 
$N=0$, 
$t=57.4$; (b) and (d) 
$N=11.1$, 
$t=53.7$; (c) and (e) 
$N=44.5$, 
$t=40.7$.
Now, we aim to elucidate the mechanisms underlying the aforementioned phenomenon. First, the distribution of the current line surrounding the liquid/vapour interface is plotted in figure 7(d,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. Figure 7(b ii,c ii) presents 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 vapour 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 (figure 7b), 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 colour scale results in figure 7(a–c). This explains the suppression of vortex formation, and such decay in vortices counteracts the pressure difference across the vapour–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 vapour converges towards the bulge, ultimately leading to the formation of a vapour jet under the influence of buoyancy. For vapour jets, bubble detachment is triggered by the amplification of instabilities and surface waves (Lin & Reitz Reference Lin and Reitz1998). Then, in the scenario with a large MF intensity (
$N=44.5$, figure 7c), the Lorentz forces entirely suppress the vortices and lead to uniformly distributed 
$u_y$, with the flow directing toward the centre of the vapour 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 vapour 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.
The temporal evolution of the space-averaged Nusselt number (
$Nu$) under different vertical MFs is shown in figure 8(a). It can be observed that, as the vertical MF intensity increases, 
$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 behaviour. Specifically, for cases with 
$N=6.3$, 
$N=11.1$ and 
$N=44.5$, the initially oscillating 
$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 figure 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 vapour phase exhibits a nearly linear distribution, the temperature distribution remains unchanged, as observed in figures 6 and 7.
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 (Reference Klimenko1981). (b) Variation of the space- and time-averaged 
$Nu$ (
$\overline {Nu}$) with vertical MF intensity in the case of single-mode film boiling. The dimensional parameter 
$B$ is used on the horizontal axis to enhance clarity in presenting the value distribution.
Figure 8(b) illustrates the space- and time-averaged Nusselt number (
$\overline {Nu}$) calculated using (2.22) for different MF intensities. In cases of constant amplitude oscillation (
$N=0$, 
$N=0.7$ and 
$N=2.8$), 
$Nu$ is averaged over time from the second trough to the seventh trough in figure 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-hand side of the dashed line represents the periodic bubble detachment area, where 
$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-hand side of the dashed line is the vapour column boiling area, where 
$Nu$ exhibits damped oscillation. The value of 
$\overline {Nu}$ sharply declines in the transition zone. After that, 
$\overline {Nu}$ gradually increases with the increase in MF intensity.
To understand the reason behind the increase in 
$\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 figure 9(a). Note that, in the same figure, the embedded panel illustrates the vapour interface in the vapour film region (
$y > 0.2$), and figure 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 vapour thickness decreases with the increase in MF intensity, leading to an increase in 
$Nu$ and 
$\overline {Nu}$. Now, the question arises: Why does the vapour thickness decrease by generating a vapour jet as the MF intensity rises? Figure 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 figure 10b–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 vapour generated through evaporation to flow quickly from the near-wall region into the vapour jet, resulting in a smaller minimum distance between the vapour interface and the overheated wall. Analysing the vertical component of the velocity contour (
$u_x$) plot (figure 10b–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 vapour 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 vapour enters the centre of the jet. The higher the bubble detachment position, the greater the vapour velocity in the jet, making vapour overflow more quickly. Thus, a higher MF intensity decreases the vapour film thickness near the overheated wall, promoting 
$Nu$ and 
$\overline {Nu}$.
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/vapour 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 MF 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$.
The abrupt cliff-like decrease in 
$\overline {Nu}$ during the transition region in figure 8(b) can be explained through the same analytical method. Figure 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 vapour 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 vapour overflow velocity, as shown by the colour scale in figure 10(e). This facilitates the rapid influx of vapour generated by evaporation into the bulge, significantly reducing the vapour film thickness during the bubble necking stage. As a result, 
$Nu$ increases, enhancing 
$\overline {Nu}$.
3.1.3. Effect of the horizontal MF
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$ 
$\textrm {(B = 0.05 T)}$, 
$N=2.8$ 
${(B = 0.1\,\textrm {T})}$, 
$N=6.3$ 
${(B = 0.15\,\textrm {T})}$, 
$N=11.1$ 
${(B = 0.2\,\textrm {T})}$ and 
$N=44.5$ 
${(B = 0.4\,\textrm {T})}$.
Figure 11 depicts the evolution of the liquid/vapour 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 2-D boiling flow forms and the vapour features as a thin sheet along the direction of the MF after complete development. Furthermore, after the vapour 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 vapour sheet. In the study of sheet instability (Hagerty & Shea Reference Hagerty and Shea1955), two modes of instability for sheets have been identified: the sinuous mode and the varicose mode. From figure 11, it is evident that the flow corresponds to the sinuous mode.
Figure 11. Temporal evolution of the vapour bubble shapes during the film-boiling flow without and with the horizontal MFs. Panels show (a) 
$N=2.8$ 
${(B = 0.1\,\textrm {T})}$, (b) 
$N=6.3$ 
${(B = 0.15\,\textrm {T})}$, (c) 
$N=11.1$ 
${(B = 0.2\,\textrm {T})}$ and (d) 
$N=44.5$ 
${(B = 0.4\,\textrm {T})}$.
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 figures 12(a) and 12(b), respectively. The corresponding distribution of Lorentz forces and vertical velocity component (
$u_x$) are presented in figure 12(c,d). Figure 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 figure 12(a,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 figure 12(c,d), explains the disappearance of the vortex perpendicular to the MF. From figure 12(c i,d i), 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 vapour. Eventually, a columnar bulge along the MF direction is formed. In contrast, for the 
$XOZ$ plane perpendicular to the MF (see figure 12(c ii,d ii), 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 vapour sheet. These conclusions align with the theoretical prediction of MHD effects (Sommeria & Moreau Reference Sommeria and Moreau1982; Davidson Reference Davidson1995), 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 figure 12(a,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 figure 12(c,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 vapour entering the bulge, making heat transfer more challenging. Consequently, the internal temperature becomes lower.
Figure 12. (a,b) Comparison of temperature field and streamline patterns in two vertical cross-sections (
$XOY$: parallel to the MF, in (a i,b i); 
$XOZ$: perpendicular to the MF, in (a ii,b ii)) under different horizontal MF 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. Panels show (a,c) 
$N=6.3$, 
$t=27.8$; (b,d,e) 
$N=11.1$, 
$t=27.8$.
Figure 13(a) presents the temporal evolution of 
$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 vapour film's bottom, leading to an increased vapour 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 vapour 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 (Reference Klimenko1981). (b) Variation of the space- and time-averaged 
$Nu$ (
$\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.
The value of 
$\overline {Nu}$ for the aforementioned cases is plotted in figure 13(b). Similar to the vertical MF, for small MF intensity cases 
$(N<6.3)$, the average of 
$Nu$ from the second trough to the seventh trough in figure 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-hand side represents the periodic detachment region of the bubble. In this region, 
$\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 
$Nu$, both attributed to the appearance of the elongated bulge in the vapour film. The right-hand side is the vapour sheet boiling region, where 
$\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 vapour sheet (figure 14a) and when the vapour sheet stably exists (figure 14b) 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 vapour 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 figure 12(c,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 vapour 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 figure 14, as depicted in figure 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 vapour film near the overheated wall is stronger, resulting in a smaller thickness of the vapour film and a larger 
$\overline {Nu}$.
Figure 14. Comparison of the bubble shapes in the vertical cross-section (
$XOZ$ plane) perpendicular to the MF 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/vapour interface under different MF intensities corresponding to figure 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.
After examining the MHD effect on the growth of vapour bubbles, it is useful to compare our findings with those of previous studies. The most relevant work is by Chandrasekhar (Reference Chandrasekhar2013), who used linear stability analysis to investigate the effects of vertical and horizontal MFs 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} k B_{0}(\sqrt {\rho _l/\rho _0}+\sqrt {\rho _v/\rho _0})+n(k^{2} B_{0}^{2}-(\rho _l-\rho _v)/(\rho _l+\rho _v) g k)=2 (\rho _l-\rho _v)/(\rho _l+\rho _v) g k^{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 vapour bubble grows more slowly under a stronger vertical MF. However, the Rayleigh–Taylor (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 
$\boldsymbol {B}$-lines are unaffected by the horizontal MF, while disturbances parallel to the 
$\boldsymbol {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 
$\boldsymbol {B}$-lines but has little impact in the direction transverse to the 
$\boldsymbol {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 (vapour column boiling for vertical MF scenario and vapour sheet boiling for horizontal MF scenario).
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 figure 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 vapour bubble. The comparison of the temporal evolution of 
$Nu$ with the MF and upon its removal is illustrated in figure 17(a). After removing the vertical MF, a noticeable decrease in 
$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 vapour column (see figure 16). As analysed in § 3.1.2, the bubble detachment behaviour inhibits the vapour influx into the columnar bulge, causing the near-wall vapour film to thicken and reducing 
$Nu$. Figure 17(b) compares the evaporation rate (here, we measure the evaporation rate using the mass flow at the interface (
$\dot {m} S_{\varGamma }$)) 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 vapour film. Note that the oscillation of 
$\dot {m} S_{\varGamma }$ after removing the MF is caused by the reappeared detachment of the vapour 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 vapour for stable columnar bulge.
Figure 16. Temporal evolution of the temperature profile (colour contours) and velocity streamline pattern (black arrow lines) in the vertical cross-section (
$XOY$) after removing the MF at 
$t=44.4$ for 
$N=44.5$. Panels show (a) 
$t=44.4$, (b) 
$t=50.0$,(c) 
$t=51.9$, (d) 
$t=54.6$, (e) 
$t=56.5$ and (f) 
$t=59.3$. It can be seen that the removal of the MF leads to vapour 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_{\varGamma }$) 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 MF removal, is due to the reappearance of vapour bubble detachment.
Now, we examine the scenario based on a horizontal MF with the intensity of 
$N=22.2$. When film boiling reaches a relatively stable vapour sheet boiling at 
$t=57.4$, the MF is suddenly removed. However, intriguingly, we observe a very interesting result: film boiling in the vapour sheet state does not appear to be affected after removing the horizontal MF, as depicted in figure 18. Furthermore, by comparing 
$Nu$ and the evaporation rate (figure 19), it is observed that, with and without the presence of a horizontal MF, both 
$Nu$ and the evaporation rate are almost identical. This indicates that the vapour sheet formed under a horizontal MF is initially stable, in contrast to the vapour column generated under a vertical MF. To explain this, we should remember that the vertical MF eliminates the vortices around the vapour bubble in an axisymmetric manner, and such a controlling effect leads to the generation of an isotropic vapour column. Therefore, the axisymmetric Lorentz forces, directing radial-wards, maintain the force balance around the vapour 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 vapour 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, figure 12(c,d) reveals that, corresponding to the stable vapour 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 2-D vapour sheet. It indicates that, for the vapour sheet scenario under horizontal MF, the Lorentz forces play little role in the force balance in the radial direction, and thus, removing it causes little influence on the vapour pattern.
Figure 18. (a–c) The temporal evolution of the liquid/vapour interface after removing the horizontal MF at 
$t=44.4$ for 
$N=44.5$. (d–f) Temperature profile (colour 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_{\varGamma }$) with the horizontal MF and after the removal of the MF at 
$t=57.4$ for 
$N=22.2$. We observe that, in (a,b), the values are almost identical before and after removing the horizontal MF.
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:
\begin{align} x=\frac{\lambda_d}{16}+0.05 / n \sum_{i=1}^{n} A(i)\left(\cos \left( \frac{2 {\rm \pi}i y}{3\lambda_d} \right)+\sin \left( \frac{2 {\rm \pi}i z}{3\lambda_d} \right)+\cos \left( \frac{2 {\rm \pi}i y}{3\lambda_d} \right)+\sin \left( \frac{2 {\rm \pi}i z}{3\lambda_d} \right) \right), \end{align}
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\, \textrm {T}}$, 
${B = 0.1\,\textrm {T}}$, 
${B = 0.2\,\textrm {T}}$, 
${B = 0.4\,\textrm {T}}$ and 
${B = 0.6 \,\textrm {T}}$. The corresponding interaction numbers (Hartmann numbers) for both vertical and horizontal MFs 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
Figure 20 depicts the evolution of the liquid/vapour interface under three different vertical MFs. In the initial stage 
$(t<60.0)$ of film boiling, as observed in the first two columns of the figure, the distribution of bubbles aligns with the theory proposed by Zuber (Reference Zuber1958) 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 vapour jet releasing small bubbles from the top. This evolution ultimately leads to vapour columnar boiling without any bubble detachment, a pattern confirmed by figure 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 figure 7 reveals that, in multi-mode boiling, the effect of the vertical MF aligns with the observations in single-mode boiling. The evolution of 
$Nu$ also mirrors the single-mode cases in the early stages. As the MF 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 figure 22.
Figure 20. Temporal evolution of the bubble shapes during multi-mode film boiling without and with the vertical MFs. Panels show (a) 
$N=0$ (
$B=0.0$ T), (b) 
$N=11.1$ (
$B=0.2$ T) and (c) 
$N=100.1$ (
$B=0.6$ T).
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. Panels show (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 (Reference Klimenko1981).
Next, let us focus on the later stages of the film-boiling process under vertical MFs. In the last two columns (
$t > 64$) of figure 20, it becomes evident that the positions of bubbles during boiling undergo displacement in both scenarios – with and without MF. Particularly, such a boiling pattern comprising vapour 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 figure 23, chosen with 
$N=11.1$ as an example since this behaviour is observed consistently regardless of the presence of the MF. In figure 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 behaviours 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, vapour 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$ variation curve shown in figure 22 exhibits small amplitude oscillations around the correlation proposed by Klimenko (Reference Klimenko1981) during the latter stage.
Figure 23. Interaction of two neighbouring 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
Figure 24 presents the evolution of the liquid/vapour interface under three different horizontal MFs. In the initial stage of boiling flow, as shown in figure 24(a i,b i,c i), 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 figure 24), boiling still occurs in the form of bubble detachment. For large MF strength (
$N=11.1$ and 100.1 in figure 24), boiling takes place in the form of a vapour sheet. Figure 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 (figure 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 figure 26, is also similar to the single-mode cases before 
$t=60.0$.
Figure 24. Temporal evolution of the bubble shapes during the multi-mode film-boiling flow without and with the horizontal MFs. Panels (a–c) show 
$N=2.8$ 
${(B=0.1\,\textrm {T})}$, 
$N=11.1$ 
${(B=0.2\,\textrm {T})}$ and 
$N=100.1$ 
${(B=0.6\,\textrm {T})}$, respectively.
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. Panels show (a) 
$N=2.8$, 
$t=42.6$ and (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 multi-mode film boiling. The dashed line represents the theoretical results of Klimenko (Reference Klimenko1981).
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 (
$t > 57.0$) figure 24(a iii,a iv,b iii,b iv,c iii,c iv). Moreover, the flow exhibits more chaotic behaviour 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 vapour sheet. As boiling proceeds, the interaction between vapour sheets intensifies compared with bubbles and vapour columns, ultimately resulting in larger-scale rupture or fusion of vapour sheets, exerting a more significant influence on the entire flow system. Figure 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 vapour sheets (A, B, C, D) becomes more pronounced. The vapour sheets are then twisted, eventually leading to the fusion of vapour sheets B, C and D at the base and fragmentation at the top, as shown in figure 27 at 
$t=48.2$. The fragments produce small bubbles, making the flow system more chaotic. After fusion, new vapour sheets E and F are generated, as shown in figure 27 at 
$t=66.7$, corresponding to the sudden increase in 
$Nu$ in figure 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 vapour sheets occur almost simultaneously, resulting in a smaller fluctuation in the distance between the interface and the overheated wall.
Figure 27. Dynamics of vapour sheets and the evolution of the streamline pattern for different time instants ranging 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, VOF methods often simplify the coalescence process (Scardovelli & Zaleski Reference Scardovelli and Zaleski1999). Additionally, recent work by Innocenti et al. (Reference Innocenti, Jaccod, Popinet and Chibbaro2021) 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 vapour bubbly flows under external MFs, 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 3-D film-boiling process under MFs have been conducted using a numerical method based on the combined MHD and phase-change model developed in our previous studies, see Zhao et al. (Reference Zhao, Zhang and Ni2022) and Zhang & Ni (Reference Zhang and Ni2018). The present study analyses 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 MF, 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 vapour 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-averaged Nusselt number (
$\overline {Nu}$) is not substantial. However, under high MF strength in the vapour column boiling area, 
$\overline {Nu}$ increases with 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 vapour 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 MF.
In the presence of a horizontal MF, film boiling gradually transitions from periodic bubble detachment to a 2-D planar vapour 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 2-D 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 vapour film during boiling in the form of a planar vapour 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 vapour sheet boiling area. This is because, during boiling in the form of a planar vapour sheet, the Lorentz force acting vertically downward on the bottom vapour film becomes more pronounced with increasing MF strength, resulting in a thinner bottom vapour film thickness. Unlike the vertical MF, it was observed that, under complete planar vapour 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 with vertical MF conditions. During this stage, bubble formation, detachment and merging occur simultaneously, causing 
$Nu$ to deviate from its original variation pattern with a significantly reduced amplitude.
Funding
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.
Appendix
Here, we demonstrate that, even when using more realistic physical properties for the liquid metal, the MHD effects on the rising characteristics of the vapour 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/vapour 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$, the influence of the vertical MF on film boiling remains consistent: the vertical MF suppresses the Rayleigh–Taylor instability, causing a delay or even complete suppression of vapour 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 vapour–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. Panels show (a i) 
$t=18.9$ and (a ii) 
$t=28.9$; (b i) 
$t=17.5$ and (b ii) 
$t=29.8$; (c i) 
$t=17.9$ and (c ii) 
$t=35.8$.
