Dynamics of the Two-Phase System in a Confined Plane Area Under Local Thermal Load Victoria B. Bekezhanova Olga N. Goncharova Institute of Computational Altai State University Modelling SB RAS Barnaul, Russia 656049 Krasnoyarsk, Russia 660036 gon@math.asu.ru vbek@icm.krasn.ru Abstract The mathematical model for simulating of the dynamics of bilayer sys- tem subjected to local thermal exposure is supposed. The model im- plies the occurrence of phase transition in the system due to evapora- tion. The diffusive type of evaporation is assumed. Boundary condi- tions imposed on the liquid – gas interface are generalized for the case of evaporative convection, using the thermodynamical relations. Nu- merical algorithm is developed to calculate basic characteristics of the liquid-gas system and the interface location at any time. Testing cal- culations are performed for one class of boundary conditions. 1 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 adjust- ment 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. Copyright © 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 2021 , Omsk, Russia, 24-04-2021, published at http://ceur-ws.org 1 2 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). Figure 1: Configuration of the two-phase system 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 + Gj , (1) ∂t Tj + ∂x Tj ∂y ψj − ∂y Tj ∂x ψj = Pr−1 −1   ∆ψj + ωj = 0, j Rej (∆Tj + δ∆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-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 + α∆T1 ). (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 , Tj present the stream function, vorticity and temperature in jth layer, G1 = (Gr1 /Re21 )∂x T1 + γ(Ga/Re21 )∂x C, G2 = (Gr2 /Re22 )∂x T2 . 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 (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 /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 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. 2 3 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: Tj x=0 = 0, Tj x=X = 0, T1 y=Y = 0, T2 y=0, x∈Q / s = 0, T2 y=0, x∈Qs = qls (t), (4) l l where Qsl is the area of the substrate where l-th heater with the temperature qls is arranged. On the part of external boundary confining the gas layer, the condition of zero vapor flux is imposed: ∂n C + α∂n T1 = 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 tempera- ture. In the dimensionless form, the relation for determining σ is given by σ(T ) = 1 − MaCa(T2 − T0 ), where Ma = (σT T∗ )/(ρ2 u∗ ν2 ) is the Marangoni number, Ca = ρ2 u∗ ν2 /σ0 is the capillary number. Here, σ0 , T0 are the reference values of the surface tension and temperature for the liquid, respectively, upon that, T0 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: ! ! 1 ∂x f ∂x f 1 s= p , p , n= −p , p . 1 + ∂x2 f 1 + ∂x2 f 1 + ∂x2 f 1 + ∂x2 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: 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 = vn n + vs s. The continuity of the temperature is postulated on the phase boundary: T1 = T2 = 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: p ∂t f + 1 + ∂x2 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. 3 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), ∂n ω2 − ρ̄ν̄∂n ω1 = F2 (t, x), (9) where F1 = Ma ∂s T + 2(1 − ρ̄ν̄) ∂s vn + vs R−1 ,  F2 (t, x) = −2 ∂s ∂n (v2 )n − ρ̄ν̄ ∂n (v1 )n + 2 1 − ρ̄ν̄ ∂s ∂x f vs R−1 + Ca−1 ∂s 1 − MaCaT R−1 −         − Gr2 Re−1 −1  p p  2 − ρ̄ν̄Gr1 Re1 T ∂x f / 1 + ∂x2 f + GaRe2 (1 − ρ̄)∂x f / 1 + ∂x2 f + Re2 (ρ̄ − 1)∂t vs + (ρ̄ − 1)vs ∂s vs + + (1 − ρ̄) ∂x f vn2 R−1 + vn (ω2 − ρ̄ω1 ) .  Here, R−1 = ∂xx f /(1 + ∂x2 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 T2 − k̄ ∂n T1 − δ k̄ ∂n C = −M, M = −E ∂n C + α∂n T1 , (10) where k̄ = k1 /k2 is the ratio of the thermal conductivity coefficients kj , M is the evaporative mass flow rate through the phase boundary, E = DLρ1 /(k2 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 + ε∗ (T1 − T0 )], ε∗ = T∗ Lµ0 /(R∗ T02 ). (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 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 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. 4 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 finite- difference 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 , Tj and position of the interface f (t, x). At t = 0 both fluids are at rest and have constant temperature (here, T0 ), 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 vn in each point of the interface is defined. 4 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 vs at Γt is calculated, it allows one to define 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 Ω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 numer- ically 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 uj , vj can be recalculated due to standard relations relating the physical variables and new required functions uj = ∂y ψj , vj = −∂x ψj , ωj = ∂x vj − ∂y uj . 9. Now, equation (8) is numerically solved and new position of the interface Γ is defined. Then, new values of vn 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. 5 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/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 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. 6 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). 5 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 References [Bek19] V. B. Bekezhanova, A. S. Ovcharova. Convection regimes induced by local boundary heating in a liquid – gas system. Journal of Fluid Mechanics, 873:441–458, 2019. [Bek20a] V. B. Bekezhanova, O. N. Goncharova. Impact of gravity on the flow pattern in a locally heated two-layer system. Microgravity Science and Technology, 32:229–243, 2020. [Bek20b] V. B. Bekezhanova, V. M. Fliagin, O. N. Goncharova, N. A. Ivanova, D. S. Klyuev. Thermocapillary deformations of a two-layer system of liquids under laser beam heating. International Journal of Multiphase Flow, 132:103429, 2020. [And12] V. K. Andreev, Yu. A. Gaponenko, O. N. Goncharova, V. V. Pukhnachov. Mathematical models of convection (de Gruyter Studies in Mathematical Physics). De Gruyter, Berlin, Boston, 2012. [Bek18] V.B. Bekezhanova, O.N. Goncharova. Problems of the Evaporative Convection (Review). Fluid Dyn., 53(Suppl.1):S69–S102, 2018. [Bek20c] V.B. Bekezhanova, O.N. Goncharova. Influence of the Dufour and Soret effects on the characteristics of evaporating liquid flows. Int. J. Heat Mass Transfer, 154:119696, 2020. 6