Mathematical and Information Technologies, MIT-2016 β€” Mathematical modeling On the Inverse Problem of a Creeping Motion in Thin Layers Victor Andreev Institute Computational Modelling SB RAS, Siberian Federal University, Krasnoyarsk, Russia andr@icm.krasn.ru Abstract. The new partially invariant solution of two-dimential mo- tions of heated viscous liquid equations is considered. For factor-system arised the initial boundary value problem is formulated. This problem is inverse one and describing of common motion of two immiscible liquids in a plane channel under the action of thermocapillary forces. As Marangoni number is small (so-called creeping flow) the problem becomes the lin- ear one. Some a priori estimates are obtained and input data conditions when solution tends to stationary one are found. In Laplace transforms the exact solution is obtained as quadratures and some numerical results of velocities behavior in layers are presented. Keywords: Thermocapillarity, a priori estimates, conjugate initial-boundary value problem, asymptotic behaviour, numerical simulation 1 Introduction It is well known that in a non-uniformly heated liquid a motion can arise. In some applications of liquid flows, a joint motion of two or more fluids with surfaces takes place. If the liquids are not soluble in each other, they form a more or less visual interfaces. The petroleum-water system is a typical example of this situation. At the present time modelling of multiphase flows taking into account different physical and chemical factors is needed for designing of cooling systems and power plants, in biomedicine, for studying the growth of crystals and films, in aerospace industry [1-4]. Nowadays, there are exact solutions of the Marangoni convection [5-7]. One of the first solutions was obtained in [8]. This is the Poiseuille stationary flow of two immiscible liquids in an inclined channel. As a rule, all such flows were considered steady and unidirectional. The stability of such flows was investigated in [9, 10]. As for non-stationary thermocapillary flows, studying of them began recently [11, 12]. Thermocapillary convection problem for two incompressible liquids separated by a closed interface in a container was investigated in [13]. Local (in time) unique solvability of the problem was obtained in Holder classes of functions. The problem of thermalcapillary 3D motion of a drop was studied in [14]. More- over, its unique solvability in Holder spaces with a power-like weight at infinity 258 Mathematical and Information Technologies, MIT-2016 β€” Mathematical modeling was established. Velocity vector field decreases at infinity in the same way as the initial data and mass forces, the temperature diverges to the constant which is the limit of the initial temperature at infinity. The present work is devoted to studying of solutions of a conjugate boundary value problem arising as a re- sult of linearization of the Navier-Stokes system supplemented with temperature equation. The description of the 2D creeping joint motion of two viscous heat conducting fluids in flat layers is also provided here. The motion arises due to thermocapillary forces imposed along two interfaces, after which the unsteady Marangoni convection begins. Such kind of convection can dominate in flows under microgravity conditions or in motions of thin liquid films. 2 Statement of the Problem The 2D motion of a viscid incompressible heat-conducting liquid in the absence of mass forces is described by the system of equations 1 𝑒1𝑑 + 𝑒1 𝑒1π‘₯ + 𝑒2 𝑒1𝑦 + 𝑝π‘₯ = 𝜈(𝑒1π‘₯π‘₯ + 𝑒1𝑦𝑦 ), (2.1) 𝜌 1 𝑒2𝑑 + 𝑒1 𝑒2π‘₯ + 𝑒2 𝑒2𝑦 + 𝑝𝑦 = 𝜈(𝑒2π‘₯π‘₯ + 𝑒2𝑦𝑦 ), (2.2) 𝜌 𝑒1π‘₯ + 𝑒2𝑦 = 0, (2.3) πœƒπ‘‘ + 𝑒1 πœƒπ‘₯ + 𝑒2 πœƒπ‘¦ = πœ’(πœƒπ‘₯π‘₯ + πœƒπ‘¦π‘¦ ). (2.4) Here 𝑒1 (π‘₯, 𝑦, 𝑑) and 𝑒2 (π‘₯, 𝑦, 𝑑) are the components of the velocity vector, 𝑝(π‘₯, 𝑦, 𝑑) is the pressure, πœƒ(π‘₯, 𝑦, 𝑑) is the temperature, 𝜌 > 0 is the density, 𝜈 > 0 is the kinematic viscosity and πœ’ > 0 is the thermal conductivity of the liquid. The quantities 𝜌 > 0, 𝜈 > 0 and πœ’ > 0 are constant. The system of equation (2.1)–(2.4) admits a four-dimential Lie subalgebra 𝐺4 = βŸ¨πœ•π‘₯ , πœ•π‘’1 + π‘‘πœ•π‘₯ , πœ•π‘ , πœ•πœƒ ⟩. Its invariants are 𝑑, 𝑦, 𝑒2 and a partially invariant solution of rank 2 and defect 3 should be sought for in the form 𝑒1 = 𝑒1 (π‘₯, 𝑦, 𝑑), 𝑒2 = 𝑣(𝑦, 𝑑), 𝑝 = 𝑝(π‘₯, 𝑦, 𝑑), πœƒ = πœƒ(π‘₯, 𝑦, 𝑑). Inserting the exact form of the solution into the equations (2.1)–(2.3) yields 𝑒1 = 𝑀(𝑦, 𝑑)π‘₯ + 𝑔(𝑦, 𝑑), 𝑀 + 𝑣𝑦 = 0, 1 𝑓 (𝑑)π‘₯2 𝑀𝑑 + 𝑣𝑀𝑦 + 𝑀2 = 𝑓 (𝑑) + πœˆπ‘€π‘¦π‘¦ , 𝑝 = 𝑑(𝑦, 𝑑) βˆ’ , (2.5) 𝜌 2 𝑑𝑦 = πœˆπ‘£π‘¦π‘¦ βˆ’ 𝑣𝑑 βˆ’ 𝑣𝑣𝑦 , 𝑔𝑑 + 𝑣𝑔𝑦 + 𝑀𝑔 = 0 with some function 𝑓 (𝑑) that is arbitrary so far. Regarding the temperature field, we assume that equation (2.4) has the so- lution of the form πœƒ = π‘Ž(𝑦, 𝑑)π‘₯2 + π‘š(𝑦, 𝑑)π‘₯ + 𝑏(𝑦, 𝑑). (2.6) 259 Mathematical and Information Technologies, MIT-2016 β€” Mathematical modeling As we see below, (2.6) is in good accord with conditions on the interface. The stationary solution of the Navier-Stokes equations in the form (2.5) for 𝑔 = 0 for pure viscous fluid was found for the first time by [15]. It describes the liquid impingement from infinity on the plane 𝑦 = 0 under the no slip condition on it. In the paper [16], this solution for the flow between two plates or for the flow in a cylindrical tube (axisymmetric analogue of solution (2.5)) was applied. It is known that the temperature dependence of the surface tension coefficient is the one of the most important factors leading to the dynamic variety of the interfacial surface. In the papers [17, 18] the stationary solutions in form (2.5), (2.6) was found at π‘Ž(𝑦, 𝑑) ≑ 0, 𝑏 = const for a flat layer with a free boundary 𝑦 = 𝑙 = const and a solid wall 𝑦 = 0. The non-uniqueness of solution depending on the physical parameters of the problem was revealed. A similar problem in the case of half space was investigated in [19]. We assume for simplicity that 𝑔(𝑦, 𝑑) ≑ 0, π‘š(𝑦, 𝑑) ≑ 0. The latter condition means that the temperature field has an extremum at π‘₯ = 0, more exactly, a maximum for π‘Ž(𝑦, 𝑑) < 0 and a minimum for π‘Ž(𝑦, 𝑑) > 0. Let us apply the solution of the form (2.5), (2.6) to described joint motion of two immiscible liquids in the flat layer 0 < 𝑦 < β„Ž considering that the wall 𝑦 = 0 and 𝑦 = β„Ž are solid and the line 𝛀 : 𝑦 = 𝑙(π‘₯, 𝑑) is their common interface, see Fig. 1. Fig. 1. Geometry of the Marangoni convection problem Introduction the index 𝑗 = 1, 2 for the liquids and using (2.5) and (2.6), we come to the conclusion that the unknowns satisfy the equations 𝑀𝑗𝑑 + 𝑣𝑗 𝑀𝑗𝑦 + 𝑀𝑗2 = πœˆπ‘— 𝑀𝑗𝑦𝑦 + 𝑓𝑗 (𝑑), 𝑀𝑗 + 𝑣𝑗𝑦 = 0, (2.7) 1 𝑓𝑗 (𝑑)π‘₯2 𝑝𝑗 = 𝑑𝑗 (𝑦, 𝑑) βˆ’ , 𝑑𝑗𝑦 = πœˆπ‘— 𝑣𝑗𝑦𝑦 βˆ’ 𝑣𝑗𝑑 βˆ’ 𝑣𝑗 𝑣𝑗𝑦 , (2.8) πœŒπ‘— 2 π‘Žπ‘—π‘‘ + 2𝑀𝑗 π‘Žπ‘— + 𝑣𝑗 π‘Žπ‘—π‘¦ = πœ’π‘— π‘Žπ‘—π‘¦π‘¦ , 𝑏𝑗𝑑 + 𝑣𝑗 𝑏𝑗𝑦 = πœ’π‘— 𝑏𝑗𝑦𝑦 + 2πœ’π‘— π‘Žπ‘— (2.9) in domain 0 < 𝑦 < 𝑙(π‘₯, 𝑑) for 𝑗 = 1 and in domain 𝑙(π‘₯, 𝑑) < 𝑦 < β„Ž for 𝑗 = 2. 260 Mathematical and Information Technologies, MIT-2016 β€” Mathematical modeling At the interface 𝑦 = 𝑙(π‘₯, 𝑑) the conditions hold [1] 𝑀1 (𝑙(π‘₯, 𝑑), 𝑑) = 𝑀2 (𝑙(π‘₯, 𝑑), 𝑑), 𝑣1 (𝑙(π‘₯, 𝑑), 𝑑) = 𝑣2 (𝑙(π‘₯, 𝑑), 𝑑), (2.10) 𝑙𝑑 + π‘₯𝑀1 (𝑙(π‘₯, 𝑑), 𝑑)𝑙π‘₯ = 𝑣1 (𝑙(π‘₯, 𝑑), 𝑑), (2.11) πœ•π‘Ž1 πœ•π‘Ž2 π‘Ž1 (𝑙(π‘₯, 𝑑), 𝑑) = π‘Ž2 (𝑙(π‘₯, 𝑑), 𝑑), π‘˜1 = π‘˜2 , πœ•π‘› πœ•π‘› πœ•π‘1 πœ•π‘2 𝑏1 (𝑙(π‘₯, 𝑑), 𝑑) = 𝑏2 (𝑙(π‘₯, 𝑑), 𝑑), π‘˜1 = π‘˜2 , πœ•π‘› πœ•π‘› π‘˜1 > 0, π‘˜2 > 0 are the heat conductivity coefficients and n = (1 + 𝑙π‘₯2 )βˆ’1/2 (βˆ’π‘™π‘₯ , 1) is the normal to the line 𝑦 = 𝑙(π‘₯, 𝑑). The dynamic condition for 𝛀 has a vector form [1] (𝑝1 βˆ’ 𝑝2 )n + 2[πœ‡2 𝐷(u2 ) βˆ’ πœ‡1 𝐷(u1 )]n = 2𝜎𝐾n + βˆ‡π›€ 𝜎, πœ‡π‘— = πœŒπ‘— πœˆπ‘— . (2.12) In (2.12) 𝐷(u) is the strain-rate tensor, 𝜎(πœƒ1 ) is the surface tension coefficient, 𝐾 is the mean curvature of the interface, whereas βˆ‡π›€ = βˆ‡ βˆ’ n(n Β· βˆ‡) on the right-hand side designates the surface gradient. For most of real liquid media the dependence 𝜎(πœƒ1 ) is approximated well by the linear function 𝜎(πœƒ1 ) = 𝜎 0 βˆ’ πœ…πœƒ1 , (2.13) where 𝜎 0 > 0 and πœ… > 0. They are assumed constant and determent by experimental methods. Projecting condition (2.12) to the tangent direction 𝜏 = (1 + 𝑙π‘₯2 )βˆ’1/2 (1, 𝑙π‘₯ ), and using (2.13), (2.6) we obtain π‘₯ 𝑙π‘₯ [πœ‡2 (𝑣2𝑦 βˆ’ 𝑀2 ) βˆ’ πœ‡1 (𝑣1𝑦 βˆ’ 𝑀1 )] + (1 βˆ’ 𝑙π‘₯2 )(πœ‡2 𝑀2𝑦 βˆ’ πœ‡1 𝑀1𝑦 ) 2 = βˆ’πœ…(πœƒ1π‘₯ + 𝑙π‘₯ πœƒ1𝑦 ) = βˆ’πœ…[2π‘Ž1 π‘₯ + 𝑙π‘₯ (π‘Ž1𝑦 π‘₯2 + 𝑏1𝑦 )]. (2.14) The projection (2.12) to the normal n yields [𝜌2 𝑓2 (𝑑) βˆ’ 𝜌1 𝑓1 (𝑑)]π‘₯2 𝜌1 𝑑1 βˆ’ 𝜌2 𝑑2 + + 2[πœ‡2 𝐷(u2 ) βˆ’ πœ‡1 𝐷(u1 )]n Β· n 2 𝑙π‘₯π‘₯ = [𝜎 0 βˆ’ πœ…(π‘Ž1 π‘₯2 + 𝑏1 )] . (2.15) (1 + 𝑙π‘₯2 )3/2 The boundary conditions on the solid walls have the form 𝑀1 (0, 𝑑) = 0, 𝑣1 (0, 𝑑) = 0, 𝑀2 (β„Ž, 𝑑) = 0, 𝑣2 (β„Ž, 𝑑) = 0, (2.16) π‘Ž1 (0, 𝑑) = π‘Ž10 (𝑑), π‘Ž2 (β„Ž, 𝑑) = π‘Ž20 (𝑑), (2.17) 𝑏1 (0, 𝑑) = 𝑏10 (𝑑), 𝑏2 (β„Ž, 𝑑) = 𝑏20 (𝑑), (2.18) with some given functions π‘Žπ‘—0 (𝑑) and 𝑏𝑗0 (𝑑). The initial conditions for the velocities are zero because of we study the properties of the solution of the problem simulating the motion only under the 261 Mathematical and Information Technologies, MIT-2016 β€” Mathematical modeling action of thermocapillary forces 𝑀𝑗 (𝑦, 0) = 0, 𝑣𝑗 (𝑦, 0) = 0. Besides, 𝑙(π‘₯, 0) = 𝑙0 (π‘₯), π‘Žπ‘— (𝑦, 0) = π‘Ž0𝑗 (𝑦), 𝑏𝑗 (𝑦, 0) = 𝑏0𝑗 (𝑦). Note several specific features of the formulated problem. This is a nonlinear and inverse one since the functions 𝑓𝑗 (𝑑) are unknowns also. It is easy to under- stand if we exclude the functions 𝑣𝑗 (𝑦, 𝑑) from the equations of mass conservation. Then the problem reduces to the conjugate problem for the functions 𝑀𝑗 (𝑦, 𝑑), π‘Žπ‘— (𝑦, 𝑑) and 𝑙(π‘₯, 𝑑). The problem for 𝑏𝑗 (𝑦, 𝑑) given 𝑣𝑗 (𝑦, 𝑑) and π‘Žπ‘— (𝑦, 𝑑) can be sep- arated. The functions 𝑑𝑗 (𝑦, 𝑑) can be recovered by quadrature from the second equation (2.8) up to a function of time. The last condition in (2.10) and the fourth from (2.16) are the additional conditions on 𝑓𝑗 (𝑑), 𝑗 = 1, 2. Let us introduce the characteristic scales of length and time as well as func- tions 𝑀𝑗 , 𝑣𝑗 , π‘Žπ‘— , 𝑑𝑗 and 𝑓𝑗 , namely, the quantities 𝑙0 , 𝑙02 /𝜈1 , πœ…π‘Ž0 𝑙0 /πœ‡1 , πœ…π‘Ž0 𝑙02 /πœ‡1 , π‘Ž0 , πœ…π‘Ž0 𝑙0 /𝜌1 , πœ…π‘Ž0 /(𝜌1 𝑙0 ), where 𝑙0 = const > 0 is the average value of thickness of the first layer of the liquid at 𝑑 = 0, π‘Ž0 = max |π‘Ž20 (𝑑) βˆ’ π‘Ž10 (𝑑)| > 0, or 𝑑>0 π‘Ž0 = max max |π‘Žπ‘—0 (𝑦)| > 0, if π‘Ž20 (𝑑) = π‘Ž10 (𝑑). In the dimensionless variables, 𝑗 𝑦 some factor appears at the nonlinear terms in (2.7), the Marangoni number M = πœ…π‘Ž0 𝑙03 /(πœ‡1 𝜈1 ). (2.19) The same applies to the kinematic condition (2.11) ¯𝑙𝑑¯ + π‘₯ Β― ¯𝑙(Β― Β―M𝑀( π‘₯, 𝑑¯), 𝑑¯)¯𝑙π‘₯Β― = MΒ― 𝑣1 (¯𝑙(Β― π‘₯, 𝑑¯), 𝑑¯). (2.20) Assume that the M β‰ͺ 1. The latter holds either in the thin layers or large viscosities. Then the nonlinear terms in the equations can be neglected and the latter become linear. In particular, the kinematic condition (2.20) has the form ¯𝑙𝑑¯ = 0, i.e. ¯𝑙 = ¯𝑙(π‘₯). Let us turn to (2.15). After transition to the dimensionless variables on the right-hand side the Weber number We = 𝜎 0 /(πœ…π‘Ž0 𝑙02 ) appears instead of 𝜎 0 . In the real conditions We ≫ 1 for the most of liquid media; for example, for the water–air system We ∼ 106 . Therefore, for these Weber numbers, (2.14) assume the form ¯𝑙π‘₯Β―π‘₯Β― = 0, i.e. ¯𝑙 = 𝛼π‘₯ + 𝑙0 . We assume later that 𝛼 = 0 and the interface is the plane 𝑦 = 𝑙0 < β„Ž parallel to the solid walls 𝑦 = 0 and 𝑦 = β„Ž; in what follows, the index 0 for 𝑙0 will be omitted. 3 A priori Estimates Let us present the so-obtained linear problem in its entirely in dimensional form 𝑀𝑗𝑑 = πœˆπ‘— 𝑀𝑗𝑦𝑦 + 𝑓𝑗 (𝑑), (3.1) 𝑀𝑗 (𝑦, 0) = 0, (3.2) 𝑀1 (0, 𝑑) = 0, 𝑀2 (β„Ž, 𝑑) = 0, (3.3) 𝑀1 (𝑙, 𝑑) = 𝑀2 (𝑙, 𝑑), (3.4) 262 Mathematical and Information Technologies, MIT-2016 β€” Mathematical modeling πœ‡2 𝑀2𝑦 (𝑙, 𝑑) βˆ’ πœ‡1 𝑀1𝑦 (𝑙, 𝑑) = βˆ’2πœ…π‘Ž1 (𝑙, 𝑑), (3.5) ∫︁ 𝑙 βˆ«οΈβ„Ž 𝑀1 (𝑧, 𝑑) 𝑑𝑧 = 0, 𝑀2 (𝑧, 𝑑) 𝑑𝑧 = 0, (3.6) 0 𝑙 where 0 < 𝑦 < 𝑙 for 𝑗 = 1 and 𝑙 < 𝑦 < β„Ž for 𝑗 = 2. The first equality in (3.6) follows from (2.10) whereas the last in the no-slip condition 𝑣2 (β„Ž, 𝑑) = 0. Let us write the problem for the functions π‘Žπ‘— (𝑦, 𝑑) π‘Žπ‘—π‘‘ = πœ’π‘— π‘Žπ‘—π‘¦π‘¦ , (3.7) π‘Žπ‘— (𝑦, 0) = π‘Ž0𝑗 (𝑦), (3.8) π‘Ž1 (0, 𝑑) = π‘Ž10 (𝑑), π‘Ž2 (β„Ž, 𝑑) = π‘Ž20 (𝑑), (3.9) π‘Ž1 (𝑙, 𝑑) = π‘Ž2 (𝑙, 𝑑), π‘˜1 π‘Ž1𝑦 (𝑙, 𝑑) = π‘˜2 π‘Ž2𝑦 (𝑙, 𝑑). (3.10) In order to obtain a priori estimates for 𝑀𝑗 (𝑦, 𝑑), 𝑓𝑗 (𝑑) of the solution of (3.1)–(3.5), it is necessary firstly to infer the estimates for the solutions of initial- boundary value problem (3.7)–(3.10). We perform the change of variables π‘Ž10 (𝑑)(𝑦 βˆ’ 𝑙)2 π‘Ž1 (𝑦, 𝑑) = π‘Ž Β―1 (𝑦, 𝑑) + , 0 6 𝑦 6 𝑙0 ≑ 𝑙, 𝑙2 (3.11) π‘Ž20 (𝑑)(𝑦 βˆ’ 𝑙)2 π‘Ž2 (𝑦, 𝑑) = π‘Ž Β―2 (𝑦, 𝑑) + , 𝑙 6 𝑦 6 β„Ž. (β„Ž βˆ’ 𝑙)2 The functions π‘Ž ¯𝑗 (𝑦, 𝑑) in their domains satisfy the equations 2πœ’1 π‘Ž10 (𝑑) π‘Žβ€²10 (𝑑)(𝑦 βˆ’ 𝑙)2 π‘Ž Β―1𝑑 = πœ’1 π‘Ž Β―1𝑦𝑦 + βˆ’ ≑ πœ’1 π‘Ž Β―1𝑦𝑦 + 𝑔1 (𝑦, 𝑑), (3.12) 𝑙2 𝑙2 2πœ’2 π‘Ž20 (𝑑) π‘Žβ€²20 (𝑑)(𝑦 βˆ’ 𝑙)2 π‘Ž Β―2𝑑 = πœ’2 π‘Ž Β―2𝑦𝑦 + βˆ’ ≑ πœ’2 π‘Ž Β―2𝑦𝑦 + 𝑔2 (𝑦, 𝑑), (3.13) (β„Ž βˆ’ 𝑙)2 (β„Ž βˆ’ 𝑙)2 where the prime denotes differentiation with respect to time. Boundary condi- tions (3.9) for π‘Ž Β―1 and π‘Ž Β―2 become homogeneous, whereas (3.10) preserve it form. Initial conditions (3.8) for π‘Ž Β―1 and π‘Ž Β―2 change π‘Ž10 (0)(𝑦 βˆ’ 𝑙)2 Β―1 (𝑦, 0) = π‘Ž01 (𝑦) βˆ’ π‘Ž Β―01 (𝑦), β‰‘π‘Ž 𝑙2 (3.14) 2 π‘Ž20 (0)(𝑦 βˆ’ 𝑙) Β―2 (𝑦, 0) = π‘Ž02 (𝑦) βˆ’ π‘Ž 2 Β―02 (𝑦). β‰‘π‘Ž 𝑙 Let us multiply (3.1), (3.2) by 𝜌1 𝑐1 π‘ŽΒ―1 and 𝜌2 𝑐2 π‘Ž Β―2 𝑐1 , 𝑐2 and integrate over the segments [0, 𝑙], [𝑙, β„Ž] taking into account (3.8) and (3.9). Then add up the result. We infer that ∫︁ 𝑙 βˆ«οΈβ„Ž ∫︁ 𝑙 βˆ«οΈβ„Ž 𝑑𝐴(𝑑) + π‘˜1 Β―21𝑦 𝑑𝑦 + π‘˜2 π‘Ž Β―22𝑦 𝑑𝑦 = 𝜌1 𝑐1 π‘Ž 𝑔1 π‘Ž Β―1 𝑑𝑦 + 𝜌2 𝑐2 𝑔2 π‘Ž Β―2 𝑑𝑦, (3.15) 𝑑𝑑 0 𝑙 0 𝑙 263 Mathematical and Information Technologies, MIT-2016 β€” Mathematical modeling ∫︁ 𝑙 βˆ«οΈβ„Ž 𝜌 1 𝑐1 𝜌 2 𝑐2 𝐴(𝑑) = Β―21 𝑑𝑦 + π‘Ž Β―22 𝑑𝑦, π‘Ž (3.16) 2 2 0 𝑙 where 𝑐𝑗 are the coefficients of the specific heat capacity. Along with (3.15) there is another identity ∫︁ 𝑙 βˆ«οΈβ„Ž [οΈ‚ ∫︁ 𝑙 βˆ«οΈβ„Ž ]οΈ‚ 1 πœ• 𝜌1 𝑐1 Β―21𝑑 𝑑𝑦 + 𝜌2 𝑐2 π‘Ž Β―22𝑑 𝑑𝑦 + π‘Ž π‘˜1 π‘Ž 2 Β―1𝑦 𝑑𝑦 + π‘˜2 π‘Ž 2 Β―2𝑦 𝑑𝑦 2 πœ•π‘‘ 0 𝑙 0 𝑙 ∫︁ 𝑙 βˆ«οΈβ„Ž = 𝜌1 𝑐1 𝑔1 π‘Ž Β―1𝑑 𝑑𝑦 + 𝜌2 𝑐2 𝑔2 π‘Ž Β―2𝑑 𝑑𝑦. (3.17) 0 𝑙 From (3.15) and (3.17) we obtain the inform estimates in 𝑦 (οΈ‚ )οΈ‚1/4 8πœ’π‘— |π‘Žπ‘— (𝑦, 𝑑)| 6 𝐹 (𝑑)𝐴(𝑑) + |π‘Žπ‘—0 (𝑑)|, (3.18) π‘˜π‘—2 where ∫︁ 𝑙 βˆ«οΈβ„Ž [οΈ‚ βˆ«οΈπ‘‘ βˆ«οΈπ‘‘ ]οΈ‚ 2π‘˜1 4πœ’1 𝑙 𝐹 (𝑑) = π‘˜1 Β―210 (𝑦) 𝑑𝑦+π‘˜2 π‘Ž Β―220 (𝑦) 𝑑𝑦+ π‘Ž π‘Ž210 (𝜏 ) π‘‘πœ + (π‘Žβ€²10 (𝜏 ))2 π‘‘πœ πœ’1 𝑙3 5 0 𝑙 0 0 [οΈ‚ βˆ«οΈπ‘‘ βˆ«οΈπ‘‘ ]οΈ‚ 2π‘˜2 4πœ’2 2 β„Žβˆ’π‘™ β€² 2 + π‘Ž20 (𝜏 ) π‘‘πœ + (π‘Ž20 (𝜏 )) π‘‘πœ ≑ 𝐹 (𝑑), (3.19) πœ’2 (β„Ž βˆ’ 𝑙)3 5 0 0 [οΈ‚ βˆšοΈƒ (οΈ‚ βˆ«οΈπ‘‘ βˆšοΈ‚ βˆ«οΈπ‘‘ )οΈ‚ βˆšοΈ€ π‘˜1 2πœ’ 𝑙 𝐴(𝑑) 6 π‘’βˆ’2𝛿𝑑 𝐴(0) + √1 π›Ώπœ 𝑒 |π‘Ž10 (𝜏 )| π‘‘πœ + 𝑒 π›Ώπœ |π‘Žβ€²10 (𝜏 )| π‘‘πœ πœ’1 𝑙3 5 0 0 βˆšοΈƒ (οΈ‚ βˆ«οΈπ‘‘ βˆšοΈ‚ βˆ«οΈπ‘‘ )οΈ‚]οΈ‚2 π‘˜2 2πœ’2 β„Žβˆ’π‘™ + βˆšοΈ€ π›Ώπœ 𝑒 |π‘Ž20 (𝜏 )| π‘‘πœ + 𝑒 π›Ώπœ |π‘Žβ€²20 (𝜏 )| π‘‘πœ . (3.20) πœ’2 (β„Ž βˆ’ 𝑙)3 5 0 0 As to functions 𝑀𝑗 (𝑦, 𝑑), 𝑓𝑗 (𝑑) the following estimates hold [οΈ‚ (οΈ‚ )οΈ‚]οΈ‚1/4 𝐸(𝑑) 4πœ…2 π‘™π‘Ž21 (𝑙, 𝑑) |𝑀1 (𝑦, 𝑑)| 6 2 𝐹 (𝑑) + , (3.21) 𝜈1 5πœ‡1 )οΈ‚1/4 (οΈ‚ 8 |𝑀2 (𝑦, 𝑑)| 6 𝐸(𝑑)𝐹2 (𝑑) , (3.22) 𝜈2 [οΈ‚ (οΈ‚ )οΈ‚]οΈ‚1/4 (οΈ‚ )οΈ‚1/4 𝐸1 (𝑑) 4πœ…2 𝑙(π‘Žβ€²1 (𝑙, 𝑑))2 12𝜈1 8𝐸(𝑑) |𝑓1 (𝑑)| 6 2 𝐹3 (𝑑) + + 2 𝐹2 (𝑑) , 𝜈1 5πœ‡1 𝑙 𝜈1 264 Mathematical and Information Technologies, MIT-2016 β€” Mathematical modeling (οΈ‚ )οΈ‚1/4 (οΈ‚ )οΈ‚1/4 8𝐸1 (𝑑) 12𝜈2 8𝐸(𝑑) |𝑓2 (𝑑)| 6 𝐹3 (𝑑) + 𝐹 2 (𝑑) , (3.23) 𝜈2 (β„Ž βˆ’ 𝑙)2 𝜈2 where βˆ«οΈπ‘‘ 𝐸(𝑑) 6 π‘’βˆ’4𝛿1 𝑑 𝐻(𝜏 )𝑒4𝛿1 𝜏 π‘‘πœ, (3.24) 0 [οΈ‚(οΈ‚ )οΈ‚1/2 ]οΈ‚ 2πœ… 8πœ’1 2 𝐻(𝑑) = 𝐹 (𝑑)𝐴(𝑑) + π‘Ž10 (𝑑) . (3.25) πœ€ π‘˜12 The functions 𝐸1 (𝑑), 𝐹1 (𝑑), 𝐹2 (𝑑) and 𝐹3 (𝑑) have the same structures as 𝐸(𝑑), 𝐹 (𝑑). 4 Stationary Flow The problem (3.1)–(3.10) has the stationary solution 𝑀𝑗𝑠 (𝑦), π‘Žπ‘ π‘— (𝑦), 𝑓𝑗𝑠 πœ…(1 βˆ’ 𝛾)π΄β„Ž(3𝑦 2 /β„Ž2 βˆ’ 2𝛾𝑦/β„Ž) 𝑀1𝑠 (𝑦) = , 2π›Ύπœ‡2 [𝛾 + πœ‡(1 βˆ’ 𝛾)] πœ…π›Ύπ΄β„Ž(3𝑦 2 /β„Ž2 βˆ’ 2(2 + 𝛾)𝑦/β„Ž + 1 + 2𝛾) 𝑀2𝑠 (𝑦) = , 2(1 βˆ’ 𝛾)πœ‡2 [𝛾 + πœ‡(1 βˆ’ 𝛾)] (π‘Žπ‘ 20 βˆ’ π‘Žπ‘ 10 ) 𝑦 π‘Žπ‘ 1 = + π‘Ž10 , (4.1) [𝛾 + π‘˜(1 βˆ’ 𝛾)] β„Ž 1 [︁ 𝑦 ]︁ π‘Žπ‘ 2 = π‘˜(π‘Žπ‘ 20 βˆ’ π‘Žπ‘ 10 ) + π‘˜π‘Žπ‘ 10 + 𝛾(1 βˆ’ π‘˜)π‘Žπ‘ 20 , 𝛾 + π‘˜(1 βˆ’ 𝛾) β„Ž 3πœ…πœˆ(1 βˆ’ 𝛾)𝐴 3πœ…π›Ύπ΄ 𝑓1𝑠 = βˆ’ , 𝑓2𝑠 = βˆ’ , π›Ύβ„ŽπœŒ2 [𝛾 + πœ‡(1 βˆ’ 𝛾)] (1 βˆ’ 𝛾)β„ŽπœŒ2 [𝛾 + πœ‡(1 βˆ’ 𝛾)] π‘Žπ‘ 1 (0) = π‘Žπ‘ 10 , π‘Žπ‘ 2 (β„Ž) = π‘Žπ‘ 20 , π‘˜ = π‘˜1 /π‘˜2 , 𝜈 = 𝜈1 /𝜈2 , 𝛾 = 𝑙/β„Ž < 1, πœ‡ = πœ‡1 /πœ‡2 , (π‘Žπ‘ 20 βˆ’ π‘Žπ‘ 10 )𝛾 𝐴= ; (4.2) 𝛾 + π‘˜(1 βˆ’ 𝛾) (οΈ‚ 3 )οΈ‚ 𝑠 πœ…(1 βˆ’ 𝛾)π΄β„Ž 𝑦 𝛾𝑦 2 𝑣1 (𝑦) = βˆ’ βˆ’ 2 , 2π›Ύπœ‡2 [𝛾 + πœ‡(1 βˆ’ 𝛾)] β„Ž3 β„Ž [οΈ‚(οΈ‚ )οΈ‚ 𝑠 πœ…π›Ύπ΄β„Ž2 𝑦3 3 (4.3) 𝑣2 (𝑦) = βˆ’ βˆ’π›Ύ 2(1 βˆ’ 𝛾)πœ‡2 [𝛾 + πœ‡(1 βˆ’ 𝛾)] β„Ž3 (οΈ‚ 2 )οΈ‚ (οΈ‚ )οΈ‚]οΈ‚ 𝑦 𝑦 βˆ’(2 + 𝛾) 2 βˆ’ 𝛾 2 + (1 + 2𝛾) βˆ’π›Ύ . β„Ž β„Ž Introducing the differences 𝑁𝑗 (𝑦, 𝑑) = π‘Žπ‘ π‘— (𝑦) βˆ’ π‘Žπ‘— (𝑦, 𝑑), 𝑀𝑗 (𝑦, 𝑑) = 𝑀𝑗𝑠 (𝑦) βˆ’ 𝑀𝑗 (𝑦, 𝑑) (4.4) 265 Mathematical and Information Technologies, MIT-2016 β€” Mathematical modeling and carrying out the calculations analogous to those in Section 2, we can prove that the solution of the nonstationary problem reaches the steady regime 𝑀𝑗𝑠 (𝑦), π‘Žπ‘ π‘— (𝑦) and 𝑓𝑗𝑠 under the conditions of convergence of the integrals ∫︁∞ ∫︁∞ ∫︁∞ 𝑒 π›Ώπœ |π‘Žπ‘ π‘—0 βˆ’ π‘Žπ‘—0 (𝜏 )| π‘‘πœ, 𝑒 π›Ώπœ |π‘Žβ€²π‘—0 (𝜏 )| π‘‘πœ, π‘’π›Ώπœ |π‘Žβ€²β€²π‘—0 (𝜏 )| π‘‘πœ. (4.5) 0 0 0 More exactly, ‖𝑀𝑗 (𝑦, 𝑑) βˆ’ 𝑀𝑗𝑠 (𝑦)β€– 6 𝑑𝑗 π‘’βˆ’π›Ώ1 𝑑 , β€–π‘Žπ‘— (𝑦, 𝑑) βˆ’ π‘Žπ‘ π‘— (𝑦)β€– 6 𝑙𝑗 π‘’βˆ’π›Ώ2 𝑑 , ‖𝑓𝑗 (𝑑) βˆ’ 𝑓𝑗𝑠 β€– 6 𝑁 π‘’βˆ’π›Ώ3 𝑑 with the positive constant 𝑑𝑗 , 𝑙𝑗 , 𝑁 , 𝛿1 , 𝛿2 , 𝛿3 depending on physical parameters of liquid and layers thicknesses. 5 Nonstationary Motion and Numerical Results To describe the nonstationary motion of two viscous thermally conducting liquids the Laplace transform will be applied to problem (3.1)–(3.10). As a result we come to boundary value problem for images π‘Ž ˆ𝑗 (𝑦, 𝑝) of functions π‘Žπ‘— (𝑦, 𝑑) 𝑝ˆ π‘Ž π‘Ž0𝑗 (𝑦) ˆ𝑦𝑦 βˆ’ π‘Ž =βˆ’ , (5.1) πœ’π‘— πœ’π‘— π‘Ž Λ†1 (0, 𝑝) = π‘Ž Λ†10 (𝑝), π‘Ž Λ†2 (β„Ž, 𝑝) = π‘Ž Λ†20 (𝑝), (5.2) π‘Ž Λ†1 (𝑙, 𝑝) = π‘Ž Λ†2 (𝑙, 𝑝), π‘˜1 π‘Ž Λ†1𝑦 (𝑙, 𝑝) = π‘˜2 π‘Ž Λ†2𝑦 (𝑙, 𝑝), (5.3) ˆ𝑗 (𝑦, 𝑝) and 𝑓ˆ𝑗 (𝑝) of functions 𝑀𝑗 (𝑦, 𝑑), 𝑓𝑗 (𝑑) and images 𝑀 𝑝 𝑓ˆ𝑗 (𝑝) ˆ𝑗𝑦𝑦 βˆ’ 𝑀 ˆ𝑗 = βˆ’ 𝑀 , (5.4) πœˆπ‘— πœˆπ‘— 𝑀 Λ†1 (0, 𝑝) = 0, 𝑀 Λ†2 (β„Ž, 𝑑) = 0, (5.5) 𝑀 Λ†1 (𝑙, 𝑝) = 𝑀 Λ†2 (𝑙, 𝑝), (5.6) Λ†2𝑦 (𝑙, 𝑝) βˆ’ πœ‡1 𝑀 πœ‡2 𝑀 Λ†1𝑦 (𝑙, 𝑝) = βˆ’2πœ…Λ† π‘Ž1 (𝑙, 𝑝), (5.7) ∫︁ 𝑙 βˆ«οΈβ„Ž 𝑀 Λ†1 (𝑦, 𝑝) 𝑑𝑦 = 0, 𝑀 Λ†2 (𝑦, 𝑝) 𝑑𝑦 = 0. (5.8) 0 𝑙 In condition (5.2) and equation (5.4) π‘Ž ˆ𝑗0 (𝑝), 𝑓ˆ𝑗 (𝑝) are images of functions π‘Žπ‘—0 (𝑑), 𝑓𝑗 (𝑑) respectively. The solutions of problem (5.1)–(5.8) can be written as βˆšοΈ‚ βˆšοΈ‚ βˆ«οΈπ‘¦ βˆšοΈ‚ 𝑝 𝑝 1 𝑝 ˆ𝑗 (𝑦, 𝑝) = 𝐢𝑗1 sh π‘Ž 𝑦+𝐢𝑗2 ch π‘¦βˆ’ √ π‘Ž0𝑗 (𝑧) sh (π‘¦βˆ’π‘§) 𝑑𝑧, (5.9) πœ’π‘— πœ’π‘— π‘πœ’π‘— πœ’π‘— 𝑦𝑗 [οΈ‚ βˆšοΈ‚ βˆšοΈ‚ ]οΈ‚ 𝑝 𝑝 𝐿𝑗 (𝑝) 𝑀 π‘Ž1 (𝑙, 𝑝) 𝐷𝑗1 sh ˆ𝑗 (𝑦, 𝑝) = βˆ’2πœ…Λ† 𝑦 + 𝐷𝑗2 ch 𝑦+ , (5.10) πœˆπ‘— πœˆπ‘— 𝑝 266 Mathematical and Information Technologies, MIT-2016 β€” Mathematical modeling where 𝑓ˆ𝑗 (𝑝) = βˆ’2πœ…Λ† π‘Ž1 (𝑙, 𝑝)𝐿𝑗 (𝑝). The values 𝐢𝑗 , 𝐢𝑗2 , 𝐷𝑗1 , 𝐷𝑗2 and 𝑓ˆ𝑗 (𝑝) determined from the boundary condi- 1 tions (5.2), (5.3), (5.5)–(5.8). Due to the cumbersome the type of these values is not given here. Let us assume that lim π‘Žπ‘—0 (𝑑) = π‘Žπ‘ π‘—0 , 𝑗 = 1, 2, using the formulas (5.9), (5.10) π‘‘β†’βˆž and presenting for the values 𝐢𝑗1 , 𝐢𝑗2 , 𝐷𝑗1 , 𝐷𝑗2 and 𝑓ˆ𝑗 (𝑝) we can prove the limit equalities lim π‘Žπ‘— (𝑦, 𝑑) = π‘Žπ‘ π‘— (𝑦), lim 𝑀𝑗 (𝑦, 𝑑) = 𝑀𝑗𝑠 (𝑦), π‘‘β†’βˆž π‘‘β†’βˆž lim 𝑓𝑗 (𝑑) = 𝑓𝑗𝑠 , π‘‘β†’βˆž where π‘Žπ‘ π‘— (𝑦), 𝑀𝑗𝑠 (𝑦), 𝑓𝑗𝑠 are determined by formulas (4.1), (4.2). Let us apply the numerical method of inversion of Laplace transformation to obtained formulas (5.9), (5.10). The graphs only for the velocities are given because the have a real physical meanings. All numerical calculations were made for the system of liquid silicon–water. Thickness of the layers is the same and equal to 1 mm. The corresponding values of the defining parameters are given in Table 1. Table 1. Physical properties of liquids Item liquid silicon water 3 𝜌, kg/m 956 998 βˆ’6 2 𝜈 Γ— 10 , m /s 10.2 1.004 π‘˜, kg Β· m/s3 Β· K 0.133 0.597 βˆ’6 2 πœ’ Γ— 10 , m /s 0.0675 0.143 Γ¦ Γ— 10βˆ’5 , kg/s2 Β· K 6.4 15.14 Figure 2–5 show the profiles of the dimensionless functions 𝑀¯𝑗 (πœ‰, 𝜏 ) = 𝑀𝑗 (𝑦, 𝑑)πœ‡2 /(πœ…π΄) (πœ‰ = 𝑦/𝑙, 𝜏 = 𝜈1 𝑑/𝑙2 are the dimensional variables) and transverse velocity 𝑣¯𝑗 (πœ‰, 𝜏 ) = 𝑣𝑗 (𝑦, 𝑑)πœ‡2 /(πœ…π΄β„Ž) with π‘Ž20 (𝑑) = 0. In particu- lar, the functions 𝑀 ¯𝑗 are negative, so reverse flows arise here. Figure 2, 3 show the results of calculations when π‘Ž10 (𝜏 ) = sin 𝜏 , π‘Ž20 (𝜏 ) = 0. That is the limit of π‘Ž10 (𝜏 ) at 𝜏 β†’ ∞ does not exist and the velocity field does not converge to a stationary one. Figure 4, 5 show an evolution of the convergence of functions 𝑀 ¯𝑗 and trans- verse velocities 𝑣¯𝑗 to stationary regime for the case π‘Ž10 (𝜏 ) = 1 + π‘’βˆ’πœ cos(10𝜏 ), π‘Ž20 (𝜏 ) = 0. These results are good agreement with the a priori estimates were obtained in Section 4. 267 Mathematical and Information Technologies, MIT-2016 β€” Mathematical modeling 0.08 wΜ„ 0.04 0 βˆ’0.04 0 0.5 1 1.5 2 ΞΎ Fig. 2. Evolution of functions 𝑀¯𝑗 for π‘Ž10 (𝜏 ) = sin 𝜏 . Total line is the stationary profiles, βˆ’βˆ’ – 𝜏 = 4, βˆ’ Β· βˆ’ – 𝜏 = 5, Β· Β· Β· – 𝜏 = 7 0.02 0.01 0 vΜ„ βˆ’0.01 βˆ’0.02 0 0.5 1 1.5 2 ΞΎ Fig. 3. Evolution of functions 𝑣¯𝑗 for π‘Ž10 (𝜏 ) = sin 𝜏 . Total line is the stationary profiles, βˆ’βˆ’ – 𝜏 = 3, βˆ’ Β· βˆ’ – 𝜏 = 6, Β· Β· Β· – 𝜏 = 8 0.08 0.04 wΜ„ 0 βˆ’0.04 0 0.5 1 1.5 2 ΞΎ Fig. 4. Evolution of functions 𝑀 ¯𝑗 for π‘Ž10 (𝜏 ) = 1 + π‘’βˆ’πœ cos(10𝜏 ). Total line is the stationary profiles, βˆ’βˆ’ – 𝜏 = 1, Β· Β· Β· – 𝜏 = 4 268 Mathematical and Information Technologies, MIT-2016 β€” Mathematical modeling 0.02 0.01 0 vΜ„ βˆ’0.01 βˆ’0.02 0 0.5 1 1.5 2 ΞΎ Fig. 5. Evolution of functions 𝑣¯𝑗 for π‘Ž10 (𝜏 ) = 1 + π‘’βˆ’πœ cos(10𝜏 ). Total line is the sta- tionary profiles, βˆ’βˆ’ – 𝜏 = 1, βˆ’ Β· βˆ’ – 𝜏 = 4, Β· Β· Β· – 𝜏 = 2 6 Conclusion The two-dimensional horizontal layer is a matter of great importance in connec- tion with the theory of convective stability applications in the design of cooling systems, in studying the growth of crystals and films, or in the aerospace in- dustry. We have presented a theoretical and numerical study of a creeping flow of two immiscible viscous heat conducting liquids in thin layers. The flow arises due to heat exchange with the localized parabolic heating of the borders and through the thermocapillary forces on the interface. The following results are obtained: (1) the exact solution describing the stationary thermocapillary con- vective flow is found; (2) a priori estimates of the initial boundary value problem are established and sufficient conditions on input data when solution tends to stationary one are obtained; (3) the solution of the non-stationary problem in the form of final analytical formulas in the Laplace representation is found and some numerical results of velocities behaviour in layers are presented. Acknowledgments. This research was supported by the Russian Foundation for Basic Research (14-01-00067). References 1. Andreev V. K., Zahvataev V. E., Ryabitskii E. A.: Thermocapillary Instability. Nauka. Sibirskoe otdelenie. Novosibirsk (2000) 2. Nepomnyashii A., Simanovskii I., Legros J.-C.: Interfacial Convection in Multilayer System. Springer. New-York (2006) 3. Narayanan R., Schwabe D.: Interfacial Fluid Gynamics and Transport Processes. Springer-Verlag. Berlin. Heidelberg. New-York (2003) 4. Zeytovnian R. Kh.: Convection in Fluids. Springer. Dordrecht. Heidelberg. London. New-York (2009) 269 Mathematical and Information Technologies, MIT-2016 β€” Mathematical modeling 5. Andreev V. K.: The Birikh Solution of Convection Equations and it Some General- ization. Preprint 1–10. Institute of Computational Modelling SB RAS. Krasnoyarsk (2010) 6. Andreev V. K., Kaptsov O. V., Pukhnachov V. V., Radionov A. A.: Applications of Group-Theoretical Methods in Hydrodynamics. Kluwer Academic Publisher. Dordnrcht-Boston-London (1998) 7. Andreev V. K., Gaponenko Y. A., Goncharova O. N., Pukhnachov V. V.: Math- ematical Models of Convection. Walter de Gruyter GmbH and Co. KG. Berlin- Boston (2012) 8. Napolitano L. G.: Plane Marangoni-Poiseuille flow two immiscible fluids. Acta Astronautica. Vol. 7. No. 4-5. 461-478 (1980) 9. Andreev V. K., Bekezhanova V. B.: Stability of Non-isothermal Fluids. Siberian Federal University (2010) 10. Andreev V. K. On a conjugate initial boundary value problem. Diff. eq. No. 5. 17 (2008) 11. Andreev V. K., Bekezhanova V.B.: Stability of Non-isothermal Fluids (Review). Jour. Appl. Mech. and Tech. Phys. Vol. 54. No. 2. 171-184 (2013) 12. Andreev V. K., Lemeshkova E.N.: Evolution of the thermocapillary motion of three liquids in a plane layer. Appl. Math. and Mech. Vol. 78. No. 4. 485-492 (2014) 13. Denisova I. V.: On the problem of thermocapillary convection for two incompress- ible fluids separated by a closed interface. Prog. Nonlinear Differ. Equ. Appl. Vol. 61. 45-64 (2005) 14. Denisova I. V.: Thermocapillary convection problem for two compressible immis- cible fluids. Microgravity Sci. Technol. Vol. 20 No. 3-4 287-291 (2008) 15. Hiemenz K.: Die Grenzschicht an einem in den gleichförmigen Flüssigkeitsstrom eingetauchten geraden Kreiszylinder. Dinglers Poliytech. J. V. 326. 321-440 (1911) 16. Brady J. F., Acrivos A.: Steady flow in a channel or tube with an accelerating surface velocity. J. Fluid Mech. V. 112. 127-150 (1981) 17. Bobkov N. N., Gupalo Yu. P.: The structure of the flow in the liquid layer and the spectrum of the boundary value problem for the nonlinear dependence of the surface tension from temperature. Appl. Math. and Mech. Vol. 60. No. 6. 1021-1028 (1996) 18. Gupalo Yu. P., Ryazantsev Yu. S.: On a thermocapillary motion of a fluid with a free surface at the nonlinear dependence of the surface tension from temperature. Fluid Dynamics. No. 5. 132-137 (1988) 19. Gupalo Yu. P., Ryazantsev Yu. S., Skvortsov A. V.: Influence of thermal capillary forces on the fluid flow with free boundary. Fluid Dynamics. No. 5. 3-7 (1989) 270