Introduction

Interest to the study of convection accompanied by the phase transition in fluidic mini-systems exposed to local thermal load is occasioned by wide application of fluidic technologies in engineering and smart-industry. Various setups and equipment using evaporating liquids as working media are used in thermophysics, microelectronics, material engineering, chemical and pharmaceutical industry, biomedicine. In the developmental stage of the technologies, the initial tune-up for whole system or certain working units can be required. For effective adjustment the knowledge of basic characteristics of the system in different conditions is indispensable. The necessary characteristics of the fluid system and parameters of phase transitions can be substantively obtained with the help of mathematical modeling. The theoretical (or numerical) results can allow one to determine the impact of particular factors on the system dynamics, to work out the control ways of appearing convective modes, as well as to specify methods to reduce undesirable effects caused by the external thermal action. In the present work the mathematical model is developed to investigate the features of convection in a multiphase system with the interface, to define the influence of liquid -gas phase transition on the phase boundary evolution and behavior of working fluids. The supposed model is a generalization of the model which was earlier developed by the authors and used for studying processes of heat-and-mass transfer in the locally heated multiphase systems without phase changes [Bek19,Bek20a,Bek20b]. Under statement of boundary conditions on the interface it is assumed that weak evaporation occurs.

Governing Equations

A bilayer system with sharp liquid -gas interface fills a plane horizontal cell being in the field of mass forces with the gravity vector g = (0, −g). External boundaries of the massif are the rigid impermeable walls. Here, heating from the substrate by thermal elements with finite size is carried out. The general scheme is sketched in Figure 1. The phase boundary is the thermocapillary surface which admits mass transfer of diffusive type due to evaporation. Working fluids (liquid and gas-vapor mixture) are considered as viscous heat-conducting media with constant thermophysical parameters. At initial time moment the system is in state of mechanical and local thermodynamical equilibrium (both fluids are quiescent, the liquid -gas surface is straight). The activating of heaters arranged on the lower wall gives rise to convective motion in both fluids, liquid evaporation and ensuing deformation of the interface. To describe the motion in each phase the Navier -Stokes equations in the Oberbeck -Boussinesq approximation are utilized. The governing equations are written in the terms of "stream function -vorticity" functions (ψ − ω). They have the following dimensionless form:

∂ t ω j + ∂ x ω j ∂ y ψ j − ∂ y ω j ∂ x ψ j = Re −1 j ∆ω j + G j , ∆ψ j + ω j = 0, ∂ t T j + ∂ x T j ∂ y ψ j − ∂ y T j ∂ x ψ j = Pr −1 j Re −1 j (∆T j + δ∆C).

(1)

The underlined term in the heat transfer equation (the last equation in (1)) is taken into account only for gas phase. Furthermore, in modelling dynamics of the gas-vapor layer, motion equations (1) is supplemented by the convection-diffusion equation describing the vapor transfer in the gas:

∂ t C + ∂ x C∂ y ψ 1 − ∂ y C∂ x ψ 1 = Pe −1 Re −1 1 (∆C + α∆T 1 ). (2)

Using equation (2), one imply that vapor is the passive admixture, which does not change any properties of the background gas. Note, that the influence of reciprocal effects of thermodiffusion and diffusive thermal conductivity occurring due to the presence of evaporated component in the gas is regarded. Here and below, the subscripts j = 1 and j = 2 corresponds to the upper (gas) and lower (liquid) layers respectively (see Figure 1); required functions ψ j , ω j , T j present the stream function, vorticity and temperature in jth layer,

