Numerical simulations on thermocapillary �ow on heated sinusoidal topography

The interaction between thermocapillary flow and substrate geometry is analyzed numerically. Taking surface tension into account, the momentum equation is derived and solved using a commercial FEM solver, Comsol Multiphysics where the effects of surface tension and surface deflection can be easily incorporated into the momentum equation. In the case of Ma ≈ Ma 𝑐 , the strong symmetric thermocapillary flow is observed when the wavelength of topography, 𝜆 𝑇 , and the wavelength of instability motion, 𝜆 , are nearly the same. This interesting phenomenon has been called flow-structure resonance. Through the numerical simulations, various flow modes, such as symmetric two-cell and four-cell modes, asymmetric two-cell mode, and oscillatory asymmetric two-cell mode are identified by changing the Marangoni number and wavelength of topography. It is clearly shown that for a certain 𝜆 𝑇 -system, the transition from oscillatory mode to steady one is possible by relaxing the previous non-deformable surface condition due to high surface tension, i.e., Cr → 0 . The present study reveals that the preferred flow mode is the complex function of the various parameters such as the Marangoni number, the Biot number, the wavelength of topography, and the capillary number.


INTRODUCTION
Surface tension is ubiquitous in natural systems and industrial processes where immiscible fluids are involved.In heat transfer systems with a free surface experiencing surface-tension variations, the fluid motion driven by the temperature-induced surface tension gradient, which has been called thermocapillary flow, is often encountered.Since thermocapillary flow plays important roles in crystal growth [1], surface patterning [2], melting [3] and coating [4], much research has been conducted theoretically and experimentally.
Sen and Davis [5] studied thermocapillary flows in a rectangular slot bounded by differently heated lateral walls.For a small aspect ratio system, Sen and Davis [5] successfully obtained approximate solutions for thermal and velocity fields.For the similar rectangular slot, many researchers extended Sen and Davis's [5] work by employing asymptotic theory and numerical simulation [6].Apart from the liquid inside the laterally heated rectangular slot, a liquid layer over a patterned substrate has been considered experimentally and theoretically.In a three-dimensional (3-D) system, where the topography has two-dimensional (2-D) regularity, Ismagilov et al. [7] experimentally consider the interaction of intrinsic hydrodynamic instability patterns and topographically imposed patterns.They showed that the pattern of Bénard-Marangoni flow can be controlled by the thickness of a fluid layer or the temperature difference between the substrate and air.Later, for the 3-D liquid layer system, Saprykin [8] showed that disjoining pressure can induce the rupture of a very thin liquid layer, where the Marangoni number is small enough.
In addition, in 2-D systems, where the topography has one-dimensional (1-D) regularity, Stroock et al. [9] controlled the direction of Bénard-Marangoni flow by using symmetric and asymmetric grooves.
They also showed that asymmetric grooves can be used as millimeter-scale micropump.Later, Alexeev et al. [10] analyzed thermocapillary flow in a silicone oil film induced by a row of identical parallel grooves experimentally and numerically.The measured flow velocity reached 1.2 mm/s for a wall temperature 40K above the ambient.Their numerical simulations were in good agreement with their experimental findings.They also showed that symmetric wall topography leads to the appearance of symmetric vortices.For the sinusoidal topography system, Kabova et al. [11] showed that the analytic approach employing a 2-D long-wave lubrication approximation is in good agreement with the numerical solution of the 2D Navier-Stokes equation.Recently, Yoo et al. [12] analyzed thermocapillary flows on heated sinusoidal topography using long-wave approximation.They also investigated the interaction between thermocapillary flow and substrate curvature through numerical simulations.
In the present study, 2-D thermocapillary flows in a thin liquid layer laid on heated sinusoidal topography are studied.To understand the interaction between thermocapillary flow and substrate curvature, systematic numerical simulations were conducted by changing the Marangoni number and the wavelength of sinusoidal topography.To identify asymmetric flows in symmetric topographies, simulations are conducted for a half-wavelength domain and one wavelength domain.Furthermore, for a given Marangoni number system, the onset condition of oscillatory motion is determined by changing the curvature of the substrate.

