STRUCTURAL CHANGES OF A RADIALLY SYMMETRIC THERMOCAPILLARY FLOW IN THE SHALLOW CAVITY

. This paper presents a numerical study of spatial transformations of a radially symmetric ﬂow in a shallow ﬂuid-ﬁlled cylindrical cavity partially covered by a solid non-deformable ﬁlm. The upper central part of the liquid has a free surface where an intense light beam is focused to produce a hot spot on the symmetry axis. The heating generates a divergent thermocapillary motion on the free surface, which causes the ﬂuid to ﬂow under the immovable solid ﬁlm. The edge of this ﬁlm induces thermal perturbations, which, at speciﬁc heat generation values, begin to increase and give rise to a gradual radial ﬂow symmetry breakdown that visually demonstrates the onset of vorticity in the azimuthal plane. Three-dimensional calculations have been carried out based on the interfacial hydrodynamics equations using Comsol Multiphysics software. The results of numerical calculations conﬁrm that the instability occurs in the area corresponding to the boundary between the free surface and the solid ﬁlm. The motion in the azimuthal direction becomes more evident with the growth of heating intensity. The vorticity in the azimuthal plane makes the ﬂow structure signiﬁcantly more complex compared to the axisymmetric radial ﬂow. Thus, there is a predominantly radial ﬂow in the area of free surface and, at the same time, the vortices with the azimuthal velocity component are observed under the ﬁlm. As in the experiment, the number of vortices is determined by the ratio of the area of free surface to the covered area. It is shown that the appearance of motion in the azimuthal direction depends on the joint action of several physical mechanisms, which have been initialized in the theoretical model by means of various piecewise thermal and mechanical conditions at the upper boundary.


