Thermal convection in a HeleShaw cell with the dependence of thermal diffusivity on temperature

The investigation of thermal convection of a fluid with the dependence of thermal diffusivity on temperature in a vertical Hele Shaw cell heated from below has been fulfilled theoretically.The expression for equilibrium temperature distribution in a cavity has been derived analytically. It has been found that the dependence of temperature on the vertical coordinate looks like a square rootlaw.The linear stability of mechanical equilibrium state against small normal perturbations has been investigated by means of Galerkin method. It has been shown that the most dangerous perturbation in a cavity under considerationis described by the mode which corresponds to the two vortexsteady flow. The numerical simulation of overcritical steady and oscillatory flows has been carried out in the approximation of plane trajectories.This simplification of theoretical modelis consistent with all experimental data on thermal convection in similar cavities. It has been shown that the inclusion of the dependence of thermal diffusivity on temperature into the mathematical model leads to the up down symmetry breakdown for the small values of over criticality.


Introduction 1.1. Thermal convection in the cavities of different shape
In experiments, the evolution of the spatial pattern of convective flows is mainly determined by the shape of the cavity and the type of heating induced. In the case of a model infinite horizontal layer, steady flow has the cellular form and the vortices size is comparable with the thickness of the layer. But in a closed tank there will be a discretisation of allowed wavelengths and flow structure is strongly limited by the shape of a cavity. The type of boundaries also play an important role: the stability of a horizontal layer with free boundaries is fundamentally different from one with rigid boundaries [1]. Also, specific flows arise in cavities due to the particular thermal conditions or the presence of any spatial symmetry.

Dependence of the kinetic coefficients on thermodynamics parameters
In fluid dynamics the coefficient of thermal conductivity and the coefficient of viscosity are usually treated in most calculations as constants. While both theoretically and experimentally, these coefficients can be shown to depend on thermodynamic parameters, the dependence of those coefficients on temperature or pressure is usually neglected as these properties are more often insignificant. For example, the dependence of viscosity on pressure may be completely ignored when dealing with an incompressible fluid. However in the presence of a significant temperature gradient, the dependence of those coefficients on temperature may become substantial. In [2,3] the temperature dependence of viscosity was taken into account in order to explain the up-down symmetry violation of steady and oscillatory convective flows in a vertical Hele-Shaw cell heated from the bottom. Thus in the presence of strong environmental conditions those dependencies can not be trivially ignored.

Dependence of thermal diffusivity on temperature
In the analysis of most of convective problems, the coefficient of thermal diffusivity χ is assumed to be constant. In reality, it can depend both on temperature and pressure from general point of view. However, experiments have shown that the dependence of thermal diffusivity on pressure is much less significant than its dependence on temperature. While extremely small value the dependence on pressure can be neglected, the extent of the effect of temperature is less clear. In particular, the dependence of thermal diffusivity on temperature could help explain the symmetry breakdown in convective flows in Hele-Shaw cell.
Free thermal convection in a vertical Hele-Shaw cell has been studied thoroughly for a long time [4]. It was shown experimentally and mathematically that the convection threshold is determined by a steady flow in a cavity heated from the below. The number of vortices is governed by the degree of the cell elongation along horizontal axis. However the effect of the temperature dependence of thermal diffusivity on the boundary of stability and non-linear convective regimes has not been examined yet. In particular the stability of these states with χ(T ) has not been studied in detail and the effect could be observable and measureable in an experiment 2 Statement of the problem 2.1 Hele-Shaw cell Let us consider rectangular cavity with specific aspect ratio when the height and length are much greater than the thickness (see Fig. 1). It means mathematically that h, l >> d. This cavity is called by the HeleShaw cell [4] The cavity is filled by a fluid with a temperature dependent thermal conductivity. The Hele-Shaw cell is heated from the bottom, Θ is the temperature difference between upper and lower heat exchangers. The approximation of plane trajectories is valid for a wide range of the governing parameters. Even chaotic regimes for large values of over-criticality visually have no transverse component of velocity (v z = 0) and the trajectories of liquid particles lie in the (x, y) plane.

