- Email: [email protected]
Multilayer system of films heated from above
International Journal of Heat and Mass Transfer 114 (2017) 992–1000
Contents lists available at ScienceDirect
International Journal of Heat and Mass Transfer journal homepage: www.elsevier.com/locate/ijhmt
Multilayer system of films heated from above A.S. Ovcharova Lavrentyev Institute of Hydrodynamics, Siberian Branch, Russian Academy of Sciences, 15, Acad. Lavrentyev Avenue, Novosibirsk 630090, Russia
a r t i c l e
i n f o
Article history: Received 11 February 2017 Received in revised form 15 June 2017 Accepted 27 June 2017
Keywords: Multilayer system Viscous nonisothermal fluid Interfaces Marangoni effect
a b s t r a c t The paper deals with the development of convective flows in a complex system of thin liquid layers occurring under the effect of thermal load. To study this process a mathematical model which is based on the decomposition of complex problem into unitary elements (modules) with some set of rules allowing their linkage with each other was constructed. Each module represents a monotypic local model. The fluid flow and heat transfer in each of these modules are described by the Navier-Stokes and thermal conductivity equations. The boundary conditions required to solve these systems of equations are represented explicitly and written in the form of conservation laws. The numerical analysis of effect of thermal load on the characteristics of the liquid layers movement depending on the Marangoni number is carried out. It is shown that physical and thermophysical properties of liquid layers play the decisive role. It is exactly they control the convective flows in layers and determine the position and shape of interfaces. The results of model problem solutions are presented. Ó 2017 Elsevier Ltd. All rights reserved.
1. Introduction Study of convective flows and interfaces in a complex system of thin liquid layers is very important due to the significance of problems, which are related to a broad industrial application of thin liquid layers (films). The continuous development of science-intensive technologies significantly expands the range of problems associated with mathematical modeling of the convection processes of viscous heat-conducting liquids. In this regard, the number of publications and even reviews on this subject is so great (and it is increasing every year) that their list is endless. Depending on the set goal, the authors use different approaches to both the description of the processes occurring in thin layers of viscous liquid, and the numerical solution of the arising tasks. Here, the correctness of the selection of an appropriate model describing the phenomenon under study is, of course, of decisive importance. Mathematical models for describing the processes of mass and heat transfer require correct conditions at the interfaces of liquids. Especially in the case of convective currents of two or more liquid mediums contacting at deformable interfaces. One of the difficulties in solving problems with deformable boundaries is that along with the unknown flow characteristics, it is also necessary to determine the position of the interfaces which change during the motion.
E-mail address: [email protected] http://dx.doi.org/10.1016/j.ijheatmasstransfer.2017.06.123 0017-9310/Ó 2017 Elsevier Ltd. All rights reserved.
System-based derivation of conditions at the interface between two immiscible liquids under nonisothermal unstable flow was first given in [1]. The modification of these considerations is given in [2]. The purpose of this work, which can be considered as a continuation of the work [3], is the construction of a mathematical model for studying the processes arising in a system consisting of any finite number of thin liquid layers with deformable interfaces. This mathematical model is formulated on the basis of known hypotheses about the physical processes provided by the implementation of conservation laws. It allows us to investigate a wide range of problems associated with the impact of temperature or surfactants on a multilayer system. All the boundary conditions for solving local tasks in each liquid layer are presented in explicit form. The structure of the model allows it to expand by adding to it new data on the physical processes ensured the implementation of conservation laws. The idea of constructing a model is based on the decomposition of complex problem into unitary elements (modules) with some set of rules allowing their linkage with each other [4]. Each module represents a monotypic local model, within which the calculation is carried out independently of the other modules of the system, while the ‘‘set of rules” defining, for example, the common boundaries between the modules and the corresponding functional relations, allow connecting these modules back to the original system. In this case the decisive role is played by physical parameters and thermophysical properties of liquid layers, which control the convective flows in the layers and determine the position and shape of interfaces.
A.S. Ovcharova / International Journal of Heat and Mass Transfer 114 (2017) 992–1000
993
Nomenclature t x,y vs,vn D P ! ! s; n d hi/h0 h0
v0
p0 g0
time Cartesian coordinates tangent and normal velocities of points lying on interface stress tensor pressure tangent and normal vectors width of thermal beam dimensionless thickness of layers length scale velocity scale pressure scale gravity scale
Dimensionless parameters Pri Prandtl number in i-th layer Capillary number Cai Mni Marangoni number Cri Crispation number G Galileo number Greek symbols w stream function x vorticity h temperature m kinematic viscosity q density
@wi ; @y
vi ¼
@wi ; @x
@ v i @ui : @x @y
2. Mathematical model
ui ¼
Consider a multilayer system consisting of thermallyconductive immiscible liquid layers as shown in Fig. 1. Here Gi are liquid layers, fj (t,x) are the interfaces between them. The bottom layer is lying on a solid substrate f J ðt; xÞ, the top layer is in contact with the air through the interface f0 (t,x); at the left and right the system of layers is bounded by two solid planes x = 0 and x = L. Each liquid layer Gi is characterized by its thickness hi, density qi, kinematic viscosity mi and surface tension coefficient ri ðTÞ. Moreover, the evaporation and condensation processes are not taken into account in this model; qi and mi are considered to be constant. Such a complex system of immiscible liquid layers can be represented as a flat ribbon graph consisting of two different types of elements: liquid layers Gi, (i = 1, . . ., I) and the interfaces between them fj (j = 1, . . ., J 1), which we will call the internal interfaces. For the unequivocal description of the graph we must determine the order of the conformity of its elements. With this purpose, we enumerate all the interfaces and all liquid layers, for example, as is done in Fig. 1. Selected numbering scheme must be fulfilled during the whole process of solving the problem. According to the chosen numbering scheme, just one liquid layer with the number i + 1 adjoins to internal interface from the bottom, while the layer with the number i adjoins from the top of the interface. Each liquid layer Gi is limited from above by the interface fj (t,x), while from below – by the interface fj+1 (t,x). Surfaces x = 0, x = L, f0 (t,x) and fJ (t,x) have certain features and are classified as a special type surfaces. The orientation of the graph in the plane is presented in Fig. 1.
The characteristic values (scale) of length x0, velocity v0 and pressure p0 should be the same for each layer of the system of layers. Here Rei = v0 hi/mi and Pri = mi/vi, where vi is the thermal diffusivity coefficient, are the Reynolds and the Prandtl numbers in i-th layer of the system, respectively. hi ¼ ðT i T 0 Þ=dT, where T 0 is the characteristic value of the temperature; dT is the characteristic temperature difference for the entire system of layers. In addition, each liquid layer is characterized by its surface tension coefficient ri ðTÞ at the interface between the liquid and air. We assume that the surface tension coefficient of each liquid is a linear function of temperature
2.1. Consider the first type of graph elements – liquid layers
xi ¼
ðiÞ ri ðTÞ ¼ r0ðiÞ ð1 rTðiÞ ðT i T 0 ÞÞ; rðiÞ 0 ¼ ri ðT 0 Þ; rT
¼
1 dri rðiÞ dT 0
ðiÞ
; rT > 0:
T¼T 0
i.e. it decreases linearly with the increase of temperature. 2.2. Interfaces. functional relations Interfaces are the second type of graph elements. Consider any internal interface between the two layers Gi and Gi+1. As we consider any internal surface, we can define it by the equation y = f (t,x). It is however important to keep in mind that the liquid with number i + 1 lies below this interface, while the liquid with the number i is over it (see Fig. 2). It is related to the fact that when constructing the model we predetermined the order of conformity of graph elements. Define the unit tangent and normal vectors at the interface f (t, x) of liquids as
The fluid flow and heat transfer in each of the liquid layer Gi, (i = 1, . . ., I) is described by the system of Navier-Stokes and thermal conductivity equations, which in the variables w, x, h (stream function, vorticity and temperature) have the form [5,6]
8 9 > > < 1 fx = ! ; s ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffi ; qffiffiffiffiffiffiffiffiffiffiffiffiffi > 2> : 1 þ f2 1 þ fx; x
8 9 > > < f = 1 ! x n ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffi ; qffiffiffiffiffiffiffiffiffiffiffiffiffi > 2> : 1 þ f2 1 þ fx; x
@ xi @ @w @ @w 1 xi i x i i ¼ Dx i ; þ @x @y Rei @t @y @x
ð1Þ
! Note that for liquid Gi+1, n is the outer normal unit vector at the interface of liquids, while for the liquid Gi for the same interface this normal unit vector is directed into the liquid.
Dwi þ xi ¼ 0;
ð2Þ
@hi @ @w @ @w 1 hi i hi i ¼ Dh i þ @y Rei Pri @t @x @y @x
ð3Þ
2.2.1. Functional relations for temperature In this paper we do not consider the evaporation (condensation) process at the liquid interface and neglect energy consumption related to its deformation (since they are negligible). Therefore,
994
A.S. Ovcharova / International Journal of Heat and Mass Transfer 114 (2017) 992–1000
Fig. 1. Sketch showing the agreement between elements of the film system placed on a solid substrate y = 0. Here, Gi are the layers, fj are the interfaces of the system; d is the width of heat beam effecting on interface f0.
the temperature conditions at the interface f (t,x) can be written in the form of temperature and heat flux continuity conditions.
hiþ1 hi ¼ 0;
@hiþ1 @h i ¼ 0; j @n @n
ð4Þ
¼ ji =jiþ1 is the ratio between heat conductivities of (i)-th where j and (i + 1)-th liquids, respectively. 2.2.2. Functional relations for the stream function From the equality of the tangential velocities of the liquids Gi and Gi+1 at their common interface f (t,x) and the volume conservation condition for each of these liquids, it follows
wiþ1 wi ¼ 0;
@wiþ1 @wi ¼ 0: @n @n
ð5Þ
2.2.3. Functional relations for the vorticity For derivation of these relations we use the kinematic and dynamic conditions at the interface [1,2]. The dynamic condition is based on the preservation of the continuity of the normal and tangent stresses on the interface, which is fulfilled if the mediums preserve their continuity. If we neglect the evaporation (condensation) process, the dynamic condition at the liquid interface f (t,x) can be represented as ! r ðTÞ ! n þ rr ðTÞ; ðDiþ1 Di Þ n ¼ R
ð6Þ
where Di, Di+1 are the stress tensors in liquids Gi and Gi+1, respectively, r ðTÞ is the interfacial tension coefficient between these liquids, R is the radius of the liquid interface curvature, rr ðTÞ is the interfacial tension gradient along the interface. Thus, it only remains to consider what is meant by the function r ðTÞ and, most importantly, how changes its distribution along the interface depending on temperature. There are several approaches for defining this function depending on physical, thermophysical and even
thermodynamic properties of liquids (for example, the polarity of the molecules). We can apply the commonly used Antonow’s rule [7], which reads as follows: interfacial tension at the interface of two immiscible liquids is equal to the difference between surface tensions of these liquids (at the interface with air or intrinsic vapor) in the conditions of mutual saturation. That is, just in our case
r ðTÞ ¼ riþ1 ðTÞ ri ðTÞ
ð7Þ
This rule is valid for many pairs of liquids, though many exceptions to this rule are known as well. It may also happen that a liquid Gi in contact with a liquid Gi+1 satisfies Antonow’s rule, while for the same liquid Gi in contact with another liquid Gi-1 this rule does not work. There are also other formulas to determine the interfacial tension between the liquids [8,9]. However, they are all based on the use of surface tension coefficients of these liquids at their contact with air or intrinsic vapor (which are assumed to be known in advance) together with empirical coefficients which are inherent of the considered liquids. In this case, we can always determine not only interfacial tension between liquids, but also its distribution along the interface depending on temperature. When solving the problem with a specific combination of liquids, the researcher determines, which formula would provide most reliable solution. In this work, when describing the model, we use Antonow’s rule. From dynamic conditions we can obtain two corollary facts. ! Carrying out scalar multiplication of (6) by vector s we get
! ! @ r ðTÞ s ðDiþ1 Di Þ n ¼ @s " # ! 2 1 f x @ 2 wiþ1 @ 2 wiþ1 @ 2 wiþ1 2f x qiþ1 miþ1 2 @y2 @x2 @y@x ! " # 2 2 2 2 1 f x @ wi @ wi @ wi qi mi 2f x 2 @y2 @x2 @y@x 2
¼
1 þ f x @ r ðTÞ : @s 2
ð8Þ
Using the kinematic condition for each of the liquids at their common interface as well as Antonow’s rule (7), the condition (8) can be reduced to the following form:
xiþ1 q mxi ¼ F 1 ðt; xÞ; where
Mni Þ m F 1 ¼ ðMniþ1 q
@h @v n v s Þ m þ 2ð1 q þ @s @s R
ð9Þ
ð9aÞ
! Carrying out scalar multiplication of (6) by vector n , we obtain ! ! r ðTÞ n ðDiþ1 Di Þ n ¼ R
Fig. 2. Layout of the layers and interface between them.
A.S. Ovcharova / International Journal of Heat and Mass Transfer 114 (2017) 992–1000
or, to give a more convenient form,
@v n f vs ðPiþ1 Pi Þ ¼ 2qiþ1 miþ1 x @n iþ1 R @v n f xv s r ðTÞ þ þ 2qi mi R @n i R
ð10Þ
Next, we consider the vector equation of motion for each of the liquids, written in the natural variables, i.e. the velocity – the pressure.
! @vi 1 ! ! ! ! þ ð v i rÞ v i ¼ rPi þ mi D v i þ g @t qi
ð11Þ
We are looking for such a solution, where vs at the interface of liquids must comply with the conditions (5) (equality of the tangential velocities of points situated at the interface). This is requirement of the physical model describing the problem under consideration. At the same time, the dynamic condition has its requirement – a pressure jump at the interface should satisfy the balance relations (6). If we take the scalar product of vector equa! tion (11) for i with the vector qi s ðt; xÞ and the equation for i + 1 ! with vector qiþ1 s ðt; xÞ and put them together, we will get the equation
ðqiþ1 qi Þ
@v s @v s @ @ þ ðqiþ1 qi Þv s ¼ ðPiþ1 Pi Þ @s @n @t @s ðqiþ1 miþ1 xiþ1 qi mi xi Þ ! ! þ ðqiþ1 qi Þ g s þ ðqiþ1 qi Þ
in which a pressure jump (Pi+1 Pi) must satisfy (6). Eliminating @ðP iþ1 Pi Þ=@s with the use of (10), after conversion to dimensionless form, we get
@ xiþ1 @ xi m q ¼ F 2 ðt; xÞ; @n @n where
ð12Þ
@ @v n @ @v n @ f xvs Þ m m þ 2 ð1 q q @s @n iþ1 @s @n i @s R 1 @ 1 ð1 Mniþ1 Caiþ1 h þ Caiþ1 @s R @ 1 Ca1 m ð1 Mn Ca h q i i i @s R 2
F 2 ðt; xÞ ¼ 2
@v s @v fx 6 1Þv s s þ ðq 1Þ G qffiffiffiffiffiffiffiffiffiffiffiffiffi 1Þ þ Reiþ1 4ðq þ ðq @t @s 2 1 þ fx # f v2 xi Þ Þ x n þ v n ðxiþ1 q ð12aÞ þð1 q R 1 f xx ; ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi R 2 3 ð1 þ f x Þ
Expression (12a) looks at first glance bulky, in fact it is quite a simple expression. Meaningful part of it consists in the first two lines, reflecting the contribution of the pressure jump. The last line is a tribute to the nonstationarity of the problem. To obtain this expression we also used Antonow’s rule (7). Recall that the relations of type (9) with the right hand part of (9a), as well as (12) with the right hand part of (12a) are valid for any internal interface between a couple of layers Gi and Gi+1, which are included to the system and have common internal interface. The position of the interface f (x,t) is determined from the kinematic condition at the interface, for example, for the lower layer. In the present case this is
ft þ
qffiffiffiffiffiffiffiffiffiffiffiffiffi 2 @wiþ1 1 þ fx ¼0 @s
q v 0 mi rðiÞ dT g h0 Cai ¼ i ðiÞ ; Mni ¼ T ; G ¼ 02 : qi v 0 mi v0 r0
Here vs, vn are tangential and normal velocity components of points lying in the interface f (x,t); R is the radius of surface curvature f (x,t); Cai is the capillary number; Mni is the Marangoni num ¼ qi =qiþ1 is the density ratio, ber; G is the Galileo number; q m ¼ mi =miþ1 is the viscosity ratio of liquids i and i + 1, respectively. plays an important role in formation of m The last parameter q the convective flows in the layers [10], especially in the cases where there is no thermal load and the convection currents are caused by other factors [11].
ð13Þ
2.3. The boundary conditions at the remaining surfaces of the system The boundary conditions for the temperature at the other boundaries take form in accordance with the requirements of the particular tasks (see numerical calculations section). At x = 0, x = L, y = fJ (x), i.e. at solid surfaces, we use the no-slip boundary conditions
wi ¼ 0;
f x v 2n þ v n ðqiþ1 xiþ1 qi xi Þ R
995
@wi ¼ 0: @n
Here we have two conditions for the stream function w while none of conditions for the vorticity x. To find the vorticity value x on a solid wall, we use the finite-difference conditions of the Thom’s type [12]. At the liquid-gas interface f0 (t,x) we neglect not only the evaporation (condensation) process but also the effects of the dynamic characteristics of the gas to the liquid motion. We can’t use functional relations for w and x, obtained for internal interfaces, as in this case we will have two conditions for x and no one condition for stream function w. Therefore, it makes sense to distinguish liquid-gas interface and write down boundary conditions on this interface according to its features. Again, we denote the interface f0 (t,x) just as f (t,x) and have the same formulations as in the case of the internal interface with just one exception. A second corollary fact of balance ratios (6) allows representation
@v n f xv s @v s v n ! ! r ðTÞ ! ¼ 2qm P ¼ 2qm nDn ¼ R @n R @s R In the derivation of functional dependencies at the internal interface we used the first representation for P
@v n f xv s P ¼ 2qm @n R Whereas to obtain the boundary conditions at the liquid-gas interface we used a second representation for P. A detailed derivation of the boundary conditions at the liquid-gas interface is available in [3], here we present only the result of this work: The position of the surface f (t,x) is determined using the kinematic condition
ft þ
qffiffiffiffiffiffiffiffiffiffiffiffiffi 2 @w 1 þ fx ¼0 @s
ð14Þ
The boundary conditions take the form:
@w ¼ v s; @n
ð15Þ
996
A.S. Ovcharova / International Journal of Heat and Mass Transfer 114 (2017) 992–1000
v @v @h x ¼ 2 s þ n þ Mn ; @s R @s where
ð16Þ
vs is the solution to the equation
Bi ¼ 1;
ðiÞ
Ai ¼ 1;
R1 ¼
SðiÞ t g þ ðf jþ1 Þt
@v s @v s @ vs þ vs ¼ 2 2 þ F; @t @s @s @x @ 1 rT dT fx @ v n F¼ 2 hÞ G qffiffiffiffiffiffiffiffiffiffiffiffiffi ð1 þ Ca1 @s R @s R @n r0 2 1 þ fx
v
ð17Þ
f
ð18Þ
2
1 þ fx
3. Numerical method Suppose we have a system consisting of I liquid layers Gi, (i = 1, . . ., I) and J + 1 interfaces fj, (j = 0, . . ., J). Each layer Gi is limited at the top by the surface fj (t,x), while at the bottom – by the surface fj+1 (t,x) as shown in Fig. 1. Map each of the area occupied by the liquid Gi in a rectangle with sides of 0 6 n 6 L, 0 6 g 6 1 using the transformation:
SðiÞ ¼ f j ðt; xÞ f jþ1 ðt; xÞ
In this case, all boundaries of each layer in the system will rest on the coordinate lines. Every Eqs. (1)–(3) for each layer i can be written in the following form
h i ðiÞ ðiÞ B12 ¼ Sn g þ ðf jþ1 Þn ;
For Ui ¼ xi
ðiÞ
B22 ¼
ð19Þ
ðiÞ2
1 þ B12 SðiÞ
R1 ¼
StðiÞ g þ ðf jþ1 Þt SðiÞ
;
ðiÞ
R2 ¼ 0:
ðiÞ
Ai ¼ 0;
R1 ¼ 0;
ðiÞ
R2 ¼ ki xi ;
where ki is the iteration parameter to solve Poisson’s equation for wi . If we take new designations:
ðiÞ
VðUi Þ ¼ B12
@ Ui @w ðiÞ @ Ui Ai Ui i þ B12 @n @g @g @ Ui @w ðiÞ @ Ui þ Ai Ui i þ B22 @n @g @n
then Eq. (19) will be
@ Ui 1
ðiÞ @ Ui ðiÞ U n ðUi Þ þ V g ðUi Þ þ R1 þ R2 ¼ @t @g Bi SðiÞ
ð21Þ
Eq. (21) for hi (temperature), xi (vorticity) and wi (stream function) has been solved by the stabilizing-correction finite-difference scheme [15,16] written in the form
Ukþ1=2 Uki i
s
¼
Ukþ1 Ukþ1=2 i i
s
1 h Bi S ¼
ðiÞ
kþ1=2 i ðiÞ @ Ui ðiÞk U kn ðUi Þ þ V kþ1=2 ðUi Þ þ R1 þ R2 g @g
1 h Bi S
ðiÞ
i k U kþ1 n ðUi Þ U n ðUi Þ
This scheme is one of the most efficient finite-difference schemes which uses fractional steps. The first fractional step gives the full approximation of the equation, while the second step is a correctional one, and is used with the aim of improving stability. The transition to new variables in the boundary conditions and finite-difference technique when solving the problem for a single liquid layer is described in detail in [3]. Here we only point out that because of connectedness of the boundary conditions at the internal interfaces, numerical solution in the y - direction should be carried out consistently, for example, moving from the bottom layer GI to the top layer G1. Next, we consider the general scheme of solution for the system of layers. Let the positions of the interfaces fj (t,x) as well as distribution of temperature hi , vorticity xi , and the stream function wi are known for time moment t k ¼ ks (s is the time step) for all liquid layers, constituting the system. Stage 1. Using Eqs. (13) and (14), we find the new position of all interfaces f j ðtkþ1 ; xÞ for time moment tkþ1 ¼ ðk þ 1Þs. Define normal velocities v n of points at these interfaces. ðiÞ
Here ðiÞ
1 ; ki
ðiÞ
@ Ui 1 @ @w ðiÞ @ Ui ðiÞ @ Ui B11 Ai Ui i ¼ þ B12 ðiÞ @t @n @g @g Bi S @n @ @ U @ U @w i i ðiÞ ðiÞ ðiÞ @ Ui ðiÞ B þ Ai Ui i þ R1 þ R2 þ B22 þ @ g 12 @n @g @n @g
ðiÞ
Ai ¼ Pri ;
UðUi Þ ¼ B11
Recall that Eqs. (14)–(18) with all the motion characteristics and the governing parameters are relevant only to the liquid-gas interface f (t,x), i.e. where we neglect the effect of the dynamic characteristics of the gas. Let summarize the previous constructions. To obtain quantitative characteristics of liquid motion and temperature distribution in the considered system of layers, we must solve Eqs. (1)–(3) for each layer. The boundary conditions required for numerical solution of these equations and determination of the interfaces are presented in an explicit form. These are conditions (4), (5), (9), (9a), (12), (12a), and (13). The boundary conditions for the liquid-gas interface, where we neglect the effect of the dynamic characteristics of the gas, are classified as a special type. These are conditions (14)–(18). It should be noted, that closing relations (i.e. boundary conditions) in this model are not only derived from conservation laws, but written in the form of the conservation laws in accordance with the assumptions and limitations, which underlie the model. Besides, these conditions are valid for any finite number of liquid layers that make up the system. In conclusion, we note that the model can be adapted for the solution of a wide range of tasks, such as calculations with variable viscosity of liquids [13], the calculations connected with the evaporation process from one layer into another one [14], etc.
B11 ¼ SðiÞ ;
Bi ¼ Pri ;
Bi ¼
R
y ¼ gSðiÞ þ f jþ1 ðt; xÞ;
ðiÞ
R2 ¼ 0:
For Ui ¼ wi
2 n
t v n ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffi
x ¼ n;
;
For Ui ¼ hi
2
þ v nx þ f x
SðiÞ
ð20Þ
ðiÞ
ðiÞ
Stage 2. Using Eq. (20) determine the coefficients B11 , B12 , B22 required to the solution of the motion and heat conduction equations in each of the system’s layers. Stage 3. At this important stage we find the boundary conditions for the solution of the motion equations: using Eqs. (9a) and (12a) define functions F 1 ðn; gÞ and F 2 ðn; gÞ, which are included in the right hand part of the boundary conditions (9) and (12); solving Eq. (17), we find the tangential velocity vs of points lying at the surface f0 (t,x) that allows finding the boundary conditions (15) and (16) at this boundary. The importance of this
997
A.S. Ovcharova / International Journal of Heat and Mass Transfer 114 (2017) 992–1000
stage is that it allows finding all boundary conditions required for solution of all equations for all liquid layers of the system. Stage 4. Finally, solving Eq. (21) with the new coefficients ðiÞ
ðiÞ
ðiÞ
B11 ðn; gÞ, B12 ðn; gÞ, B22 ðn; gÞ consequently for temperature hi , vorticity xi , and stream function wi for all layers of the system, we obtain a new distribution of these functions at the new time step for the iteration l = 0. The iteration process is necessary due to the nonlinearity of Eq. (21). Iteration continues until the following condition
lþ1 l Ui ðn; gÞ Ui ðn; gÞ
x1 ¼ L=2 mh0 ;
x2 ¼ L=2 þ mh0
We consider the concentrated thermal load acting on film free surface of such type:
(
hðt; xÞ ¼
h ðxÞ; if x2 6 x 6 x3 @h ; @n
otherwise
Here, L is the length of all layers of the system; m is positive number, which determines a width of the heat beam affecting the film free surface. At solid walls (x = 0, x = L, y = 0) the temperature h = 0. In all calculations m = 3, the width of the thermal beam is d = 6h0, h⁄ = 1. In problems of this type there are a large number of governing parameters, even if considering just one layer. When considering a larger number of liquid layers, the number of governing parameters increases multiple times. We focused on the study of the influence of Marangoni effect on the behavior of the system of liquid
layers and interfaces between them. In fact, exactly this effect causes the convection motions in the layers and, as the consequence, the deformation of interfaces. Figs. 3–6 present the calculation results, carried out for the system consisting of two thin layers. Liquid G1 is situated between the surfaces f1 and f0, while liquid G2 is confined by the surfaces f2 and f1. Layers thickness, governing parameters and physical properties of liquids are chosen in such a way as to demonstrate essential deformation of interfaces up to the rupture of the layer. The calculations are performed for following parameters: L/h1 = 90, q1 = 780 kg/m3, Ca1 = 0.025
m1 ¼ 1:2 106 m2 =s,
Pr1 = 16,
L/h2 = 90, q2 = 1980 kg/m3, m2 ¼ 1 106 m2 =s, Pr2 = 7, Ca2 = 0.02 the ratio of the thermal conductivities of the first and second ¼ 0:5. liquids is j Fig. 3a presents the general pattern of convection motions in layers at Mn1 = 15 and Mn2 = 5 at the time moment corresponding to the rupture of the upper layer. Convection flow in the lower layer is accompanied by the two vortexes, which constrain the motion of the lower liquid G2 within the domain of impact of thermal beam. The liquid at the top is bordered by air, whose dynamic characteristics are neglected in this work. Because Mn1 > Mn2, the rupture occurs at the lateral parts of the interface f1, rather than in its lowest points (Fig. 3b), because as a result of heating, the surface tension of the liquid in upper layer changes significantly faster than the interfacial tension between liquids (see formulas (9), (9a) and (12), (12a)). Consider the situation where Mn1 = 5 and Mn2 = 15, while all the other characteristics of both liquids remain unchanged. In this case, the convective motion of the lower liquid G2 also will be accompanied by two vortex formations, restricting this motion within certain limits. The intensity of convective liquid motion G2 will rise slightly, though it rises essentially in the upper liquid G1 (Fig. 4a). Convection motions, caused by the action of the Marangoni effect along the surface f0 are superposed by the perturbations caused by the same effect at the surface f1. Together, they significantly increase the intensity of convective flows in liquid G1. The surface f0 cannot contact the lateral walls of the depression, as it was shown in Fig. 3b, because the Marangoni number of the liquid G1 is much less than the difference between the Marangoni numbers of these liquids. Therefore, the rapture of the upper layer occurs at the periphery d of the width of the thermal beam affecting the surface f0. At that, ‘‘corners,” formed by the surfaces f0 and f1 destroy the vortex formations in the liquid G2. We have carried out numerical experiments with an even greater difference between the Marangoni numbers of both liquids, while all other physical characteristics remained the same. The results of these studies are presented in Fig. 5a and b (Mn1 = 30, Mn2 = 5) and Fig. 6a and b (Mn1 = 5, Mn2 = 30). As we see, this resulted in change of just the quantitative characteristics of liquid flow, while the general pattern of the liquid motion and the nature of rupture of the top liquid layer remain unchanged. We have also considered the case where the Marangoni number in the top liquid layer Mn1 = 0. This means that the liquid surface tension does not depend on temperature. For greater effect, between two layers of liquids we placed another layer with the same property Mn2 = 0. The calculations were performed for the following parameters: L/h1 = 180, q1 = 780 kg/m3, Ca1 = 0.025, Mn1 = 0,
m1 ¼ 2 106 m2 =s,
Pr1 = 16,
L/h2 = 90, q2 = 998 kg/m3, Ca2 = 0.0125, Mn2 = 0,
m2 ¼ 1 106 m2 =s,
Pr2 = 7,
998
A.S. Ovcharova / International Journal of Heat and Mass Transfer 114 (2017) 992–1000
Fig. 3. (a) The general pattern of convection flows in system consisting of two layers for Marangoni numbers Mn1 = 15 and Mn2 = 5 at the time moment corresponding to the rupture of the upper layer t* = 22.03. (b) Location of the interfaces f0 and f1 in area of acting of the thermal beam.
Fig. 4. (a) The general pattern of convection flows in system consisting of two layers for Marangoni numbers Mn1 = 5 and Mn2 = 15 at the time moment corresponding to the rupture of the upper layer t* = 15.58. (b) Location of the interfaces f0 and f1 in area of acting of the thermal beam.
Fig. 5. (a) The general pattern of convection flows in system consisting of two layers for Marangoni numbers Mn1 = 30 and Mn2 = 5 at the time moment corresponding to the rupture of the upper layer t* = 7.82. (b) Location of the interfaces f0 and f1 in area of acting of the thermal beam.
Fig. 6. (a) The general pattern of convection flows in system consisting of two layers for Marangoni numbers Mn1 = 5 and Mn2 = 30 at the time moment corresponding to the rupture of the upper layer t* = 8.87. (b) Location of the interfaces f0 and f1 in area of acting of the thermal beam.
L/h3 = 90, q3 = 1200 kg/m3, m3 ¼ 1:2 106 m2 =s, Pr3 = 15, Ca3 = 0.015, Mn3 = 15. the ratio between thermal conductivities of the first and second 1 ¼ 0:5; liquids j the ratio between thermal conductivities of the second and 2 ¼ 2. third liquids j Figs. 7 and 8 represent sequentially the appearance and development of flows in all liquid layers, as well as the position of the interfaces. At the beginning of the system’s heating, as would be expected, nothing happens until the beginning of liquid heating G3. Tensions arising at the surface f2 due to the Marangoni effect, generate convection flow in the liquid G3, whereas in the liquids
G2 and G1, as the stream function isolines show, only the shear flows take place, which are caused by deformations of interfaces (Fig. 7). At that, the deformations of all interfaces happens simultaneously, because they are induced by action of hydrodynamic effect, but not the thermal one, i.e. the pressure jump. Here, the Marangoni stresses play a role of start mechanism creating this jump on the interface f2 (see formulas (9), (9a), (12), (12a)). Fig. 8 shows liquid motion pattern close to the steady-state flow in the considered system of liquid layers. The liquid motion G3 within the area of thermal beam is accompanied by two fairly large vortexes, which occupy more than 2/3 of the layer thickness G3 at a standstill. Vortex formations do not allow the surface f2 to expand further to the ends of the domain in x - direction. The surface f2
A.S. Ovcharova / International Journal of Heat and Mass Transfer 114 (2017) 992–1000
999
Fig. 7. The general pattern of convection flows and location of the interfaces f0, f1 and f2 in system consisting of three layers for Marangoni numbers Mn1 = 0, Mn2 = 0 and Mn3 = 15 at the time moment t* = 8.7.
Fig. 8. The general pattern of convection flows and location of the interfaces f0, f1 and f2 in system consisting of three layers for Marangoni numbers Mn1 = 0, Mn2 = 0 and Mn3 = 15 at the time moment t* = 130.
was constructed which is based on the decomposition of complex problem into unitary elements (modules) with some set of rules allowing their linkage with each other. The liquid flow and heat transfer in each of these modules is described by the system of Navier-Stokes and thermal conductivity equations. The boundary conditions needed to solve these systems are presented explicitly and written in the form of conservation laws. We conducted numerical analysis of the heat load effect on the liquid layers flow characteristics depending on the Marangoni number. It is shown that physical parameters and thermophysical properties of liquid layers play the decisive role in the process, since exactly these properties control the convective flows in the layers as well as determine the position and shape of interfaces.
Conflict of interest Fig. 9. The temporal change of normal velocities of all interfaces fj (x, t), (j = 0,2) in the point x = L/2. Here, 1 (2,3) is the graphics of the normal velocities of interfaces f0 (L/2, t), f1 (L/2,t), f2 (L/2,t), respectively.
cannot move down as well, because of the established thermal equilibrium. In this case the reason for any significant changes in position of interfaces f0 and f1 disappears. This situation is reflected in Fig. 9. This figure presents the temporal change of normal velocities of all interfaces fj (x, t), (j = 0,2) in the point x = L/2. The figure shows that all 3 velocities tend to zero, at that, the normal velocity of the surface f2 tends to zero faster than the other two velocities. We would like to draw attention to the fact, as the viscosities of the liquids control the shape of the interfaces and their normal velocities. 5. Conclusions We have considered the development of convective flows occurring in the complex system of thin liquid layers under the effect of thermal load. To study this process, a mathematical model
The authors declared that there is no conflict of interest.
References [1] L.G. Napolitano, Thermodynamics and dynamics of surface phases, Acta Astronaut. 6 (9) (1979) 1093–1112. [2] V.K. Andreev, YuA. Gaponenko, O.N. Goncharova, V.V. Pukhnachev, Mathematical Models of Convection, De Gryuter Publ, Berlin-Boston, 2012. [3] A.S. Ovcharova, Rupture of liquid film, placed over deep fluid, under action of thermal load, Int. J. Heat Mass Transf. 78 (2014) 294–301. [4] A.F. Voevodin, S.M. Shugrin, Methods of Solution of One-dimensional Evolutional Systems, Nauka, Novosibirsk, 1993 (in Russian). [5] N.E. Kochin, I.A. Kibel, N.V. Rose, Theor. Hydrodyn., vol. 2, Fizmatlit, Moscow, 1963 [in Russian]. [6] L.M. Milne-Thompson, Theoretical Hydrodynamics, London Macmillan and Co. LTD, New York, 1960. [7] N.K. Adam, The Physics and Chemistry of Surfaces, Oxford University Press, Humphrey Milford, London, 1941. [8] V.A. Volkov, Colloid Chemistry, Moscow State Textile University, Moscow, 2001 (in Russian). [9] A.W. Adamson, A.P. Gast, Physical Chemistry of Surfaces, A Wiley-Interscience Publication John Wiley and Sons, Inc., New York, Chichester, Weinheim, Brisbane, Singapore, Toronto, 1999.
1000
A.S. Ovcharova / International Journal of Heat and Mass Transfer 114 (2017) 992–1000
[10] O.N. Goncharova, O.A. Kabov, V.V. Pukhnachov, Solutions of special type describing the three dimensional thermocapillary flows with an interface, Int. J. Heat Mass Transf. 55 (2012) 715–725. [11] R. Kushnir, V. Segal, A. Ullmann, N. Brauner, Inclined two-layered stratified channel flows: Long wave stability analysis of multiple solution regions, Int. J. Multiph. Flow 62 (2014) 17–29. [12] P.J. Roache, Computational Fluid Dynamics, Hermosa Publishers Albuquerque, New Mexico, 1976. [13] A.S. Ovcharova, Leveling a capillary ridge generated by substrate geometry, Comput. Math. Math. Phys. 46 (2) (2006) 305–314.
[14] V.V. Kuznetsov, M.V. Bartashevich, O.A. Kabov, Interfacial balance equations for diffusion evaporation and exact solution for weightless drop, Microgravity Sci. Technol. 24 (2012) 17–31, http://dx.doi.org/10.1007/s12217-011-9285-2. [15] Jim Douglas Jr., H.H. Rachford Jr., On the numerical solution of heat conduction problems in two and three space variables, Trans. Am. Math. Soc. 82 (2) (1956) 421–439. [16] N.N. Yanenko, The Method of Fractional Steps (The Solution of Problems of Mathematical Physics in Several Variables), Springer Verlag, Berlin, Heidelberg, New York, 1971 (Engl. transl. eg. by M. Holt).
Contents lists available at ScienceDirect
International Journal of Heat and Mass Transfer journal homepage: www.elsevier.com/locate/ijhmt
Multilayer system of films heated from above A.S. Ovcharova Lavrentyev Institute of Hydrodynamics, Siberian Branch, Russian Academy of Sciences, 15, Acad. Lavrentyev Avenue, Novosibirsk 630090, Russia
a r t i c l e
i n f o
Article history: Received 11 February 2017 Received in revised form 15 June 2017 Accepted 27 June 2017
Keywords: Multilayer system Viscous nonisothermal fluid Interfaces Marangoni effect
a b s t r a c t The paper deals with the development of convective flows in a complex system of thin liquid layers occurring under the effect of thermal load. To study this process a mathematical model which is based on the decomposition of complex problem into unitary elements (modules) with some set of rules allowing their linkage with each other was constructed. Each module represents a monotypic local model. The fluid flow and heat transfer in each of these modules are described by the Navier-Stokes and thermal conductivity equations. The boundary conditions required to solve these systems of equations are represented explicitly and written in the form of conservation laws. The numerical analysis of effect of thermal load on the characteristics of the liquid layers movement depending on the Marangoni number is carried out. It is shown that physical and thermophysical properties of liquid layers play the decisive role. It is exactly they control the convective flows in layers and determine the position and shape of interfaces. The results of model problem solutions are presented. Ó 2017 Elsevier Ltd. All rights reserved.
1. Introduction Study of convective flows and interfaces in a complex system of thin liquid layers is very important due to the significance of problems, which are related to a broad industrial application of thin liquid layers (films). The continuous development of science-intensive technologies significantly expands the range of problems associated with mathematical modeling of the convection processes of viscous heat-conducting liquids. In this regard, the number of publications and even reviews on this subject is so great (and it is increasing every year) that their list is endless. Depending on the set goal, the authors use different approaches to both the description of the processes occurring in thin layers of viscous liquid, and the numerical solution of the arising tasks. Here, the correctness of the selection of an appropriate model describing the phenomenon under study is, of course, of decisive importance. Mathematical models for describing the processes of mass and heat transfer require correct conditions at the interfaces of liquids. Especially in the case of convective currents of two or more liquid mediums contacting at deformable interfaces. One of the difficulties in solving problems with deformable boundaries is that along with the unknown flow characteristics, it is also necessary to determine the position of the interfaces which change during the motion.
E-mail address: [email protected] http://dx.doi.org/10.1016/j.ijheatmasstransfer.2017.06.123 0017-9310/Ó 2017 Elsevier Ltd. All rights reserved.
System-based derivation of conditions at the interface between two immiscible liquids under nonisothermal unstable flow was first given in [1]. The modification of these considerations is given in [2]. The purpose of this work, which can be considered as a continuation of the work [3], is the construction of a mathematical model for studying the processes arising in a system consisting of any finite number of thin liquid layers with deformable interfaces. This mathematical model is formulated on the basis of known hypotheses about the physical processes provided by the implementation of conservation laws. It allows us to investigate a wide range of problems associated with the impact of temperature or surfactants on a multilayer system. All the boundary conditions for solving local tasks in each liquid layer are presented in explicit form. The structure of the model allows it to expand by adding to it new data on the physical processes ensured the implementation of conservation laws. The idea of constructing a model is based on the decomposition of complex problem into unitary elements (modules) with some set of rules allowing their linkage with each other [4]. Each module represents a monotypic local model, within which the calculation is carried out independently of the other modules of the system, while the ‘‘set of rules” defining, for example, the common boundaries between the modules and the corresponding functional relations, allow connecting these modules back to the original system. In this case the decisive role is played by physical parameters and thermophysical properties of liquid layers, which control the convective flows in the layers and determine the position and shape of interfaces.
A.S. Ovcharova / International Journal of Heat and Mass Transfer 114 (2017) 992–1000
993
Nomenclature t x,y vs,vn D P ! ! s; n d hi/h0 h0
v0
p0 g0
time Cartesian coordinates tangent and normal velocities of points lying on interface stress tensor pressure tangent and normal vectors width of thermal beam dimensionless thickness of layers length scale velocity scale pressure scale gravity scale
Dimensionless parameters Pri Prandtl number in i-th layer Capillary number Cai Mni Marangoni number Cri Crispation number G Galileo number Greek symbols w stream function x vorticity h temperature m kinematic viscosity q density
@wi ; @y
vi ¼
@wi ; @x
@ v i @ui : @x @y
2. Mathematical model
ui ¼
Consider a multilayer system consisting of thermallyconductive immiscible liquid layers as shown in Fig. 1. Here Gi are liquid layers, fj (t,x) are the interfaces between them. The bottom layer is lying on a solid substrate f J ðt; xÞ, the top layer is in contact with the air through the interface f0 (t,x); at the left and right the system of layers is bounded by two solid planes x = 0 and x = L. Each liquid layer Gi is characterized by its thickness hi, density qi, kinematic viscosity mi and surface tension coefficient ri ðTÞ. Moreover, the evaporation and condensation processes are not taken into account in this model; qi and mi are considered to be constant. Such a complex system of immiscible liquid layers can be represented as a flat ribbon graph consisting of two different types of elements: liquid layers Gi, (i = 1, . . ., I) and the interfaces between them fj (j = 1, . . ., J 1), which we will call the internal interfaces. For the unequivocal description of the graph we must determine the order of the conformity of its elements. With this purpose, we enumerate all the interfaces and all liquid layers, for example, as is done in Fig. 1. Selected numbering scheme must be fulfilled during the whole process of solving the problem. According to the chosen numbering scheme, just one liquid layer with the number i + 1 adjoins to internal interface from the bottom, while the layer with the number i adjoins from the top of the interface. Each liquid layer Gi is limited from above by the interface fj (t,x), while from below – by the interface fj+1 (t,x). Surfaces x = 0, x = L, f0 (t,x) and fJ (t,x) have certain features and are classified as a special type surfaces. The orientation of the graph in the plane is presented in Fig. 1.
The characteristic values (scale) of length x0, velocity v0 and pressure p0 should be the same for each layer of the system of layers. Here Rei = v0 hi/mi and Pri = mi/vi, where vi is the thermal diffusivity coefficient, are the Reynolds and the Prandtl numbers in i-th layer of the system, respectively. hi ¼ ðT i T 0 Þ=dT, where T 0 is the characteristic value of the temperature; dT is the characteristic temperature difference for the entire system of layers. In addition, each liquid layer is characterized by its surface tension coefficient ri ðTÞ at the interface between the liquid and air. We assume that the surface tension coefficient of each liquid is a linear function of temperature
2.1. Consider the first type of graph elements – liquid layers
xi ¼
ðiÞ ri ðTÞ ¼ r0ðiÞ ð1 rTðiÞ ðT i T 0 ÞÞ; rðiÞ 0 ¼ ri ðT 0 Þ; rT
¼
1 dri rðiÞ dT 0
ðiÞ
; rT > 0:
T¼T 0
i.e. it decreases linearly with the increase of temperature. 2.2. Interfaces. functional relations Interfaces are the second type of graph elements. Consider any internal interface between the two layers Gi and Gi+1. As we consider any internal surface, we can define it by the equation y = f (t,x). It is however important to keep in mind that the liquid with number i + 1 lies below this interface, while the liquid with the number i is over it (see Fig. 2). It is related to the fact that when constructing the model we predetermined the order of conformity of graph elements. Define the unit tangent and normal vectors at the interface f (t, x) of liquids as
The fluid flow and heat transfer in each of the liquid layer Gi, (i = 1, . . ., I) is described by the system of Navier-Stokes and thermal conductivity equations, which in the variables w, x, h (stream function, vorticity and temperature) have the form [5,6]
8 9 > > < 1 fx = ! ; s ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffi ; qffiffiffiffiffiffiffiffiffiffiffiffiffi > 2> : 1 þ f2 1 þ fx; x
8 9 > > < f = 1 ! x n ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffi ; qffiffiffiffiffiffiffiffiffiffiffiffiffi > 2> : 1 þ f2 1 þ fx; x
@ xi @ @w @ @w 1 xi i x i i ¼ Dx i ; þ @x @y Rei @t @y @x
ð1Þ
! Note that for liquid Gi+1, n is the outer normal unit vector at the interface of liquids, while for the liquid Gi for the same interface this normal unit vector is directed into the liquid.
Dwi þ xi ¼ 0;
ð2Þ
@hi @ @w @ @w 1 hi i hi i ¼ Dh i þ @y Rei Pri @t @x @y @x
ð3Þ
2.2.1. Functional relations for temperature In this paper we do not consider the evaporation (condensation) process at the liquid interface and neglect energy consumption related to its deformation (since they are negligible). Therefore,
994
A.S. Ovcharova / International Journal of Heat and Mass Transfer 114 (2017) 992–1000
Fig. 1. Sketch showing the agreement between elements of the film system placed on a solid substrate y = 0. Here, Gi are the layers, fj are the interfaces of the system; d is the width of heat beam effecting on interface f0.
the temperature conditions at the interface f (t,x) can be written in the form of temperature and heat flux continuity conditions.
hiþ1 hi ¼ 0;
@hiþ1 @h i ¼ 0; j @n @n
ð4Þ
¼ ji =jiþ1 is the ratio between heat conductivities of (i)-th where j and (i + 1)-th liquids, respectively. 2.2.2. Functional relations for the stream function From the equality of the tangential velocities of the liquids Gi and Gi+1 at their common interface f (t,x) and the volume conservation condition for each of these liquids, it follows
wiþ1 wi ¼ 0;
@wiþ1 @wi ¼ 0: @n @n
ð5Þ
2.2.3. Functional relations for the vorticity For derivation of these relations we use the kinematic and dynamic conditions at the interface [1,2]. The dynamic condition is based on the preservation of the continuity of the normal and tangent stresses on the interface, which is fulfilled if the mediums preserve their continuity. If we neglect the evaporation (condensation) process, the dynamic condition at the liquid interface f (t,x) can be represented as ! r ðTÞ ! n þ rr ðTÞ; ðDiþ1 Di Þ n ¼ R
ð6Þ
where Di, Di+1 are the stress tensors in liquids Gi and Gi+1, respectively, r ðTÞ is the interfacial tension coefficient between these liquids, R is the radius of the liquid interface curvature, rr ðTÞ is the interfacial tension gradient along the interface. Thus, it only remains to consider what is meant by the function r ðTÞ and, most importantly, how changes its distribution along the interface depending on temperature. There are several approaches for defining this function depending on physical, thermophysical and even
thermodynamic properties of liquids (for example, the polarity of the molecules). We can apply the commonly used Antonow’s rule [7], which reads as follows: interfacial tension at the interface of two immiscible liquids is equal to the difference between surface tensions of these liquids (at the interface with air or intrinsic vapor) in the conditions of mutual saturation. That is, just in our case
r ðTÞ ¼ riþ1 ðTÞ ri ðTÞ
ð7Þ
This rule is valid for many pairs of liquids, though many exceptions to this rule are known as well. It may also happen that a liquid Gi in contact with a liquid Gi+1 satisfies Antonow’s rule, while for the same liquid Gi in contact with another liquid Gi-1 this rule does not work. There are also other formulas to determine the interfacial tension between the liquids [8,9]. However, they are all based on the use of surface tension coefficients of these liquids at their contact with air or intrinsic vapor (which are assumed to be known in advance) together with empirical coefficients which are inherent of the considered liquids. In this case, we can always determine not only interfacial tension between liquids, but also its distribution along the interface depending on temperature. When solving the problem with a specific combination of liquids, the researcher determines, which formula would provide most reliable solution. In this work, when describing the model, we use Antonow’s rule. From dynamic conditions we can obtain two corollary facts. ! Carrying out scalar multiplication of (6) by vector s we get
! ! @ r ðTÞ s ðDiþ1 Di Þ n ¼ @s " # ! 2 1 f x @ 2 wiþ1 @ 2 wiþ1 @ 2 wiþ1 2f x qiþ1 miþ1 2 @y2 @x2 @y@x ! " # 2 2 2 2 1 f x @ wi @ wi @ wi qi mi 2f x 2 @y2 @x2 @y@x 2
¼
1 þ f x @ r ðTÞ : @s 2
ð8Þ
Using the kinematic condition for each of the liquids at their common interface as well as Antonow’s rule (7), the condition (8) can be reduced to the following form:
xiþ1 q mxi ¼ F 1 ðt; xÞ; where
Mni Þ m F 1 ¼ ðMniþ1 q
@h @v n v s Þ m þ 2ð1 q þ @s @s R
ð9Þ
ð9aÞ
! Carrying out scalar multiplication of (6) by vector n , we obtain ! ! r ðTÞ n ðDiþ1 Di Þ n ¼ R
Fig. 2. Layout of the layers and interface between them.
A.S. Ovcharova / International Journal of Heat and Mass Transfer 114 (2017) 992–1000
or, to give a more convenient form,
@v n f vs ðPiþ1 Pi Þ ¼ 2qiþ1 miþ1 x @n iþ1 R @v n f xv s r ðTÞ þ þ 2qi mi R @n i R
ð10Þ
Next, we consider the vector equation of motion for each of the liquids, written in the natural variables, i.e. the velocity – the pressure.
! @vi 1 ! ! ! ! þ ð v i rÞ v i ¼ rPi þ mi D v i þ g @t qi
ð11Þ
We are looking for such a solution, where vs at the interface of liquids must comply with the conditions (5) (equality of the tangential velocities of points situated at the interface). This is requirement of the physical model describing the problem under consideration. At the same time, the dynamic condition has its requirement – a pressure jump at the interface should satisfy the balance relations (6). If we take the scalar product of vector equa! tion (11) for i with the vector qi s ðt; xÞ and the equation for i + 1 ! with vector qiþ1 s ðt; xÞ and put them together, we will get the equation
ðqiþ1 qi Þ
@v s @v s @ @ þ ðqiþ1 qi Þv s ¼ ðPiþ1 Pi Þ @s @n @t @s ðqiþ1 miþ1 xiþ1 qi mi xi Þ ! ! þ ðqiþ1 qi Þ g s þ ðqiþ1 qi Þ
in which a pressure jump (Pi+1 Pi) must satisfy (6). Eliminating @ðP iþ1 Pi Þ=@s with the use of (10), after conversion to dimensionless form, we get
@ xiþ1 @ xi m q ¼ F 2 ðt; xÞ; @n @n where
ð12Þ
@ @v n @ @v n @ f xvs Þ m m þ 2 ð1 q q @s @n iþ1 @s @n i @s R 1 @ 1 ð1 Mniþ1 Caiþ1 h þ Caiþ1 @s R @ 1 Ca1 m ð1 Mn Ca h q i i i @s R 2
F 2 ðt; xÞ ¼ 2
@v s @v fx 6 1Þv s s þ ðq 1Þ G qffiffiffiffiffiffiffiffiffiffiffiffiffi 1Þ þ Reiþ1 4ðq þ ðq @t @s 2 1 þ fx # f v2 xi Þ Þ x n þ v n ðxiþ1 q ð12aÞ þð1 q R 1 f xx ; ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi R 2 3 ð1 þ f x Þ
Expression (12a) looks at first glance bulky, in fact it is quite a simple expression. Meaningful part of it consists in the first two lines, reflecting the contribution of the pressure jump. The last line is a tribute to the nonstationarity of the problem. To obtain this expression we also used Antonow’s rule (7). Recall that the relations of type (9) with the right hand part of (9a), as well as (12) with the right hand part of (12a) are valid for any internal interface between a couple of layers Gi and Gi+1, which are included to the system and have common internal interface. The position of the interface f (x,t) is determined from the kinematic condition at the interface, for example, for the lower layer. In the present case this is
ft þ
qffiffiffiffiffiffiffiffiffiffiffiffiffi 2 @wiþ1 1 þ fx ¼0 @s
q v 0 mi rðiÞ dT g h0 Cai ¼ i ðiÞ ; Mni ¼ T ; G ¼ 02 : qi v 0 mi v0 r0
Here vs, vn are tangential and normal velocity components of points lying in the interface f (x,t); R is the radius of surface curvature f (x,t); Cai is the capillary number; Mni is the Marangoni num ¼ qi =qiþ1 is the density ratio, ber; G is the Galileo number; q m ¼ mi =miþ1 is the viscosity ratio of liquids i and i + 1, respectively. plays an important role in formation of m The last parameter q the convective flows in the layers [10], especially in the cases where there is no thermal load and the convection currents are caused by other factors [11].
ð13Þ
2.3. The boundary conditions at the remaining surfaces of the system The boundary conditions for the temperature at the other boundaries take form in accordance with the requirements of the particular tasks (see numerical calculations section). At x = 0, x = L, y = fJ (x), i.e. at solid surfaces, we use the no-slip boundary conditions
wi ¼ 0;
f x v 2n þ v n ðqiþ1 xiþ1 qi xi Þ R
995
@wi ¼ 0: @n
Here we have two conditions for the stream function w while none of conditions for the vorticity x. To find the vorticity value x on a solid wall, we use the finite-difference conditions of the Thom’s type [12]. At the liquid-gas interface f0 (t,x) we neglect not only the evaporation (condensation) process but also the effects of the dynamic characteristics of the gas to the liquid motion. We can’t use functional relations for w and x, obtained for internal interfaces, as in this case we will have two conditions for x and no one condition for stream function w. Therefore, it makes sense to distinguish liquid-gas interface and write down boundary conditions on this interface according to its features. Again, we denote the interface f0 (t,x) just as f (t,x) and have the same formulations as in the case of the internal interface with just one exception. A second corollary fact of balance ratios (6) allows representation
@v n f xv s @v s v n ! ! r ðTÞ ! ¼ 2qm P ¼ 2qm nDn ¼ R @n R @s R In the derivation of functional dependencies at the internal interface we used the first representation for P
@v n f xv s P ¼ 2qm @n R Whereas to obtain the boundary conditions at the liquid-gas interface we used a second representation for P. A detailed derivation of the boundary conditions at the liquid-gas interface is available in [3], here we present only the result of this work: The position of the surface f (t,x) is determined using the kinematic condition
ft þ
qffiffiffiffiffiffiffiffiffiffiffiffiffi 2 @w 1 þ fx ¼0 @s
ð14Þ
The boundary conditions take the form:
@w ¼ v s; @n
ð15Þ
996
A.S. Ovcharova / International Journal of Heat and Mass Transfer 114 (2017) 992–1000
v @v @h x ¼ 2 s þ n þ Mn ; @s R @s where
ð16Þ
vs is the solution to the equation
Bi ¼ 1;
ðiÞ
Ai ¼ 1;
R1 ¼
SðiÞ t g þ ðf jþ1 Þt
@v s @v s @ vs þ vs ¼ 2 2 þ F; @t @s @s @x @ 1 rT dT fx @ v n F¼ 2 hÞ G qffiffiffiffiffiffiffiffiffiffiffiffiffi ð1 þ Ca1 @s R @s R @n r0 2 1 þ fx
v
ð17Þ
f
ð18Þ
2
1 þ fx
3. Numerical method Suppose we have a system consisting of I liquid layers Gi, (i = 1, . . ., I) and J + 1 interfaces fj, (j = 0, . . ., J). Each layer Gi is limited at the top by the surface fj (t,x), while at the bottom – by the surface fj+1 (t,x) as shown in Fig. 1. Map each of the area occupied by the liquid Gi in a rectangle with sides of 0 6 n 6 L, 0 6 g 6 1 using the transformation:
SðiÞ ¼ f j ðt; xÞ f jþ1 ðt; xÞ
In this case, all boundaries of each layer in the system will rest on the coordinate lines. Every Eqs. (1)–(3) for each layer i can be written in the following form
h i ðiÞ ðiÞ B12 ¼ Sn g þ ðf jþ1 Þn ;
For Ui ¼ xi
ðiÞ
B22 ¼
ð19Þ
ðiÞ2
1 þ B12 SðiÞ
R1 ¼
StðiÞ g þ ðf jþ1 Þt SðiÞ
;
ðiÞ
R2 ¼ 0:
ðiÞ
Ai ¼ 0;
R1 ¼ 0;
ðiÞ
R2 ¼ ki xi ;
where ki is the iteration parameter to solve Poisson’s equation for wi . If we take new designations:
ðiÞ
VðUi Þ ¼ B12
@ Ui @w ðiÞ @ Ui Ai Ui i þ B12 @n @g @g @ Ui @w ðiÞ @ Ui þ Ai Ui i þ B22 @n @g @n
then Eq. (19) will be
@ Ui 1
ðiÞ @ Ui ðiÞ U n ðUi Þ þ V g ðUi Þ þ R1 þ R2 ¼ @t @g Bi SðiÞ
ð21Þ
Eq. (21) for hi (temperature), xi (vorticity) and wi (stream function) has been solved by the stabilizing-correction finite-difference scheme [15,16] written in the form
Ukþ1=2 Uki i
s
¼
Ukþ1 Ukþ1=2 i i
s
1 h Bi S ¼
ðiÞ
kþ1=2 i ðiÞ @ Ui ðiÞk U kn ðUi Þ þ V kþ1=2 ðUi Þ þ R1 þ R2 g @g
1 h Bi S
ðiÞ
i k U kþ1 n ðUi Þ U n ðUi Þ
This scheme is one of the most efficient finite-difference schemes which uses fractional steps. The first fractional step gives the full approximation of the equation, while the second step is a correctional one, and is used with the aim of improving stability. The transition to new variables in the boundary conditions and finite-difference technique when solving the problem for a single liquid layer is described in detail in [3]. Here we only point out that because of connectedness of the boundary conditions at the internal interfaces, numerical solution in the y - direction should be carried out consistently, for example, moving from the bottom layer GI to the top layer G1. Next, we consider the general scheme of solution for the system of layers. Let the positions of the interfaces fj (t,x) as well as distribution of temperature hi , vorticity xi , and the stream function wi are known for time moment t k ¼ ks (s is the time step) for all liquid layers, constituting the system. Stage 1. Using Eqs. (13) and (14), we find the new position of all interfaces f j ðtkþ1 ; xÞ for time moment tkþ1 ¼ ðk þ 1Þs. Define normal velocities v n of points at these interfaces. ðiÞ
Here ðiÞ
1 ; ki
ðiÞ
@ Ui 1 @ @w ðiÞ @ Ui ðiÞ @ Ui B11 Ai Ui i ¼ þ B12 ðiÞ @t @n @g @g Bi S @n @ @ U @ U @w i i ðiÞ ðiÞ ðiÞ @ Ui ðiÞ B þ Ai Ui i þ R1 þ R2 þ B22 þ @ g 12 @n @g @n @g
ðiÞ
Ai ¼ Pri ;
UðUi Þ ¼ B11
Recall that Eqs. (14)–(18) with all the motion characteristics and the governing parameters are relevant only to the liquid-gas interface f (t,x), i.e. where we neglect the effect of the dynamic characteristics of the gas. Let summarize the previous constructions. To obtain quantitative characteristics of liquid motion and temperature distribution in the considered system of layers, we must solve Eqs. (1)–(3) for each layer. The boundary conditions required for numerical solution of these equations and determination of the interfaces are presented in an explicit form. These are conditions (4), (5), (9), (9a), (12), (12a), and (13). The boundary conditions for the liquid-gas interface, where we neglect the effect of the dynamic characteristics of the gas, are classified as a special type. These are conditions (14)–(18). It should be noted, that closing relations (i.e. boundary conditions) in this model are not only derived from conservation laws, but written in the form of the conservation laws in accordance with the assumptions and limitations, which underlie the model. Besides, these conditions are valid for any finite number of liquid layers that make up the system. In conclusion, we note that the model can be adapted for the solution of a wide range of tasks, such as calculations with variable viscosity of liquids [13], the calculations connected with the evaporation process from one layer into another one [14], etc.
B11 ¼ SðiÞ ;
Bi ¼ Pri ;
Bi ¼
R
y ¼ gSðiÞ þ f jþ1 ðt; xÞ;
ðiÞ
R2 ¼ 0:
For Ui ¼ wi
2 n
t v n ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffi
x ¼ n;
;
For Ui ¼ hi
2
þ v nx þ f x
SðiÞ
ð20Þ
ðiÞ
ðiÞ
Stage 2. Using Eq. (20) determine the coefficients B11 , B12 , B22 required to the solution of the motion and heat conduction equations in each of the system’s layers. Stage 3. At this important stage we find the boundary conditions for the solution of the motion equations: using Eqs. (9a) and (12a) define functions F 1 ðn; gÞ and F 2 ðn; gÞ, which are included in the right hand part of the boundary conditions (9) and (12); solving Eq. (17), we find the tangential velocity vs of points lying at the surface f0 (t,x) that allows finding the boundary conditions (15) and (16) at this boundary. The importance of this
997
A.S. Ovcharova / International Journal of Heat and Mass Transfer 114 (2017) 992–1000
stage is that it allows finding all boundary conditions required for solution of all equations for all liquid layers of the system. Stage 4. Finally, solving Eq. (21) with the new coefficients ðiÞ
ðiÞ
ðiÞ
B11 ðn; gÞ, B12 ðn; gÞ, B22 ðn; gÞ consequently for temperature hi , vorticity xi , and stream function wi for all layers of the system, we obtain a new distribution of these functions at the new time step for the iteration l = 0. The iteration process is necessary due to the nonlinearity of Eq. (21). Iteration continues until the following condition
lþ1 l Ui ðn; gÞ Ui ðn; gÞ
x1 ¼ L=2 mh0 ;
x2 ¼ L=2 þ mh0
We consider the concentrated thermal load acting on film free surface of such type:
(
hðt; xÞ ¼
h ðxÞ; if x2 6 x 6 x3 @h ; @n
otherwise
Here, L is the length of all layers of the system; m is positive number, which determines a width of the heat beam affecting the film free surface. At solid walls (x = 0, x = L, y = 0) the temperature h = 0. In all calculations m = 3, the width of the thermal beam is d = 6h0, h⁄ = 1. In problems of this type there are a large number of governing parameters, even if considering just one layer. When considering a larger number of liquid layers, the number of governing parameters increases multiple times. We focused on the study of the influence of Marangoni effect on the behavior of the system of liquid
layers and interfaces between them. In fact, exactly this effect causes the convection motions in the layers and, as the consequence, the deformation of interfaces. Figs. 3–6 present the calculation results, carried out for the system consisting of two thin layers. Liquid G1 is situated between the surfaces f1 and f0, while liquid G2 is confined by the surfaces f2 and f1. Layers thickness, governing parameters and physical properties of liquids are chosen in such a way as to demonstrate essential deformation of interfaces up to the rupture of the layer. The calculations are performed for following parameters: L/h1 = 90, q1 = 780 kg/m3, Ca1 = 0.025
m1 ¼ 1:2 106 m2 =s,
Pr1 = 16,
L/h2 = 90, q2 = 1980 kg/m3, m2 ¼ 1 106 m2 =s, Pr2 = 7, Ca2 = 0.02 the ratio of the thermal conductivities of the first and second ¼ 0:5. liquids is j Fig. 3a presents the general pattern of convection motions in layers at Mn1 = 15 and Mn2 = 5 at the time moment corresponding to the rupture of the upper layer. Convection flow in the lower layer is accompanied by the two vortexes, which constrain the motion of the lower liquid G2 within the domain of impact of thermal beam. The liquid at the top is bordered by air, whose dynamic characteristics are neglected in this work. Because Mn1 > Mn2, the rupture occurs at the lateral parts of the interface f1, rather than in its lowest points (Fig. 3b), because as a result of heating, the surface tension of the liquid in upper layer changes significantly faster than the interfacial tension between liquids (see formulas (9), (9a) and (12), (12a)). Consider the situation where Mn1 = 5 and Mn2 = 15, while all the other characteristics of both liquids remain unchanged. In this case, the convective motion of the lower liquid G2 also will be accompanied by two vortex formations, restricting this motion within certain limits. The intensity of convective liquid motion G2 will rise slightly, though it rises essentially in the upper liquid G1 (Fig. 4a). Convection motions, caused by the action of the Marangoni effect along the surface f0 are superposed by the perturbations caused by the same effect at the surface f1. Together, they significantly increase the intensity of convective flows in liquid G1. The surface f0 cannot contact the lateral walls of the depression, as it was shown in Fig. 3b, because the Marangoni number of the liquid G1 is much less than the difference between the Marangoni numbers of these liquids. Therefore, the rapture of the upper layer occurs at the periphery d of the width of the thermal beam affecting the surface f0. At that, ‘‘corners,” formed by the surfaces f0 and f1 destroy the vortex formations in the liquid G2. We have carried out numerical experiments with an even greater difference between the Marangoni numbers of both liquids, while all other physical characteristics remained the same. The results of these studies are presented in Fig. 5a and b (Mn1 = 30, Mn2 = 5) and Fig. 6a and b (Mn1 = 5, Mn2 = 30). As we see, this resulted in change of just the quantitative characteristics of liquid flow, while the general pattern of the liquid motion and the nature of rupture of the top liquid layer remain unchanged. We have also considered the case where the Marangoni number in the top liquid layer Mn1 = 0. This means that the liquid surface tension does not depend on temperature. For greater effect, between two layers of liquids we placed another layer with the same property Mn2 = 0. The calculations were performed for the following parameters: L/h1 = 180, q1 = 780 kg/m3, Ca1 = 0.025, Mn1 = 0,
m1 ¼ 2 106 m2 =s,
Pr1 = 16,
L/h2 = 90, q2 = 998 kg/m3, Ca2 = 0.0125, Mn2 = 0,
m2 ¼ 1 106 m2 =s,
Pr2 = 7,
998
A.S. Ovcharova / International Journal of Heat and Mass Transfer 114 (2017) 992–1000
Fig. 3. (a) The general pattern of convection flows in system consisting of two layers for Marangoni numbers Mn1 = 15 and Mn2 = 5 at the time moment corresponding to the rupture of the upper layer t* = 22.03. (b) Location of the interfaces f0 and f1 in area of acting of the thermal beam.
Fig. 4. (a) The general pattern of convection flows in system consisting of two layers for Marangoni numbers Mn1 = 5 and Mn2 = 15 at the time moment corresponding to the rupture of the upper layer t* = 15.58. (b) Location of the interfaces f0 and f1 in area of acting of the thermal beam.
Fig. 5. (a) The general pattern of convection flows in system consisting of two layers for Marangoni numbers Mn1 = 30 and Mn2 = 5 at the time moment corresponding to the rupture of the upper layer t* = 7.82. (b) Location of the interfaces f0 and f1 in area of acting of the thermal beam.
Fig. 6. (a) The general pattern of convection flows in system consisting of two layers for Marangoni numbers Mn1 = 5 and Mn2 = 30 at the time moment corresponding to the rupture of the upper layer t* = 8.87. (b) Location of the interfaces f0 and f1 in area of acting of the thermal beam.
L/h3 = 90, q3 = 1200 kg/m3, m3 ¼ 1:2 106 m2 =s, Pr3 = 15, Ca3 = 0.015, Mn3 = 15. the ratio between thermal conductivities of the first and second 1 ¼ 0:5; liquids j the ratio between thermal conductivities of the second and 2 ¼ 2. third liquids j Figs. 7 and 8 represent sequentially the appearance and development of flows in all liquid layers, as well as the position of the interfaces. At the beginning of the system’s heating, as would be expected, nothing happens until the beginning of liquid heating G3. Tensions arising at the surface f2 due to the Marangoni effect, generate convection flow in the liquid G3, whereas in the liquids
G2 and G1, as the stream function isolines show, only the shear flows take place, which are caused by deformations of interfaces (Fig. 7). At that, the deformations of all interfaces happens simultaneously, because they are induced by action of hydrodynamic effect, but not the thermal one, i.e. the pressure jump. Here, the Marangoni stresses play a role of start mechanism creating this jump on the interface f2 (see formulas (9), (9a), (12), (12a)). Fig. 8 shows liquid motion pattern close to the steady-state flow in the considered system of liquid layers. The liquid motion G3 within the area of thermal beam is accompanied by two fairly large vortexes, which occupy more than 2/3 of the layer thickness G3 at a standstill. Vortex formations do not allow the surface f2 to expand further to the ends of the domain in x - direction. The surface f2
A.S. Ovcharova / International Journal of Heat and Mass Transfer 114 (2017) 992–1000
999
Fig. 7. The general pattern of convection flows and location of the interfaces f0, f1 and f2 in system consisting of three layers for Marangoni numbers Mn1 = 0, Mn2 = 0 and Mn3 = 15 at the time moment t* = 8.7.
Fig. 8. The general pattern of convection flows and location of the interfaces f0, f1 and f2 in system consisting of three layers for Marangoni numbers Mn1 = 0, Mn2 = 0 and Mn3 = 15 at the time moment t* = 130.
was constructed which is based on the decomposition of complex problem into unitary elements (modules) with some set of rules allowing their linkage with each other. The liquid flow and heat transfer in each of these modules is described by the system of Navier-Stokes and thermal conductivity equations. The boundary conditions needed to solve these systems are presented explicitly and written in the form of conservation laws. We conducted numerical analysis of the heat load effect on the liquid layers flow characteristics depending on the Marangoni number. It is shown that physical parameters and thermophysical properties of liquid layers play the decisive role in the process, since exactly these properties control the convective flows in the layers as well as determine the position and shape of interfaces.
Conflict of interest Fig. 9. The temporal change of normal velocities of all interfaces fj (x, t), (j = 0,2) in the point x = L/2. Here, 1 (2,3) is the graphics of the normal velocities of interfaces f0 (L/2, t), f1 (L/2,t), f2 (L/2,t), respectively.
cannot move down as well, because of the established thermal equilibrium. In this case the reason for any significant changes in position of interfaces f0 and f1 disappears. This situation is reflected in Fig. 9. This figure presents the temporal change of normal velocities of all interfaces fj (x, t), (j = 0,2) in the point x = L/2. The figure shows that all 3 velocities tend to zero, at that, the normal velocity of the surface f2 tends to zero faster than the other two velocities. We would like to draw attention to the fact, as the viscosities of the liquids control the shape of the interfaces and their normal velocities. 5. Conclusions We have considered the development of convective flows occurring in the complex system of thin liquid layers under the effect of thermal load. To study this process, a mathematical model
The authors declared that there is no conflict of interest.
References [1] L.G. Napolitano, Thermodynamics and dynamics of surface phases, Acta Astronaut. 6 (9) (1979) 1093–1112. [2] V.K. Andreev, YuA. Gaponenko, O.N. Goncharova, V.V. Pukhnachev, Mathematical Models of Convection, De Gryuter Publ, Berlin-Boston, 2012. [3] A.S. Ovcharova, Rupture of liquid film, placed over deep fluid, under action of thermal load, Int. J. Heat Mass Transf. 78 (2014) 294–301. [4] A.F. Voevodin, S.M. Shugrin, Methods of Solution of One-dimensional Evolutional Systems, Nauka, Novosibirsk, 1993 (in Russian). [5] N.E. Kochin, I.A. Kibel, N.V. Rose, Theor. Hydrodyn., vol. 2, Fizmatlit, Moscow, 1963 [in Russian]. [6] L.M. Milne-Thompson, Theoretical Hydrodynamics, London Macmillan and Co. LTD, New York, 1960. [7] N.K. Adam, The Physics and Chemistry of Surfaces, Oxford University Press, Humphrey Milford, London, 1941. [8] V.A. Volkov, Colloid Chemistry, Moscow State Textile University, Moscow, 2001 (in Russian). [9] A.W. Adamson, A.P. Gast, Physical Chemistry of Surfaces, A Wiley-Interscience Publication John Wiley and Sons, Inc., New York, Chichester, Weinheim, Brisbane, Singapore, Toronto, 1999.
1000
A.S. Ovcharova / International Journal of Heat and Mass Transfer 114 (2017) 992–1000
[10] O.N. Goncharova, O.A. Kabov, V.V. Pukhnachov, Solutions of special type describing the three dimensional thermocapillary flows with an interface, Int. J. Heat Mass Transf. 55 (2012) 715–725. [11] R. Kushnir, V. Segal, A. Ullmann, N. Brauner, Inclined two-layered stratified channel flows: Long wave stability analysis of multiple solution regions, Int. J. Multiph. Flow 62 (2014) 17–29. [12] P.J. Roache, Computational Fluid Dynamics, Hermosa Publishers Albuquerque, New Mexico, 1976. [13] A.S. Ovcharova, Leveling a capillary ridge generated by substrate geometry, Comput. Math. Math. Phys. 46 (2) (2006) 305–314.
[14] V.V. Kuznetsov, M.V. Bartashevich, O.A. Kabov, Interfacial balance equations for diffusion evaporation and exact solution for weightless drop, Microgravity Sci. Technol. 24 (2012) 17–31, http://dx.doi.org/10.1007/s12217-011-9285-2. [15] Jim Douglas Jr., H.H. Rachford Jr., On the numerical solution of heat conduction problems in two and three space variables, Trans. Am. Math. Soc. 82 (2) (1956) 421–439. [16] N.N. Yanenko, The Method of Fractional Steps (The Solution of Problems of Mathematical Physics in Several Variables), Springer Verlag, Berlin, Heidelberg, New York, 1971 (Engl. transl. eg. by M. Holt).