Introduction
Interest in the fundamental problems of physics of interfacial phenomena has increased significantly over the last few decades.Moreover, the study of heat and mass transfer generated by surface effects is of great practical importance.Successful manufacturing applications in this field are directly related to the synthesis of different sciences on micro-and macroscopic levels.The necessity for gaining deeper theoretical insights into interfacial processes is dictated by active investigations in physical chemistry [26], nanotechnology and condensed matter science [12,20], continuous media mechanics [3,15] and thermal physics [10].Acquiring a new knowledge base will allow finding new solutions to existing technological problems.Therefore, the application of interfacial phenomena requires formulation of the problems in terms of fundamental physics, which ensures better understanding of the basic features of the processes.
For this purpose, some specific model objects are considered, among which an infinite fluid layer in interfacial hydrodynamics occupies a special place because of its sufficient simplicity and geometry commonality.This somewhat artificial statement leads to the situation when the system is studied by means of particular mathematical tools of continuous media mechanics using the theory of differential equations and the method of asymptotic expansions.Such an approach permits solving governing equations in the main approximation.The famous papers of Pearson [18], Nield [17] and Birikh [1] are devoted to the investigation of interfacial convective instability in the infinite horizontal layer where the basic properties of heat and mass transfer are considered in relation to the thermocapillary effect.The research of Sen and Davis [21], a logical continuation of the above mentioned approach, focuses on finding a solution of the thermocapillary problem for a closed horizontal layer.Nowadays, there are lots of applications that use analytical methods to describe the interfacial hydrodynamic phenomena with such additional complicating factors as evaporation [13], surface contamination [6] or electroconvection effects [2].
It is worth noting here that the application of asymptotic methods enables one to determine a basic solution, analyze its stability and evaluate the non-linear effects in relation to the expansion order.Moreover, these techniques can be used to describe the Benard-Marangoni convection in the case of a deformable surface [19].There is a wide range of issues on instability control via a modulation of a parameter [16].Mikishev and Nepomnyashchy [14] demonstrated that, for most of these problems, the solution gives the boundaries of various stability regimes and makes it possible to visualize the competition between the complex spatial patterns induced by nonlinear effects.
Sometimes in practice we have to deal with different types of liquid-gas interface contamination.Even a small amount of foreign substance can produce a significant impact on the macroscopic behavior of the system.On the other hand, the introduction of an equation for surface concentration and a practical set of media parameters into the mathematical model makes the problem very complicated.Therefore, it is not sufficient to use only analytical methods for investigation of the coupled dynamics of a surfactant and a parent liquid.One similar phenomenon, studied by the method of asymptotic expansion, is described in the papers of Homsy and co-authors [5,11], where the film of insoluble surfactant partially covers the surface of a liquid in a shallow two-dimensional rectangular cavity with a movable upper boundary.Because of the non-uniformity of temperature distribution over the surface, there occurs a thermocapillary effect that forces the film to shift towards the cold area of the cavity.This results in the non-uniformity of a surfactant concentration at the interface, which gives rise to the concentration-capillary mechanism, which resists the initial flow propulsion.It has been found that certain conditions arise when the interface is divided into two distinct regions (with and without surfactant).The point separating the surface into these two parts is called a stagnant point.Shmyrov and Mizev [24,25] investigated experimentally this phenomenon using the PIV method.The dependence of the stagnant point position on the elasticity number (the parameter determining the ratio between the thermo-and concentration-capillary forces) is the main characteristic used to make a comparison between experimental and theoretical results.The disagreement between the results of these experiments and Homsy's theory is due to the fact the theory does not take into account the temperature profile nonlinearity occurred during convection and the possibility of a phase transition in a surfactant film.The mathematical model described in [25] includes all these factors.The numerical solution by Demin and Petukhov [25] gave the results similar to those obtained experimentally and provided an opportunity to describe the behavior of a surfactant in a wide range of governing parameters.In [7,8], the dynamics of such a surfactant was studied during the relaxation process performed to achieve a uniform distribution after the heating stops.It was demonstrated that the surfactant film, propagating along the free surface, behaves itself in an appreciably nonlinear manner.First, it accelerates, but then, because of the complicated temperature distribution at the interface, it begins to move downwards.Only after this deceleration, the free window sharply collapses.At the final stage, the surfactant spreads uniformly over the whole surface.This effect can be identified as a rapidly damped oscillation.It is important to note that this phenomenon, initially observed experimentally, has been described in detail using only numerical methods.
In [6], the effect of solubility of a surfactant in a parent liquid on its concentration distribution at the interface was investigated in the experiments conducted in a vertical Hele-Shaw cell heated non-uniformly from above.It was shown that a concentration "tongue" may occur in the course of surfactant ablation by the steady large-scale flow.
One more complexity of the problem is to consider the process of surfactant redistribution as a two dimensional distribution with two degrees of freedom.Shmyrova investigated this issue experimentally in [22].Unlike other studies, the working cavity had a cylindrical shape.The point heat source was in the geometrical center of the free surface.Initially, the surfactant covered the entire free surface of the liquid.After the heater intensity reached a certain threshold, the radial thermocapillary flow shifted the film to the cold external boundary of the cavity, purifying the central hot region which formed the sharp stagnant line.Shmyrova analyzed the volumetric flow shape occurred due to its interaction with a mobile film.At the sight from above, the flow under the free surface had a uniform radial structure, indicating the fluid motion towards the cavity periphery.When the thermocapillary flow reached the stagnant line, it moved under the film edge, and the velocity field unexpectedly got an azimuthal component.Thus, the flow was split into a series of vortices in the horizontal plane.Note that the same mechanism of vortex generation was observed on the surface of a bubble located in an incoming flow [23].Noticeably different statements demonstrate general features of a hydrodynamic instability in the fluid flow.As far as we know, there are no theoretical studies involving a realistic and quantitative 3D description of the behavior of a fluid under surfactant.In this paper, we use numerical simulations to find a possible explanation for the occurrence of secondary vortices under the surfactant film after the basic radial thermocapillary flow breaks up.

