Numerical Simulation of Thermocapillary Convection with Evaporation Induced by Boundary Heating

Nomenclature

1  Introduction

Thermocapillary convection is a subject of extensive research due to the great importance of the results which can be applied in developing various technologies in the fields of materials science and microfluidic systems. The mechanism of occurrence of convective motion is related to changes in the surface tension of a liquid. The changes can result from an external thermal action [1,2], concentration effects [3–5] or evaporation [6–8] on the liquid surface. If the combined action of these effects occurs, they can be amplified or they can compete with each other, causing the alteration of flow regimes and instabilities of different nature. The control of thermally induced and/or concentration-induced motion and the dynamical operation of transport processes in the volume of a liquid with a free surface or in multicomponent systems with interfaces are crucial for the development of new innovative approaches in material engineering and various smart technologies. These include coating technologies connected with obtaining coversheets with given functional characteristics [4,9–11], methods applied upon surface purification [5,12,13], bioengineering and medical technologies, including microfluidic sorting [14,15], targeted drug delivery [16] and local elevation/reduction of components in fluid solutions [2,17], as well as various approaches directed towards the enhancement of heat transfer characteristics in electronic devices [18].

Fluidic systems with working media being in different aggregate states and filling domains with outer solid boundaries as well as limited volumes of the liquid contacting with a rigid substrate have a three-phase contact line near which the convection characteristics can also significantly vary [19,20]. The development of mathematical models based on the conservation laws of mass, momentum and energy for describing the fluid dynamics in domains with solid walls and internal interfaces is fraught with difficulties. The main challenges are in choosing a way to describe the three-phase line motion (or three-phase points for two-dimensional statements), correct determining the dynamic contact angle, and substantiating the kinematic compatibility of the boundary conditions of different types for the velocity function [21–24]. The problem still has no well-posed closure. For the first time, this inconsistency between the conditions on the free surface and the no-slip ones in the contact point neighborhood was pointed out in [21] and was mathematically proved in [22] with minimal assumptions for the velocity and free surface smoothness. In the common case, the energetic norm of the velocity field will become infinite which means that the energy becomes unbounded. There are various approaches to correctly formulate mathematical models for describing flows of a viscous incompressible liquid with a moving contact line or a moving contact point in the two-dimensional statement. The approaches allow one to replace the no-slip conditions with the slip conditions on the solid wall sections near the contact line, or the asymptotic approach can be used; in addition, one can also consider the suggestion about the contact angle being equal to π in the leakage regime or to “zero” in the runoff one as well as other statements. The solvability of the initial boundary-value problems was studied for a number of mathematical models of fluid flows with dynamic contact angles [23–25]. The well-posedness of the statement with the free boundary and dynamic contact angle was proved in [24,25] in the case when the no-slip condition on the solid boundary was changed to the condition of proportionality of tangential stresses to tangential velocities. The existence of a solution and asymptotic expansion in the neighborhood of the contact point were proved, and the regularity of the free boundary was shown in the statement with a uniformly moving contact point. Two-dimensional and axis-symmetric stationary problems were studied in [23] in the case of capillary numbers being small. Theoretical results of the study of apparent contact angles in a problem of an evaporating liquid placed on a solid substrate are presented in [26]. In [27], a moving contact line problem for two-phase incompressible flows was investigated based on the kinematic approach using an evolution equation for the contact angle in terms of the transporting velocity field; this approach and other common models including the gradients of surface tension and mass transfer across the fluid–fluid interface were analyzed. The behavior of the contact angle relative to the contact point velocity was investigated in physical experiments (see [28] and papers cited there). In [28], the numerical investigations of isothermal fluid flows with the dynamic contact angle were carried out for a problem with the proportionality conditions of tangential stresses and velocities on solid walls; here, rather good agreement with experimental results was confirmed and the monotonous character of the dependence of the contact angle on the contact point velocity was shown. In [29], this problem was studied in terms of the influence of the thermal boundary regimes on the contact point velocity. The thorough consideration of the statements in terms of contact lines is due to the problems of heat transfer processes in fluids, especially, in the case of evaporation [20,30]. The evaporation rate and heat flux density in the vicinity of the contact lines are very high [31]. The value of the contact angle depends on the temperature of the heated solid surface and on the evaporation rate as well as on the contact line velocity. One of the effects observed experimentally is the contact angle hysteresis. The effect consists in the existing difference between the advancing and receding contact angles [32] (see review in [20] and references in [20,32]). This phenomenon is studied based on different approaches in a significant number of papers. An ambiguity appears even in the definition of the static value of the contact angle. Experiments on wetting dynamics allow one to investigate various problems of the contact line motion concerning the surface heterogeneities. On a small scale, they lead to contact line deformation while on a macroscopic scale, they result in contact angle hysteresis.

In the present paper, we study the structure and characteristics of convective regimes in a liquid–gas system under diffusive-type evaporation based on numerical modeling. We investigate the influence of different types of thermal action on the flow structure, behavior of the phase boundary and dynamic contact angle. We consider the cases of local heating from the substrate and distributed lateral thermal impact. The simulation is carried out in the framework of a two-dimensional mathematical model whose implementation admits the correct agreement between the no-slip boundary condition and the moving contact points [21,33].

2  Problem Statement

2.1 Basic Assumptions and Description of the Physical System

