SOLUTION OF THE LINEAR PROBLEM OF THERMAL CONVECTION IN LIQUID ROTATING LAYER Текст научной статьи по специальности «Математика»

DOI: 10.17516/1997-1397-2022-15-3-273-280 УДК 536.252

Solution of the Linear Problem of Thermal Convection in Liquid Rotating Layer

Victor K. Andreev*

Institute of Computational Modeling SB RAS Krasnoyarsk, Russian Federation

Liliya I. Latonova^

Siberian Federal University Krasnoyarsk, Russian Federation

Received 29.10.2021, received in revised form 19.12.2021, accepted 10.02.2022 Abstract. The initial boundary problem arising in the modeling of viscous fluid creeping rotational motion in a flat layer was solved. A stationary solution was found. The quadrature solution in images was obtained using the Laplace transform method. The time convergence of the the non-stationary problem solution to the established stationary solution was proved under certain conditions on the temperature distribution on the walls.

Keywords: thermal convection, Laplace transform, stationary solution.

Citation: L.I. Latonova, V.K. Andreev, Solution of the Linear Problem of Thermal Convection in Liquid Rotating Layer, J. Sib. Fed. Univ. Math. Phys., 2022, 15(3), 273-280. DOI: 10.17516/1997-1397-2022-15-3-273-280.

1. Problem statement

Let us assume that the fields of pressure and temperature velocities are rotationally symmetrical. Then, their values depend only on r = \Jx2 + y2, z, and time t in a cylindrical coordinate system. Moreover, we suppose that the only external force acting on the fluid is the centrifugal force. Then [1], the momentum, continuity, and energy equations can be written as

ut + uur + wuz — 2wv--=--pr + v --— ш2 ßrQ,

r p V r2 /

uv / v \

vt + uvr + wvz + 2wu +--= v Д v--7Г ,

r V r2)

r

wt + uwz + wwz = — pz + vAw, (1.1)

P

u

ur +---+ wz = 0,

r

&t + u&r + w&z = xA©,

where A = d2 fdr2 + r-1d/dr + d2¡dz2 is the axisymmetric part of Laplace operator.

Equations (1.1) are written in the rotating coordinate system with constant angular velocity w relatively to the original inertial system. Its rotation axis and the z axis of the cylindrical

*andr@icm.krasn.ru tllatonova@sfu-kras.ru © Siberian Federal University. All rights reserved

coordinate system r, y, z are coincide. The radial and axial components of the velocity are denoted as u and w, respectively, and v is the deviation of the rotational velocity component from the solid rotation velocity wr. The quantity p characterizes the pressure deviation from equilibrium pressure: pw2r2/2; and O is the temperature deviation from the mean value. The positive constants p, v, x, P are the physical liquid characteristics: density, kinematic viscosity, thermal diffusivity, and volumetric expansion coefficient.

The solution of a system (1.1) is sought in the form [2]

u = rf (z,t), v = rg(z,t), w = w(z,t),

