Introduction

The problems of convection under phase transitions is of interest due to wide application of fluid media in different technological equipments. Bright examples of fluidic technologies are liquid cooling in various systems of thermal control and thermal coating. The experimental testing of any set-up or technique is preceded by theoretical study of regularities of heat and mass transfer processes, behavior of working fluids under different conditions. The main problem under theoretical investigation of evaporative convection is the choice of a way of description of considered phenomena [Bek18a]. The most widely used mathematical model to research convection accompanied by phase changes is built upon the use of the Navier -Stokes equations or their approximations. To consider the temperature effects in whole two-phase system and concentration effects in a gas phase appearing due to the liquid vaporization, the equations of heat and mass transfer are invoked additionally to the motion equations. The resulting system of governing equation possesses natural properties of symmetry with respect to space and time, and as a consequence of this, it provides the implementation of the fundamental conservation laws. Except that the system of thermal concentration convection equations is nonlinear one and has high order, it is irrelevant to certain type. All of these heightens the value of exact solutions which take the basic symmetry properties peculiar to the original system of equations. It implies physical plausibility of results obtained with the help of the solutions.

There are a few available exact solution to describe convection in the two-phase systems [Bek18a]. Examples of meaningful solutions for studying steady flows with phase transitions are the Ostroumov -Birikh type solutions, which can be effectively applied to define basic characteristics of the convective regimes for various conditions. Two-dimensional solutions of such a class allowed one to specify the features of the flow regimes with evaporation/condensation depending on the different parameters of the problem [Bek16,Bek18b,Bek18c,She18,She19,Bek20a,Bek20b]. In all the cited works the theoretical results were obtained on the assumption on uniform character of evaporation. In the present work, the generalization of the Ostroumov -Birikh type solution is derived for the case of inhomogeneous evaporation when the evaporation rate is not constant and changes along the channel.

Governing Equations And Functional Form Of Partially Invariant Solution

The joint flow of evaporating liquid and co-current gas-vapor flux in a plane infinite horizontal channel is considered. The system is oriented so that the mass force vector g has the coordinates g = (0, −g) in the Cartesian coordinate system (x, y). The lower and upper boundaries of the channel y = −l and y = h are the solid impermeable walls. The fluids contact along the common thermocapillary surface Γ admitting phase transition due to evaporation/condensation. The interface remains to be flat, it is given by the equation y = 0.

The Navier -Stokes equations in the Oberbeck -Boussinesq approximation is utilized for description of the flow in both phases. In plane case the governing equations have the following form:

u ∂u ∂x + v ∂u ∂y = − 1 ρ ∂p ∂x + ν ∂ 2 u ∂x 2 + ∂ 2 u ∂y 2 ,(1)u ∂v ∂x + v ∂v ∂y = − 1 ρ ∂p ∂y + ν ∂ 2 v ∂x 2 + ∂ 2 v ∂y 2 + g (βT + γC) ,(2)∂u ∂x + ∂v ∂y = 0,(3)u ∂T ∂x + v ∂T ∂y = χ ∂ 2 T ∂x 2 + ∂ 2 T ∂y 2 + δ ∂ 2 C ∂x 2 + ∂ 2 C ∂y 2 . (4)

The terms containing function C are taken into account in modelling the flows in the gas phase. The vapor transfer in the carrier gas is described by the diffusion equation [Lan1987]:

u ∂C ∂x + v ∂C ∂y = D ∂ 2 C ∂x 2 + ∂ 2 C ∂y 2 + α ∂ 2 T ∂x 2 + ∂ 2 T ∂y 2 . (5)