We investigate joint non-stationary flows in a two-layer system filling a plane rectangular cuvette (Fig. 1). The domains Ω1 and Ω2 are occupied by gas and liquid media, respectively, and divided by a common interface which is a thermocapillary surface admitting the mass transfer of the diffusive type due to evaporation/condensation. The phase boundary is considered as a smooth curve embedded in two-dimensional space. The curve is defined by the equation y=f(t,x) . The flow domain Ω=Ω1∩Ω2 has an external boundary ∂Ω=∂Ω1∩∂Ω2,∂Ω={x=0,x=X,y=0,y=Y}; here, the segments composing the outer edge are solid impermeable walls. The points where the interface between two immiscible fluids joins the lateral walls x=0 and x=X are three-phase contact points. The angles between the surface and the rigid boundaries are contact angles.

Figure 1: The scheme of a two-phase flow in the Cartesian coordinates

Initially, both fluids are at rest, and the interface is straight so that the left and right contact angles φ are equal to 90°. The media are in the state of intermutual saturation with the average concentration vapor in the gas C0 . The value is taken as an equilibrium concentration corresponding to the temperature T0 , which characterizes the initial uniform temperature field in the whole two-phase system. The motion in the bilayer system is induced by an external thermal load applied either from the lower boundary y=0 by activating the heaters with the finite size arranged on the substrate or from the lateral walls due to uniform distributed heating. The temperature effect results in convective motion in both fluids, increase in the evaporation rate and interface deformations, causing changes in the contact angles. The fluid motion occurs from the given initial state which is the state of mechanical and local thermodynamical equilibrium.

We assume that the system is in the terrestrial gravity field g=(0,−g) , g = 9.81 m/s2, and both fluids are incompressible viscous heat-conducting media. Then, the basic mechanisms governing the convective flows with evaporation in the bilayer system are buoyancy and thermocapillary forces caused by external heating of the cell walls. For describing the dynamics of the system we use a two-sided model of evaporative convection on the basis of the Oberbeck–Boussinesq approximation proposed in [33]. The formulation of the problem within the frame of the mathematical model implies the transition from true physical variables of the velocity and pressure to stream and vorticity functions (ψ and ω) and the introduction of additional auxiliary functions allowing one to set all the boundary conditions in the form of conservation laws in the explicit form.

2.2 Governing Equations

To model the convective heat and mass transfer in the domains Ωj, we use the governing equations written in ψ−ω terms in the dimensionless form [33]:

∂ωj∂t+∂∂x(ωj∂ψj∂y)−∂∂y(ωj∂ψj∂x)=1RejΔωj+GrjRej2∂Tj∂x+δj1γGaRe12∂C∂x,

Δψj+ωj=0, (1)

∂Tj∂t+∂∂x(Tj∂ψj∂y)−∂∂y(Tj∂ψj∂x)=1PrjRej(ΔTj+δj1αCΔC).

Assuming that the vapor is a passive admixture which does not change any properties of the background gas, we add the convection-diffusion equation

∂C∂t+∂∂x(C∂ψ1∂y)−∂∂y(C∂ψ1∂x)=1PeRe1(ΔC+αTΔT1). (2)

to model the process of vapor transfer in the gas. The subscripts j=1 and j=2 correspond to the upper Ω1 and lower Ω2 domains, respectively (see Fig. 1), and δj1 is the Kronecker delta symbol. The stream function ψj ( uj=∂ψj∂y,vj=−∂ψj∂x ), vorticity ωj , temperature Tj in the jth layer and vapor content C in the upper layer are the required functions. Here, uj,vj are the velocity components, and Tj defines the temperature deviation of the jth fluid from the equilibrium temperature T0. The problem contains the following similarity criteria: Reynolds numbers Rej=u∗h∗/νj , Grashof numbers Grj=βjT∗gh∗3/νj2, Prandtl numbers Prj=νj/χj , Galilei number Ga=gh∗3/ν12, and diffusive Peclet number Pe=u∗h∗/D. The presented similarity criteria are calculated using the physical parameters of the media: βj (coefficient of thermal expansion), νj (kinematic viscosity), χj (thermal diffusivity), D (diffusion coefficient) as well as the characteristic scales of the problem such as temperature drop T∗ , velocity u∗=ν1/h∗ (relaxation velocity of viscous stresses in the gas phase), and length h∗=Y (height of the cuvette). In contrast to the classical Oberbeck–Boussinesq and molecular transport equations, relations (1) and (2) take into account the effect of the concentration density expansion and thermodiffusion and diffusive thermal conductivity effects which appear in the gas phase due to the presence of a vaporized component. Here, the coefficients αT and αC are non-dimensional analogues of the Soret and Dufour parameters αT∗ and αC,∗ , respectively ( αT=αT∗T∗ and αC=αC∗/T∗ ); γ is the coefficient of concentration expansion.

2.3 Initial and Boundary Conditions

The following conditions define the initial state of the system at the time moment t=0 :

ψj=0,Tj=0,C=C0. (3)

The initial location of the interface is set by the equation y=f0=y0/h∗.

On the outer boundary ∂Ω the conditions for the stream functions ψj are derived as consequences of the no-slip conditions for the velocity vectors, and the corresponding temperature regimes are prescribed. Furthermore, we impose the zero vapor flux condition on the part of the outer boundary contacting with the gas layer

ψj|∂Ωj =0,∂ψj∂n|∂Ωj=0,