Problem statement
In this study, we apply the same design scheme for the experiment as before [22].So, we consider a cylindrical shallow cavity with partially open upper boundary.The working liquid is water characterized by a sufficiently strong temperature dependence of the surface tension coefficient.The cavity of radius R and height h has isothermal boundaries, see Figure 1.There is an optical heating in the geometrical center of the free surface so that the resulting hot spot has a characteristic radius the value of which is dependent on the energy distribution of the light beam.
In [7,8,24,25], analysis of the one-dimensional motion of the film on the surface between the wide boundaries of the vertical Hele-Shaw cell made it possible to solve the problem with a mobile insoluble surfactant in full statement, taking into account both thermo-and concentration capillary mechanisms at the interface.This approach provided the opportunity to study the time evolution of the longitudinal deformation of a film.It was shown that the joint action of thermo-and concentration capillary forces in all these cases led to a stationary state with a complex nonlinear distribution of surfactant concentration.Even a phase transition from the gas state to the liquid phase was described sufficiently well.This indicates that there is a good agreement between the experiment and the simulation results.
In the case of a cylindrical cavity, the motion of a surfactant is considered as a two-dimensional motion, while the volumetric flow of a parent liquid as a three-dimensional flow.Besides, the results of numerical simulation show that the magnitude of the longitudinal velocity component is much higher at the free boundary than under the surfactant film.The extremely high velocity gradients exist near the stagnant line and under the film.Under these circumstances, if we want to get any productive result, the problem statement should be simplified.Accordingly, the limiting case of a solid non-compressible immovable film with a no-slip condition for the flow velocity is considered.
This approach gives a possibility to describe the flow of this hydrodynamic system using the classical equations of thermal gravitational convection [9] with addition of a thermocapillary force at the free surface.Namely, under the Boussinesq approximation, we can write the Navier-Stokes equation in dimensional form as: Then, to describe the convective system, we introduce the equation for heat transfer: and the mass conservation equation (known as a continuity equation) for incompressible fluid: Here, v, p and T are the dimensional velocity and pressure fields and the deviation of temperature from the reference value, respectively; ν and χ are the kinematic viscosity and thermal diffusivity coefficients, ρ is the mean density of fluid, β is the thermal expansion coefficient, g is the gravity acceleration vector counterdirected to the z-axis.According to the above mentioned assumptions, the surfactant film is immovable and, hence, there are no governing equation for it.
Let us now specify strict boundary conditions for the system of equations (2.1)-(2.1).The solid film has the form of a ring (Fig. 1), which is placed at the periphery of the upper boundary and touches the cavity wall.The condition of the surfactant film, when it is shifted by the flow from the centre to the border, can be considered as its final state.Thus, the region of surfactant-free interface has a constant radius R * .The nonuniformity of surface tension temperature initiates thermocapillary flows on the free surface.In the experiments, the characteristic temperature difference was estimated as 1-2 K, which allows using a linear approximation for surface tension: where σ 0 is the surface tension at the reference temperature, σ T is the coefficient of temperature dependence of surface tension.This law yields the following tangential stress balance on the free surface: In this boundary condition, η is the coefficient of dynamic viscosity.
In the experiment, the heating was implemented using the optical contactless method, and therefore the value of energy flux across the surface has to diminish smoothly as a function of the radial coordinate.Thus, the temperature distribution is determined by the heat flux condition according to the formula: Using (2.5), the heating can be realized so that the thermal flux gradually decreases with time, see Figure 2.
The amplitude A determines the maximum value of heating intensity.The parameters k 1 and k 2 refer to the characteristic radius of a hot spot and the duration of a transient heat stage, respectively.We present two types of computational domains, of which one comprises the whole cylindrical cavity (general case) and the other corresponds to the segment with an apex angle 90 • .The first domain is considered to demonstrate the rough effect of fragmentation under the film on the vortices with a definite wave number.The problem shown in the second domain is solved to quantitatively analyze the inclination of a vortex plane with the greatest possible accuracy.The standard no-slip condition is imposed on the velocity at all solid boundaries of the cavity in both cases (Fig. 1).It is assumed that the film, bottom and vertical wall of the cylindrical cavity are isothermal, and that their temperature is kept constant.Thus we can write As one can see (the segment in Fig. 1b), there are additional conditions for non-penetration of velocity on the vertical boundaries I and II (non-zero tangential velocity components): ) This indicates the absence of the flow across these segments for the tangential momentum component.Different variants of possible temperature conditions are given below in Section 3, separately for each segment boundary (I and II ).