Equations of thermal convection
Now let us take the well known generalized equations of the thermal convection in Boussinesq approximation to describe the behaviour of the incom-pressible fluid in a Hele-Shaw cell. This system includes the Navier-Stokes equation, generalized equation of heat transfer and mass conservation law: Here v, p and T are dimensional fields of velocity, pressure and temperature. Parameters v and χ are the coefficients of kinematic viscosity and thermal diffusivity, respectively, β is the coefficient of thermal expansion, ρ is the average density of the fluid, g is the gravitational acceleration, and γ is the unit vector oriented vertically upward. We assume the viscosity to be constant to simplify our problem. On the other hand, thermal diffusivity is determined by where κ and C p are the coefficient of thermal conductivity and heat capacity, respectively. The dependence of χ on temperature like in (2.2) comes from treating the density and thermal capacity as constants.
Let us summarize experimental data [5,6] and suppose that the dependence of thermal diffusivity on temperature can be expressed by the simple formula: Here T is the deviation of temperature on conventional convective zero, α is dimensional material coefficient which characterizes the dependence of thermal diffusivity on temperature. Theoretically it can be negative as well as positive. If α is negative, the full temperature additive in formula (2.3) must be sufficiently small and have to be less than the unit. Note that thermal diffusivity is essentially a positive parameter.
The sides of the cavity are supposed to be rigid, therefore the velocity vanishes on the boundaries. As a result of the heating from below the specific temperature distribution takes place on the boundary. Namely, the unknown velocity and temperature fields in equations (2.1), (2.2) must satisfy to the following boundary conditions: These relations account for the existence of no slip condition on all sides of the cell. Furthermore, the cavity has boundaries with high thermal conductivity.
Let us fulfil the numerical analysis of the problem (2.1), (2.2) in the terms of non-dimensional variables. We shall use following set of units during the simulation: We use the ordinary conventions in thermal convection: the viscous and thermal conductive units to measure the time and velocity, respectively. The equations of thermal convection in non-dimensional form can be written as There are three governing parameters in equations system (2.4) − (2.6) First and second parameters are the Rayleigh and Prandtl numbers, correspondingly. Non-dimensional parameter ε describes the dependence of thermal diffusivity on temperature.
In our statement the conditions on the upper and lower boundaries of the Hele-Shaw cell for the nondimensional fields of velocity and temperature have the form: Note that in experiment the variation of Rayleigh number implies usually the change of temperature difference on the upper and lower boundaries. But parameter ε also depends on Θ. Thus the growth of Rayleigh number for the fixed ε must be interpreted as the increase of the buoyancy force in comparison with the effect of thermal diffusivity dependence on temperature for the constant heating.
In addition it is important to estimate the value of parameter ε for any typical case. Thermal diffusivity and kinematic viscosity for metal melts can be evaluated as χ ∼ 10 −6 m 2 /s, v ∼ 10 −7 m 2 /s Pr ∼ 10 −2 and β ∼ 10 −4 1/K. For the value of Rayleigh number 3 · 10 2 (over-criticality is equal to 1.07 in Hele-Shaw cell with aspect ratio for wide boundaries 20:40) and thickness 3 mm(d = 1.5 mm) the character temperature difference on the heat exchangers is estimated as 4.4 K. Thus, for α ∼ 0.051/K [5] nondimensional parameter ε is evaluated as 0.22.

Mechanical equilibrium state
The mechanical equilibrium state takes place for small values of Rayleigh number. Mechanical equilibrium is characterized by the following conditions: As a result, the equations system (2.4), (2.5) turns into The exact solution of these equations has to be found taking into account the expression (2.6) for thermal diffusivity and boundary conditions for tem-perature. For different signs of the ε the solution can be written in following form: where H is the non-dimensional height. The lower mathematical signs correspond to negative values of ε i upper signs describe the temperature distribution for positive ε.
The curves characterizing the temperature distribution in dependence on vertical coordinate are presented in Fig. 2 for different values of ε.
In spite of linearity of the expression (2.6) for χ(T ) the profile of the temperature is appreciably nonlinear in the state of mechanical equilibrium. Nevertheless, it is easy to see that the dependence becomes linear in the limiting case ε = 0.