In equations (1) -(5) the following notations are used: u, v are the longitudinal and transversal components of the velocity vector, correspondingly, p is the modified pressure function which defines the deviation of real pressure from the hydrostatic one, T is the fluid temperature, C is the vapor concentration in background gas, ρ is the fluid density (for each phase it is some reference value of the density), ν, χ, β are the coefficients of kinematic viscosity, heat diffusivity and thermal expansion, respectively, D is the coefficient of vapor diffusion in the gas, γ is the concentration coefficient of density, the parameters δ and α are the Dufour and Soret coefficients respectively, characterizing the effects of diffusive thermal conductivity and thermodiffusion. System (1) -( 5) admits an exact solution of a special type [Bek16]:

u i = u i (y), v i = 0, T j = T j (x, y) = a i 1 + a i 2 y x + ϑ j (y), C = C(x, y) = (b 1 + b 2 y) x + φ(y), p j = p j (x, y). (6)

Here, index i denotes characteristics of working fluid: i = 1 and i = 2 regard to the liquid in the lower layer and the gas-vapor mixture in the upper one, respectively. Explicit expressions for p j functions and additive terms ϑ j (y), φ(y) as well as exact values of the parameters a i 2 , b j (i = 1, 2; j = 1, 2) which are the constants, are determined with the help of the boundary conditions.

Boundary Conditions And Closing Relations

The boundary conditions are formulated taking into account the character of dependence of the unknown functions on spatial variables defined by (6). The no-slip conditions on the rigid walls are valid for the velocity functions:

u 1 (−l) = 0, u 2 (h) = 0,(7)

It is supposed that distributed thermal load is applied on the substrate

T 1 (x, −l) = A 1 x + ϑ − ,(8)

whereas the upper wall is thermally-insulated:

T 2y + δC y y=h = 0. (9)

The condition of zero concentration of vapor is prescribed on the upper wall:

C (x, h) = 0. (10)

The last condition defines property of the boundary to instantaneously absorb the vapor. In the experimental setup it is provided by a freezing out of the vapor on the wall. The applicability of the first kind boundary condition for the vapor concentration function is discussed in [Bek18d] in the framework of three-dimensional statement of the evaporative convection problem.

Setting the boundary conditions on the interface Γ, one should take into account that it is not deformed and given by the equation y = 0. The kinematic condition is fulfilled identically in view of the exact solution form. The continuity conditions for the longitudinal component of velocity vector and temperature at the interface are valid:

u 1 = u 2 , T 1 = T 2 . (11)

The continuity of the transversal component of the velocity vector follows from the kinematic condition. Projecting the dynamic condition on the unit tangential and normal vectors to the interface, there result two scalar relations:

ρ 1 ν 1 u 1y = ρ 2 ν 2 u 2y − σ T T x y=0 , p 1 = p 2 .(12)

Here, the impact of the surface forces is characterized by σ T > 0 which gives the temperature coefficient of the surface tension σ, σ = σ 0 − σ T (T − T 0 ); σ 0 > 0, T 0 are the reference values of the surface tension and the liquid temperature. The heat balance condition is written considering the thermodiffusion effects and mass balance equation:

κ 1 T 1y − κ 2 T 2y − δκ 2 C y y=0 = −λM, M = −Dρ 2 (C y + αT 2y | y=0 ) .(13)

Here, λ is the latent heat of evaporation, M defines the evaporation mass flow rate of the liquid at the interface.

In the present paper, the exact solution is derived under assumption, that the evaporation rate is not constant and changes along the channel according to the linear law:

M = M (x) = M 0 + M x x.(14)

The evaporation mass flow rate M is one of the important qualitative characteristics of the convective regimes with phase change. The positive (negative) values of M correspond to evaporation of the liquid into the gas flow (vapor condensation from the gas phase).

The condition for C function on the phase boundary sets the saturated vapor concentration and presents the linearized form of the equation which is the consequence of the Clapeyron -Clausius and Mendeleev -Clapeyron equations:

C y=0 = C * [1 + ε(T 2 y=0 − T 0 )].(15)

Here, ε = λµ 0 /(RT 2 0 ), µ 0 is the molar mass of the evaporating liquid, R is the universal gas constant, C * is the saturated vapor concentration at T 2 = T 0 .

