UDC 551.465.11

GEOPHYSICS

A. S. BLOKHIN, N. S. BLOKHINA

ONSET OF CONVECTION IN A LIQUID NEAR THE TEMPERATURE OF DENSITY INVERSION

(Presented by Academician V. V. Shuleikin on 11 XII 1969)

Natural convection in various bodies of water during the period of autumn cooling and spring warming is substantially affected by the presence, in water, of a density maximum at \(4^\circ\). This density maximum also determines the features of the hydrological regime of bodies of water in the winter period, when their surface is covered with ice.

Recently, problems of free convection of a heavy liquid have been studied by a number of authors using numerical methods \((^{1-7})\). In these investigations it was assumed that the density of the liquid is a linear function of temperature

\[ \rho(T)=\rho(T_0)-\beta(T-T_0), \tag{1} \]

however, for water near \(4^\circ\) the following relation is valid

\[ \rho(T)=\rho(4^\circ)-\gamma(T-T_0)^2. \tag{2} \]

Such a dependence, usually called density inversion, should lead to a pattern of motion substantially different from the pattern of convection in a liquid with a linear dependence of density on temperature (1). If a liquid is heated from below and inside its mass there exists the isotherm \(4^\circ\), then below this isotherm the liquid layer has density stratification that is Rayleigh-unstable, while above it the stratification is stable. The interaction of the unstable and stable layers determines the qualitative peculiarity of water convection in the vicinity of \(4^\circ\). Under laboratory conditions, interesting experimental studies of convection in a layer of water formed by melting ice were carried out by Y.-C. Yen and F. Galea \((^{8,9})\). In a paper by C. Tien \((^{10})\), an attempt was made at a theoretical analysis of the influence of density inversion on the onset of convective motion.

In the present paper we give the results of a numerical study of free convection in a heavy viscous liquid with a quadratic dependence of density on temperature (2). In view of the limited capabilities of computing machines, we considered a two-dimensional problem of stationary convective flows inside a square cavity with side \(a\). On the lower wall the temperature was assumed constant, \(T_1=8^\circ\), on the upper wall \(T_2=0^\circ\), and on the lateral walls—varying according to a linear law. As the origin we chose the lower left corner of the square; the \(X\)-axis is directed along the lower side of the square to the right, and the \(Y\)-axis along the side upward.

We write the system of equations in dimensionless form, choosing the following quantities as units of measurement: for the unit of distance, the side of the square \(a\); of time, \(a^2/\nu\); of temperature, \((T_1-T_2)\); of the stream function, the kinematic viscosity of the liquid \(\nu\):

\[ \begin{aligned} \frac{\partial \varphi}{\partial t} +\left( \frac{\partial \Psi}{\partial y}\frac{\partial \varphi}{\partial x} - \frac{\partial \Psi}{\partial x}\frac{\partial \varphi}{\partial y} \right) &= \Delta\varphi + G\,\frac{\partial (T-0.5)^2}{\partial x}, \\ \frac{\partial T}{\partial t} +\left( \frac{\partial \Psi}{\partial y}\frac{\partial T}{\partial x} - \frac{\partial \Psi}{\partial x}\frac{\partial T}{\partial y} \right) &= \frac{1}{P}\Delta T. \end{aligned} \tag{3} \]

Here the field of the stream function is \(\Psi(x,y)\); the temperature field is \(T(x,y)\); the vorticity field is \(\varphi(x,y)=-\Delta\Psi(x,y)\) \((\Delta=\partial^2/\partial x^2+\partial^2/\partial y^2)\); the Grashof number is \(G=[g\gamma(T_1-T_2)^2/\nu^2]a^3\) \((g\) is the acceleration of gravity); the Prandtl number is \(P=\nu/\chi\) \((\chi\) is the thermal diffusivity).

We write the boundary conditions in the form

\[ \begin{gathered} T(x,0)=1,\qquad T(x,1)=0,\qquad 0\le x\le 1,\\ T(0,y)=T(1,y)=1-y,\qquad 0\le y\le 1,\\ \Psi(x,0)=\Psi(x,1)=\partial\Psi(x,0)/\partial y=\partial\Psi(x,1)/\partial y=0,\\ 0\le x\le 1,\\ \Psi(0,y)=\Psi(1,y)=\partial\Psi(0,y)/\partial x=\partial\Psi(1,y)/\partial x=0,\qquad 0\le y\le 1. \end{gathered} \tag{4} \]

The system of equations (3) and the boundary conditions (4) were replaced by an explicit difference scheme and boundary conditions described in work \((^3)\). In the present work they are not written out because of their bulkiness. The calculations were carried out on a BESM-4 computer using a \(17\times17\) grid. The time step was selected experimentally from the condition of stability of the computation. As initial conditions we took:

\[ T_0(x,y)=1-y \tag{5} \]