Linear stability problem
Let us divide the full fields of velocity, temperature and pressure into the two parts which correspond to basic state (index 0) and small non-stationary perturbations (prime): Equilibrium state in Hele-Shaw cell is considered as a basic one therefore v o = 0. Let us substitute (3.1) into the starting equations (2.4), (2.5) and linearize them. The resulting evolutionary equations for small disturbances have the form: Here v, p and T are the perturbations of velocity, pressure and temperature, correspondingly. The lateral vertical sides of the cavity are isothermal.
In addition on upper and lower heat exchangers the boundary conditions for unknown fields of the velocity and deviation of temperature from equilibrium profile have to be written as Preliminarily, it is convenient to rewrite initial differential equations T, ψ ∼ (ϑ, σ)e −λt e inπx/L cos πz 2 Here ϑ, σ are the amplitudes of the respective perturbations; L is nondimensional length, λ is an exponential decay rate, n is the integer value that characterizes the periodicity of the perturbation along the x axis. The stream function is connected with the components of velocity by means of relations v x = ∂ψ/∂y, v y = −∂ψ/∂x.
There are two distinct cases: ε is positive or negative. Later we will only examine the case with ε > 0 in our calculations that means the growth of thermal diffusivity with the temperature.
In the following stage we assume the decrement to be equal zero. The perturbations of this type do not grow and do not decrease. These perturbations are called by the neutral one. In this case Rayleigh number Ra plays the role of the eigenvalue in the spectral amplitude problem. It depends on all parameters like H, L, n and ε. It should be emphasized that the linear problem for small disturbances does not contain the Prandtl number. So the boundary of stability does not depend on this parameter.

Limiting case of ε = 0
At first let us consider the limiting case of ε = 0 and find analytically the expression for neutral curves i.e. the function of critical Rayleigh number in dependence on the parameters of the problem. Setting the exponential decay rate λ to zero, we shall be interested in neutral disturbances. If dependence of thermal diffusivity on temperature is absent, the temperature distribution along vertical is linear and the equations system for amplitudes ξ and θ of small perturbations can be represented in the following form [3]: Here ∆ 1 is the plane Laplacian. We shall search the solution of these equations in the form of simple harmonics The analysis of this formula for L = 20, H = 40 shows that the three lower values of the Rayleigh number in our spectrum correspond to the modes (2, 1), (3, 1) and (1, 1) (here the indexes describe the periodicity of the flow along axes x, y and conform to the different values of n and m, respectively).

General case of arbitrary ε
In general case for arbitrary values of parameter ε the calculations were fulfilled by means of Galerkin method [8,9]. The coefficients in equations (3.2), (3.3) depend on only y. So the solution on x can be found in the form of simple harmonics (3.6). As a result the system of two ordinary differential equations with variable coefficients for amplitudes σ and ϑ takes place In a limiting case of ε = 0 the analysis of (3.8), (3.9) leads to important result that the three most dangerous mode for the cavity with the aspect ratio 2:20:40 are (2, 1), (3, 1) and (1, 1).
Thus, let us apply the Galerkin's procedure for both unknown functions σ and θ with only one basis function: (σ, ϑ) ∼ sin(πy/H).
Termwise integration of (3.11), (3.12) with the weight was carried out numerically. The resulting data for critical Rayleigh numbers in dependence on parameter ε are presented in Fig. 3 in the case of a cavity with L = 20, H = 40. As parameter ε increases, the critical Rayleigh number increased too: when ε = 0, Ra c = 281 and when ε = 0.8, Ra c = 417.

Non-linear regimes of convection
The geometry of the problem permits to use the approximation of plane trajectories. As a result it is convenient to execute numerical simulation on the basis of the full non-linear equations written in terms of stream function and vorticity. Let us exclude the pressure and introduce stream function by the same way as in the course of the solution of the linear problem. Finally the standard system of non-linear equations in terms of the stream function ψ, vorticity ϕ and temperature T has the form Boundary conditions for equations system (4.1), (4.2) have to be written in the form: x = 0, L : ψ = ψy = 0, T = T 0 (y); y = 0, H : ψ = ψ x = 0, T = 1, 0.
It is convenient to divide the full temperature field on equilibrium part T o (y) and deviation T 1 (x, y, t). The small thickness of the cavity permits to impose the certain form of the solution depending on the transversal coordinate z and apply the Galerkin procedure. Let us separate the variables and search the solution as ψ = σ(x, y, t) cos(πz/2), where the σ, ϑ are the amplitudes of stream function and temperature which depend on time and coordinates in the plane of wide boundaries. After the averaging on z with the corresponding weight the equations for amplitudes can be written as ∂ϑ ∂y x = 0, L : σ = σ y = 0, ϑ = 0; y = 0, H : σ = σ x = 0, ϑ = 0.