Numerical procedure
3), describing the heat-and mass transfer of the convective system, together with boundary conditions (2.4)-(2.9),were solved numerically using Comsol Multiphysics software.The calculations were performed on the computing cluster Triton at the Institute of Continuous Media Mechanics UB RAS.
Numerical results showing good convergence were obtained for the chosen number of cells.All simulations were carried out with MUMPS (multifrontal massively parallel sparse direct solver).In the full computational domain, the total number of cells varied in the range N = 1 40 300 ÷ 1 52 000 for different values of the free zone radius R * = 3 ÷ 5 cm.In the segment domain, the number of cells was in the range N = 2 51 200 ÷ 2 61 700.
The stability of the calculation scheme was verified by evaluating the smallness of the Courant number.In physical terms, the Courant condition means that, at definite temporal discretization, a fluid particle should not move more than one spatial step over the time interval.The Courant number is computed for all cells of the calculation domain at each instant of time.
The analytical study of Homsy [5] shows that the velocity field is continuous in the nearest vicinity of the film edge.Indeed, it is possible to get numerical instability near this zone during calculations.To avoid it, we create a non-uniform tetrahedral mesh with crowding near the centre of the interface and near the film edge (Fig. 3) and use the moderate values of governing parameters.Thus, the velocity gradients in the nearest vicinity of a stagnant line are not too high, and the streamlines are continuous in this region.
In the process of numerical simulation, the initial conditions (zero velocity and temperature fields) are not related to the final structure of vortices.At the same time, there occur convective patterns with definite azimuthal periodicity, which show time-dependent stability.Our simulations indicate that the secondary flow with vortices in a horizontal plane is stable for the chosen set of parameters.

Limiting cases
We start by considering two special situations in the framework of our statement.These limiting cases represent the liquid, completely covered by the film, and the clear free surface.For the completely covered surface, boundary condition (2.5) was modified, i.e., it was applied to the entire upper boundary.The structure of these flows, obtained via numerical simulation, is given in Figure 4.By analyzing the shape of streamline  projections, one can conclude that the radial flow with only poloidal vorticity occurs in the velocity field, i.e., no azimuthal component exists in the flow.As seen from above (top view), there is no qualitative difference between these two limiting cases.The thermal gravitational and thermocapillary flows (Fig. 4a, b) spread in the radial direction from the hot centre of the cavity to the cold solid boundary at its periphery.The side view shows that the streamline projections slightly differ from each other by only the vortex center position.In the cavity covered by the film, it is located at the point with the radial coordinate r ≈ R/2, and, in the open cavity, it is shifted to the wall.Almost the same effect was obtained in [25] for the flow in a Hele-Shaw cell.It should be noted that, in the experiment described in [22], the azimuthal instability developed at the final stage of the transient process.At this stage, the thermocapillary flow compressed the surfactant film until the stationary state was reached.Obviously, the description of the effects associated with the initial mobility of the film is not included in the problem statement.It can be supposed that our simplification has no impact on the final state of the convective system regardless of different character of the flow interaction with the obstacle on its way.Hence, we can conclude that, in the limiting cases under consideration, the radial flow is the only consequence of the obstacle absence at the surface, i.e., there is no any other reason why the instability appears in the flow.The numerical simulation for the full cylindrical cavity (Fig. 1a) and for segment (Fig. 1b) shown below confirms this deduction.