Tj|x=0 =TL(t),Tj|x=X =TR(t),T1|y=Y =0,T2|y=0, x∉Qls =0,T2|y=0, x∈Qls =qls(t), (4)

(∂C∂n+αT∂T1∂n)∂Ω1=0.

Here, TL , TR are the constants giving the heating temperature of the lateral walls, Qls ( l=1,2 ) is the area of the substrate where the lth heater with the temperature qls is arranged (see Fig. 1).

The interface conditions are derived on the basis of the conservation laws for mass, momentum and energy, while other ones are the result of using additional suppositions with regard to the character of the processes under study [7,34]. In the problem being discussed, the basic assumptions are the following:

(i) the thermocapillary surface y=f(x,t) is characterized by the surface tension σ which is a linear function on the temperature, σ(T)=1−MaCa(T2−T0). Here, Ma=σTT∗/(ρ2u∗ν2) is the Marangoni number, Ca=ρ2u∗ν2/σ0 is the capillary number; σ0 , T0 are the reference values of the surface tension and temperature for the liquid, respectively, and σT is the temperature coefficient of the surface tension (both constants σ0 and σT are positive);

(ii) we consider that T0 is the local thermodynamic equilibrium temperature, and there occurs only the diffusive mass transfer due to evaporation or condensation through the interface (the convective mass flux is not taken into account).

Then, complementing the kinematic, dynamic and energy conditions with the standard continuity conditions for the velocity and temperature and the relation for the vapor concentration on the phase boundary, we can write the interface conditions at y=f(x,t) in terms of the ψ−ω functions in the following form [33]:

∂f∂t+1+(∂f∂t)2∂ψ2∂s=0,

ω2−ρνω1=F1(t,x),∂ω2∂n−ρν∂ω1∂n=F2(t,x),

∂T2∂n−k∂T1∂n−αCk∂C∂n=−EM, (5)

T1=T2,ψ1=ψ2,∂ψ2∂n−∂ψ1∂n=0,

C=C0(1+εT1).

Here, ∂/∂s and ∂/∂n are the differential operators for derivatives in the direction of the tangent and normal vectors to the interface s and n, respectively (see Fig. 1), ρ=ρ1/ρ2 , ν=ν1/ν2 , k=k1/k2 is the ratio of the thermal conductivity coefficients kj , E=LDρ1/(k2T∗) is the evaporation number, L is the latent heat of vaporization, ε=T∗Lμ0/(RgT02) , μ0 is the molar mass of the evaporating liquid, and Rg is the universal gas constant. The functions F1 and F2 are expressed with the help of the auxiliary functions vs and vn which have the meaning of the tangent and normal components of the velocity vector for any fluid particle lying on the interface (the velocity at any point on the surface can be presented as v=vnn+vss) :

F1(t,x)=Ma∂T∂s+2(1−ρν)(∂vn∂s+vsR),

F2(t,x)=−2[∂∂s(∂(v2)n∂n−ρν∂(v1)n∂n)]+2[(1−ρν)∂∂s(∂f∂xvsR)]+

+1Ca∂∂s[(1−MaCaT1)R−1]−(Gr2Re2−ρνGr1Re1)T1∂f∂x[1+(∂f∂x)2]−1/2+

+GaRe2(1−ρ+γρC)∂f∂x[1+(∂f∂x)2]−1/2+

+Re2[(ρ−1)∂vs∂t+(ρ−1)vs∂vs∂s+(1−ρ)∂f∂xvn2R+vn(ω2−ρω1)],

M=−(∂C∂n+αT∂T1∂n). (6)

The sign of M determines the regime of the phase transitions. If M>0 , then the liquid evaporates into the gas, if M<0 , then the vapor condensation occurs. Given the chosen way of nondimensionalization, the characteristic value of the evaporative mass flow rate is M∗=Dρ1/h∗.

A detailed substantiation of the problem statement presented is given in [33]. Therein, the consistent interpretation of the three-phase contact points in the physical system is explained in accordance with the conception imposed in [21], where the idea of the kinematic compatibility of the no-slip condition and moving three-phase contact point was convincingly justified in terms of the displacing and displaced fluids.

3  Results of Numerical Modeling

3.1 Brief Remarks on the Numerical Algorithm and Calculation Parameters

We use the numerical algorithm proposed in [33] in order to solve problem (1)–(5). A distinctive peculiarity of the method is that in its implementation there is no need to explicitly keep track of the location of the contact points. The algorithm provides calculations in the canonic computation region Ω^={[0,1]×[0,1]} instead of the physical domains Ωj with a curvilinear interface. At each time step, new spatial variables ξ and ζ are introduced: x=ξ,y=ζ(Y−f(t,ξ))+f(t,ξ) for Ω1 and x=ξ,y=ζf(t,ξ) for Ω2. All the equations and boundary relations are rewritten in the new variables ξ and ζ . The physical contact points are corner points in Ω^ , which simultaneously belong to the side and lower (for Ω1 ) or upper (for Ω2 ) boundaries. Moreover, the contact angles (which are the angles approaching the interface to the sidewalls) can be directly determined with the help of the numerical algorithm. The interpretation of the three-phase contact point within the framework of the numerical model is dictated by the mathematical model, enabling the departure of the fluid material point from the contact point into both liquid and gas volume phases. Such a treatment of the three-phase contact point is consistent with the presupposition on the slow velocity of the point motion [21]. The complete description of the numerical algorithm used, key aspects, validation and testing procedures and formulas for the calculation of the contact angles are presented in detail in [33]. The numerical algorithm is based on the use of an unconditionally stable finite-difference scheme and formally has the second order of accuracy. Numerical experiments on a sequence of refined grids were carried out in order to test the computational method and to calculate the experimental order of convergence that was established to be close to the second one. Note that there are no experimental data for comparison which would completely correspond to the considered problem. The straightforward verification and validation of the numerical algorithm were performed in [1], where the method used was modified for a more complicated bilayer system with a deformable interface and a free surface. Therein, a good qualitative agreement with the experimental data was demonstrated, the calculated and measured values of different characteristics were presented and the reasons for quantitative differences related to inherent errors were discussed.