To close the problem statement, the gas flow rate is preset:

Q = h 0 ρ 2 u 2 (y) dy.(16)

Thus, we have two sets of equations, each of them governs the heat and mass transfer in a certain fluid, and interface conditions which relate these equations. The conditions imposed on the liquid -gas boundary allow one to correctly describe the phase transition of diffusive type. The problem statement used is substantiated in [Bek16,Bek20b]. Solving the problem, one can define the basic characteristics of the bilayer system: the velocity u i , pressure p i , temperature T i in ith phase and vapor concentration C in gas layer.

Explicit Expressions For Required Functions

The fulfilment of condition of temperature continuity at the interface (the second equality in (11)) results in the following equality: a i 1 = A (i = 1, 2). The value determines the interface longitudinal temperature gradient characterizing the intensity thermal effects on the phase boundary, and as consequence, the intensity of evaporation and thermocapillary convection.

Substitution of (6) in equations allows one to derive the explicit expressions for all the required functions:

u i (y) = c i 3 + c i 2 y + c i 1 y 2 2 + L i 3 y 3 6 + L i 4 y 4 24 , p i (x, y) = d i 1 + d i 2 y + d i 3 y 2 2 x + c i 8 + K i 1 y + K i 2 y 2 2 + K i 3 y 3 3 + K i 4 y 4 4 + +K i 5 y 5 5 + K i 6 y 6 6 + K i 7 y 7 7 + K i 8 y 8 8 , T i (x, y) = A + a i 2 y x + c i 5 + c i 4 y + N i 2 y 2 2 + N i 3 y 3 6 + N i 4 y 4 24 + N i 5 y 5 120 + N i 6 y 6 720 + N i 7 y 7 1008 , C (x, y) = (b 1 + b 2 y) x + c 7 + c 6 y +

Coefficients c i j (i = 1, 2; j = 1, ..., 5; 8), c 6 , c 7 are the integration constants; the common scheme for finding of c i j will be outlined below. Coefficients L i 4 , L i 3 , S j , K i m (i = 1, 2; j = 2, ..., 7; m = 1, ..., 8) are calculated via physical parameters of the problem g, β i , ν i , χ i , ρ i , D, γ, coefficients defining the longitudinal gradients for the temperature and vapor concentration A, a i 2 , b i (i = 1, 2), and integration constants c i j .

Algorithm For Calculating The Integration Constants

Let the gas flow rate (16) and values of the longitudinal temperature gradients A 1 and ϑ − (see ( 9)) be preset.

The fulfilment of boundary conditions (7) -(15) and relationship (16) results to a system of equations for finding the integration constants c i j (i = 1, 2; j = 1, ..., 5; 8), c 6 , c 7 . Determining these constants, the velocity and temperature profiles, the pressure distributions for both phases and the vapor concentration in the gas layer are calculated based on formulas (17). Below, the algorithm for finding all the unknown parameters and constants is proposed. (v) Temperature distribution (8) determines the parameter a 1 2 = (A − A 1 )/l, and heat transfer condition (13) at the interface with the help of the expressions from (iii) leads to a relationship for the longitudinal temperature gradients A and A 1 :

A 1 = A 1 + (l/h)Dρ 2 (1 − αδ)λC * ε/κ 1 .

(vi) Two linear algebraic equations for the integration constants follow from the conditions defining the thermal boundary regime (8) and ( 13):

ϑ 1 (−l) = ϑ − , κ 1 c 1 4 − κ 2 c 2 4 − δκ 2 c 6 = −λM 0 . (vii)

Values of the unknown integration constants {c i 1 , c i 2 , c i 3 } (i = 1, 2) are calculated by solving the linear algebraic equation system that is a consequence of no-slip conditions (7), balance relations (12), condition of velocity continuity in (11) and equality (16) giving the gas flow rate.