Full cylindrical cavity
Now, we analyze the intermediate case when the immovable film covers the surface partially.A relatively intense thermocapillary motion along the free zone is always observed even for the total weak heat flux.This is due to a relatively small radius of the cavity and the intense localized heating which produces high temperature gradients at the open surface.First, we consider vortex generation in the azimuthal plane under the film and then compare the total number vortices with those obtained experimentally.The vortex generation is caused in this case by the self-organization system in the full computational domain.Due to the formation of convective structures in the cylindrical cavity, we get approximately the same number of vortices for different values of the film width (Fig. 5).The radial spreading of a surfactant is seen along the free surface for R < R * .The noticeable difference viewed from above appears when the liquid flows under the film.The radial flow structure significantly changes because of the instability arising after the interaction of the flow with the film edge, which leads to the distortion of trajectories of the liquid particles.One can see that this instability is responsible for the occurrence of azimuthal components in the velocity field.Thus, the flow with almost periodic structure in the azimuthal plane develops under the film (Fig. 6).Analysis of the numerical results obtained in this study has revealed that the vortices tend to take the form of a natural circle with a characteristic diameter equal to the width of the covered zone.The flow separation in the form of vortices in the horizontal plane occurs as a result of a self-organizing process.The boundaries between vortices occur due to the radial direct and re-entrant flows of the fluid under the film.The azimuthal periodicity of these flows leads to the formation of vortices, the area of which depends on the position of a stagnant line.
In the limiting case of a radial flow, the change in the direction of a liquid particle was realized in a definitely poloidal section.For partially covered surface, the reverse current occurs in the inclined plane, and therefore the projection of motion can be seen on the horizontal plane.This allows the visualization of the vorticity in the azimuthal plane (Fig. 7a).The transition of the reverse flow to the horizontal section is energetically preferable for such system; due to this, the vortices are formed under the film with regular periodicity.
In the experiment, the value R * is determined by taking into account such parameters as cavity radius, heating intensity and insoluble surfactant concentration.When the equilibrium condition for thermal-and concentration-capillary forces is established, then the system reaches the steady state, promoting thus the formation of a definite integer number of vortex pairs.At the same time, one can observe the correspondence between the position of a stagnant line and the number of vortices.
In the simulation, R * is the chosen parameter, and therefore it is difficult to find at once its precise value for definite heat intensity needed to adequately simulate the flows similar to those obtained experimentally.Thus, only the final steady stages can be compared for the intermediate case 0 < R * < R.
In the above section, the existence of two limiting cases with symmetrical poloidal flows in the cavity was demonstrated.Note that in the calculations for intermediate values R * there must be one more continuous transition.This permits us to describe the transformation of the inclined vortex with the azimuthal components of the velocity field in the poloidal flow in relation to the gradually reduced thermal non-uniformity.The characteristic values of the inclination angle of the vortex plane with respect to the horizontal section as a function of the heat flux amplitude is given in Figure 7b.According to the results of numerical simulation, the growth of heat intensity A leads to a decrease in the inclination angle ϕ, which is measured from the horizon.In the limiting case A → 0, the angle ϕ → 90 • and the slow radial flow take place in the poloidal plane.For smaller values of R * it is easy to obtain the significantly inclined vortices even for the relatively weak thermal fluxes A. Due to the growth of R * , there occurs a necessity to increase the characteristic value A so that the observable inclination of the vortex plane can be kept.Finally, if R * tends to the R, then, for the radius bigger than the definite value, no inclination of the vortex plane can be detected.
A separate study is required to analyze the stability of the radial flow with respect to small perturbations; the full non-linear problem is already solved numerically.This approach also permits one to describe the transitions between different convective regimes.In our case, for the fixed free zone radius, the radial flow symmetry breakdown begins to be observed even for small values of the heat gradient, and the inclination of the vortex plane is directly proportional to the increase of heat intensity (Fig. 7b).Thus, in a certain sense, the heat gradient plays a role of any order parameter.An increase in this parameter leads to the symmetry breakdown of the basic radially symmetric flow.In the frame of our hydrodynamic problem, Figure 7b is an analogue of the bifurcation diagram where the heat gradient is a bifurcation parameter.
The article [22] is focused on the experimental investigation of the effect of vortices origination under the mobile surfactant film.The authors of [22] suppose that all vortices appear only in the nearest vicinity under the surface with negligible penetration into the volume.Therefore, the laser knife of PIV facility was oriented horizontally to illuminate the space of liquid extremely close to interface.The figure for solid mask, which partially covers the liquid, shows intensive reflection of laser knife from the inner part of this film.It is obvious from physical point of view, that there is no motion in this region because of no slip condition on the solid mask.Also, this figure demonstrates that initially radially-symmetric streamlines in the free region slightly deviate in horizontal section after interaction with stagnant line.In other words, we detect only fragments of these streamlines.The streamlines outgoing into the volume are not seen in this figure.Evidently, it is possible to identify vortices entirely by mean of laser illumination in whole volume.Unfortunately, in this article there is no investigation of the flow in whole space.It does not give us any quantitative data to check it in numerical simulation.