The dynamics of convective heat and mass transfer in a closed cuvette filled with the liquid and gas media divided by an interface under various heating conditions is studied on the example of a benzene–air system. The physicochemical parameters of the working fluids are given in Table 1 [35,36]. The geometrical characteristics of the system are chosen as follows: the cuvette length and height are X=0.04 m and Y=0.002 m, respectively. At the initial state, at t=0 the layer thicknesses are equal to 10−3 m. All the dimensionless parameters and characteristic quantities are calculated according to the corresponding formulas and given in Table 2.

We consider the cases of non-stationary heating from below (Heating mode I) and from the side of the lateral walls (Heating mode II). In Heating mode I, the temperature of the lateral walls does not change, TL=TR=0 , and the thermal load is produced by two heaters arranged on the substrate. The elements have the same size |Qjs|=0.004m and operate in different regimes. From the power up time the first thermal element (Heater 1, see Fig. 1) has the constant temperature q1s(0)=0.1 (the value corresponds to the real temperature 11∘C ), and operates for 100 s and then, it is switched off so that q1s(100)=0 . The second heater functions in the commutated mode at varied temperature: q2s(0)=0.1 , q2s(10)=0.15,q2s(20)=0.2, q2s(30)=0.25, q2s(40)=0.3, q2s(50)=0.35, q2s(60)=0.4, q2s(70)=0.45. After 70 s, this heater continues to work at constant temperature q2s=0.45 (the corresponding dimensional value is 14.5∘C so that the maximum temperature drop by 4.5 degrees is provided). In Heating mode II, the substrate temperature remains to be constant, T2=0 , and time-dependent heating of both lateral walls is realized in the same switching regime. The temperature changes throughout the length of the side boundaries according to the rule: TL(0)=TR(0)=0.1,TL(10)=TR(10)=0.15,TL(20)=TR(20)=0.2,TL(30)=TR(30)=0.25,TL(40)=TR(40)=0.3,TL(50)=TR(50)=0.35,TL(60)=TR(60)=0.4,TL(70)=TR(70)=0.45. After 70 s, the temperature of the walls is kept constant. We analyze the structure of the flows and variations of the evaporative mass flow rate induced by the external thermal load, as well as the behavior of the interface and contact angles in each case.

3.2 Heating Mode I

Figs. 2–6 demonstrate the evolution of the basic characteristics of the bilayer system upon local non-stationary heating. We present the parameters of the system at certain time moments in order to show the qualitative alteration of the flow pattern accompanied by the formation of typical transition regimes (Figs. 2–5) as well as the ultimate mode (Fig. 6) which are realized in the system heated from below.

Figure 2: The characteristics of the benzine–air system 2.5 s after switching on the bottom heaters: temperature distributions and isotherms (a); temperature distributions near the interface and interface shape, the formation of thermocapillary flexures in the zones of thermal impact (b); velocity field, the appearance of double-vortex motion above the heaters (c); tangential velocity vs at the interface (d); vapor concentration field in the gas phase and concentration isolines, the appearance of specific strips with elevated vapor content (e); evaporation mass flow rate M (f)

Figure 3: The characteristics of the benzine–air system 5 s after switching on the bottom heaters: temperature distributions and isotherms (a); temperature distributions near the interface and interface shape, the thermocapillary-driven dimples (b); velocity field, the transition to the multicellular regime upon heating (c); tangential velocity vs at the interface, the augmentation of the motion intensity (d); vapor concentration field in the gas phase and concentration isolines, the forked concentration patterns above the thermal elements (e); evaporation mass flow rate M (f)

Figure 4: The characteristics of the benzine–air system 15 s after switching on the bottom heaters: temperature distributions and isotherms (a); temperature distributions near the interface and interface shape, marked changes in the interface curvature accompanied by the wave behavior of the surface (b); velocity field, the reverse transition to double-vortex motion in the domains of local heating (c); tangential velocity vs at the interface (d); vapor concentration field in the gas phase and concentration isolines, the global “fork” concentration pattern in the gas layer (e); evaporation mass flow rate M (f)

Figure 5: The characteristics of the benzine–air system 64 s after switching on the bottom heaters: temperature distributions and isotherms (a); temperature distributions near the interface and interface shape, the decay of the thermocapillary dimple above “weak” Heater 1 (b); velocity field, the degeneration of the multicellular regime (c); tangential velocity vs at the interface (d); vapor concentration field in the gas phase and concentration isolines, “club-shaped” concentration structure (e); evaporation mass flow rate M, the intensification of evaporation above Heater 2 (f)