\[ \begin{aligned} \Psi_0(x,y)={}&(1+\cos 2\pi (x-0.5))(1+\cos 2\pi (y-0.5))\\ &\times(1+0.076\cos 2\pi (x-0.5)-0.132\cos 2\pi (y-0.5))^* . \end{aligned} \tag{6} \]

Fig. 1. \(a\) — field \(T\); \(b\) — field \(\psi\)

When 4–5 decimal places were established in successive printouts of the hydrodynamic fields (the printouts were made every 0.25 units of dimensionless time), it was considered that the flow had reached a stationary regime.

Variants were calculated for \(G=2000;\ 3000;\ 3250;\ 3500;\ 3750;\ 4000;\ 6000;\ 7000;\ 8000\). The Prandtl number was taken equal to \(P=11.59\). This value corresponds to the temperature of density inversion of water.

For \(G=2000;\ 3000\) the initial motion decayed: the field \(\Psi(x,y)\to0\), and the temperature field \(T(x,y)\to1-y\), i.e., stationary convective motion was absent.

For \(G=3250\), a stationary convective flow was established (Fig. 1). To regions with stratification Rayleigh-unstable there corresponds its own independent circulation system; to regions with stable stratification, their own, with the circulation of the lower unstable zone penetrating into the stable zone \((T>0.5)\). The flow is asymmetric with respect to the vertical line passing through the middle of the cavity. In the zone with stable stratification, the isotherms differ little or even practically not at all from the isotherms of the temperature field (5), which is established in the absence of convective motion. In the zone of unstable stratification, the convective currents deform the isotherms of field (5): the rising fluid carries the horizontal isotherms of field (5) upward, and the descending fluid carries them downward.

\[ \text{* See formula (1.10) of work }(^4). \]

For \(G = 3500;\ 3750;\ 4000\) the general character of the motion is preserved; however, the velocities of the liquid flow increase with increasing Grashof number, the small vortex in the lower right-hand corner decreases in size and disappears, while the vortex in the lower left-hand corner grows. The isotherms in the zone of unstable stratification are distorted more strongly, while in the stable zone they remain practically unchanged.

For \(G = 6000\) (see Fig. 1) the steady flow assumed a different character. The circulation both in the lower and in the upper zone split into two vortices each. The hydrodynamic fields became symmetric with respect to the vertical line passing through the middle of the cavity. The liquid rises upward near the walls and descends downward in the middle of the cavity. It should be noted that in this case the isotherms in the zone of stable stratification differ little from the horizontal isotherms (5). For \(G = 7000\) and \(G = 8000\) the picture is the same as for \(G = 6000\), only the motion has become more intense.

It is interesting to compare the symmetric branch of the flow investigated by us with the corresponding branch of the flow of a liquid with a linear dependence of density on temperature. In [4] it was established that in this case the symmetric branch of the flow is metastable; we obtained a stable symmetric branch. From Fig. 10 of paper [4] it is seen that the liquid rises upward in the middle of the cavity and descends downward along the edges; therefore the isotherms are carried upward by the rising liquid in the middle of the cavity and descend downward near the walls. This picture is the reverse of the flow pattern presented by us in Fig. 1 for \(G = 6000\).

The upper independent circulation isolates the lower circulation of warmer liquid from the cold upper wall. Indeed, the velocities in the upper zone of stable stratification are 20–30 times smaller than the velocities in the lower, unstable zone. The isotherms near the upper wall practically do not differ from the horizontal isotherms of field (5). This means that the temperature gradients and, consequently, the heat fluxes in the zone of stable stratification are the same as in the case of the absence of convective fluxes. Thus, the cooling of a reservoir slows down after a layer of liquid with temperature below \(4^\circ\) forms in its upper part, which agrees with numerous observations under natural conditions.

Moscow State University
named after M. V. Lomonosov

Received
30 XI 1969

REFERENCES

1. W. A. Deardorff, J. Atmosph. Sci., 21, No. 4 (1964).
2. J. E. Fromm, Phys. Fluids, 8, No. 10, 1757 (1965).
3. G. Z. Gershuni, E. M. Zhukhovitskii, E. L. Tarunin, Izv. AN SSSR, Mekh. zhidk. i gaz., No. 5, 56 (1966).
4. G. Z. Gershuni, E. M. Zhukhovitskii, E. L. Tarunin, ibid., No. 6, 93 (1966).
5. M. R. Samuels, S. V. Churchill, Am. Inst. Chem. Eng. J., 13, No. 1, 77 (1967).
6. V. I. Polezhaev, Izv. AN SSSR, Mekh. zhidk. i gaz., No. 2, 103 (1967).
7. V. I. Polezhaev, ibid., No. 5, 124 (1968).
8. Y.-C. Yen, Phys. Fluids, 11, No. 6, 1263 (1968).
9. Y.-C. Yen, F. Galea, ibid., 12, No. 3, 509 (1969).
10. C. Tien, Am. Inst. Chem. Eng. J., 14, 652 (1968).