ON A CREEPING 3D CONVECTIVE MOTION OF FLUIDSWITH AN ISOTHERMAL INTERFACE Текст научной статьи по специальности «Математика»

Journal of Siberian Federal University. Mathematics & Physics 2020, 13(6), 661-669

DOI: 10.17516/1997-1397-2020-13-6-661-669 УДК 517.977.55:536.25

On a Creeping 3D Convective Motion of Fluids with an Isothermal Interface

Viktor K. Andreev*

Institute of Computational Modelling SB RAS Krasnoyarsk, Russian Federation Siberian Federal University Krasnoyarsk, Russian Federation

Received 22.06.2020, received in revised form 02.07.2020, accepted 20.09.2020 Abstract. In the work the 3D two-layer motion of liquids, the velocity field of which has a special form, is considered. The arising conjugate initial boundary value problem for the Oberbek-Boussinesq model is reduced to a system of ten integrodifferential equations with full conditions on a flat interface. It is shown that for small Marangoni numbers the stationary problem can have up to two solutions. The case when the stationary flow arises due to a change in the internal interphase energy is analyzed separately.

Keywords: Oberbek-Boussinesq model, interphase energy, creeping flow, inverse problem. Citation: V.K. Andreev, On a Creeping 3D Convective Motion of Fluids with an Isothermal Interface, J. Sib. Fed. Univ. Math. Phys., 2020, 13(6), 661-669. DOI: 10.17516/1997-1397-2020-13-6-661-669.

1. Statement of the problem and basic equations

Suppose that two viscous heat-conducting fluids with a common interface 2 = l1 <l2 move in a layer \x\ < to, \y\ < to, 0 < z < l2, lj are constants. The fluid " 1" occupies the region 0 < z <l1 and fluid " 2" occupies the region li < z <l2. The planes z = 0 and z = l2 are solid fixed walls, the force of gravity is directed perpendicular to the layers. Oberbeck-Boussinesq equations are used as a mathematical model of fluid motion. Solutions are sought in a special way

uj(x,y) = (fj(z,t) + hj(z,t))x, Vj(x,y) = {fj(z,t) - hj(z,t))y, Wj = -2 fj(1)

J z 0

—Pj = bj (z,t)x2 + dj (z,t)y2 + qj (z,t), (2)

pj

Tj = aj (z,t)x2 + cj (z,t)y2 + 6j (z,t), (3)

where Uj (x,y,z,t), Vj (x,y,z,t), Wj (x,y,z,t) are projections of velocity vectors on the x, y, z axis, respectively; pj (x,y,z,t) are pressures; pj are constants of density; Tj (x,y,z,t) are absolute temperatures, j = 1, 2. The functions fj, hj, bj, dj, qj, aj, Cj, Qj are new unknown function.

Substitution of the formulas (1)-(3) in the systems of Oberbeck-Boussinesq equations leads to the following systems

fjt + j + h2 - 2fjzi fj (£,t) + gj f (aj (£,t) + Cj (£,t)) = Vj fjzz + nji(t), (4)

z0 z0

*andr@icm.krasn.ru © Siberian Federal University. All rights reserved

hjt + 2fjhj - 2hjZ f fj (Ç, t) dÇ + gfa f (aj(Ç,t) - cj(Ç,t)) dÇ = Vjhjzz + nj2(t), (5) Jz0 Jz0

ajt + 2aj (fj + hj ) - 2ajz fj (Ç, t) dÇ = Xjajzz, (6)

z0

cjt + 2cj (fj - hj ) - 2cjz fj (Ç, t) dÇ = XjCjzz, (7)

z0

Ojt - 29jj fj (Ç, t) dÇ = XjOjzz + 2xj (aj + cj ). (8)

z0

Here Vj > 0, Xj > 0, Pj > 0 are constants of kinematic viscosities, thermal diffusivities and thermal expansion coefficients of liquids; nj1(t), nj2(t) are arbitrary functions of time. By the known functions aj, cj the functions bj, dj are determined by quadratures

bj (z,t) = gj f aj (Ç,t) dÇ - nji (t), dj (z,t) = gpj f cj (Ç,t) dÇ - nj2(t). (9)

z0 z0

In the integral terms, the constant z0 is equal to "0" for the first fluid (j = 1) and l1 for the second fluid (j = 2). It can be verified that pressures in liquids are determined as follows

1

* =

gßj aj(Ç,t) dÇ - nji(t) X2 + gßj Cj(Ç,t) dÇ - nj2(t)