Figure 6: The characteristics of the ultimate convective regime in the benzine–air system 200 s after switching on the bottom heaters: temperature distributions and isotherms (a); temperature distributions near the interface and interface shape, the stable thermocapillary flexure in the heating zone (b); velocity field, the eventual dual-vortex regime (c); tangential velocity vs at the interface (d); vapor concentration field in the gas phase and concentration isolines (e); evaporation mass flow rate M (f)

We note that the convective mechanism completely defines the flow structure and character of changes in the basic characteristics. The heat transfer from the thermal elements triggers cascade responses affecting the system dynamics. The onset of motion is caused by the buoyancy force; here, the topological pattern of the flow is determined by the number and size of the heaters. The quantity of vortices forming typical convective cells at the first stage is equal to 2n, where n is the number of thermal elements. All the swirls have pairwise opposite circulation. Ascendant currents transferring heat to the interface are formed in the domains of thermal impact directly above the thermal elements. The velocity of rising convective flows is much higher than the velocity of the diffusive transfer in any direction (this is indicated by the structures of temperature and velocity fields in Figs. 2a, 2c). Once the heat reaches the interface, there appear surface effects. The surface tension of the liquid changes and the interface is distorted (Fig. 2b). The motion of the liquid along the interface occurs due to the action of the Marangoni force (the tangential velocity vs in Fig. 2d characterizes the intensity of thermocapillary spreading, and the dimensional values of the velocity are calculated by multiplying by u∗ from Table 2). Also, we note that the profile of vs at the interface allows one to determine the rotation direction within the appearing vortices: the positive and negative values of vs correspond to the clockwise and counterclockwise motion, respectively. We observe an increase in the evaporation rate above the regions where the ascendant flows transfer the hot liquid to the surface; here, vapor condensation occurs in the domains of down-flows (Fig. 2f). In the gas phase, vapor diffusion results in the formation of peculiar concentration structures (Fig. 2e). At the initial stage, the concentration strips appear above each of the heaters in the gas layer. The transverse size of these strips depends on the size of the thermal element. The larger the heating zone is, the wider is the strip. Here, heavier vapor diffuses from the warmed interface to the cold upper wall due to the thermodiffusion effect. Therefore, the maximum concentration gradients are vertically directed.

The initial convective regime is replaced by a multicellular mode whose characteristics are shown in Fig. 3. The new state also represents a transient mode which is characterized by a more complex flow structure with a higher velocity of motion (compare the values of vs in Figs. 2 and 3). The action of the Marangoni forces results in the alteration of the temperature field in the liquid (Fig. 3a), and there follows the transformation of the topological pattern of the flow (Fig. 3c). There appears a near-surface boundary layer in the heating regions, which will be retained further. The permanent heat supply from the heaters enhances the interface deformations and evaporation rate. We observe pronounced thermocapillary flexures above the heating zones and convex menisci (Fig. 3b). The number of molecules transferred from the heater to the liquid surface possessing energy to disrupt intermolecular bonds increases, therefore evaporation becomes more intensive (Fig. 3f) and the vapor content in the gas phase grows (Fig. 3e). All the characteristics of the bilayer system increase monotonically until the change-over of Heater 2.

The jump-like change in the temperature of Heater 2 leads to the development of oscillatory regimes with the behavior of the system characterized by the oscillations of the basic parameters near certain values. The wave behavior of the interface distinguishes these regimes from the first two modes (Fig. 4). The liquid–gas boundary undergoes oscillations, accompanied by periodical irregular variations of the surface curvature over the entire length (Fig. 4b). We observe a gradual transformation of the vortices caused by the temperature difference of the heaters. The quadruple-vortex motion above each thermal element degenerates into a two-fold one; here, the cells above the hotter heater have several deformed cores with more intensive motion (Figs. 4c, 4d). It explains the higher evaporation rate in the action zone of the second heater (Fig. 4f). The “fork” structure of the concentration field in the gas phase is typical for this mode (Fig. 4e). This oscillatory regime is also a transient one. With time it evolves into the next mode in which the oscillatory behavior of all the characteristics persists. Strong heating under the exposure to Heater 2 attenuates the surface tension of the interface and it becomes cambered (Fig. 5b) and continues to oscillate. Every change-over of Heater 2 results in a decrease in the liquid thickness in the zone of the thermal impact. Therefore, the velocity of convective rising flows grows. It entails an increase in the velocity of thermocapillary spreading and intensification of evaporation (Figs. 5d, 5f). The “club-shaped” concentration structure followed by the “fork” pattern is preserved in the gas phase (Fig. 5e).

After switching off Heater 1, the double-vortex motion in each fluid is settled. The flow occupies almost the entire cuvette. There is one zone with an ascendant flow above Heater 2 (Fig. 6c), where the maximum evaporation rate is reached (Fig. 6f). The typical thermocapillary flexure persists here (Fig. 6b). We note that the most intensive vapor condensation occurs in the vicinity of the lateral wall near the operating heater. The whole interface continues to oscillate at a certain finite frequency, but the pronounced surface waves propagating from the heating zone to the periphery vanish. The regime is ultimate; here, the maximum values of the flow velocity and evaporation rate also oscillate near certain values. The oscillatory character of convection is induced by the first abrupt change in the temperature of Heater 2, and it is retained even when the intensity of thermal load stops changing. Thus, the principal mechanism which governs the characteristics of the appearing transient and eventual modes is a convective one.

3.3 Heating Mode II