G 1 = (Gr 1 /Re 2 1 )∂ x T 1 + γ(Ga/Re 2 1 )∂ x C, G 2 = (Gr 2 /Re 2 2 )∂ x T 2 . Parameters Gr j = β j T * gh 3 1 /ν 2 j (Grashof number), Re j = u * h 1 /ν j (Reynolds number), Pr j = ν j /χ j (Prandtl number), Ga = gh 3 1 /ν 2 1 (Galilei number), Pe = u * h 1 /D (diffusive Peclet

number) are the basic similarity criteria. They are assumed to be given and defined through the fluid physical parameters (coefficients of thermal expansion β j , kinematic viscosity ν j , thermal diffusivity χ j , diffusion D) as well as the characteristic scales of basic parameters (temperature drop T * , velocity of viscous stresses relaxation u * = ν 1 /h 1 in the gas phase, height of the upper layer h 1 in the unperturbed (initial) state). Parameters δ and α are the non-dimensional analogues of the Dufour and Soret coefficients respectively, γ is the coefficient of concentration expansion. Symbols ∂ t , ∂ x , ∂ y denote the partial-differential operators with respect to independent variables t, x, y correspondingly.

Boundary Conditions

On the outer solid boundaries of the cuvette (x = 0, x = X, y = 0, y = Y ) the no-slip conditions are fulfilled. In terms of stream functions ψ j the relations take form:

ψ j = 0, ∂ n ψ j = 0,(3)

where ∂ n denotes the normal derivative to corresponding boundary. For the temperature functions, the first kind conditions are set on all the external walls. Here, the presence of thermal elements on the substrate (y = 0) is taken into account:

T j x=0 = 0, T j x=X = 0, T 1 y=Y = 0, T 2 y=0, x / ∈Q s l = 0, T 2 y=0, x∈Q s l = q s l (t),(4)

where Q s l is the area of the substrate where l-th heater with the temperature q s l is arranged. On the part of external boundary confining the gas layer, the condition of zero vapor flux is imposed:

∂ n C + α∂ n T 1 = 0.

(5)

The common sharp interface Γ t dividing the liquid and gas phase is interpreted as the Gibbs surface which has zero thickness and mass. It is defined by the equation y = f (t, x). The tangential forces act lengthwise of the interface, at this, the surface tension of liquid σ is the function linearly depending on the temperature. In the dimensionless form, the relation for determining σ is given by σ(T ) = 1 − MaCa(T 2 − T 0 ), where Ma = (σ T T * )/(ρ 2 u * ν 2 ) is the Marangoni number, Ca = ρ 2 u * ν 2 /σ 0 is the capillary number. Here, σ 0 , T 0 are the reference values of the surface tension and temperature for the liquid, respectively, upon that, T 0 is the temperature of the local thermodynamical equilibrium; σ T is the temperature coefficient, σ 0 , σ T > 0. At the initial instant t = 0 both fluids are at rest, and interface Γ t has zero curvature. The actuation of heaters of finite size placed on the lower walls of the cuvette results in onset of convective motion in the layers and changes in the surface tension inducing the liquid motion along the interface and ensuing interface deformations.

Boundary conditions on Γ t are derived based on conservation laws for mass, momentum and energy and some additional assumptions [And12]. To formulate the conditions in ψ − ω variables, the unit vectors of tangent and normal lines to Γ t are introduced:

s = 1 1 + ∂ 2 x f , ∂ x f 1 + ∂ 2 x f , n = − ∂ x f 1 + ∂ 2 x f , 1 1 + ∂ 2 x f

. Such a form for n is obtained taking into account the orientation of the normal vector, n is the vector of external normal to Γ t for lower fluid (Figure 1). For all the points which lie on the interface, one can determine the normal and tangent components of the velocity: v n = −∂ s ψ, v s = ∂ n ψ; here, symbols ∂ s , ∂ n denote the derivatives in the tangential and normal directions. Therefore, the velocity of any point on the interface can be written as

v = v n n + v s s.

The continuity of the temperature is postulated on the phase boundary:

T 1 = T 2 = T,(6)

so that T is the common value of the temperature for both media on the interface. Analogous condition for the total vector velocity leads to the following relations for the stream functions on Γ t :

ψ 1 = ψ 2 , ∂ n ψ 2 − ∂ n ψ 1 = 0.(7)

In deriving these conditions, the volume conservation requirement for each of the fluids is additionally taken into account.

The kinematic conditions representing the material character of the phase boundary is written in the form:

∂ t f + 1 + ∂ 2 x f ∂ s ψ 2 = 0. (8)