J Zo \ L Zo

y2 - 2vj fj - gz+

jj

(10)

+ gPj i % + 2 f(z - OfAU) d£ + 2((Z f (U) dA + qjo(t),

J z 0 J z 0 \J zo J

with arbitrary functions qjo(t).

Remark 1. The velocity field (1), proposed in [1] is a special case of the velocity field for the Navier-Stokes equations [2].

2. Boundary and initial conditions

On solid boundaries, the sticking conditions for the velocities are satisfied, which implies equalities

W 2

fi(0,t) = h1(0,t)=0, f2(l2,t) = h2(l2,t) = h(i,t) di = 0 (11)

(12)

And the temperature is set

ai(0,t) = a0(t), c1(0,t) = c0(t), 91(0, t) = 91(t), a2(l2,t) = a2 (t), c2(l2,t) = c2(t), 92(l2,t) = 92(t).

The top wall can also be thermally insulated

a2z (l2,t) = c2z (l2,t) = 92z (l2,t) = 0, (13)

To formulate the conditions on the undeformed interface z = li , we assume that the surface tension depends linearly on temperature

a(T)= ao - œ(T - To), (14)

whereao, œ, T0 are given positive constants, T(x,y,li;t) is temperature on this border.

On the interface 2 = li there are equalities of velocities and temperatures. Taking into account the representation (1), (3) we get [3]

fi(li,t) = f2(h,t), h1(l1,t) = h2(h,t), a1(l1,t) = a2(lut), ci(li,t) = C2(h,t), 0i(h,t) = 02(h,t). Tangential stresses are reduced to two relations

№f2z(h,t) - l^ifiz(li,t) = -œ(ai(li,t) + c1(l1,t)), №h2z(li,t) - mhiz(li,t) = -œ(ai(li,t) - ci(li,t)),

(15)

(16)

where Hj = pjVj are dynamic viscosity of liquids.

The kinematic condition for a fixed and non-deformable interface (w1(l1,t) = w2(l1,t) = 0) is equivalent to the integral equality

r fi(^t) d£ = 0. (17)

Jo

The energy condition [3], taking into account the assumptions (8), can be written as

k2T2z(x,y,li,t) - kiTiz(x,y,li,t) = '<eT(x,y,li,t)divru. (18)

where kj are constant coefficients of thermal conductivity of liquids; divru is surface divergence of the velocity vector; T(x, y, l1,t) = T1(x, y, l1,t) = T2(x, y,l1,t). Since in our case divru = ux+vy, then using the formulas (1), (3) from (18) we derive the relations

k2a2z(l1,t) - k1a1z(h,t) = 2&a1(h,t)f1(h,t),

k2C2z (l1,t) - k1C1z (h,t) = 2lbC1(l1,t)f1(l1,t), (19)

k2Q2z(l1,t) - k1Q1z(l1,t) = 2*Q1(l1,t)f1(l1 ,t).

The relation order of equation right-hand side (18) to the first terms of its left-hand side is estimated by the parameter E = '■s2Q*/^2k2 (for the second term ¡j1k1), where Q* is the characteristic temperature on the interface [3]. These parameters for ordinary liquid media are small and instead of (18) the equality of heat fluxes is used. However, for low-viscosity liquids and small kj the right-hand side in (18) (right-hand sides in (19)) must be taken into account, for example, for cryogenic media [3].

At the initial moment of time, all functions are set

fj (z, 0)= fjo(z), hj (z, 0)= hjo(z), aj (z, 0)= ajo(z), Cj (z, 0)= Cjo(z), Qj (z, 0)= Qj 0(z), (20)

that satisfy the conditions of agreement with (12), (13), (15)-(17), (19). For example, f10(l1) = = f2o(h) etc.

Remark 2. The formulated initial-boundary value problem (4)-(9), (11)-(17), (19), (20) is the inverse, since the functions nj1(t), nj2(t) must be found along with its solution. For a complete statement of this problem, two more conditions must be set

pi 1 pi 2

/ h^t) d£ = 0, / h2(Z,t) d£ = 0, (21)

o i1

which together with the integral equalities (11), (17) mean closedness of motion.

3. Dimensionless variables

We introduce dimensionless variables and parameters

t = ft, e = x = ^ Pj = j m = * =1 = t < i

ll l2 X2 Xj M2 l2

j xï j XI V2X1 j xim'

Hj (t,T ) = XÏÏM hj (z,t), Aj (£,t ) = jM, Cj ) = ,