In contrast to the case of local heating from below, the dynamical characteristics as well as the behavior of the liquid–gas interface do not qualitatively change with time upon lateral heating. Figs. 7–9 allow one to track the dynamics of the basic characteristics when the bilayer system undergoes non-stationary lateral heating. One can see that boundary thermal and concentration layers (images (a) and (e) in Figs. 7–9) develop close to the side walls. The structure of the corresponding fields is explained by a slight diffusive transfer in any directions in comparison with the convective transfer. This directly follows from the heat and molecular transport equations. Given the values of the Reynolds and Prandtl numbers for the liquid layer (at j=2 ), one can see that the term in the right side of the energy equation (the third relation in (1)) related to the heat diffusion has the value which is by an order of magnitude lower than that of the convective terms in the left side. In the gas phase (at j=1 ), the Re and Pr numbers are comparable. Thus, in the absence of essential density stratification caused by temperature drops in the vertical direction, heat is transferred due to heat conductivity more slowly in comparison with the case of heating from below. This effect is also typical for vapor diffusion in the gas layer, resulting in the formation of peculiar wall concentration structures with the elevated vapor content near the upper boundary (images (e) in Figs. 7–9). The occurrence of the patterns is driven by the Soret effect, causing the diffusion of heavy benzine vapor from the hot domains into the regions with a lower temperature.

Figure 7: The characteristics of the benzine–air system at t=2.5 s in lateral heating: temperature distributions and isotherms (a); temperature distributions near the interface and interface shape, the formation of thermocapillary bump (b); velocity field, the appearance of the boundary layers near the side walls (c); tangential velocity vs at the interface (d); vapor concentration field in the gas phase and concentration isolines (e); evaporation mass flow rate M (f)

Figure 8: The characteristics of the benzine–air system at t=64 s upon lateral heating: temperature distributions and isotherms (a); temperature distributions near the interface and interface shape (b); velocity field (c); tangential velocity vs at the interface (d); vapor concentration field in the gas phase and concentration isolines (e); evaporation mass flow rate M (f)

Figure 9: The characteristics of the ultimate convective regime in the benzine–air system at t=200 s upon lateral heating: temperature distributions and isotherms (a); temperature distributions near the interface and interface shape (b); velocity field (c); tangential velocity vs at the interface (d); vapor concentration field in the gas phase and concentration isolines (e); evaporation mass flow rate M (f)

The specific temperature distribution entails other behavior of the interface. The phase boundary has a convex upward form in response to lateral heating. It results from the combined action of the surface tension forces and arising pressure gradients. The lateral heating attenuates the interface surface tension. At the same time, the pressure from the gas side near the cuvette edges is higher than that from the liquid, and therefore, the phase boundary hangs down therein. Out of the thermal “spots”, the pressure drop vanishes and the interface is forced out upwards. It leads to the formation of a thermocapillary bump. The higher the intensity of the external thermal load, the larger the pressure drop and the larger the height of the interface rise (compare images (b) in Figs. 7–9).

As to the velocity field pattern, we note that the wall swirls near both side boundaries evolve into the vortex structures with a complex symmetry. The eventual flow regime is characterized by the coexistence of pronounced boundary layers with a near-wall core and a slow-convecting zone. Here, the velocity of motion induced by the distributed lateral heating is well above those in the case of the local thermal load from the substrate (compare the values of vs in Figs. 6d and 9d). The transverse size of the arising structures is about one-third of the cuvette width so that motionless domains persist within the central part of the fluids. The evolution of the velocity field is exemplified by the images (c) in Figs. 7–9. Thus, the typical pattern of the convective mode occurring upon the heat supply from the side walls settles immediately after the start of heating and develops monotonously.

3.4 Behavior of the Contact Angle

The behavior of the dynamic contact angle is illustrated by an example of the angle φ on the right lateral wall. One chooses the three-phase contact point closest to operating Heater 2, since the interface undergoes the largest deformations in the zone of the temperature action produced by the thermal element. We see that the local heating from below causes quite little variations of the contact angle (see the zoomed curve reflecting the changes in the contact angle with time in Fig. 10a, and dashed line in Fig. 10b). The enhancement of the thermal load leads to a gradual decrease in the contact angle. Here, there is a certain time lag when the system responds to changes in the intensity of thermal head. Even when the heater temperature q2s stops changing, the contact angle continues to vary up to a certain limiting value. Then, the interface relaxes monotonically. It is accompanied by an increase in the contact angle and by the establishment of a new equilibrium value of φ which is lower than the initial angle φ=90∘. Thus, a concave meniscus is formed in the ultimate mode under local heating from below. Small variations of φ are explained by small changes in the surface energies on the interfaces between the phases determining the equilibrium contact angle. The character of changes in the surface tension of the liquid–gas boundary σ has the largest impact. The surface tension drastically attenuates in the action area of the heater (it is clear from both the formula for the σ function and the interface shape), whereas it is hardly varied near the edges. It explains the observed behavior of the contact angle upon the local thermal load from below.

Figure 10: The dynamic contact angle φ under local heating from below (a) and distributed thermal load from the side walls (solid line) and local heating from below (dashed line) (b)

Another character of the contact angle alterations is observed under lateral heating. The intensive supply of heat from the side walls results in a rapid significant decrease in the φ angle (Fig. 10b). The process is a sort of “negative” feedback of the system to changes in the balance between the convective heat transfer and the interfacial thermal effects caused by the applied thermal load. It discontinues due to gradual equalization of the surface temperature near the cuvette edges, and there occurs a subsequent increase in the contact angle up to the moment of stabilization of all the system parameters after setting the constant temperature TR on the lateral wall. The convex meniscus with φ>90∘ persists in the eventual regime under Heating mode II.