SYSTEM AND GOVERNING EQYATIONS
The system considered here is the liquid film coated over the substrate having sinusoidal topography.
The average thickness of the liquid layer and initial temperature are  and   , respectively.For time  ≥ 0, the substrate is isothermally heated with   and, through the liquid-air interface, the liquid layer is cooled by the convective heat transfer.The schematic diagram of the basic system of pure diffusion is shown in Fig. 1.As shown in this figure, the liquid-air interface temperature is higher at the locations of thin film (over the topography crests) than at the locations of the thick film (topography troughs).This temperature gradient can induce thermocapillary flow near the liquid-are interface region.
For the thin liquid film, by neglecting the gravity effect, the governing equations for mass, momentum, and thermal energy are [5,12] ∇ •  = 0, Here, , , ,  are velocity vectors, time, pressure, and temperature, respectively.The liquid has density, ρ , viscosity,  , and thermal diffusivity,  .In the present study, we didn't consider the disjoining pressure which plays an important role in the rupture of a very thin liquid layer [10].In addition, the surface tension of the liquid, σ, has the following linear relation with the temperature: where  =   ⁄ |   .The proper initial and boundary conditions are  =  and  =   at  = 0, (5a)  =  and  =   at  = (), where () is the sinusoidal function describing the shape of substrate topography.At the air-liquid interface, the following conditions have been applied: where  is the thermal conductivity of liquid and q is the heat transfer coefficient at the interface,  is the unit normal vector from the liquid side to the air side,   and (= − + (∇ + ∇  )) are the total bulk stress tensor in the air and liquid phases, respectively.In addition, ∇  = ∇ − ( • ∇) is the surface gradient operator, and the term −(∇ • ) represents the curvature , i.e.,  = −(∇ • ).
Furthermore, at the air-liquid interface,  = (, , ), the following kinematic condition should be satisfied [13]: In the present study, the sinusoidal topography is expressed as where  0 and   are the amplitude and the wavelength of the topography, respectively.The periodicity of the substrate is considered by the symmetricity of domains summarized in Fig. 2. As shown in this figure because the domain   can be obtained by translating the domain  −1 by   in the positive -direction and domains  0/2 and  1/2 are symmetrical about the line  =   2 ⁄ , the present half-wavelength domain corresponds to the periodic system having a period   .Similarly, the present one-wavelength domain represents the periodic system having a period 2  .
Using  2  ⁄ ,  ,   ⁄ ,  3  ⁄ and Δ(=   −   ) time, length, velocity, pressure, and temperature scales, Eqs. ( 1)-( 3) can be nondimensionalized as The proper initial and boundary conditions are  =  and θ = 0 at  = 0, (12a) Then, the kinematic condition ( 7) can be reduced as where the unit outward normal vector is Here, the important parameters, the dimensionless amplitude of topography  0 , the aspect ratio , the Prandtl number Pr, the Biot number Bi, the Marangoni number Ma, and the capillary number Ca are defined as [14,15]  0 = , where ν(=   ⁄ ) is the kinematic viscosity of the liquid.Material properties of water and silicone oil for calculating the above dimensionless groups are summarized by Alexeev et al. [9].