1 Nj (t) = j, Qj ) = .

jV ' XÏM jK ' a*l\M

(22)

Here Pj are Prandtl numbers, Gj are Grashof numbers, M is Marangoni number. It is further believed that a* = max la^t) > 0 and the characteristic temperature at the interface is

t^o

0* = a*l\.

In the new variables, the system (4)-(8) will be rewritten as follows

ri i r£

FjT + M

F2 + H2 - 2Fj« Fj(Z, t) dÇ + Gj I (Aj (Z, t) + Cj (Z, t)) dZ

HjT - 2M

zo/l2 J Jz0/Î2 (23)

= Pj l2£j Fj«« + Nji(T ),

« ] r«

' Z0/I2

Fj Hj - 2Hj« Fj (Z, t ) dZ + Gj (Aj (Z, t ) - Cj (Z, t )) dZ =

zo/12 (24)

= Pj Hj«« + jT ),

r«

AjT + 2MAj (Fj + Hj ) - 2MAj« f Fj (Z, t) dZ = l2£j Aj««, (25)

J Z0/I2

CjT + 2MCj (Fj - Hj ) - 2MCjtf Fj (Z, t) dZ = l2SjCj««, (26)

J Z0/I2

QjT - 2MQj« Fj (Z, t) dZ = l2ejQj«« + 2ej (Aj + Cj ). (27)

J Zo/l2

In integral expressions for j = 1 the z0 = 0 and at j = 2 we have z0 = li, so that 0 < £ < l in the first layer and l < £ < 1 in the second layer.

The boundary conditions (11)-(13), (15)-(17), (19), (21) are rewritten as

Fi(0,r) = H1(0,r)=0, Fi(1,t) = H2(1,r) = l F2(£,r) d£ = 0, (28)

Ai(0,r ) = aI(t ), cI(0,t ) = Ci (t ), qI(0,t ) = Qi(r ), ai(1,t ) = Ai(r ), C2 (1,t ) = Ci (t ), Qi(1,r ) = Qi(r ),

Aiç(1, t) = Cie(1, t) = Qi6(1, t) = 0, (30)

Fi(l,r ) = Fi(l,T ), Hi(l,r ) = Hi(l,r ), Ai (l,r )= Ai(l,r ), Ci(l,r ) = Ci(l,r), Qi(l,r ) = Qi (l,r ),

Fie(l,r) - MFie(l,r) = -M(Ai(l,r) + Ci(l,r)), Hie(l,r) - MHie(l,r) = -M(Ai(l,r) - C\(l,r)),

(29)

(30)

(31)

f fi(z,t) dt = 0, Jo

a2?(1,t) - kAie(1,t) = 2EA1(l,r)F1 (1,t), C2i(1,t) - kC1i(1,t) = 2ECi(l, t)F1(l, t), Q2t(1,t) - kQie(1,t) = 2EQi(l, t)Fi(1,t),

hi(£,t) = 0,

h2 (z,t) dt = 0.

(33)

(34)

(35)

(36)

The initial data (20) will be of the form

Fj(Z, 0) = Fjo(S), Hj(S, 0) = Hjo(S), Aj(S, 0) = jS), Cj (S, 0) = Cjo(S), Qj (S, 0) = Qjo(Z).

4. Stationary creeping flow with an isothermal interface

In this case, the right-hand sides (32) must be zero. It means that Af(l,T) = Cf(l,T) = 0 and the task set above will be redefined. Here we consider the creeping motion (M C 1). It is necessary to assume that the initial initial data are of the order M. Let M ^ 0, then the equations (23)-(27) will be linear and the right-hand sides of the boundary conditions are equal to zero. However, the relations (34) remain nonlinear.

Remark 3. If, assume that Aj(Z,t) = 0, Cj (Z,t) = 0, then the interface will be isothermal: Ti(x, y,l,T)= T2(x,y,l,T) = 01(1,t) = 02(1,t)=0.

In this paragraph, we assume that the upper plane is thermally insulated and conditions (30) are satisfied on it; initial data (36) are omitted. Let Af, Cf, Qf are specified stationary values of boundary conditions (29). Not complicated, but rather long calculations lead to representations

