The mathematical model for simulating of the dynamics of bilayer system subjected to local thermal exposure is supposed. The model implies the occurrence of phase transition in the system due to evaporation. The di usive type of evaporation is assumed. Boundary conditions imposed on the liquid { gas interface are generalized for the case of evaporative convection, using the thermodynamical relations. Numerical algorithm is developed to calculate basic characteristics of the liquid-gas system and the interface location at any time. Testing calculations are performed for one class of boundary conditions.
Copyright © by the paper's authors. Use permitted under Creative Commons License Attribution 4.0 International (CC BY 4.0).
A bilayer system with sharp liquid { gas interface lls a plane horizontal cell being in the eld 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 nite size is carried out. The general scheme is sketched in Figure 1. The phase boundary is the thermocapillary surface which admits mass transfer of di usive type due to evaporation. Working uids (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 uids are quiescent, the liquid { gas surface is straight).
The activating of heaters arranged on the lower wall gives rise to convective motion in both uids, 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:
!j + Gj ; j + !j = 0;
C): 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-di usion equation describing the vapor transfer in the gas:
T1): Using equation (2), one imply that vapor is the passive admixture, which does not change any properties of the background gas. Note, that the in uence of reciprocal e ects of thermodi usion and di usive 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 , Tj present the stream function, vorticity and temperature in jth layer, G1 = (Gr1=Re12)@xT1 + (Ga=Re12)@xC, G2 = (Gr2=Re22)@xT2. Parameters Grj = j T gh31= j2 (Grashof number), Rej = u h1= j (Reynolds number), Prj = j = j (Prandtl number), Ga = gh31= 12 (Galilei number), Pe = u h1=D (di usive Peclet number) are the basic similarity criteria. They are assumed to be given and de ned through the uid physical parameters (coe cients of thermal expansion j , kinematic viscosity j , thermal di usivity j , di usion D) as well as the characteristic scales of basic parameters (temperature drop T , velocity of viscous stresses relaxation u = 1=h1 in the gas phase, height of the upper layer h1 in the unperturbed (initial) state). Parameters and are the non-dimensional analogues of the Dufour and Soret coe cients respectively, is the coe cient of concentration expansion. Symbols @t; @x; @y denote the partial-di erential operators with respect to independent variables t; x; y correspondingly. (1) (2)
On the outer solid boundaries of the cuvette (x = 0, x = X, y = 0, y = Y ) the no-slip conditions are ful lled. In terms of stream functions j the relations take form: where @n denotes the normal derivative to corresponding boundary.
For the temperature functions, the rst kind conditions are set on all the external walls. Here, the presence of thermal elements on the substrate (y = 0) is taken into account:
Tj x=0 = 0;
Tj x=X = 0;
T1 y=Y = 0;
T2 y=0; x2=Qls = 0;
T2 y=0; x2Qls = qls(t); where Qls is the area of the substrate where l-th heater with the temperature qls is arranged.
On the part of external boundary con ning the gas layer, the condition of zero vapor ux is imposed: s = ! ; n = ! : 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 uid (Figure 1). For all the points which lie on the interface, one can determine the normal and tangent components of the velocity: vn = @s , vs = @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 = vnn + vss.
The continuity of the temperature is postulated on the phase boundary: (3) (4) (5) (6) (7) (8) 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: In deriving these conditions, the volume conservation requirement for each of the uids is additionally taken into account.
The kinematic conditions representing the material character of the phase boundary is written in the form:
T1 = T2 = T; 1 = 2;
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 = F1(t; x); (9) where
F2(t; x) = + 2 1
MaCaT R 1 + (1 !1) : Here, R 1 = @xxf =(1 + @x2f )3=2, = 1= 2 and = 1= 2 are the ratios of uids 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 di usive mass transfer due to evaporation:
M;
M = C = C [1 + " (T1
T0)]; " = T L 0=(R T02): 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 T1 = T0). 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 thermodi usion e ects 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 uids, vapor content in the gas layer, evaporative mass ow rate and pro le of the interface at every instant. 4
To solve the problem the numerical algorithm supposed in [Bek20a] is modi ed 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 nitedi erence scheme of stabilizing correction which is unconditionally stable one and has formally the second order of accuracy. As a nal 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 , Tj and position of the interface f (t; x). At t = 0 both uids are at rest and have constant temperature (here, T0), the surface t is at, f (t; x) = const. With known basic characteristics, equation (8) is numerically solved, thereby new position of t is found, and the normal velocity vn in each point of the interface is de ned. (10) (11) 2. At each time step, new spatial variables are introduced according to laws: x = , y = (Y in 1 and x = , y = f ( ; t) in 2. 3. Sweep procedure is used to compute the nite-di erence motion, heat- and mass transfer equations in each layer. 4. The tangential velocity vs at t is calculated, it allows one to de ne right-hand sides F1 and F2 in conditions (9). 5. The temperature functions Tj 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 conditions (5), (11). 1 on the basis of equation (2) and boundary 7. Equations of momentum transfer (the rst equation in (1)) with boundary conditions (9) on t are numerically solved with known functions T , C and f . Thus, !j are de ned. At this, on the xed 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 uj , vj can be recalculated due to standard relations relating the physical variables and new required functions uj = @y j , vj = @x j , !j = @xvj @yuj . 9. Now, equation (8) is numerically solved and new position of the interface vn are computed. is de ned. Then, new values of 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. 5
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 lls the vessel with length X = 0:2 m and height Y = 0:004 m. The terrestrial conditions with g = 9:81 m/s2 are considered. In the unperturbed state the thickness of each layer hj was taken to be equal to 2 10 3 m. The system characteristics were numerically investigated under di erent operating modes of the heaters of various size arranged on the substrate. The rst 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 elds as well as vapor concentration in the gas layer were investigated. The behavior of the liquid { gas interface and variations of the evaporative mass ow 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 uidic mini-systems. 6
The rst author is grateful for nancial supports to the Krasnoyarsk Mathematical Center nanced 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- 20211384) (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).