P = 1K(t)r2 + Appw2r2(ln r - 1U h(z,t), (1.2)

2 w 2 \ a 2,

r

O = A In - + T(z,t), a

where A and a is the constant dimensions of temperature and length correspondingly. The substitution of (1.2) in (1.1) results in the system

ft + wfz - 2wg + f2 - g2 = -1 K(t) + vfzz - w2pT,

p

gt + wgz + 2wf + 2fg = vgzz, 2f + wz =0, (1.3)

Tt + wTz + Af = xTzz, wt + ww z = -1 hz + vwzz ■

p

The solution of (1.2) may be interpreted as the following. A viscous heat-conductive liquid fills the layer between flat walls z = ±a rotating with angular velocity w = const around the z axis. The no-slip condition u(r, ±a,t) = 0, v(r, ±a,t) = 0, w(r, ±a,t) =0 is satisfied on them. At the initial instant the velocity and temperature distributions are specified consistent with (1.2) formulas. On the rotation axis r = 0 sinks or sources of heat are distributed with constant linear density 2nAk (k > 0 is the constant liquid thermal diffusivity coefficient). The solid walls (planes) bounding the liquid are ideally heat conductive. All the assumptions above lead to the formulation of an initial boundary value problem for the system (1.3)

f = -2woz(z), g = go(z), w = w0(z), T = T0(z), \z\ < a, t = 0; (L4)

f = g = 0, w = 0, T = Th2(t), z = ±a, t> 0; (1.5)

with the specified functions w0(z), g0(z), T0(z), Tlj2(t). The conditions of thermal insulation of one (or both) walls can be used instead of the last in (1.5), for instance T(-a,t) = T\(t), Tz(a,t) = 0. Note, that for smooth solutions the agreement conditions should be satisfied

wo(±a) = 0, wo z (±a) = 0, go(±a) = 0, To(±a)= Tia(0) (To(-a)= Tx(0), To z (a)=0).

(1.6)

Let us introduce the dimensionless variables by

a2- - _

t = — t, z = az, f = ujR2 f, g = wRg w = a^R2w, T = RAT,

v _ _ 2 (1.7)

K = pw2RK, h = pu2 a2 Rh R = —, P = e = ßA,

v X

where R, P, e are the Reynolds, Prandtl, and Boussinesq numbers correspondingly. Since d2/dt = va-2d2/dt, d2/dz = a-1 d2/dz, we obtain the following system by substituting (1.7) into (1.3) and omitting the upper bars

ft + R3wfz - 2g + R3f2 - Rg2 = fzz - K(t) - eT,

gt + R3wgz + 2R2f = gzz, 2f + Wz = 0, (18)

Tt + R2wTz + R2 f = PTzz, Wt + R3WWz = -hz + Wzz, |z| < 1, t> 0.

The conditions (1.4), (1.5), (1.6) remain unchangeable, it is just needed to take into account that Izl < 1. In addition, w0(z) = wR2w0(z), g0(z) = uRg0(z), T0(z) = RAT0(z), Ti^(t) = = RAT 1,2(t) in the initial data.

2. Linear initial boundary value problem

Let be R << 1; such movements are called creeping. In practice they arise due to the high kinematic viscosity, cross-sectional layer size fineness or small angular velocity w. Assuming that

f = fo + Rfi +----, g = go + Rgi +----, w = w0 + Rw1 +----,

T = To + RTi + ■■■ , K = Ko + RKi + ■■■ , (2.:)

and substituting it into (1.8) we obtain the initial boundary value problem in the zero approximation (the subscript "0" is omitted)

ft - 2g = fzz - K(t) - eT,

gt = gzz, 2f + wz = 0, (2.2)

Tt = pTzz, wt = Wzz - hz, \z\ < 1, t > 0;

f (z, 0) = -2woz(z), g(z, 0) = go(z), T(z, 0) = To(z), w(z, 0) = wo(z), \z\ ^ 1;

f(±1,t)=0, g(±1,t)=0, T(±1,t)= T12(t), w(±1,t) = 0, t > 0. (2.4)

Note, that

r-1

J f (z,t) dz = 0, (2.5)

I-1

what follows from the third equation in (2.2) and non-slip condition (2.4): w(±1,t) = 0. The integral equality (2.5) is correct also for the general problem (1.3), (1.4), (1.5). This additional condition is used to compute the part of radial pressure "gradient", which is the function K(t), see (1.2). Thus, the problem under consideration is an inverse problem.

Let us find the stationary solution of system (2.2)-(2.5). It is denoted as f s(z), gs(z), ws(z), Ts(z), Ks, hs(z) and corresponds to the data TS2 = const. Simple calculations lead to the next

formulas

1

gs(z) = 0, Ts(z) = 2 ((TS - TS) z + TS + TS),

f s(z) = ^ (TS - TS) (z3 - z)

12

-

Ks = - + TS), (2.6)

wS(z) = 24(TS -TS) z -1)2,

hs(z) = h0 + - (Tf - T2!) z (z2 - 1) , h0 = const. 6

The real fields of velocities us(r,z), vs(r,z), ws(z), pressure ps(r,z), and temperature OS(r,z) are given by (1.2).

The solution of inverse problem (2.2)-(2.5) can be obtained using the partition method in the form of Fourier series. First, the functions g(z,t), T(z,t) are to be found as solutions of the first classical initial boundary value problems for the heat conduction equations [3]. After that f(z, t) and K(t) should be determined taking into account the overloading condition (2.5). The function w(z, t) can be recovered by quadrature from the third equation of the system (2.2), and h(z,t) can be found by the latter from (2.2). This solution procedure is rather cumbersome. Here we use the Laplace transform method to find a solution [4]. Let

u(z,s)= u(z,t)e-st dt

Jo

be the Laplace transform for the function u. Since

d

ut(z,s) = Su(z, s) - uo(z), uzz = the problem for f(z, s), g(z, s), T(z, s), K(s) takes the form

fzz - s f = eT - 2g + K - fo (z) gzz - sg = -go(z), Tzz - PsT = -PTo(z), \z\ < 1,

(2.7)

where Tlj2(s) is the Laplace transform of the specified functions Tlj2(t). Moreover, the next conditions are satisfied

