A Multiplicity Of The Evaporative Convection Regimes Modeled On The Basis Of Exact Solutions Olga N. Goncharova Victoria B. Bekezhanova Altai State University Institute of Computational Barnaul, Russia 656049 Modelling SB RAS gon@math.asu.ru Krasnoyarsk, Russia 660036 vbek@icm.krasn.ru Abstract Problems of heat and mass transfer in liquid media with evapora- tion/condensation at interfaces require the development of mathemat- ical models and numerical methods, the study of solvability ques- tions and construction of exact solutions. The paper provides a short overview of results for description of the convective flows with phase transition obtained on the basis of modelling with the help of new exact solutions of the classical convection equations. Classification of three- dimensional two-layer flows is carried out; flow control methods at dif- ferent temperature and vapor concentration conditions at the channel boundaries are investigated. 1 Introduction The problems of heat and mass transfer in the liquid media with evaporation/condensation at the interface remain one of the most difficult problems in the convection theory. They require development or clarification of mathematical models, the study of solvability of the initial-boundary value problems, construction of the exact solutions of governing equations or elaboration of the effective numerical algorithms. The importance of theoretical and experimental study of the convective flows with interfaces, including phase transition boundaries, is explained by the need to optimize and improve applied developments in the field of fluidic technologies, such as liquid cooling, liquid-layer information recording systems and techniques, which are associated with the operation of heat exchangers, evaporators, heat pipes. Physical experiments performed in the frame of international projects at the universities of Belgium, France and Russia (at the Institute of Thermophysics SB RAS) allow us to study the various characteristics of the flows of evaporating liquids and co-current gas fluxes under conditions of gravitational fields of different intensities [Cel08, Ior11]. The structure of such flows is determined by various factors: natural and thermocapillary convection, additional shear stresses induced by a gas flow, mass transfer through the interface due to evaporation or condensation etc. 2 Governing Equations Based on the Navier – Stokes equations in the Oberbek – Boussinesq approximation, a mathematical model has been developed that takes into account the effects of thermal diffusion and diffusion thermal conductivity. Copyright ⃝ c by the paper’s authors. Use permitted under Creative Commons License Attribution 4.0 International (CC BY 4.0). In: Sergei S. Goncharov, Yuri G. Evtushenko (eds.): Proceedings of the Workshop on Applied Mathematics and Fundamental Computer Science 2020 , Omsk, Russia, 30-04-2020, published at http://ceur-ws.org 1 Figure 1: Flow Geometry An evaporative liquid and a gas-vapor mixture fill an infinite channel with interface (see Fig. 1). The bound- aries of the domains Ω1 , Ω2 Ω1 = {(x, y, z) : −x0 < x < 0, 0 < y < 1, −∞ < z < ∞}, Ω2 = {(x, y, z) : 0 < x < x0 , 0 < y < 1, −∞ < z < ∞} are the impermeable immovable rigid walls defined by the equations x = −x0 , x = x0 , y = 0, y = 1. The thermocapillary forces act along the interface Γ defined by the equation x = 0. The surface tension σ linearly depends on the temperature T : σ = σ0 − σT (T − T0 ); σ0 , T0 are the reference values of the surface tension and liquid temperature, respectively, σT > 0 is the temperature coefficient of the surface tension. The Cartesian coordinate system (x, y, z) is chosen in such a way that the gravitational vector is directed opposite to the Ox axis (g = (−g, 0, 0)). The width of the channel h is the characteristic linear size; u∗ is the characteristic velocity, T∗ is the characteristic temperature drop, p∗ = ρ1 u2∗ is the characteristic pressure. The velocity vectors vi = (ui , vi , wi ), the pressure deviation from the hydrostatic one pi , temperature Ti and vapor concentration C functions satisfy the convection equations written in the non-dimensional form as follows: ηiv (vi · ∇)vi = −ηip ∇pi + ∆vi + Gi , div vi = 0, (1) Re ηiT vi · ∇Ti = (∆Ti + δ∆C), (2) Re Pr 1 v2 · ∇C = (∆C + α∆T2 ). (3) Pe Last equation (3) and the underlined term in (2) are used to model the gas-vapor motion in Ω2 . The parameters and functions indicated with i = 1 and i = 2 are related to the media, filling the layers Ω1 ((liquid) and Ω2 (gas-vapor), respectively. The following notations are used: G1 = −i(Gr/Re2 )T1 , G2 = −i β(Gr/Re2 )T2 + ) γ(Ga/Re2 )Cs , i is the unit vector of the Ox axis. The arising dimensionless parameters are determined here: η1p = 1, η2p = 1/ρ, η1v = 1, η2v = ν, η1T = 1, η2T = χ; Re = u∗ h/ν1 is the Reynolds number, Pr = ν1 /χ1 is the Prandtl number, Gr = β1 T∗ gh3 /ν12 is the Grashof number, Ga = gh3 /ν12 is the Galilei number, Pe = u∗ h/D is the diffusion Peclet number; g = |g|; ρ = ρ2 /ρ1 , ν = ν2 /ν1 , χ = χ2 /χ1 , β = β2 /β1 are the ratios of the densities, coefficients of the kinematic viscosity, thermal diffusivity, and thermal expansion, respectively; γ is the concentration coefficient of the gas density, D is the coefficient of vapor diffusion; α = α∗ T∗ , δ = δ∗ /T∗ ; the dimensional coefficients α∗ and δ∗ characterize the Soret and Dufour effects in the gas phase. 2.1 Exact Solutions Of The Convection Equations In the system with evaporation/condenstation at the thermocapillary interface the phase transition leads to formation of a temperature gradient, which arises due to a decrease in the average kinetic energy of the liquid volume. For the first time, an exact solution describing convective flows in the presence of an arbitrarily oriented temperature gradient was presented by Ostroumov [Ost52] and later re-obtained by Birikh for an infinite layer with a free boundary [Bir66]. Upon that, a temperature gradient linear in longitudinal coordinate is applied at the boundaries of the flow domain. Solutions of the convection equations having a special form, which allows us to talk about the presence of a longitudinal temperature gradient, are commonly called the Ostroumov – Birikh type solutions. 2 Here, the exact solution of the governing equations (1) – (3) is considered, which is a generalization of the Ostroumov – Birikh solution in the three-dimensional case. It describes three-dimensional flows in an infinite channel under the influence of applied longitudinal (horizontal) temperature gradient, and is interpreted as a solution describing a two-layer flow under the conditions of the “liquid – gas” phase transition on the working section of a sufficiently long vessel. It is worth pointing out that the steady motion is characterized by a resulting temperature gradient, in general case it is non-uniform in the vertical direction and formed due to action of combination of different factors, namely, thermocapillary and thermodiffusion effects, phase transition etc. The constructed solution of the convection equations (1) – (3) has the following form: ui = ui (x, y), vi = vi (x, y), wi = wi (x, y), (4) Gr Gr Ga p1 = −A xz + q1 (x, y), p2 = −Aρβ 2 xz + Bργ xz + q2 (x, y), (5) Re2 Re Re2 Ti = −Az + Θi (x, y), (6) C = Bz + Φ(x, y). (7) It is assumed that the interface Γ between the layers Ω1 and Ω2 is determined by the equation x = 0. The coefficients A and B (see (6), (7)) specify the constant longitudinal gradients of temperature and concentration along the interface. If A∗ and B∗ are dimensional gradients, then the non-dimensional ones are equal to A = A∗ h/T∗ , B = B∗ h, respectively. Due to the form of the solution solution (4) – (7) a splitting of the original 3D problem statement into a set of the two-dimensional ones is possible (for details of the reduction we refer to [Bek18a]). As a result we model the stationary 3D convective two-layer flows with evaporation at the interface when solving the chain of 2D problems. It should be noted that solution (4) – (7) is considered to be “exact” one in the extended (group) sense [Bek18a]. 2.2 Boundary Conditions Solution (4) – (7) satisfies the equations (1) – (3) and the following boundary conditions [Bek18a, Bek18b]. On the solid channel walls the no-slip conditions for the velocity vectors and the conditions of thermal insulation of the lateral walls are set. The case of the absence of the vapor flux on the upper and lateral rigid boundaries is considered. At the thermocapillary interface Γ the kinematic and dynamic conditions are fulfilled. The heat transfer condition and the mass balance equation are taken into consideration. The concentration of saturated vapor is set at the phase boundary. This condition is the linearized form of the relation derived as a consequence of the Clapeyron – Clausius equation and Mendeleev – Clapeyron equation for ideal gas, and it is used for the vapor concentration function at the interface. Additionally the fulfillment of the continuity conditions of the velocity and temperature is required at the interface. Complete statement of the problem under study is presented and substantiated in [Bek18a]. Numerical investigations of the set of 2D problems for computation of the input functions ui (x, y), vi (x, y), wi (x, y) and Θi (x, y), Φ(x, y) are carried out with the help of own code developed on the basis of a finite- difference scheme of second approximation order and being unconditionally stable, known as the method of alternative direction. Detailed description of the numerical algorithm is presented in [Bek20]. 3 Results And Discussion Mathematical modeling of two-layer flows on the basis of an exact solution makes it possible to quickly and efficiently study the influence of the thermophysical properties of the working media, the gravitational field, the thermal load at the boundaries, and the geometric parameters of the system on the characteristics of the flows and the parameters of the arising regimes. The constructed solution allows us to describe the formation of lengthwise thermocapillary structures observed in thermophysical experiments and ensembles of convective cells that are elements of the spatially periodic pattern of flows [Kab11, Bek16]. The spatial dimensions and shape of the cells can vary depending on the geometry of the system (thickness of the layers of the working media), gravity, the type of liquid coolant and the intensity of external thermal impact. The solution predicts the modes in which thermal shafts, thermal rolls with defects (structures like “thermal horns” or “thermal torches”) occur. In all the case, the temperature field completely determines the characteristics of the vapor content in the gas. 3 The effect of thermal diffusion effects in a gas-vapor mixture on the parameters of the mass evaporation rate is revealed. We present results of numerical investigations, which describe topology, thermal characteristics and vapor distributions for the liquid-gas systems like ethanol – nitrogen, HFE-7100 – nitrogen and FC-72 – nitrogen char- acterized by certain physico-chemical properties [Bek18a, Bek20]. The values of the dimensionless parameters are computed, when the characteristic temperature drop and characteristic length are equal to T∗ = 10 K and h = 10−2 m. The characteristic velocity is chosen equal to the velocity of viscous stresses relaxation ν/h so that u∗ = {0.15 · 10−5 , 0.38 · 10−4 , 0.47 · 10−4 } m/s for (1) ethanol, (2) HFE-7100 and (3) FC-72. The thermal load intensity is determined by the coefficient of the A characterizing the longitudinal temperature gradient; A has the values A = 0.01, A = 0.1 and A = 0.3 corresponding to the dimensional ones equal to A∗ = 10 K/m, A∗ = 100 K/m and A∗ = 300 K/m, respectively. The calculations were carried out for fluid layers with the thickness of the liquid layer hl = 0.1 or hl = 0.25 and thickness of the gas-vapor layer hg = 0.5 for all the considered cases. The appropriate dimensional values of the layer thicknesses are h∗l = 0.1 or h∗l = 0.25 cm and h∗g = 0.5 cm. Values of the non-dimensional complexes (see notations to equations (1) – (3)) are presented in Table 1. Some important parameters are also specified in Table 1 because they are included to the boundary conditions. They are the Marangoni number M a = σT T∗ h/(ρ1 ν1 χ1 ), capillary number Ca = ρ1 ν1 u∗ /σ0 , dimensionless latent heat evaporation parameter L = λDρ2 /(κ1 T∗ ), λ is the latent heat of vaporization. The influence of intensity of the gravitational field and interface thermal regime under change of geometry of the flow domain, especially of thickness of the liquid layer, has been studied. 3.1 Formation Of Thermal “Horns” And Appearance Of Hot Interfacial Film With change in heat load determined by the value of the longitudinal temperature gradient A the alteration of the temperature field and topological flow structure as a whole is observed. In the case, when A∗ = 300 K/m, the formation of thermal “horns” near the lateral walls can occurs in the ethanol – nitrogen system (see figure 2, 3). These “horns” appear due to the near-wall condensation which is typical for this system. Reconstruction of the temperature field can also be caused by change in gravity intensity determined by g value. Figures 2, 3 present the difference in the ethanol – nitrogen flow structure appearing under weak and terrestrial gravity. Here, the case when the liquid layer thickness is half the gas layer thickness is considered. For sufficiently large values of A∗ the spatial structure of the hydrodynamic field is becoming more complicated. The figures show trajectories of the fluid particles and projections of the stream tubes in the cross sections of the channel. Figure 2: Flow Topology (Streamlines and Trajectories) and Thermal Field in Ethanol – Nitrogen System; A∗ = 300 K/m; hl = 2.5 mm; hg = 5 mm; g = g0 · 10−2 m/s2 In the case, when the thickness of the liquid layer is equal to the thickness of the gas layer, in [Bek18a] it was elucidated that, while maintaining structures such as thermal horns in the near-wall regions, the hot shaft in the central part of the channel observed under microgravity is replaced by the cold one in terrestrial conditions. Such thermal field alteration is explained by the fact, that the thermocapillary effect is suppressed by gravity in the normal gravitational field, that prevents the transfer of liquid from the hot region to the cold one. Note, that with a decrease in the thickness of the liquid layer, the thermal picture in the ethanol – nitrogen system becomes completely different: a hot film appears on the surface of the liquid; and the thermal shaft (and, consequently concentration ones) split into two smaller rolls, see figure 4. 4 Figure 3: Flow Topology (Streamlines and Trajectories) and Thermal Field in Ethanol – Nitrogen System; A∗ = 300 K/m; hl = 2.5 mm; hg = 5 mm; g = g0 = 9.81 m/s2 Figure 4: Flow Topology (Streamlines and Trajectories) and Thermal Field in Ethanol – Nitrogen System; A∗ = 300 K/m; hl = 1 mm; hg = 5 mm; g = g0 = 9.81 m/s2 3.2 Formation Of Thermal And Concentration Shafts Modelling two-layers convective flows on the basis of exact solutions, we can predict the formation of the lon- gitudinal structures observed during physical experiments. The constructed solution allows us to describe the appearance of both thermocapillary longitudinal rolls and thermal (and concentration) shafts. The figures 5, 6 show the flow topology and temperature distribution for the HFE-7100 – nitrogen and ethanol – nitrogen systems, respectively, in the case when the liquid layer thickness is 2.5 mm and g value correspond to the normal gravity; here, the value of the longitudinal temperature gradient A∗ is equal to 100 K/m. With a change in gravity, the intensity of fluid flows changes (the values of the velocity modulus were compared). Spatial size of thermal and concentration shafts also changes. It is seen, that depending on the liquid coolant type, the topological flow pattern changes (figures 5, 6; left). A six-vortex flow with characteristic small angular vortices forms in the HFE-7100 liquid (figure 5; left). The size of the heat shaft also varies significantly compared to the ethanol liquid (figures 5, 6; right). These various flow characteristics is due to differences in the viscous and transport properties, including the heat conductivity and thermal diffusivity of liquids, as well as to significant difference in the latent heat of vaporization, which is significantly lower in the hydrofluoroethers than in ethanol. Finally, we demonstrate differences in the flow structure of the presented coolants from the liquid FC-72, see figures 7, 8. Here, the conditions of normal gravity are fixed, the same values of the liquid layer thickness (hl = 2.5 mm) and longitudinal temperature gradient (A∗ = 300 K/m and A∗ = 100 K/m) are chosen as for ethanol – nitrogen and HFE-7100 – nitrogen systems, respectively. Differences in the flow topology and tempera- ture distribution in the channel are clearly manifested. For FC-72 – nitrogen system it is characteristic formation of zones with counter motion near the bottom wall and change in shape of vortices (figures 7, 8, left). As for 5 Figure 5: Flow Topology (Streamlines and Trajectories) and Thermal Field in HFE-7100 – Nitrogen System; A∗ = 100 K/m; hl = 2.5 mm; hg = 5 mm; g = g0 = 9.81 m/s2 Figure 6: Flow Topology (Streamlines and Trajectories) and Thermal Field in Ethanol – Nitrogen System; A∗ = 100 K/m; hl = 2.5 mm; hg = 5 mm; g = g0 = 9.81 m/s2 the pattern of the temperature field, the occurrence of longitudinal thermal shaft on the interface (as for HFE- 7100 – nitrogen system) and cold thermocline within liquid layer (figures 7, 8, right) is predicted by the exact solution under study. In the case of thin liquid layers, qualitative differences have been established in the flow topology and in for- mation of thermal shafts or thermoclines, both while maintaining the thermophysical properties of the working liquid and when they change. A different number of vortex structures (stream tubes) can be observed under nor- mal and reduced gravity. Thus, parameters and characteristics of evaporative convection regimes are determined by a host of factors including both properties of working media and external actions. 4 Acknowledgements The work was financially supported by the Russian Foundation for Basic Research, the Government of the Kras- noyarsk region, and the Krasnoyarsk Region Science Foundation in the frame of the scientific project “Theoretical and experimental investigation of heat and mass transfer processes in the two-phase system of thermal control” (grant No. 18-41-242005). References [Bek18a] V.B. Bekezhanova, O.N. Goncharova. Modeling of three dimensional thermocapillary flows with evap- oration at the interface based on the solutions of a special type of the convection equations. Appl. Math. Model., 62:145–162, 2018. 6 Figure 7: Flow Topology (Streamlines and Trajectories) and Thermal Field in FC-72 – Nitrogen System; A∗ = 300 K/m; hl = 2.5 mm; hg = 5 mm; g = g0 = 9.81 m/s2 Figure 8: Flow Topology (Streamlines and Trajectories) and Thermal Field in FC-72 – Nitrogen System; A∗ = 100 K/m; hl = 2.5 mm; hg = 5 mm; g = g0 = 9.81 m/s2 [Bek20] V.B. Bekezhanova, O.N. Goncharova. Numerical study of the evaporative convection regimes in a three-dimensional channel for different types of liquid-phase coolant. Int. J. Thermal Sci., 156:106491, 2020. [Bek18b] V.B. Bekezhanova, O.N. Goncharova. Problems of the Evaporative Convection (Review). Fluid Dyn., 53(1):S69–S102, 2018. [Bek16] V.B. Bekezhanova, O.A. Kabov. Influence of internal energy variations of the interface on the stability of film flow. Interfacial Phenomena and Heat Transfer, 4(2-3):133-156, 2016. [Bir66] R.V. Birikh. Thermocapillary convection in a horizontal layer of liquid. J. Appl. Mech. Tech. Phys., 3:43–45, 1966. [Cel08] G.P. Celata, C. Colin, P. Colinet, P. Di Marco, T. Gambaryan-Roisman, O. Kabov. Bubbles, Drops, Films: Transferring Heat in Space. Europhysics News, 39(4):23-25, 2008. [Ior11] C.S. Iorio, O.N. Goncharova, O.A. Kabov. Heat and mass transfer control by evaporative thermal pattering of thin liquid layers. Computational Thermal Sci., 3(4):333-342, 2011. [Kab11] O.A. Kabov, D.V. Zaitsev, V.V. Cheverda, A. Bar-Cohen. Evaporation and flow dynamics of thin, shear-driven liquid films in microgap channels. Exp. Therm. Fluid Sci., 35(5):825-831, 2011. [Ost52] G.A. Ostroumov. Free convection under the conditions of an internal problem. Gostekhizdat Press, Moscow – Leningrad, 1952 [in Russian]. 7 Table 1. Values of dimensionless parameters Dimensionless parameter \ Ethanol HFE-7100 FC-72 Medium Density ratio ρ 0.15 · 10−2 0.8 · 10−3 0.7 · 10−3 Kinematic viscosity ratio ν 10 39.5 32 Thermal expansion ratio β 3.4 2 2.41 Thermal conductivity ratio k 0.16 0.39 0.48 Thermal diffusivity ratio χ 337 750 920 Reynolds number Re 1 1 1 Prandtl number Pr 17 9.5 14.31 Pecklet number Pe 0.1 0.54 · 10−1 0.67 · 10−1 Grashof number Gr at g0 4.7 · 104 1.22 · 106 6.84 · 105 and g0 · 10−2 4.7 · 102 1.22 · 104 6.84 · 103 Galilei number Ga at g0 4.36 · 106 6.8 · 107 4.5 · 107 and g0 · 10−2 4.36 · 104 6.8 · 105 4.5 · 105 Density concentration coefficient γ −0.62 −0.5 −0.5 Marangoni number Ma 7.6 · 104 5 · 105 1.9 · 105 Capillary number Ca 0.81 · 10−5 0.16 · 10−5 0.31 · 10−5 Latent heat evaporation parameter L 8.28 1.33 1.41 Saturated vapor concentration parameters C∗ 0.1 0.45 0.6 and ϵC 0.1 0.04 0.04 Dufour coefficient δ 10−5 10−5 10−5 Soret effect α 10−5 10−3 10−3 8