(viii) Formulas determining the thermal and vapor concentration boundary regimes (9), (10), conditions at interface setting the temperature continuity (11) and saturated vapor concentration (15), and heat balance equation (13) lead to the system of the linear algebraic equations for calculation {c i 4 , c i 5 , c 6 , c 7 }, i = 1, 2. (ix) The value of M 0 will be computed with the help of obtained values c 2 4 and c 6 .

Realizing this algorithm, all the integration constants {c i 1 , c i 2 , c i 3 } and {c i 4 , c i 5 , c 6 , c 7 } are found. Hence, the mass flow rate of evaporation M at the interface can be determined in the form (14). It should be clear that a part of parameters contained in the exact solution is associated with the some compatibility conditions dictated by the boundary regimes.

Influence Of Gas Flow Rate On Characteristics Of Two-Phase Flows

The liquid -gas two-layer system is considered. The HFE-7100 fluid (HydroFluoroEther) and nitrogen are chosen as working liquid and gas respectively. The HFE-7100 liquid volatilizes from the lower layer, so that a mixture of nitrogen and HFE-7100-liquid vapor is in the upper layer. The physical parameters of the working media are given below in the order {HFE-7100, nitrogen} (or only HFE-7100): ρ = {1.5 • 10 3 , 1.2} kg/m 3 ; ν = {0.38 • 10 −6 , 0.15

• 10 −4 } m 2 /s; β = {1.8 • 10 −3 , 3.67 • 10 −3 } K −1 ; χ = {0.4 • 10 −7 , 0.3 • 10 −4 } m 2 /s; κ = {0.07, 0.02717} W/(m•K); σ T = 1.14 • 10 −4 N/(m•K); D = 0.7 • 10 −5 m 2 /s; L = 1.11 • 10 5 W•s/kg; C * = 0.45 (correspond to equilibrium temperature T 0 = 293.15 K); µ 0 = 0.25 kg/mol; γ = −0.5; ε * = 0.04 K −1 ; δ = 10 −5 K; α = 5 • 10 −3 K −1 .

The influence of the gas flow rate Q on the flow patterns which arise in the bilayer system is investigated. The typical distributions of the basic characteristics of the convective regimes with evaporation obtained on the basis of the exact solutions under study are presented in Figs. 1, 2. The parameters defining the external thermal action applied on the substrate are the longitudinal temperature gradient A 1 and value ϑ − (see condition (8)). In the calculation, the values of A 1 equal to −10 or 10 K/m and ϑ − = 293.15 K are set. The negative (positive) values of A 1 correspond to cooling (heating) of the lower boundary of the flow domain in the direction of the longitudinal axes. Values of the gas flow rate Q are assumed to be 0.1R 0 , R 0 and 10R 0 , where R 0 = 9.6 • 10 −6 kg/(m 2 •s). The following parameters are fixed for all the cases under consideration: thicknesses of the liquid and gas layers l = 2 mm and h = 6 mm, respectively, and gravity acceleration g = 9.81 m 2 /s.

Changes in gas flow rate can lead to significant alteration of the flow topology. If the lower wall cools along the longitudinal axes x with intensity A 1 = −10 K/m, then negative thermal gradient A is formed on the phase boundary (see corresponding value of A in Table 1). It induces the thermocapillary effect action causing the liquid motion from hot domain in zone with lower temperature along the interface (see velocity profiles shown in Figs. 1(a,d,g)). In this case, action of the Marangoni effect is co-directional to shear stresses generated by gas flux. At small and moderate Q, the peculiar layering of the velocity profile occurs, which is appeared by the formation of a wedge near the interface (Fig. 1(a,d)). One can tell about dominant impact of the thermocapillary effect in these cases, since the maximum velocity is achieved in the liquid layer. Further increase in gas flow rate leads to change of dominant mechanism, the velocity field is governed by large shear stresses. Upon that, the maximum velocity is observed in the gas layer (Fig. 1(g)).