NUMERICAL SIMULATIONS
In the two-dimensional (, )-domain, the above equations ( 9)-( 13) were solved by using Finite Element Method (FEM) package Comsol Multiphysics  [16] which has been used in similar problems [17][18][19].The "Laminar Flow" module of the COMSOL was used to solve the continuity equation ( 9) and the Navier-Stokes equation (10).The heat transfer equation ( 11) was implemented by using the "Heat Transfer in Fluids" module, where the flow fields were adopted from the Navier-Stokes equation.
To implement the condition at the air-liquid interface, the "Deformed Geometry" module was employed and the interface was allowed to move according to the kinematic condition (13).
To incorporate the boundary condition (12d), first, we applied the following open surface boundary condition at the liquid-air interface: then, we added the surface tension effect by the following weak form: where  ̃ is the test function for the velocity vector and Γ is the air-liquid interface.Using the surface divergence theorem, the second integral can be [20] where ∂Γ is the contour bounding the free surface, (=  × ) is the outward binormal vector, which is tangential to the fluid surface and normal to the edge of the surface [21].As discussed by Walkley et al. [22] the last contour integral of Eq. ( 16) can be neglected if the Dirihilet and periodic conditions are imposed on ∂Γ [12].However, based on Alexeev et al. 's [9] experimental finding: that symmetric wall topography leads to the appearance of symmetric vortices, we used the following symmetry conditions: where   is an outward pointing normal to the symmetry plane, and  is a tangent vector in the symmetry line.In case the contour integral line ∂Γ lies on the symmetry plane,   corresponds to the binormal vector , then  ̃•  = 0 should be satisfied on the symmetry condition (17a).For these typical cases, i.e., the Dirichlet and the symmetric conditions on ∂Γ, the following terms should be added to the open boundary condition imposed on the air-solution interface as a weak contribution: Even though for the simplicity of the analysis, the above second term has been neglected under the assumption of the perfect flat interface [12], here, we kept the effect of the capillary number on the thermocapillary flows.
The COMSOL Multiphysics has a wide variety of options in choosing the time-dependent solver.In the present study, 1 st or 2 nd order, variable step size, and Backward Differentiation Formulae (BDF) were used.The order of the BDF solver is determined by the degree of the interpolating polynomial.
At each time step, the system of nonlinear algebraic equations is linearized employing the Newton method, and the resulting linearized system is solved by PARDISO direct solver which is fast, robust, and multi-core capable.The scaled absolute tolerance factors of 0.05 and 1 were set for the concentration and the pressure, respectively.The relative tolerance was 1 × 10 −4 .Any critical differences were not found by changing these tolerances.In the present study, the intensity of instability motion is characterized by the maximum difference in interfacial temperature, ∆  , and the maximum of the velocity magnitude ||  defined as