Results for the segment
The above section discusses the general features of structural changes in the radial thermocapillary flow, which are responsible for the occurrence of vorticity in the azimuthal plane.This effect has to be considered as a secondary convective behavior because of the interaction of the thermocapillary flow with the solid film.Now it is important to investigate in detail the form of these vortices through considering the small segment and by taking fine grating for this computation.Thus, the calculation domain is the segment with the apex angle of 90 • and with specific boundary conditions (2.8), (2.9) on the lateral sides.To avoid the need to solve the singularity problem on the symmetry axis, the angular part around this line is truncated for the working domain.Thermal properties on the lateral sides I and II correspond to the artificial conditions of two types: heat-insulated ( ∂T ∂n = 0) or isothermal boundaries (T = T 0 ), as for the film and solid walls.The results of the numerical simulation show an absolutely same scenario for the formation of azimuthal vortices.After the initial  radial divergent motion along the free surface the liquid reaches the obstacle (the edge of an immovable film).Then the trajectory starts to distort, thus giving rise to the azimuthal component in the velocity field.As in all previous cases, these vortices exist only under the film.The view from above of the complex three-dimensional flow is given in Figure 8a for the heat-insulated lateral boundaries.As one can see, there are four vortices with opposite signs of rotation.The vortices are slightly different in size as in the experiment in which the convective flows demonstrate the variety of structures and sensitivity to initial conditions.Certainly, despite the non-roughness of the system, the flow passing through the segment was a steady flow.Figure 9 presents the time dependence of the maximum value of velocity during the obvious transition to the stationary regime.Exactly the same dependencies are typical for all calculations described in this article, but it is interesting to note that the case of isothermal lateral sides shows an absolutely symmetric pair of vortices (Fig. 8b).Thermal perturbations disappear on these boundaries, which makes the flow more stable.Hence it follows that the symmetric convective pattern can be obtained without adjusting the heat flux and the open surface radius during calculations.

Conclusions
The calculations, performed in current investigation, allow to visualize the origination of azimuthal vortices under the film of insoluble surfactant in cylindrical cavity heating from above by the point source.As in experiment with the real film of insoluble surfactant our numerical simulation for the model case of solid up boundary reproduces these azimuthal vortices that can appear only under the covered part of the surface.It has been shown that paramount role for the vortices formation plays the collision of radial thermocapillary flow with the edge of the film.The second important result of this work is in the coincidence of the vortices number at horizontal plane in the course of calculations and experiment for different intermediate values of the ring width.Satisfactory conformity in the case of films with average width is explained by the ability to get the reliable data for the area of surfactant film in experiment.In addition, it permits to be sure in uniformity of physical properties of the surface active matter film.Thus, the transitional stage of surface purification by radial thermocapillary flow is not critically important in the process of azimuthal vortices formation under the film near the cavity periphery.
There is an opinion that the origination of azimuthal vortices is connected with a weak mobility of surfactant film that is also observed in experiments.However, we think that it is wrong because the mathematical continuity of solution demands in limiting case of stationary film to impose the no-slip condition on solid surface.In experiment the film of insoluble surfactant is practically immovable.Therefore, the factor of mobility is also must be negligibly small in describable phenomenon.Explanation of this effect has some analogy with the problem of radial spreading of the liquid in divergent diffusor.During the radially symmetric spreading of incompressible fluid the reverse flow in a part of diffusor can appear in a cross section.In addition, the up layer is hotter than lower liquid therefore the buoyancy force supports divergent current and impedes the motion of liquid particle at once downward with following return along the bottom into the centre of the cavity.Thus, the inclination of the motion plane for a fluid element is energetically favorable and the velocity field obtains the azimuthal component.The analysis of trajectories of liquid particles shows that the azimuthal vortices are localized in the up part of the cavity under the film as in experiment.The calculations demonstrate the dependence of vortices plane inclination into the horizontal section on the heating intensity and radius of the free surface.Continuous transition of the numerical solution is received in both limiting cases for completely covered and open cavities.This result also corresponds to experiment.

Figure 2 .
Figure 2. Heating intensity I heat [Wt/m 2 ] as a function of the radial coordinate r [m] for different instants of time: lines 1-5 correspond to t = 0.5, 2, 3, 5 and 6 s.

Figure 3 .
Figure 3. Working mesh: general case of cylindrical cavity (a), particular case of segment (b).