The relationship is the equation for determining state of the interface at every instant. The use of such a form of the kinematic condition implies that weak evaporation occurs in the bilayer system under study, i. e. the convective mass flux through the interface is not taken into account.

Let us formulate conditions setting the force balance on the interface. In the scalar form the condition is written as complex of two relations which are the analogues of the tangential and normal components of the dynamic condition. In the ψ − ω terms they are transformed in the following form:

ω 2 − ρνω 1 = F 1 (t, x), ∂ n ω 2 − ρν∂ n ω 1 = F 2 (t, x),(9)

where

F 1 = Ma ∂ s T + 2(1 − ρν) ∂ s v n + v s R −1 , F 2 (t, x) = −2 ∂ s ∂ n (v 2 ) n − ρν ∂ n (v 1 ) n + 2 1 − ρν ∂ s ∂ x f v s R −1 + Ca −1 ∂ s 1 − MaCaT R −1 − − Gr 2 Re −1 2 − ρνGr 1 Re −1 1 T ∂ x f / 1 + ∂ 2 x f + GaRe 2 (1 − ρ)∂ x f / 1 + ∂ 2 x f + Re 2 (ρ − 1)∂ t v s + (ρ − 1)v s ∂ s v s + + (1 − ρ) ∂ x f v 2 n R −1 + v n (ω 2 − ρω 1 ) . Here, R −1 = ∂ xx f /(1 + ∂ 2

x f ) 3/2 , ρ = ρ 1 /ρ 2 and ν = ν 1 /ν 2 are the ratios of fluids densities and kinematic viscosities. It should be noted, that the conditions do not consider the dynamical action of evaporant on the liquid. Detailed derivation of relations ( 9) is presented in [Bek19].

The heat balance condition takes into account the diffusive mass transfer due to evaporation:

∂ n T 2 − k ∂ n T 1 − δ k ∂ n C = −M, M = −E ∂ n C + α∂ n T 1 ,(10)

where k = k 1 /k 2 is the ratio of the thermal conductivity coefficients k j , M is the evaporative mass flow rate through the phase boundary, E = DLρ 1 /(k 2 T * ) is the evaporation number, L is the latent heat of vaporization.

The expression for M function is obtained with the help of the mass balance equation on Γ t . Value M is the qualitative parameter allowing one to specify the regimes of phase transition. If M > 0 then liquid evaporation into gas occurs in the system, if M < 0 then vapor condensation takes place.

To close the problem, one need to state a boundary condition for the vapor concentration on the interface. One of the most frequently used approaches is to use the Clapeyron -Clausius and Mendeleev -Clapeyron equations [Bek18]. In the present work the linearized relation being the result of these two equations is used:

C = C * [1 + ε * (T 1 − T 0 )], ε * = T * Lµ 0 /(R * T 2 0 ).(11)

Here, µ 0 is the molar mass of the evaporating liquid, R is the universal gas constant, C * is the equilibrium concentration of saturated vapor (concentration at T 1 = T 0 ). It should be noted, that one can use the full analogue of the condition such a form. In both cases this condition presents one of conditions of the local thermodynamical equilibrium. Remark 1. The relevance of taking into account of the thermodiffusion effects in the governing equations and boundary conditions for similar problems is substantiated in [Bek20c].

Solving stated initial boundary value problem allows one to obtain all the basic characteristics of the multiphase system: values of the velocity and temperature for both fluids, vapor content in the gas layer, evaporative mass flow rate and profile of the interface at every instant.

General Scheme Of Numerical Method

To solve the problem the numerical algorithm supposed in [Bek20a] is modified for the case under study, when mass transfer due to evaporation and changes in vapor content in gas phase are considered. The outline of the algorithm are given below.