RESULTS AND DISCUSSIONS
For the limiting case of  0 → 0 , i.e., flat substrate, Nield [15] suggested the critical Marangoni number and corresponding wavenumber as where  is the wavenumber and  is the wavelength of instability motion.In the case of  0 = 0.05, Pr = 100, Bi = 0.1, Ma = 10 3 and Cr = 10 −6 , we tested the interaction between the wavelength of instability motion, , and that of topography,   , which have been called intrinsic and topographies periodicities, respectively [7,9].As shown in Figs. 3 and 4, for the cases of small  0 and Ma ≈ Ma  , we can expect the strong instability motion around   ≈ (= 2 2.028 ⁄ ), i.e.,  =  1.014 ⁄ .This typical wavelength of topography has been called the resonance wavelength [8,23].Kelly and Pal [23] reported similar results for the buoyancy-driven Rayleigh-Benard system.They showed that at the resonant frequency, i.e.,   = (= 2   ⁄ ), the heat transfer rate enhancement is possible even for Ra < Ra  .For the rigid-rigid boundaries system, the critical Rayleigh number and corresponding wave number are [24] Ra  = 1708 and   = 3.117.
In their Nu − 1 vs. (Ra Ra c ⁄ ) plot of small  0 the system, they showed that convective motion can exist even for (Ra Ra c ⁄ ) < 1 , which was called a quasi-conduction region.They also showed the transition from a convex quai-conduction region to a concave critical region, where (Ra Ra c ⁄ ) > 1 and the vertical temperature gradient plays a certain role.As shown in Fig. 4, where ||  vs.
(Ma Ma c ⁄ ) the plot is given, a similar convex-concave transition is occurred near This means that in the case of (Ma Ma c ⁄ ) < 1, the thermocapillary flows are driven by the horizontal temperature gradient due to the inhomogeneous topography of the substrate.
In addition, the effect of the calculation domain on the flow fields is compared in Fig. 5.It should be noted that regardless of aspect ratio, flow fields are symmetric about the line  =  2 ⁄ , even if we imposed symmetry conditions at  = 0 and  = .This means that the periodic unit is the domain ) which corresponds to the domain  0 .For this small Ma -system, any asymmetric motion isn't observed, and the transition from two-cell mode instability motion proceeded around ~4 (see Fig. 4).For a symmetric grooves system, Stroock et al. [8] observed the transition from period-matched two-cell mode convection to period doubled one-cell mode one by increasing fluid thickness  =  2 ⁄ to  = .In Stroock et al.'s [8] experiment, the Maragoni number is doubled and the dimensionless amplitude of topography is halved from  0 = 0.5 to  0 = 0.25 , direct comparison with the present system, i.e, fixed Ma(= 10 3 ) and  0 (= 0.05), is not possible.However, as shown in Fig. 5, their mode transitioning phenomenon can be seen in the present analysis.
If we increase the Marangoni number to Ma = 2 × 10 3 , the flow field is slightly different from that for Ma = 10 3 .Figure 5 shows that the transition from two-cell symmetric mode to two-cell asymmetric mode flow proceeded ~3.Further, an increase in the aspect ratio to  = 5 enables us to observe oscillatory thermocapillary flow more clearly at Ma = 5 × 10 3 .As shown in Fig. 9, for the given Marangoni number, we can control the flow mode by changing the wavelength of topography,   .This controllability of the flow mode has never been studied.For some characteristic time marked in Fig. 9(a), flow fields are summarized in Fig. 9(c) which shows an oscillatory two-cell mode is the preferred one for the present oscillatory system.However, this oscillatory mode motion cannot be survived in the slightly larger or smaller aspect ratio system.This means that the preferred flow mode is a complex function of the Marangoni number and the aspect ratio.
Until now, the interaction between the Marangoni number and the wavelength of topography is focused by fixing Ca = 10 −6 .However, the capillary number effect cannot be ignored for the thin liquid layer because Ca~ −1 .Here, we studied the effect of the capillary number on the flow mode of the previous oscillatory system, i.e.,  0 = 0.05 , Pr = 100 , Bi = 0.1 and Ma = 5 × 10 3 .As shown in Fig. 10, for the case of Ca ≤ 10 −4 the oscillatory motion is a dominant thermocapillary mode and the flow mode transition explained in Fig. 9 is proceeded.This kind of oscillatory motion can induce the temporal variation of the liquid thickness (see Fig. 10(b)).It is interesting that for the oscillatory flow regime, i.e., Ca ≤ 10 −4 , the capillary number enhances the variation of liquid layer thickness, however, a deformable surface reduces the maximum variation of liquid layer thickness for the steady flow regime.Furthermore, as shown in Fig. 10(c) deformable surfaces reduce the maximum surface temperature difference, which corresponds to the driving force of the thermocapillary convection.Figure 10(d) explains the effect of the capillary number on the flow field.For the case of Ca = 5 × 10 −4 , we can expect steady four-cell mode thermocapillary flow which cannot be observed for the systems of Ca ≤ 10 −4 .For the steady four-cell mode system, i.e., Ca = 5 × 10 −4 , flow fields for the transition from two-cell mode to four-cell mode are summarized in Fig. 10(e).The above findings mean that non-deformable surface due to high surface tension, Cr → 0 prevents flow mode transition from steady two-cell symmetric mode to more stable steady four-cell symmetric mode.This kind of interaction among the capillary number, the Marangoni number, and the geometry of topography have never been considered.

CONCLUSIONS
The interaction between thermocapillary flow and substrate curvature was analyzed systematically.
For a weak flow condition, i.e., the case of Ma ≈ Ma  , the strong thermocapillary motion is observed when the wavelength of topography,   , and the wavelength of instability motion, , are nearly the same.Furthermore, the calculation in the half-wavelength domain is quite enough to explain the weak thermocapillary motion.However, the thermocapillary flow structure deviates from the symmetry of the half-wavelength domain as Ma increases.For the system of Ma = 2 × 10 3 , Bi = 0.1 and Ca → 0, this asymmetric effect is remarkable for a certain range of   .To understand the resonance between the surface tension-driven motion and the wavelength of topography more fully, more refined work to consider the effects of the amplitude of topography,  0 , the Biot number, Bi, and the capillary number, Ca, are strongly needed.

5 .Fig. 8 ,
Fig. 8, where the effects of the Marangoni number on the thermocapillary flow field are summarized, we can oscillatory mode thermocapillary motion for the range of Ma = 1.5 × 10 3 .It seems that for the case of  = 4.0, as increasing the Marangoni number, we can observe the flow mode transition from steady symmetric four-cell mode flow to steady asymmetric two-cell mode one through transient asymmetric two-cell mode flow.