F2(0

H2(0

Fi(t)

X

' P212

Ai(t) Ci([ 1 "

ait + A{, A2(0 --HZ + cs, C2(t)

a2 = ail + A\, Y2 = Yil + CS,

~ Pil2 G2(a2 + Y2) ^

GJ ai + Yi + Aj + CS A - Niit

24 t3 - 1 6

6 * j 2

- 2 (z2 -d) - N2i (t2 -1)

2

+ Di t, + D2(t -1),

(37)

Hi (t) = X

1

Gi

P212

Pil2 _ G2(a2 - Y2) ^

ai - Yi t4 + Asi - CS - Nig!

24 t3 - 1 6

- 2 (t2 -1)) - N2 (t2 -1)

+ D3t,

+ D4(t - 1).

(38)

The constants Df,... ,D4 are found from the integral equalities (28), (33), (35):

1

Di = D2 =

D3 =

D4 =

3Pil

X

P2l2 1

3Pil X

P2l2

k - Gl( (ai + Yi)P + (Ai + CS)l)

_ ii i v 20 4 j

N2i(l + 2) + G2(a2 + Y2)(l2 +2l - 1)"

3

4

N G ((ai - Yi)l2 + (Ai-CSr

Ni2 - Gi{-20-+-4

~N22(l + 2) + G2(a2 - Y2)(l2 +2l - 1)"

i

0

2

By virtue of (37),

a-2 + Y2 = (ai + Yi)l + A\ + Cf, a2 - Y2 = (ai - Yi)l + A\ - Cf.

(40)

To determine the remaining unknowns ai,y1, Nii, Nu, Nu, N^, there are relations Fi(l) = Fi(l), Fie (l) = MFie (l), Hi(l) = Hi(l), Hie = MHie,

a2«(l) - kAi«(l) = 2l-2EAi(l)Fi(l), C2«(l) - kCi«(l) = 2l ECi(l)Fi(l)

where

Further,

Fi(l) = Pi

1 n + G f(ai + yi)l2 + (a±ch - 6 nii + gi{-40-+-12—

Ai(l) = ail + Ai, Ci(l)= y il + Cf.

F2(l) = -

Hi(l)=2Pi

X(l -1)2 6P2l2

G2(a2 + Y2)(l - 1)

2

+ N2i

Gi(ai - Yi)l2 , Gi(Aï - Cf)l Ni2 - + - — -

20 6 3

H2(l) = -

X(l -1)2 6P2l2

G2(a2 - Y2)(l - 1)

2

+ N22

Fi«(l) = PÏI

2 3 -G

- 3 Ni i + w G i(a i + Yi)l2 + (AÏ + C{ )l

F2«(l) = -

X(l -1) P2l2

2 N + G2(a2 + Y2)(l - 1)

3 N21 +-4-

Hi «(l) = PÏÎ

2 3

- 2 N12 + - G 1 (a 1 - Yi )l2 + -G Ai - Cf)l

20

12

H2« (l) = -

X(l -1) P2l2

2N + G2(a2 - Y2)(l - 1) 3 22 + 4

(41)

(42)

(43)

Now from the first two equalities (41) we find N11 and N21; from the last two equalities (41) we find N12 and N22; from the last two equalities, taking into account the formulas (40), we define a1 + y1, a1 — y1, and therefore a1, y1. Below we find the indicated values for A1 = Cf. This is the case of radial heating of the substrate. Here a2 +y2 = (a1 +Y1)1+2Af, a2 — y2 = (a1 — y1)1. Let's consider the simplest option: a1 = y1 (A^£) = C(£)). Then a2 + y2 = 2(a1 + A1), a2 = y2 and the formulas (37)-(43) are greatly simplified. Unknown will be a1, N11, N21, N12, N22. Calculations show that in the general case

N12 = N22 =0, N11 = Kiai + K2AÏ, N21 = Ksai + KAi,

(44)

where

Ki

K2

Gil2

K3=

1

20 6(l + ¡(1 - l))

Gil 1

6

— (1 + l <' - ") + GV << -1)3

10

4

KA =

6(l + ¡(1 - l)) Pl2

(1 - l)(l + ¡,(1 -l)) Pl2[

Gi l2{ 1+5i (1 -1))+Gf (l -1)3_

3Gil2 , G2(l - 1)2

20

+

Pl

4+¡(1 -1)

(1 - l)(l+¡(1 -1))

Gil + G2(l - 1)2 (3l + \

—+ pl U+¡(1 -1))

P=

p1 P2 .

The constant a\ is the solution of the quadratic equation

EK— + ^^ + EA\(K2 + ai + EKAi-1 = 0.

If the quantity

Ö =

ikP1l

+ EAlK + Kil-1)) - 4EzKiK2A{l

1

(46)

(47)

is positive, then there are two solutions of the equation (46), which means that there are two stationary solutions to the two-layer system. For S = 0 there is one stationary solution, and for S < 0 there are no solutions.

Remark 4. Fpr l = j(1 + j)-1 we get N22 = 1, N12 = (1 — l)(pl)-1, and the formulas (44), (45) retain their form with the replacement of j by j = l(1 — l)-1. As for the functions Qj (£), they are determined by the formulas

Qi(Z) = Qi + aZ - -lf + Aie Q2(0 = b + 2(ail + Ai)(2Z - Z2),

where

a = l

k +

2EFi(l) l2

i

(k + l - l2)(ail + Ai) + 2

, 0 < z < l,

l < Z < 1,

2EFi(l) f ail

(48)

l

(f + Ai)

(49)

b = Qi + ail +3(l - 6)ai - 2A\(l + 2)1-1.

In (48), (49) QI is the dimensionless temperature on the substrate at the origin of coordinates, and F\(l) is given by the equality (42) at a\ = yi, a2 = y2 = a.\l + A\, = Of, and a\ is a solution to the equation (46).

Fig. 1. Dependence 3(E) for various Grashof numbers Gi ; Ai = 0.1

Figs. 1-3 shows the dependences S(E) for various values of dimensionless parameters. All calculations are given for the transformer oil-formic acid system. The dimensionless parameters

2

4

of the physical system are as follows: p = 0.74, v = 15.41, x = 0.71, к = 0.41, в = Pi/P-1 = 1.46, Pi = 308.2, P2 = 14.2. Fig. 1 represent the dependence S(E) or various Grashof numbers G1, G2 = вG1. It can be seen that as G1 grows, the region of existence of two solutions decreases.

Fig. 2 illustrates the dependence S(E) for various values of the dimensionless parameter A®. Here, for certain values of the parameter E, as Ai grows, the region where there are no solutions increases. In the case when Ai < 0 there are always two solutions. Fig. 3 shows the dependence S(E) from the geometric parameter l = l1l2-1 < 1. In this case, with an increase in the thickness of the lower layer, the region of existence of two solutions increases.

Fig. 2. Dependence S(E) for various parameter values Ai

Fig. 3. Dependence S(E) for various values of the geometric parameter I

Conclusion

In the article, the problem of three-dimensional two-layer motion with a special velocity field is reduced to the inverse conjugate problem for a system of one-dimensional integro-differential equations. In the case of a stationary flow at low Marangoni numbers, the solution is obtained in the analytical form. It is shown that, depending on the physical and geometric parameters, two stationary modes can exist. For the transformer oil - formic acid system, the effect of changes in interfacial internal energy on the number of stationary solutions has been studied.

This research was supported by the Russian Foundation for Basic Research (20-01-00234) and Krasnoyarsk Mathematical Center and financed by the Ministry of Science and Higher Education of the Russian Federation in the framework of the establishment and development of regional Centers for Mathematics Research and Education (Agreement no. 075-02-2020-1631).

References

[1] V.K.Andreev, Yu.A.Gaponenko, O.N.Goncharova, V.V.Pukhnachev, Mathematical Models of Convection, Berlin, Boston, De Gruyter, 2020.

[2] N.Aristov, D.V.Knyazev, A.D.Polyanin, Exact solutions of the Navier-Stokes equations with the linear dependence of velocity components on two space variables, Theoretical Foundations of Chemical Engineering, 43(2009), no. 5, 642-662. DOI: 10.1134/S0040579509050066

[3] V.K.Andreev, V.E.Zahvataev, E.A.Ryabitskii, Thermocapillary Instability, Nauka, Siberian brunch, Novosibirsk, 2000 (in Russian).

Об одном ползущем трехмерном конвективном движении жидкостей с изотермической границей раздела

Виктор К. Андреев

Институт вычислительного моделирования СО РАН Красноярск, Российская Федерация

Аннотация. В работе рассматривается двухслойное трехмерное движение жидкостей, поле скоростей которых имеет специальный вид. Возникающая сопряжённая начально-краевая задача для модели Обербека-Буссинеска сведена к системе десяти интегродифференциальных уравнений с полными условиями на плоской поверхности раздела. Показано, что для малых чисел Марангони её стационарный аналог может иметь до двух решений, которые находятся в явном виде. Отдельно проанализирован случай, когда стационарное течение возникает за счет изменения внутренней межфазной энергии.

Ключевые слова: модель Обербека-Буссинеска, межфазная энергия, ползущее течение, обратная задача.