4  Concluding Remarks

The dynamics of the two-phase system with a deformable interface admitting the mass transfer due to evaporation in the box-shaped cuvette subjected to various heating conditions is investigated at an example of the system of working media like benzine–air. The study is carried out in the framework of a full two-sided model based on the Boussinesq approximation. The original approach used to state the problem in terms of new auxiliary functions allows one to derive the corresponding numerical model that enables one not only to calculate all the system characteristics, including the mass evaporation rate, but also to scrutinize the behavior of the dynamic contact angles, depending on the heating mode.

Based on the numerical simulation, considerable differences in the behavior of the bilayer fluidic system are elucidated. An eventual regime under the local heating from below is established through the transient modes, whereas the inertia of the system behavior is typical for the case of the distributed lateral heating, giving rise to monotonous changes in the hydrodynamic, temperature and concentration characteristics. Here, the amplitudes of the interface deformations under the thermal load applied on the side walls is lower than those under the other type of heating. Two various scenarios of the interface evolution are realized, depending on the type of the external thermal action. A thermocapillary flexure and a thermocapillary bump appear upon heating from the substrate and side walls, respectively. In the first case, convective heat and mass exchange occurs in the whole volume of the fluids, and equilibrium concave menisci are formed on the edges in the ultimate steady regime. The second heating regime leads to the appearance of a boundary layer and near-wall temperature and concentration patterns; here, convex menisci are observed on the heated lateral walls. In both cases, the contact angle variation in time is described by the numerical model.

The results obtained contribute to a better understanding of the interaction of the basic mechanisms governing the convective regimes in systems with internal interfaces and points of three-phase contact. The results of the present study can aid in designing equipment for thermal coating or drying and working sections employed in control tools for various applications as well as in developing methods for the formation of patterns of the required structure and morphology on solid surfaces.

Acknowledgement: The authors express gratitude to Dr. O.M. Lavrenteva for useful discussions.

Funding Statement: The work of O.N. Goncharova was carried out in accordance with the State Assignment of the Russian Ministry of Science and Higher Education entitled “Modern Models of Hydrodynamics for Environmental Management, Industrial Systems and Polar Mechanics”, Govt. Contract Code: FZMW-2024-0003, https://minobrnauki.gov.ru/.

Author Contributions: The authors confirm equal contribution to the paper: study conception and design, data collection, analysis and interpretation of results, draft manuscript preparation: V.B. Bekezhanova. O.N. Goncharova. All authors reviewed the results and approved the final version of the manuscript.

Availability of Data and Materials: All the data are based on publicly available documents, whose data sources were cited in the manuscript.

Conflicts of Interest: The authors declare that they have no conflicts of interest to report regarding the present study.

References

1. Bekezhanova VB, Fliagin VM, Goncharova ON, Ivanova NA, Klyuev DS. Thermocapillary deformations of a two-layer system of liquids under laser beam heating. Int J Multiphas Flow. 2020;132:103429. [Google Scholar]

2. Al-Muzaiqer MA, Kolegov KS, Ivanova NA, Fliagin VM. Nonuniform heating of a substrate in evaporative lithography. Phys Fluids. 2021;33:092101. [Google Scholar]

3. Zuev AL, Kostarev KG. Certain pecularities of solutocapillary convection. Phys-Uspekhi. 2008;51(10):1027–45. [Google Scholar]

4. Farzeena C, Vinay TV, Lekshmi BS, Ragisha CM, Varanakkottu SN. Innovations in exploiting photo-controlled Marangoni flows for soft matter actuations. Soft Matter. 2023;19(28):5223–43. [Google Scholar] [PubMed]

5. Feldmann D, Maduar S, Santer M, Lomadze N, Vinogradova O, Santer S. Manipulation of small particles at solid liquid interface: light driven diffusioosmosis. Sci Rep. 2016;6(1):36443. [Google Scholar] [PubMed]

6. Routh AF, Russel WB. Horizontal drying fronts during solvent evaporation from latex films. American Institute of Chemical Engineers Journal. 1998;44:2088–98. [Google Scholar]

7. Bekezhanova VB, Goncharova ON. Problems of the evaporative convection (review). Fluid Dyn. 2018;53(1):S69–102. [Google Scholar]

8. Kolegov KS, Barash LY. Applying droplets and films in evaporative lithography. Adv Colloid Interfac. 2020;285:102271. [Google Scholar]

9. Prevo BG, Hon EW, Velev OD. Assembly and characterization of colloid-based antireflective coatings on multicrystalline silicon solar cells. J Mater Chem. 2007;17(8):791–9. [Google Scholar]

10. Hatton B, Mishchenko L, Davis S, Sandhage KH, Aizenberg J. Assembly of large-area, highly ordered, crack-free inverse opal films. Proc Natl Acad Sci. 2010;107(23):10354–9. [Google Scholar] [PubMed]

11. Hammond PT. Building biomedical materials layer-by-layer. Mater Today. 2012;15(5):196–206. [Google Scholar]