T±1,s)=0, T(±l,s)=0, T(±l,s)= Th2(s),

r1 N (2.8)

J f(z, s) dz = 0.

Thus, we obtain the boundary value problem (2.7), (2.8) in Laplace images for ODE systems.

Remark 1. The functions T\ 2(t) can have a finite number of the discontinuities of the first kind [4].

After simple calculations, we obtain a quadrature representation of the solution to the prob-

lem (2.7), (2.8)

9(z,s) =

T(z,s)

go(y)sh [Vs(1 - y)\ dy sh [^(z + -

1 sh(^Vs) - ,

1 iz

J go(y)sh[Mz - y)\ dy,

sh(2\PS) I

T^s)sh [vpS(i - z)] + T2(s)sh [VpS(z +1)

+

/P C1 r i r

- T0(y)sh VPS(1 - y) dy sh VPS(z +1) s J_i L J L

-J Pi To(y)sh \VP~s(z - y)] dy,

Vs J_ i L J

hz,s) = M ( - 1Ï - 1

> —

(2.9)

s \ ch^/s J y/s sh(^y's) x J F (y, s) sh [ yfs(1 - y)] dy sh [ Vs(z +1)] J F (y, s) sh [ V^(z - y)] dy,

where

F(z,s)= eT(z,s) - 2T(z,s) - fo(z).

(2.10)

Now, from equality (2.8) and representation f(z, s) (2.9) we obtain

K (s)

3

(1 - ch(2v/s))

J F(y,s)sh [- y)] dy+

2s cth^L %/s sh (2yfs) J_i

+ J J F (y, s)sh [ ^fs(z - y)] dydz

(2.11)

with F(z, s) defining by (2.9), (2.10).

The functions w(z,s), h(z,s) are determined from (2.2) taking into account differentiation properties of the Laplace transform by the following formulas

W(z,s) = -2 f(y,s) dy,

1

h(z, s) = ho(s) + Wz(z, s) - sW(z, s) + wo(z) =

(2.12)

= ho(s) - 2f(z, s) + 2s J 1 f(y, s) dy + wo(z)

with an arbitrary function ho(s) and the function f(z, s) determined by (2.9).

Under the assumptions that the Laplace transform T1j2(s), dT1y2/dt exists and that there is the limit T1t2 (t) = T(2 = const the following holds because of the property of limit

relations for the Laplace transform (see [4])

lim sTh2(s) = lim Ti2(t)= T3^-

s—ï o t^œ '

(2.13)

Let us demonstrate that lims^o sK(s) = Ks, where Ks is given by (2.6), i. e. that limt^TO K(t) = Ks. It is obviously that sjj(z,s) « 0, s ^ 0. Now we proceed to consider

l

1

z

z

the first approximation of the function sT(z, s) using the Taylor series expansion of hyperbolic functions:

sT(z, s)

1

2\fP~s

TlVP~s(1 - z) + TSVP~s(z + 1)

(2.14)

= 2 T + T2s + T - Ts) z] = Ts(z). Taking into account (2.6) and (2.13) we can obtain provided s ^ 0

sK (s)

j eTS^~s(1 - y) dy + J J eTS^~s(z - y) dydz

3e ~4

f 1 r z

■ J ((TS - TS)y + TS + TS)(1 - y) dy+ + j \fm - TS) y + TS + TS)(z - y) dydz

(2.15)

3e ~4

2 2 4

3(TS - TS) - 2(TS + TS) - -3(TS - TS) + 3(TS + TS)

= T (-¡TS + TS))= kS.

Here, the Taylor series expansions of the following functions were taken into account with the retention of the main terms

3

3

3

2s cth

v^ 2s (j; + ^ + •••) 2^~s (1 + o(s)) 2^~s

1 - ch (2^fs) yfs sh(2^s)

1-1- #-•••

(2.16)

(2V~s +

(2 VS)3

+

-1.

Now consider the limit limS^0 sf(s, z). Since

KS ( ch^z)

s \ 2 ch yfs

1

KS ( 1 + {,fsz)2/2 + o(s2) - 1 - (^~s)2/2 - o(s2)

1 + (Vs)2/2 + o(s2)

k s---

2 2

the following can be derived

sf(z,s) « KS

+

e

1) _ e

2) 4s

z

1

[(TS - TS) y + TS + TS]V~s(1 - y) dyyß(z + 1)+

2V~sJ-i

e = 2

(TS - TS) y + TS + TS

,3

zz

(TS - TS) y - (TS - TS)-

V~s(z - y) dy = fs.

