Stability analysis of wall‐attached Bénard–Marangoni convection in a vertical magnetic field

The threshold for the onset of thermocapillary flow in a planar liquid layer heated from below is increased by a vertical magnetic field when the liquid is a good electric conductor. The magnetic damping effect is reduced when the induced eddy currents are blocked by insulating side walls. Neutral conditions for this specific Bénard–Marangoni stability problem with a vertical field and side walls are obtained numerically for three‐dimensional perturbations assumed periodic in one horizontal direction. The domain is bounded by a free‐slip wall at the bottom, a free surface at the top and two free‐slip lateral walls in the other horizontal direction. Buoyancy forces and surface deformations are neglected and a constant heat flux is imposed on the free surface. Upon increasing the magnetic induction, the least stable modes become localized near the side walls and the convective threshold increases at a lower rate than for the least stable bulk mode.


INTRODUCTION
Thermal convection in a liquid layer can appear as the result of an instability if the layer is heated from below.In addition to the buoyancy generated by density variation with temperature, the thermocapillary effect can contribute as a forcing mechanism when the layer has a free surface or interface.This specific setup is known as Bénard-Marangoni convection (BMC).It is one of the generic model systems where a flow organizes into cellular patterns (e.g., hexagonal Bénard cells).Many works have explored the various BMC patterns and their stability in common liquids, for example, silicone oils [1,2].Thermocapillary flows also occur in electrically conducting liquids, specifically in molten semiconductors during crystal growth or during electron beam evaporation of metals for coating processes.While the plane layer BMC is not directly relevant for such applications, it can nonetheless be interesting to investigate it also for these types of liquids.Their electrical conductivity allows one to use magnetic fields to influence the flow through electromagnetic induction.
The buoyancy-driven analog of BMC, that is, Rayleigh-Bénard convection (RBC), is of considerable interest in magnetohydrodynamics (MHD) because astrophysical flows in the Sun or the liquid cores of planets are driven by convection and accompanied by magnetic fields [3].Magnetic fields tend to suppress the convection.Recent works on magnetoconvection in the RBC system have revealed that this magnetic damping effect is less pronounced near a side wall, that is, the flow can persist near the side walls while it is already completely suppressed in the bulk by magnetic F I G U R E 1 Sketch of the geometry.The liquid container is rectangular with electrically insulating side and bottom walls.The free surface is flat and electrically insulating with a prescribed constant heat flux density .Width and depth of the liquid layer are denoted by   and .The applied uniform magnetic field  is vertical.Walls are assumed to be stress-free.
damping [4][5][6][7][8].This provides the main motivation for the present work, which is concerned with the effect of side walls on the linear instability of BMC in a magnetic field.Although side walls were considered in thermocapillary convection with lateral heating, BMC in a magnetic field with side walls has not been studied to the knowledge of the author.There are, however, previous works based on normal mode analysis that were concerned with the primary instability of the basic quiescent state in the bulk of the layer [9][10][11].These as well as the present study can be of interest, for example, in the context of the Czochralski method for crystal growth.In this process, temperature fluctuations in the melt are reduced by the application of static and/or time-dependent magnetic fields in order to avoid imperfections in the crystal [12].
The paper is divided into three sections devoted to the mathematical model and governing equations, the numerical method, and to the stability results, respectively.Brief conclusions are provided at the end.

