УДК 517.9
Two-dimensional Plane Thermocapillary Flow of Two Immiscible Liquids
Elena N. Lemeshkova*
Institute of computational modelling SB RAS Akademgorodok, 50/44, Krasnoyarsk, 660036
Russia
Received 23.01.2019, received in revised form 29.03.2019, accepted 02.04.2019 The problem of two-dimensional stationary flow of two immiscible liquids in a plane channel with rigid walls is considered. A temperature distribution is specified on one of the walls and another wall is heat-insulated. The interfacial energy change is taken into account on the common interface. The temperature in liquids is distributed according to a quadratic law. It agrees with velocities field of the Hiemenz type. The corresponding conjugate boundary value problem is nonlinear and inverse with respect to pressure gradients along the channel. The Tau method is used for the solution of the problem . Three different solutions are obtained. It is established numerically that obtained solutions converge to the solutions of the slow flow problem with decreasing the Marangoni number. For each of the solutions the characteristic flow structures are constructed.
Keywords: interface, thermocapillary, inverse problem, Tau method. DOI: 10.17516/1997-1397-2019-12-3-310-316.
For many liquids the surface tension is well approximated by the linear function. Then energy equality is simplified and has the form [1]
i. I. nrV (1)
k2 dn - kl dn = divr U. (1)
where = —da/30, a(0) is the surface tension coefficient. In equality (1) kj are coefficients of thermal conductivity and 0j are temperature of liquids, j = 1,2; 0 = 0i = 0i and u = ui = u2 are temperature and velocity vector at the interface r, n is the normal to r directed to the second liquid.
The ratio of right-hand side of equation (1) to the terms in left-hand side of equation (1) is estimated by the parameter E = x0*/^j kj (for the first term it is necessary to assume j = 1 and for the second term j = 2), ^j are dynamic viscosities, 0* is the characteristic temperature on the interface. Parameter E determines the influence of interphase energy on the motion of liquids inside the layers. For ordinary liquids at room temperature this parameter is small. For example, it was obtained in experiments that E ~ 5 • 10~4 for the air - ethyl alcohol system at 0* = 15°C. Therefore, the right-hand side of equation (1) is often omitted and we have the equalities of the heat flux across the interface. However, for low-viscosity liquids and some cryogenic liquids (for example,liquid C02) the influence of interfacial energy must be taken into account. It is known that viscosity decreases rapidly with increasing temperature. Calculations of motion of bubbles in various liquids [2] showed that value E = 0(1) is achieved at sufficiently high temperatures. The maximum values of E are attained near the critical points. So, for water we have E ~ 0.02
*[email protected] © Siberian Federal University. All rights reserved
at e = 303.15 K; E - 0.6 at d = 573.15 K; E - 0.7 at d = 623.15 K (critical point for water is 0Kp = 647.30 K).
In the paper we study a two-dimensional stationary flow of two immiscible liquids with a common interface at which the total energy condition is satisfied. The influence of dimensionless parameters on the structures of emerging flows is studied.
1. Problem formulation
To describe a two-layer stationary flow of viscous heat-conducting fluids in layers bounded by solid walls y = 0, y = h with a common interface y = l < h we assume the following form of velocity and temperature fields
\{x,y)= Wj {y)x, u2j{x,y)= Vj (y), (x,y) = aj (y)x2 + bj (y),
(2)
where 0 < y < l for j = 1, l < y < h for j = 2; u1, u2 are velocity vector components, e is temperature. Such representation of the velocity field corresponds to the well-known Hiemenz solution [3]. The substitution of expression (2) into the equations of motion and heat transfer results in the nonlinear systems of equations for functions wj (y), vj (y), aj (y) and bj (y)
vj wjy + wj = vj wjyy + fj , Wj + Vj
2wj aj + vj ajy = Xj ajyy, vj jbjy = Xj jbjyy + 2Xj aj,
0,
(3)
where Vj are the kinematic viscosities, pj = Hj/vj are the densities, fj are the constants. The pressures in liquids are distributed as
1 vj2 x2
—Pj(x, y) = vj vjy - ~f - fjY + d°j ' d°j = const.
The values of fj characterize pressure gradients along the x axis.
Conditions on the rigid walls are
wi(0)= vi (0)=0, ai(0) = aio, bi(0) = bio, W2(h) = V2(h)=0, a2(h) = a2o, b2(h) = b2o.
The following relations are satisfied on the interface y = l
wi(l) = W2 (l), vi(l)= V2(l) = 0, ai (l)= a2(l), bi(l) = 62(1), \i2W2y(l) - Miwiy(l) = -2œai(l), k2a'2y (l) — kiaiy (l) = œai(l)wi(l), k2b2y (l) — kibiy (l) = œbi(l)wi(l).
(4)
(5)
(6)
The first four conditions in (6) follow from the continuity of the velocity and temperature fields at the interface, and the fifth condition is the dynamic condition for tangential stresses. The last two conditions were obtained by taking into account relation (1) and the linear temperature dependence of the surface tension coefficient. From the condition for normal stresses we obtain that the interface remains flat. This assumption can be fulfilled, for example, under the action of sufficiently large capillary pressure [4].
u
Taking into account the no-slip conditions on the walls, vertical velocities Vj (y) are excluded from the continuity equations
vi(y) = - wi(z) dz, 0 < y < l; V2(y) = - W2(z) dz, l < y < h. (7)
./0 Jy
Let us introduce dimensionless functions and parameters
Wj (e)
h2 Mxi
wj (y), Aj (0
aj (y)
aio
Fj
h4
œai0h3 „
M = -— , E :
X1M2
Mxi œ aio
2 fj, Ç
y h
„ v,
Pj = —,
j Xj
2 h2
where Pj are the Prandtl numbers, M is the Marangoni number. Parameters M, E can be either positive or negative. Then we have the following inverse adjoint boundary value problem in dimensionless variables
W.« = M
AiÇÇ = M
W-2 - WiÇ /" Wi(z) dz
0
2AiWi - Aie Wi(z) dz
0
W
2«
XM
P2
A2ÇÇ = xM
W22 - W2Ç Ç W2(z) dz 2A2W2 - A2Ç i W2(z) dz
Fi
Pi, 0 <e <7,
xF2
P2,
(8)
(9)
, 7 <£< 1,
Boundary conditions for this problem follow from (5), (6) and (7)
Wi(0)=0, W2(1)=0, A1(0) = 1, A2(1) = — = S,
«10
Wi(7)= W2(7), Ai(7)= A2(Y),
W2?(7) - MWi5(7) = -2AI(Y), A25(7) - A(7) = EAi(Y)WI(Y),
r W1(z) dz = 0, i W2(z) dz = 0,
JO J-y
where 7 = l/h < 1, n = f^2, x = x1/x2, k = k1/k2. Integral conditions in (10) allow us to find unknown constants (pressure gradients along the layers) Fj, j = 1, 2. Functions bj(y) are found after solving problem (8)-(10) and they do not affect the velocity field in the layers. Remark. Suppose that |M| ^ 1 and we seek the solution of problem (8)-(10) in the form Wj = W0 + M Wj + ..., Fj = F0+MF1 + ..., Aj = A0+MA1 + .... Then, zeroth approximation has the form [5]
WO(0 =v(i.-T7)2 F20 (-ae2 + 27^), A?(e) = 1 + ce 0 < e < 7,
(10)
6y2P1
XF0
w0(0 = jPr He2 + 2(y + 2)e -1 - 27), A°(O = C2d -1) + s, y < e < 1,
6P2
h
y
1 + yCi - 5
C2 — -:-, F1
Y - 1
F0 — VÎI-Y2 F0, F0
3yPiz
Y2
v(y - 1)(y + M(1 - Y))'
(11)
and parameter Z = 1 + jC1 is the solution of the quadratic equation
y2(1 -y)2E ^
2(y + M(1 - Y))
Z2 + (y + k(1 - y))Z - (k(1 - y) + y5) — 0.
Expressions for dimensionless vertical speeds are determined with the use of (7). Simple calculations show that when E = 0 (no interfacial energy effect) there is unique solution of problem (8)-(10) for small Marangoni numbers.
2. Numerical method and calculation results
To solve problem (8)-(10) the Tau method is used. This method is a modification of the Galerkin method [6]. Let us introduce new variables: = for j = 1 and = (1 — £)/(l — y) for j = 2. Then problem (8)-(10) can be rewritten in the form (primes are omitted)
Li(Wi,Fi) = W^ -
Y2 M
~P7
f 7Î
W2 - WW W1(z) dz 10
N1(W1,A1) = A1ÇÇ - y2M
x(1- y)2M L2(W2,F2 ) = W2ÇÇ - X( - Y)
P2
N2(W2, A2) = A2ÇÇ - x(1 - y)2M
0
0,
2A1W1 - AW W1(z) dz 0
1-(1-y)« W22 - W2Ç / W2(z) dz
0
2A2W2 - A2Ç j W2(z) dz
Y2F1 P1
0, 0<e< 1,
x(1 - y)2 F2 —0 P2 '
— 0, 0<e< 1,
W1(0) — W2(0) —0, A1(0) —1, A2Ç (0) —0, W1(1) — W2(1), A1(1) — A2(1),
i W1(z) dz = 0, i W2(z) dz = 0. Jo Jo
An approximate solution of problem (12)-(16) is sought in the following form
n n
Wjn(€) = E WkR№, Ajn(€) = £ AkRk(0,
k=1
k=0
(12)
(13)
(14)
1 ,, 1 k
-W2Ç (1) + ^W^ (1) —2A1(1), --A2Ç (1) + -A1Ç (1) — -EA1(1)W1(1), (15)
1 - Y Y 1 - Y Y
(16)
(17)
where Rk (z) = Pk (2z — 1) are shifted Legendre polynomials, z g [0,1], Pk (z) are ordinary Legendre polynomials [7]. Taken into account orthogonality of the Legendre polynomials Rk (z) on the interval [0, 1]
f1 1 \
Rk(z)Rm(z) dz — 5kmhm, hm — --—- , 5km — <
Jo 2m +1 [
1, k — m,
0, k — m,
we obtain from integral conditions (16) that W0 = W20 = 0. Coefficients Wjf, A^ and constants Fi, F2 are determined from the system of Galerkin approximations
f Lj (Wjn, Fj )Rm(£) = 0, f Nj (Wjn,Ajn)Rm(0 d£ = 0, m = 0,...,n - 2, j = 1, 2 (18)
00
and transformed boundary conditions (14)-(15)
k=1
k=0
]T(-1)k Wk = 0, ]T (-l)kW2k =0, ]T (-l)k Ak = l^Ak R'k (0) = 0,
k=1 k=0
n n n n
Ewk = £Wk, Y,Ak = £Ak,
k=1
k=1
k=0
k=0
^ E W2k R'k (i) + Rk (i) = ^ Ak,
Y k= 1 Y k = 1 k=0 .. n J n n n
Y^£A2k R'k (i) + -J2Ak Rk (i) = -E £ A^ Wk.
Y k=0 Y k=0 k=0 k = 1
(19)
We also taken into account that Rk(1) = 1, Rk(0) = ( — 1)k. Thus, equations (18), (19) form a closed system of algebraic nonlinear equations for coefficients Wk, Ak and constants Fj, j = 1,2.
Calculations were carried out for parameters of water (j = 1) and water vapor (j = 2) system on the saturation line at temperature 300°C: p\ = 712 kg/m3, p2 = 46.8 kg/m3, vi = = 1.3 • 10-7 m2/s, v2 = 4.2 • 10-7 m2/s, \\ = 1.4 • 10~7 m2/s, X2 = 2.5 • 10~7 m2/s, k1 = = 0.54 W/m•◦ C, k2 = 0.0719 W/m^°C, « = 1.440"4 N/m•◦ C, P1 = 0.96, P2 = 1.69, E = 0.6, M = 16.5, 5 = 1(a1o = a2o), h =1 • 10~9 m, l = 0.5 • 10~9 m and n = 15. Three different values of dimensionless constants F1, F2 were obtained: {F/ = —0.998, F2 = —3.4}, {F2 = 26.241, F22 = 82.04} and {Fj3 = —68.86, F23 = —352.18} (superscript indicates the solution number). The obtained values at n = 15 and at n = 16 differ by 10~2°, 10~12 and 10~7 for F1, F2 and F3 respectively. It means good convergence of the Tau method in solving this boundary value problem. Fig. 1 shows the profiles of dimensionless functions Wj (£) and transverse velocities Vj (£) for values F1 and F2, respectively. Here functions W (£) and V (£) coincide with functions Wj (£) and Vj (£), j = 1, 2 on their domains of definition. Profiles W(£) and V(£) qualitatively coincide
for F3 and F1. However, the flow is more intense for F3 because max IV(£, F3)| = 0.365,
«£[0,1]
max \W(£, F3)| = 6.04, and max |V(£,F 1)| =0.0032, max \W(£,F1 )| = 0.044.
«£[0,1]' V " «£[0,1]
Fig. 2 shows the relationship between dimensionless functions Wj (£) and transverse velocities Vj (£) for F1 and dimensionless parameter 5. It is seen that with decreasing 5 the flow intensity decreases. For F2 and F3 the situation is similar.
One should note that dimensionless parameters E, M, P1, P2 affect only the intensity of flows. For example, for F1 with decreasing number E the flow rate increases and with increasing Prandtl number P1 the flow rate decreases. It was also established that with decreasing Marangoni number M the obtained solutions F1 and F2 tend to solutions of the model problem Fj01 =
-1.036,
F201
-3.506 and F0
02
654.63,
F202
2215.12, respectively (see (11)). For
example, at M = 0.001 we obtain F01 - F/
I0-9 and ¡F02 - F2Y
I0-3, j = 1, 2.
This research was supported by the Russian Foundation for Basic Research (grant 17-0100229).
Fig. 1. Profiles of dimensionless functions W(£) (—) and transverse velocities V(£) (—) for F1 (a) and F2 (b)
Fig. 2. Profiles of dimensionless functions W(£) and transverse velocities V(£) for F1
References
[1] V.K.Andreev, V.E.Zahvataev, E.A.Ryabitskii, Thermocapillary Instability, Nauka, Sibirskoe otdelenie, Novosibirsk, 2000 (in Russian).
[2] F.E.Torres, E.Helborzheimer, Temperature gradients and drag effects produced by convection of interfacial internal energy around bubbles, Phys. Fluids A, 5(1993), no. 3, 537-549.
[3] K.Hiemenz, Die Grenzschicht an einem in den gleichförmigen Flüssigkeitsstrom eingetauchten geraden Kreiszylinder, Dinglers Poliytech. J, 326(1911), 321-440.
[4] R.Kh.Zeytounian,The Benard-Marangoni thermocapillary-instability problem, 41(1998), no. 3, 241-267.
[5] V.K.Andreev, The solutions properties of conjugated nonlinear boundary value problem describing a stationary flow of a two liquids in the channel, Some actual problems of modern mathematics and mathematical education. Herzen readings, Scientific conference materials, 2018, 18-23.
[6] C.A.J.Fletcher, Computational Galerkin method, Springer-Verlag, 1984.
[7] Gabor Szego, Orthogonal polynomials American Mathematical Society Colloquium Publications, Vol. 23, Revised ed. American Mathematical Society, Providence, R.I., 1959.
Двумерное плоское термокапиллярное течение двух несмешивающихся жидкостей
Елена Н. Лемешкова
Институт вычислительного моделирования СО РАН Академгородок, 50/44, Красноярск, 660036
Россия
Изучается задача о двумерном стационарном течении двух несмешивающихся жидкостей в плоском канале. На твердых стенках канал поддерживает заданное распределение температуры. Жидкости контактируют через общую поверхность раздела, на которой учитываются затраты энергии на ее деформацию. Температура в жидкостях распределена по квадратичному закону, что согласуется с полем скоростей типа Хименца. Математический анализ такого течения приводит к возникновению сопряженной краевой задачи, которая является нелинейной и обратной относительно градиентов давлений вдоль канала. Применение к ней тау-метода показывает, что она имеет три различных решения. Численно установлено, что полученные решения с уменьшением числа Марангони сходятся к решениям задачи о ползущем течении. Для каждого из решений построены характерные структуры течения.
Ключевые слова: граница раздела, термокапиллярность, обратная задача, тау-метод.