Description of numerical procedure
The method of finite differences was applied to calculate over-critical regimes. Numerical code was written in programming language Fortran-90 . The problem was considered in terms of the vorticity and stream function i.e. the so-called two-fields method of the solution was used. Explicit scheme was realized to simulate the dynamics of convective system. The basic mesh contained 37:31 nodes. The first order coordinate derivatives were approximated by the central differences with the second order accuracy. The time derivative was expressed by the one-sided difference with the first order accuracy. The Laplace operator was factorized on the base of the three point scheme and had the second order accuracy. The values of vorticity on the lateral boundaries were found over the formula of Thom and Aplte [10]. The step size of the time was calculated in accordance with the necessity of the stability of numerical procedure over the formula τ = min h 2 x , h 2 y /4δ, where h x , h y are the coordinate steps on x and y, δ is the empirical parameter greater than unit. The Poisson equation for the stream function ψ was solved by the method of simple iterations [10]. Over the numerical simulation pseudoviscosity method was used to get the snap fields of the temperature T and the stream function ψ. To analyze the oscillatory regimes maximum and minimum values of stream function ψ max and ψ min were calculated.

Results and discussion
The numerical calculations with isothermal wide boundaries of the Hele-Shaw cell yielded the following results. The first over-critical regime with the lowest value of Ra was found to correspond to the characteristic two vortex flow. This result was independent of the value of ε. For ε = 0.2, the relation between the flow regime and Ra was investigated in detail.
At low Rayleigh numbers below Ra = 330, a stationary equilibrium state flow was observed. As Ra was increased, we found a two vortex flow regime with left-right symmetry but this flow had no up-down symmetry. The isolines of stream function and temperature are shown in Figs. 5,6 for Ra = 360, Pr = 7, and ε = 0.5. The displacement of vortices to the upper part of the cavity is observed. There is a simple explanation of this effect.
The bowed equilibrium profile of temperature can be replaced approximately by the broken line with two straight pieces (dashed lines, Fig. 2). There are two regions in a cavity at the plane of wide boundaries with different characteristic values of temperature gradient. For small values of supercriticality the higher part of temperature distribution with large derivative induces convection but the lower part of the broken line may correspond to the small temperature gradient which is not sufficient to produce convective motion. Therefore the stagnant zone has been formed in lower part of the Hele-Shaw cell. The space with the slow movement can be observed from the fields of stream function and temperature (Figs. 5, 6).
At Ra = 455 and beyond, we found an oscillatory four vortex flow with reunification of corner vortices. Initially the oscillations in this regime were periodic and contained only one frequency at Ra = 455. However in the range 460 < Ra < 500, the dependence of ψ max on time clearly contained more frequencies and a Fourier analysis is required to understand these complexities. We also looked at the dependence of ψ max on Ra with different values of parameter ε for flows over the threshold of convection. It can be seen that for all values of ε, ψ max increases with the Rayleigh number.
The x-intercepts of these graph had to be extrapolated from the graphs as calculating the value of ψ max close to zero would have large errors induced by our numerical approach for big values of ε. Furthermore, the x-intercepts obtained here are close to the projections that can be obtained from Fig. 3. The difference in the two values obtained is due to the error introduced by the using of only one basis function in Ga− lerkin's method.
If we introduce point perturbation at the initial moment of time, the two-vortex steady regime is characterised by the down flow along lateral sides and upward movement along centerline. This flow gives the thermal spot in central part of the cavity (Fig. 6).  Fig. 7 for Ra = 500. The Fig. 8 demonstrates that the oscillations become more non-linear and their amplitude increases with the growth of governing parameter Ra. At the same time the period of oscillations has the tendency to be smaller with the increase of Rayleigh number.
The numerical values of period can be found in the table. For large values of super-criticality this reunification of corner vortices becomes nonperiodical.

Conclusion
Convective flows in a Hele-Shaw cell have been investigated theoretically when the dependence of thermal diffusivity on temperature is taken into account. It has been shown for the cavity with aspect ratio 2 : 20 : 40 that the inclusion of this factor in our model leads to the symmetry breakdown of the steady twovortex flow for small values of over-criticality. Visually it looks as the displacement of vortices to upper part of the cavity and formation of the stagnation zone near the lower boundary. The different oscillatory flows originate for bigger values of Rayleigh number. It was found that these oscillatory regimes correspond to various deviations on the well-known fourvortex flow with alternating reunification of opposite corner vortices.