12. Du F, Felts JR, Xie X, Song J, Li Y, Rosenberger MR, et al. Laser-induced nanoscale thermocapillary flow for purification of aligned arrays of single-walled carbon nanotubes. Am Chem Soc Nano. 2014;8(12):12641–9. [Google Scholar]

13. Ivanova N, Starov VM, Trybala A, Flyagin VM. Removal of micrometer size particles from surfaces using laser-induced thermocapillary flow: experimental results. J Colloid Interface Sci. 2016;473:120–5. [Google Scholar] [PubMed]

14. Chhasatia VH, Sun Y. Interaction of bi-dispersed particles with contact line in an evaporating colloidal drop. Soft Matter. 2011;7(21):10135–43. [Google Scholar]

15. Wyatt Shields CIV, Reyes CD, Lopez GP. Microfluidic cell sorting: a review of the advances in the separation of cells from debulking to rare cell isolation. Lab Chip. 2015;15(5):1230–49. [Google Scholar] [PubMed]

16. Takhistov P, Chang HC. Complex stain morphologies. Ind Eng Chem Res. 2002;41(25):6256–69. [Google Scholar]

17. Al-Muzaiqer MA, Ivanova NA, Fliagin VM, Lebedev-Stepanov PV. Transport and assembling microparticles via Marangoni flows in heating and cooling modes. Colloids Surf A: Physicochem Eng Asp. 2021;621(5):126550. [Google Scholar]

18. Chen L, Wang Y, Qi C, Tang Z, Tian Z. A review of methods based on nanofluids and biomimetic structures for the optimization of heat transfer in electronic devices. Fluid Dyn Mater Process. 2022;18(5):1205–27. doi:https://doi.org/10.32604/fdmp.2022.021200. [Google Scholar] [CrossRef]

19. Helseth LE, Fischer TM. Particle interactions near the contact line in liquid drops. Phys Rev E. 2003;68(4):042601. [Google Scholar]

20. Ajaev VS, Kabov OA. Heat and mass transfer near contact lines on heated surfaces. Int J Heat Mass Tran. 2017;108:918–32. [Google Scholar]

21. Dussan EBV, Davis SH. On the motion of a fluid-fluid interface along a solid surface. J Fluid Mech. 1974;65:71–95. [Google Scholar]

22. Pukhnachov VV, Solonnikov VA. On the problem of dynamic contact angle. J Appl Math Mech. 1982;46(6):771–9. [Google Scholar]

23. Baiocchi C, Pukhnachov VV. Problems with unilateral constrains for the Navier-Stokes equations and the problem of dynamic contact angle. J Appl Mech Tech Phy. 1990;31(2):185–97. [Google Scholar]

24. Kroener D. The flow of a fluid with a free boundary and dynamic contact angle. ZAMM J Appl Mat Mech. 1987;67(5):T304–6. [Google Scholar]

25. Kroener D. Asymptotic expansion for a flow with a dynamic contact angle. In: Heywood J, Masuda K, Rautmann R, Solonnikov VA, editors. The Navier–Stokes Equations: Theory and Numerical Methods. Berlin: Springer; 1990. p. 49–59. [Google Scholar]

26. Rednikov AY, Colinet P. Contact angles for perfectly wetting pure liquids evaporating into air: between de Gennes-type and other classical models. Phys Rev Fluids. 2020;5:114007. [Google Scholar]

27. Fricke M, Marić T, Bothe D. A kinematic evolution equation for the dynamic contact angle and some consequences. Phys D: Nonlinear Phenom. 2019;394:26–43. [Google Scholar]

28. Doerfler W, Goncharova O, Kroener D. Fluid flow with dynamic contactangle: numerical simulation. ZAMM J Appl Math Mech. 2002;82(3):167–76. [Google Scholar]

29. Goncharova ON, Marchuk IV, Zakurdaeva AV. Investigation of behavior of the dynamic contact angle in a problem of convective flows. Eurasian J Math Comput Appl. 2017;5(4):27–42. [Google Scholar]

30. Ajaev VS, Gambaryan-Roisman T, Stephan P. Static and dynamic contact angles of evaporating liquids on heated surfaces. J Colloid Interf Sci. 2010;342:550–8. [Google Scholar]

31. Karchevsky AL, Cheverda VV, Marchuk IV, Gigola TG, Sulyaeva VS, Kabov OA. Heat flux density evaluation in the region of contact line of drop on a sapphire surface using infrared thermography measurements. Microgravity Sci Technol. 2021;33:53. [Google Scholar]

32. Nikolayev VS, Gavrilyuk SL, Gouin H. Modeling of the moving deformed triple contact line: influence of the fluid inertia. J Colloid Interf Sci. 2006;302:605–12. [Google Scholar]

33. Bekezhanova VB, Goncharova ON. Thermodiffusion effects in a two-phase system with the thermocapillary deformable interface exposed to local heating. Int J Multiphas Flow. 2022;152:103429. [Google Scholar]

34. Andreev VK, Gaponenko YA, Goncharova ON, Pukhnachov VV. Mathematical models of convection. Berlin: De Gruyter; 2020. p. 1–432. [Google Scholar]

35. Shliomis MI, Yakushin VI. Convection in a binary system with evaporation. Uchenye Zapiski Permskogo Universiteta. Seriya: Gidrodinamika. 1972;4:129–40 [In Russian]. [Google Scholar]

36. Vargaftik NB. Handbook of thermophysical properties of gases and liquids [In Russian]. Moscow: Nauka; 1972. [Google Scholar]