If the substrate heats with temperature gradient A 1 = 10 K/m, then positive interface gradient A is set at the interface (see Table 2). Due to the Marangoni effect the reverse flow appears in the liquid layer (Fig. 2), and this effect is dominant mechanism defining the flow pattern under small gas flow rate Q = 0.1R 0 . In this case the maximum values of velocity is achieved in the liquid phase (Fig. 2(a)). At moderate values of Q, there occurs the competition of the influence of the dynamical effect arising due to the shear stresses and thermocapillary one. The intensity of motion in the gas and liquid layers is almost the same (Fig. 2(d)). At large Q the competition is maintained. The global counterflow in the liquid persists, although significant increase of the velocity in gas phase is observed (Fig. 2(g)).

For both types of external thermal load, the variations of the gas flow rate do not give rise to alteration of the fields of temperature and vapor concentration as well as to significant quantitative changes in thermal and concentration characteristics. Some quantitative parameters of the convective regimes are presented in Tables 1, 2. According to point (v) of the aforesaid algorithm (see Section 5), the interface temperature gradient A is determined by the values of A 1 for certain pair of the working media and system geometry. One can see, that the variations of Q does not result to the change in the interface thermal gradient. In the tables, the following characteristics are given:

T defines the temperature drop in the bilayer system, T max and |u| max are the maximum values of the temperature and velocity modulus in whole system respectively, C max is the maximum value of the vapor concentration in the gas phase, M specifies the percentage change of the evaporative mass flow rate on the working segment.

Table 1: Parameters of the HFE-7100 -nitrogen system in the channel with the thermally-insulated upper wall with changes in Q at fixed One can see from the data given in the tables, the impact of the gas flow rate on the evaporation rate as well as on the vapor content in the gas layer is rather weak. Two order of magnitude increase in Q does practically not change of evaporation rate and its variations along the working segment being characterized by M . It is provided by thermal regime on the upper boundary.

A 1 = −10 K/m Q, A, T , T max , C max |u| max ×, M 0 • 10 4 , M x • 10 4 , M , kg/(m 2 •s) K/m K K ×10 −3 ,Q at fixed A 1 = 10 K/m Q, A, T , T max , C max |u| max ×, M 0 • 10 4 , M x • 10 4 , M , kg/(m 2 •s) K/m K K ×10 −3 , m/s kg/(m 2 •s) kg/(m•

Summary

In the present work an exact solution of the problem of evaporative convection problem is constructed. The physical interpretation of the solution is the following: it allows one to describe a flow of the evaporating liquid driven by the gas flux on the working segment subjected by distributed thermal load from the substrate in a quite long minichannel with thermally-insulated upper boundary. The new solution takes into account the inhomogeneous character of phase transition. The explicit expressions for all the unknown functions are derived.

The impact of the flow rate in the gas layer is studied at an example of the bilayer system of certain working fluids. It is established that changes in gas flow rate slightly influence the thermal characteristics of the system and the parameters of vapor content in the considered conditions (i. e. for chosen the geometrical configuration, intensity of external thermal load and gravity action) despite the significant changes in the flow topology. One can assume, that due to thermal isolation of the upper boundary one can control the evaporation rate and vapor concentration in the carrier gas in quite large range of variations of the gas flow rate. industrial systems and polar mechanics" (Govt. contract code: FZMW-2020-0008) (the problem statement; construction of the presented type of the Ostroumov -Birikh solution; idea to apply this solution for description of the evaporative convection; obtaining the system of the linear algebraic equations for the integration constants). The second author is grateful for financial supports the Krasnoyarsk Mathematical Center financed by the Ministry of Science and Higher Education of the Russian Federation in the framework of the establishment and development of regional Centers for Mathematics Research and Education (Agreement No. 075-02-2021-1384) (verification of the constructed exact solution; explanation of the group-analytical properties; visualization of the Ostroumov -Birikh solution; explanation of the physical mechanisms of the convective flows described by presented problem statement).