MATHEMATICAL MODEL
The thermocapillary effect gives rise to a tangential stress on the liquid surface on account of surface tension changes with temperature.A complete description of BMC, therefore, needs to take the adjacent fluid phase into account.In the case of an adjacent gas phase, one can represent its effect on the liquid through simplified boundary conditions.This so-called one-layer model [13] serves as the starting point for the present problem.For reasons of convenience, the flow domain will be assumed periodic in one horizontal direction.The liquid is, therefore, only confined by two parallel side walls in addition to the bottom wall.The geometry and coordinate system is shown in Figure 1.The free surface is assumed to remain flat.This assumption holds when the surface tension  becomes very large, which corresponds to a vanishing crispation number.The long-wave mode of thermocapillary instability is thereby excluded [2].The kinematic boundary condition for the velocity field  on the flat surface is   = 0.The remaining mechanical boundary condition is the tangential stress balance, that is,     = ∇  , where  denotes the dynamic viscosity and the subscript  the projection of the vector on the surface.In addition, a fixed heat flux density  is prescribed on the surface.This is a special case of Newton's law of cooling and corresponds to a vanishing Biot number.It was applied in a previous work [14].Pressure, velocity, and temperature fields satisfy the incompressible Navier-Stokes and energy equations.
The electromagnetic induction problem is treated in the framework of the quasistatic approximation [3].It holds provided that the magnetic Reynolds number is small.In this case, the induced magnetic field is negligible compared to the externally imposed magnetic field  =   .Moreover, the induced electric field  turns out to be curl-free.The induced current density  can, therefore, be computed from Ohm's law for a moving conductor and the charge conservation condition ∇ ⋅  = 0.The electric field is represented through the negative electric potential gradient, that is,  = −∇.The domain boundaries are assumed to be electrically insulating, that is, the normal component of  is zero.
For pure heat conduction, the temperature distribution is linear in .Its slope follows from the prescribed heat flux, and the temperature difference between bottom and top boundaries is Δ 0 = ∕, where  is the thermal conductivity.It serves as unit of temperature for nondimensionalization.The units of length, velocity, and time are chosen as the layer depth , ∕, and  2 ∕, where  is the thermal diffusivity.The units of the magnetic field and the electric potential are  and .With this choice of units, the governing equations take the following nondimensional form: ∇ ⋅  = 0. ( The quantity  is the difference between the temperature and the linear conduction profile, that is,  =  −  + 1∕2.The last term in the momentum equation ( 1) is the Lorentz force density.The buoyancy force is neglected.Equation ( 4) is Ohm's law.The condition (5) determines the electric potential , which satisfies a Poisson equation.Further symbols are the Prandtl number  = ∕ and the Hartmann number which characterizes the magnitude of the Lorentz forces relative to viscous forces.The symbol  is the electrical conductivity of the liquid.The boundary conditions at the free surface are The tangential stress boundary conditions contain the Marangoni number which is proportional to the derivative of the surface tension  with respect to the temperature .The minus sign is introduced since  typically decreases with .

NUMERICAL METHOD
In the present work, it is assumed that the instability is of stationary type.The nonlinear terms and the time derivatives in Equations ( 1) and ( 2) are, therefore, neglected.The momentum and continuity equations reduce to the Stokes problem with additional Lorentz force.To avoid complications stemming from the coupling between pressure and velocity, one can choose a representation for  that satisfies the continuity equation automatically.By that, the pressure can be eliminated.This approach was employed for linear stability analysis of MHD duct flow [15] and is adopted in the present work.The velocity field is written as the curl of a vector streamfunction , that is, and the gauge condition ∇ ⋅  = 0 is imposed to determine  uniquely.The dependence on  is represented by the normal mode ansatz with wavenumber  for all fields, for example, (, ) exp() for the temperature perturbation.
The velocity components then read Equations for   and   are obtained by taking the curl of the definition ( 9) and the momentum equation (1).They are The quantities   and   are the and -components of the vorticity field ∇ × .Equations for the remaining quantities are The last equation ( 17) is obtained by substitution of Ohm's law (4) into Equation ( 5).Combined with boundary conditions (7) on the free surface and further boundary conditions on the bottom and side walls specified below, Equations ( 12)-( 17) represent a linear eigenvalue problem for the Marangoni number  that must be solved numerically.A suitable discretization of this problem is obtained by a spectral collocation method with Chebyshev polynomials   () = cos( arccos()).The scalar fields such as  are expanded as where   =   ∕, −  ∕2 ≤  ≤   ∕2 and −1∕2 ≤  ≤ 1∕2.The Poisson equations ( 14)-( 17) and boundary conditions are imposed pointwise at the Gauss-Lobatto collocation points where   + 1 and   + 1 are the number of expansion terms with respect to  and .
The boundary conditions for the vector stream function and vorticity components are determined with the help of equations (11).Zero normal velocity on the boundaries with  = .requires   = 0 and     = 0. On the boundaries with  = .the corresponding conditions are   = 0 and     = 0.The zero stress boundary conditions on  = .are satisfied if  2     = 0 and  3    = 0. On the bottom wall, these conditions are  2    = 0 and  3    = 0.The tangential stress (Marangoni) boundary conditions at  = 1∕2 take the form    2     +  3    = −     ,   =    .
The remaining boundary conditions for Equations ( 16) and ( 17) are obvious.Since the boundary conditions (zero normal velocity) for Equations ( 12) and ( 13) only involve   and   , respectively, one can represent the expansion coefficients of   and   by linear invertible maps through those of   and   (assuming the latter are augmented by the zero boundary values to be imposed on either   ,   or its normal derivatives).The expansions for   and   , therefore, contain   + 3 and   + 3 terms, respectively.The values of   ,   or its derivatives in Equations ( 14)-( 17) (and associated boundary conditions) at the collocation points are represented through expansion coefficients of   and   via these linear invertible maps.As a result of the collocation approximation, one obtains a vector  of unknowns containing the expansion coefficients of   ,   , , and  with a size of 4(  + 1)(  + 1) and a generalized linear eigenvalue problem The method was implemented in Matlab [16] using the default double precision.Notice that   ,   , and  are real variables.According to Equations ( 11) and ( 16), both   and  would be imaginary quantities.They are considered as real variables in the code and below.Problem (21) was solved with Matlab's eig routine to find all eigenvalues and eigenvectors.The routine also works with a matrix  whose rank is smaller than the rank of  (as it is the case for Equation (21).It associates spurious solutions with infinite eigenvalues.In the following discussion, it is assumed that the set of finite eigenvalues is sorted in ascending order, that is, the smallest one will be referred to as the first eigenvalue, and so forth.

RESULTS
The parameters in problem (21) are the Hartmann number , the wavenumber , and the aspect ratio   .To keep the computational cost manageable,   was set to relatively small values.Most computations were carried out for   = 2.The wavenumber was varied over a relatively large range at a given  and   .The value of  was systematically increased, which required increasing the numbers   and   of Chebyshev modes.The limits in the computations were   ≤ 60 and   ≤ 50.Sufficient resolution was verified by comparing results obtained with different   and   .Typically, between 5 and 10, modes were added per direction for this purpose.The variation of  with  is shown in Figure 2A-C for  ≤ 20.Also included are the neutral curves for the unbounded plane layer that were computed in Boeck [14].The corresponding eigenmodes are purely two-dimensional, that is, without variation in .Thanks to the chosen boundary conditions with zero heat flux and free slip on the side walls, these modes for the plane layer are also compatible with the present configuration with side walls in the case  = 0.This is apparent in Figure 2A where the red curve always coincides with an eigenvalue of problem (21) and matches the first eigenvalue for  > 1.5.This result serves as a verification of the Matlab code for the case without magnetic field.
When the magnetic field is present, the agreement between the plane layer and the present configuration is lost.The first and second eigenvalue become smaller for the case with side walls.For  = 20, the larger eigenvalues approach the curve for the plane layer.The first and second eigenvalues approach one another for  = 10 (Figure 2B) and essentially coincide for  = 20 (Figure 2C) except at very small .This applies equally for all larger values of  ≤ 60 considered in this work.In addition, the gap between the first and second eigenvalues and the next ones increases with .The qualitative behavior is analogous to Figure 2C.The first and second eigenvalues are also insensitive to   provided the latter is large enough.Figure 2D demonstrates this for the first eigenvalue at  = 20, which remains effectively unchanged between   = 2 and   = 3 (except for  → 0).The spatial structure of the first and second modes is shown in Figures 3 and 4 for  = 40 and  = 20.One can see that at  = 40, there is already a significant separation between the regions of significant velocity magnitude, which are localized close to the side walls.For the temperature perturbation, this separation is not as pronounced since it is smoother than   (cf.Equation 16).For  = 20, the separation (central gap) is not as strong, and regions of large |  | are broader.It is also apparent that the first and second eigenmode have different symmetries with respect to .The temperature perturbation (, ) and   (, ),   (, ) are even functions of  and   (, ) is an odd function of  for the first eigenmode (Figure 3).These features are reversed for the second eigenmode, that is,  and   ,   are odd and   is (A) (B) F I G U R E 5 Dependence of the critical value   (A) and corresponding wavenumber   (B) on  for   = 2 (identified as wall mode) and for the unbounded plane layer [14].
an even function of  (Figure 4).As expected by the coinciding first and second eigenvalues, the two modes themselves become indistinguishable either in the positive or negative half of the domain ( > 0 or  < 0) as  increases, that is, one can transform one into the other by symmetry operations.Notice also that the representations of the eigenmodes in Figures 3 and 4 do not represent their structure in one and the same plane  = .The -dependence differs for the different components of velocity and .It can be inferred from Equations ( 11) and ( 16).
The minimum of the first eigenvalue with respect to  represents the stability limit, that is, the critical value   for the onset of convection.Figure 5 compares   and the corresponding wavenumber   obtained with   = 2 (wall mode) with the results for the unbounded plane layer.A clear power law scaling on  is not seen in Figure 5A although   seems to increase roughly as  3∕2 at the upper limit of .The wavenumber   shown in Figure 5B does not show a power law dependence.Instead it might approach a finite limit.

CONCLUSIONS
The numerical study confirms that the onset of BMC with a vertical magnetic field occurs near the side walls as it is the case for RBC.The separation in the critical Marangoni numbers   for onset near the wall and in the also increases with  as it is the case for the critical Rayleigh numbers   in RBC.A finite limit of the wavenumber   for  → ∞ as for RBC could not be demonstrated convincingly but seems possible.An asymptotic dependence   ∼  3∕2 would match the prediction   ∼  3∕2 obtained by Busse [5].Exploring these limits requires a different approach since the numerical method does not allow one to increase the resolution by much.The free-slip walls in the present work were chosen as a compromise to keep the mathematical model simple.It seems clear that a finite value of the temperature perturbation should occur at the top corners of the domain (cf.Figures 3A and  4A).The latter of the equations (20) implies a finite value of   at these corners.However, a no-slip condition on the side walls imposes   = 0, that   would be discontinuous at the corner.A significant enhancement of the present model seems to be required in order to resolve this inconsistency and to take further effects into account that would be present in actual experiments with liquid metals [17].

A C K N O W L E D G M E N T S
The author is grateful to C. Karcher for helpful comments on the manuscript.
Open access funding enabled and organized by Projekt DEAL.

O R C I D
The gauge condition allows one to express the -component of  by   (, ) = −    (, ) −     (, ).

2
Dependence of  on the wavenumber  for   = 2 and  = 0 (A),  = 10 (B), and  = 20 (C).The first four eigenvalues are shown.The neutral stability curve for the plane layer with free-slip bottom is included for comparison.(D) Influence of   on the first eigenvalue at  = 20.