(2.17)

By direct substitution it is easy to show that

sw(z, s) « -2 fz f Sdy = -e (TS - TS)(^ - ^

1

42

= 24(TS - TS) (z4 - 2z2 + 1)=■ sh(z,s) « h0 - 2fS = hS,

(2.18)

6

2

s

2

z

2

z

a

where hS = lims^0 sh0(s). We have proven the

Theorem. Under conditions (2.13), fo(z), g0(z), T0(z) € C[— 1,1] the solution of a nonstation-ary inverse initial boundary value problem (2.2)-(2.5) converges to the stationary solution (2.6) with t ^ m.

Note, that initial values of function K(t) can be found directly from the problem (2.2)-(2.5). The solution formulas (2.9) obtained in the images can be transformed into Fourier series. To show it for the function g(z,t) we will use the first formula for g(z, s) from (2.9). Note, that g(z, s) cannot be translated directly into the original space since the second term does not tend to zero at s ^ m. It can be seen that

g(z,s)=J ^ G(z, y)go(y) dy, (2.19)

where

n( , 1 jsh -s(y + 1)sh -s(1 - z), -1 < y < z;

G(z,y, s) = —=——.—< (2.20)

v 7 yfssh(2^s) [shvTs(z +1)sh^s(1 - y), z < y < 1

is the Green's function for the operator d?/d,z2 — s with zero first-type boundary conditions at z = ±1. It is clear that G(z, y,s) ^ 0 at s ^ m for any z,y € [—1; 1].

Now we can use the result from [5], p. 273, formula No. 188, namely that the image of the function G(z, y, s) corresponds to the original

-n n t

nnz n n ye

n=1

therefore g(z, s) corresponds to the Fourier series

r-l œ „ 1

E2 21

sinnnz sinnnye-n n t = T(z,y,t), (2.21)

/1 w p i

r(z,y,t)go(y) dy = go(y)sin nny dy sin nnze-nl+1. (2.22)

1 n=iJ-1

It is easy to verify that the series (2.21) are the solution to the initial boundary value problem for g(z,t). It is classical provided there is the agreement condition g0( — 1) = g0(1) = 0 and

g'o(y) € l2(-1, 1)

1 r 1

f1 1

gon = go(y)sin nnydy =--

J-1 nn

go(y)cos nny

g'o (y)cos nnydy

1 -1

1 f1 1

— g'o(y)cos nnydy = -ß(n). nn J _1 n

(2.23)

whence it follows that \g0n\ < + 1 /2(n). Then

2 n2 2

w w 11 w 1

\g(z,t)\ ^Y, \gon\ < E 21 + 12 132(n) < m, (2.24)

2 n 2

n=1 n=1 n=1

1

as f g'0(y)cos nny dy ^ 0, n ^ m. The convergence to zero velocity for the function g(z,t) is -1

determined from the inequality

w w

\g(z,t)\ < e-nHY. \gon\e-+(n-1)t < e-nHY. \gon\ = Ce-n(2.25)

since the series \gon\ converges as noted above.

n=1

This work is supported by the 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-2022-876).

References

[1] G.Z.Gershuni, E.M.Zhukhovitsky, Convective stability of an incompressible fluid, Moscow, Nauka, 1972 (in Russian).

[2] V.V.Pukhnachev, Exact Solutions of the Hydrodynamic Equations Derived from Partially Invariant Solutions, Applied Mechanics and Technical Physics, 44(2003), no. 3, 317-323. DOI: 10.1023/A:1023472921305

[3] A.D.Polyanin, Handbook. Linear equations of mathematical physics, Moscow, Fizmatlit, 2001 (in Russian).

[4] M.A.Lavrentyev, B.V.Shabat, Methods of Theory of Complex Variable Functions, Moscow, Nauka, 1973 (in Russian).

[5] G.Doetsch, Introduction to the Theory and Applications of the Laplace Transformation, Springer Verlag, Berlin, 1974.

Решение линейной задачи тепловой конвекции во вращающемся слое жидкости

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

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

Лилия И. Латонова

Siberian Federal University Красноярск, Российская Федерация

Аннотация. Решена начально-краевая задача, возникающая при моделировании ползущего вращательного движения вязкой жидкости в плоском слое. Найдено стационарное решение. С помощью метода преобразования Лапласа решение в изображениях получено в квадратурах. Доказано, что при некоторых условиях на распределение температуры на стенках решение нестационарной задачи сходится с ростом времени к найденному стационарному решению.

Ключевые слова: тепловая конвекция, преобразование Лапласа, стационарное решение.