The computational scheme assumes a transformation of physical domains Ω j (see Figure 1) with curvilinear boundaries into the canonic computation regions with straight limiting lines. The algorithm uses the finitedifference scheme of stabilizing correction which is unconditionally stable one and has formally the second order of accuracy. As a final of approximation, the system of linear algebraic equations is resulted; its solution is stably calculated with the help of variants of the Gaussian elimination: sweep method (or the Thomas algorithm in the spatial variable directions) and original sweep method with parameters. 1. Let the system is some state characterized by known distributions ψ j , ω j , T j and position of the interface f (t, x). At t = 0 both fluids are at rest and have constant temperature (here, T 0 ), the surface Γ t is flat, f (t, x) = const. With known basic characteristics, equation ( 8) is numerically solved, thereby new position of Γ t is found, and the normal velocity v n in each point of the interface is defined.

2. At each time step, new spatial variables are introduced according to laws:

x = ξ, y = ζ(Y − f (ξ, t)) + f (ξ, t)

in Ω 1 and x = ξ, y = ζf (ξ, t) in Ω 2 .

3. Sweep procedure is used to compute the finite-difference motion, heat-and mass transfer equations in each layer.

4. The tangential velocity v s at Γ t is calculated, it allows one to define right-hand sides F 1 and F 2 in conditions (9).

5. The temperature functions T j are found numerically as solution of the corresponding boundary problem including the energy equations (the last equation in system (1)) and boundary conditions (4), ( 6), (10).

6. The vapor concentration function C is calculated in domain Ω 1 on the basis of equation ( 2) and boundary conditions (5), (11).

7. Equations of momentum transfer (the first equation in (1)) with boundary conditions (9) on Γ t are numerically solved with known functions T , C and f . Thus, ω j are defined. At this, on the fixed boundaries the Thom conditions which are resulted from (3) are utilized.

8. To calculate the unknown functions ψ j , in each time step we introduce the convergent iteration processes to solve problem including the Poisson equation (the second equation in (1)) and boundary conditions (3), ( 7), (8). The velocity vector components u j , v j can be recalculated due to standard relations relating the physical variables and new required functions

u j = ∂ y ψ j , v j = −∂ x ψ j , ω j = ∂ x v j − ∂ y u j .

9. Now, equation ( 8) is numerically solved and new position of the interface Γ is defined. Then, new values of v n are computed.

10. Transition to the step 2 is carried out.

This numerical algorithm was realized as author's FORTRAN code. The calculations were carried out on compute cluster in the Institute of Computational Modelling SB RAS.

Results Of Testing Calculations

In the framework of proposed statement the numerical simulation of the dynamics of two-phase benzine -air system exposed to local thermal load was performed with the help of the developed numerical algorithm. The system fills the vessel with length X = 0.2 m and height Y = 0.004 m. The terrestrial conditions with g = 9.81 m/s 2 are considered. In the unperturbed state the thickness of each layer h j was taken to be equal to 2 • 10 −3 m. The system characteristics were numerically investigated under different operating modes of the heaters of various size arranged on the substrate. The first thermal element has size 0.0275 m and heats with constant temperature. The second heater has length equal to 0.02 m and works in commutated mode, when its temperature can abruptly be varied. The character of changes in thermal and hydrodynamical fields as well as vapor concentration in the gas layer were investigated. The behavior of the liquid -gas interface and variations of the evaporative mass flow rate were analyzed. The visualization tools allowing one to keep track of all characteristics with time were developed. Working window of the programm is presented in Figure 2. Each window of the working display shows a certain parameter of the system in each time moment. Results demonstrate that the proposed approach gives reliable description of the changes in basic parameters of two-phase system under study subjected to local heating and can be applied to modelling real physical fluidic mini-systems.

Acknowledgements

The first author is grateful for financial supports to 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) (testing the numerical algorithm). The work of the second author was carried out in accordance with the State Assignment of the Russian Ministry of Science and Higher Education entitled "Modern methods of hydrodynamics for environmental management, industrial systems and polar mechanics" (Govt. contract code: FZMW-2020-0008) (the problem statement).

Figure 2: Visualization of the two-phase system parameters: the distribution of the vapor concentration is shown in the window 1 ; the temperature field near the interface and shape of Γ t (black solid curve) are demonstrated in the window 2 ; superposed fields of the temperature and velocity are given in the window 3 ; the tangential velocity and evaporative mass flow rate at all points of the interface are displayed in the windows 4 and 5 respectively