Low-dimensional model of turbulent Rayleigh-Bénard convection in a Cartesian cell with square domain

Low-dimensional model of turbulent Rayleigh-Bénard convection in a Cartesian cell with square domain

Jorge Bailon-Cuba and Jörg Schumacher Institut für Thermo- und Fluiddynamik,
Technische Universität Ilmenau,
Postfach 100565, D-98693 Ilmenau, Germany
July 3, 2019

A low-dimensional model (LDM) for turbulent Rayleigh-Bénard convection in a Cartesian cell with square domain, based on the Galerkin projection of the Boussinesq equations onto a finite set of empirical eigenfunctions, is presented. The empirical eigenfunctions are obtained from a joint Proper Orthogonal Decomposition (POD) of the velocity and temperature fields using the Snapshot Method on the basis of a direct numerical simulation (DNS). The resulting LDM is a quadratic inhomogeneous system of coupled ordinary differential equations which we use to describe the long-time temporal evolution of the large-scale mode amplitudes for a Rayleigh number of and a Prandtl number of 0.7. The truncation to a finite number of degrees of freedom, that does not exceed a number of 310 for the present, requires the additional implementation of an eddy viscosity-diffusivity to capture the missing dissipation of the small-scale modes. The magnitude of this additional dissipation mechanism is determined by requiring statistical stationarity and a total dissipation that corresponds with the original DNS data. We compare the performance of two models, a constant so-called Heisenberg viscosity–diffusivity and a mode-dependent or modal one. The latter viscosity–diffusivity model turns out to reproduce the large-scale properties of the turbulent convection qualitatively well, even for a model with only a few hundred POD modes.

44.25.+f, 47.27.ed

I Introduction

For most turbulent flows in nature and technology, it is impossible to resolve all relevant degrees of freedom. Systematic methods to derive models with a reduced number of degrees of freedom from the full set of nonlinear fluid equations are thus necessary. Low-dimensional modeling of transient and turbulent flows using Galerkin projection onto the empirical basis functions which are obtained from a proper orthogonal decomposition (POD) is one such well established method. Lumley1971 (); berkooz (); moehlis_2 () POD and the development of Galerkin models based on POD modes has been applied to a number of fundamental hydrodynamic flow problems, including simple wall-bounded shear flows, aubry (); berkooz_2 (); moehlis_3 (); moehlis_1 () flows over cavities cazemier (); kalb () or in the wake of a cylinder.deane_1 (); ma_1 (); Noack2003 (); Rowley2006 () The development of low-dimensional models (LDM) based on POD modes has also been extended in several directions such as to the balanced POD method Rowley2005 (); Rowley2008 () or to unsteady flow problems Noack2010 () for which fast and slow flow modes are separated. Most of these cases have been studied for laminar or transitional flows at lower or moderate Reynolds numbers.

With increasing Reynolds number the flows become turbulent, the number of degrees of freedom grows rapidly and their nonlinear couplings are increasingly relevant. The truncation of the set of nonlinear ordinary differential equations (ODE) which follows from Galerkin projection introduces always a cut-off of these mode interactions and removes couplings between the degrees of freedom which are necessary for the transfer of kinetic energy from large to small scales. An additional dissipation mechanism has to be implemented in the low-dimensional model to account for the dominant dissipation by the truncated degrees of freedom. The particular way of truncation can then alter the dynamics in the LDM.

Several approaches to this problem have been suggested in the past. Aubry et al. aubry () used directly the energy transfer between resolved modes and unresolved modes at smaller scales to formulate a spectral closure in the truncated system of ODEs. Moehlis et al. moehlis_3 () presented streamwise-invariant truncations for a plane Couette flow, and showed that very low-dimensional models with up to ten degrees of freedom can reproduce transient flow phenomena in low-Reynolds-number shear flows. They also found that the detailed behavior of their LDM depends in a subtle manner on the modes included and that a proper account of the symmetries of the system is crucial. Later Smith et al. moehlis_1 (); moehlis_2 () included streamwise variations in their model. They also introduced a linear damping term, but only when the particular POD mode expansion coefficient is significantly anti-correlated with its time derivative, . Cazemier et al. cazemier () constructed a LDM for driven cavity flows, consisting of the 80 most energetic POD modes computed from 700 snapshots of a direct numerical simulation (DNS). To study the time evolution of the truncated ODE system, a slightly different linear damping term is introduced in their model. This term is calculated from the requirement that the energy of the ODE system is conserved in a statistically stationary sense.

A few attempts to derive LDMs are reported for Rayleigh-Bénard (RB) convection, despite being one of the most comprehensively studied flows.Ahlers2009 (); Lohse2010 () Studies of RB convection in a finite box, based on the POD procedure, have been mostly done by Sirovich and co-workers. sirovich_4 (); sirovich_5 (); sirovich_6 (); sirovich_2 (); sirovich_3 () Sirovich and Park sirovich_5 (); sirovich_6 () discussed the importance of the discrete symmetries describing the velocity–temperature fluctuations field. Deane and Sirovich sirovich_2 () made a parametric study of the POD mode spectra for small Rayleigh numbers .

Only recently, a snapshot method has been applied to turbulent RB convection in a closed cylindrical cell for Rayleigh numbers up to and cell aspect ratios between one half and three.Bailon2010 () In this work, emphasis was given to relating the first POD modes to the large-scale flow circulation which is always present in a closed turbulent convection cell. Ahlers2009 () The disentanglement of the temperature and velocity fields into POD modes allowed the authors to quantify the amount of heat which is transported by the particular POD modes through the convection cell. A change of the large-scale flow from a one-roll to a two-roll pattern, which is observed when the aspect ratio is increased beyond one at a fixed Rayleigh number, was in line with a decrease of transported heat by the primary mode compared to the secondary POD mode.

As a correspondence of the few POD studies of RB convection, only a few works exist with an emphasis on developing a LDM by a Galerkin projection of the Boussinesq equations onto the most energetic POD modes. Tarman tarman_1 () derived a model from POD modes which have been however separately extracted from the velocity and temperature fields. In a second work he proposed an algorithm which incorporates the lost dissipation due to truncation.tarman_2 () Besides the cutoff index based on the energy (mode index ), a second index based on the dissipation () was considered. The time dependence of the modes with indexes was expressed as the quotient of the corresponding nonlinear and dissipation coefficient. No closed forms for the constant coefficients in the ODE system were however obtained in any of these works.

In the present work, we want to extend these studies of RB convection in several directions. First, we construct a LDM for the evolution of the POD mode coefficients in the case of turbulent Rayleigh-Bénard convection in a Cartesian cell with periodic side walls and isothermal free-slip top and bottom square planes. It is essential to use POD modes of the combined four-vector velocity-temperature field.Bailon2010 () In this derivation, it turns out that a cubic term due to the interaction of the velocity with the mean temperature field (denoted as ) becomes linear as a consequence of the orthogonality of the POD modes. The other terms which arise in the Galerkin projection are a linear production term , a linear dissipation term , a quadratic nonlinear term , and a constant term corresponding to the dissipation due to the mean temperature field. Second, our studies will extend previous works tarman_2 (); sirovich_2 (); sirovich_6 () in terms of the magnitude of the Rayleigh number of convection. A case with is considered for which RB convection is turbulent and a DNS data record exists. Third, we are interested in the long-time behavior of the dynamics in the LDM. With a view to more complex convection flows in the future, we are seeking for the least set of POD modes that can reproduce characteristic dynamics of turbulent convection.

A solution which includes the additional dissipation due to the neglected less energetic POD modes has to be considered by an additional eddy viscosity–diffusivity, . First, we present the so–called Heisenberg model with a constant which exerts the same fraction of dissipation on all POD modes. As will be shown, this closure requires at least a minimum number of degrees of freedom, in particular with respect to the vertical direction, for a qualitatively correct description of the flow. As a consequence, two LDMs with 210 and 310 degrees of freedom, respectively, are chosen. They are taken from a set of 15708 modes (see Sec. III B). As will be seen, this model fails to reproduce the large-scale evolution of convection. For the larger of the two sets of modes, the model relaxes to a statistically stationary state which contains too much energy. Second, we refine this model and include a mode–dependent (or modal) eddy viscosity–diffusivity. The magnitude of both eddy viscosity–diffusivity contributions has to be estimated. In order to do so, we will follow a procedure that has been suggested by Cazemier et al. cazemier (). The second model yields much more realistic large-scale variations of the most energetic modes, also reproducing with reasonable accuracy the energy spectrum and the turbulence statistics. Therefore, a significant part of the present work discusses the impact of both types of eddy viscosity-diffusivity on the dynamics of the LDM with different number of degrees of freedom and how it compares to DNS.

The outline of the paper is as follows. The equations of motion, the basic idea of POD – in particular for the method of snapshots – is discussed in the next section. The construction of the LDM by Galerkin projection onto POD modes of RB convection follows in Sec. III. In this section, the results of the time integration of the LDM with both eddy viscosity–diffusivity schemes, and the agreement with the DNS are also discussed. We conclude with a summary and give an outlook.

Ii Methods

ii.1 Equations of motion and numerical scheme

Turbulent Rayleigh-Bénard convection is governed by the Boussinesq equations. They are brought into a dimensionless form by rescaling with the domain height , the diffusive time scale with being the thermal diffusivity, the temperature difference and follow to


where is the velocity field, the departure from the linear conduction temperature profile, and is the kinematic pressure. Dimensionless parameters are the Rayleigh number and the Prandtl number . Besides diffusivity , they contain the kinematic viscosity , the gravitational acceleration , and the thermal expansion coefficient . Note that the total temperature field is given in our notation by (see also Ref. sirovich_5 ())


The vector is the direction in which buoyancy and gravity work and in which the mean temperature gradient is established. The dimensions of the cell are , where from now on are the horizontal and the vertical dimensionless coordinates. The aspect ratio is fixed to . For convenience, the origin of the coordinate system is in the center of the cell. Therefore, , , and with and . The -, - and -components of the velocity field will be denoted by , and , respectively. The boundary conditions are periodic in and , and free-slip in . This means that at the hot bottom plane at and the cold top plane at the following conditions hold:


For the present boundary conditions, the flow can be decomposed in


where the mean components, e.g. are ensemble averages obtained by averaging over the horizontal - plane and time, i.e., a sequence of statistically independent snapshots. The ensemble average for the velocity component is thus given by


Figure 1 shows the mean vertical profiles of and together with the linear thermal conduction profile.

Figure 1: Mean vertical profiles , and the linear conduction profile (see Eq. (4)). Data are obtained from the DNS data record.

Considering the fact that for our problem, due to symmetry considerations, , we can take the four-vector field


for the POD analysis and LDM derivation. Our analysis is in the statistically stationary regime of convective turbulence. The ensemble average of (3) yields then an expression linking the mean and the fluctuating components of the flow in the form


Following Sirovich et al. sirovich_4 () temperatures and velocities are rescaled with and , respectively. Here is the Nusselt number. The quantity is obtained by demanding that the turbulent heat flux in the center of the cell must be equal to that due to diffusion at the boundary.

The highly-resolved data record is obtained by a pseudo-spectral DNS which uses fast Fourier transformations.Schumacher2008 (); Schumacher2009 () Time stepping is done by a second-order Runge Kutta scheme. For most of the work, we consider a data set with full three-dimensional turbulence snapshots which are separated by four convective time units from each other. The computational grid consists of points in -, -, and -directions, respectively. The spectral resolution is given by . Here, and the Kolmogorov dissipation length. The Rayleigh number is and the Prandtl number .

ii.2 Proper Orthogonal Decomposition (POD)

The POD is a model reduction technique that extracts the most energetic modes from a set of realizations or snapshots of the flow. These POD modes are used as a basis for Galerkin projections of the full set of nonlinear equations thus reducing the infinite-dimensional space of solutions to a finite-dimensional system.berkooz (); sirovich () The two-point correlation tensor or covariance matrix of the four-vector field is defined by


where the asterisk denotes the complex conjugate, the time average and . For Rayleigh-Bénard convection in Cartesian domains with two homogeneous (invariant with respect to translations) directions Eq. (10) takes the form Holmes_1 ()


For a kernel (11), the eigenfunctions have the form


where are integers for the and directions, respectively. The superscript denotes a particular POD mode. The determination of follows then from


where is the Fourier transform of with respect to the homogeneous directions and . The kernel is calculated from the numerical data set by first taking the discrete Fourier transform of each realization in the horizontal plane,


and then averaging the correlation over the entire ensemble of data,


Thus the kernel is Hermitian, non-negative and on physical grounds square integrable, such that the existence of a complete set of vector eigenfunctions given by (13) is assured. Complex conjugation of Eq. (13) and use of Eq. (12) implies that


Due to the reality of the physical space fields Eq. (14) implies that


The associated expansion of the velocity field in terms of the modes is given as


and again reality of the four-vector field implies that . The index runs over the POD modes. The coefficients are calculated by the scalar product in L,

Next, the discrete Fourier transform of is introduced together with the ansatz of in (12) and the orthogonality of the complex exponentials. This gives


ii.3 The method of snapshots

The snapshot method is one way to obtain the POD modes, particularly when the computational grid size becomes large. The time coordinate in the equations above has to be substituted now by an index that runs over the sequence of snapshots. It is based on the fact that (15) is a degenerate kernel.sirovich_4 () Consequently, an eigenfunction of the kernel can be represented as


where as an element of the symmetry group as discussed in Appendix A. The explicit use of symmetries in the problem at hand enlarges the data record with originally snapshots to a total number of snapshots thus improving the convergence. For DNS snapshots we thus end up with samples that can be used to evaluate the POD modes. Replacing the kernel in (13) and using Eq. (21) results to


where represent any two snapshots (including all possible symmetries). Then (22) is the matrix problem which yields the eigenvalues and eigenfunctions . It is clear that it determines just of the empirical eigenfunctions for a fixed tupel . The eigenvalue of the matrix is the total energy (kinetic energy plus temperature variance) of the th POD mode for . Recall that in the present case the four-velocity field is expanded into Fourier modes with respect to and which is characterized by wavenumbers and , respectively. This results theoretically in an infinite set of POD modes.

Iii Results

iii.1 Galerkin Projection of the Boussinesq equations onto the POD modes

Given the full set of nonlinear Boussinesq equations and the POD modes extracted by a snapshot method from the DNS data, we can proceed to derive the LDM. This requires first a Galerkin projection step. Using the dimensionless units introduced in section II.1, Eqns. (2) and (3) can be rewritten together in a four-vector notation with respect to ,


where . Here, correspond to the three spatial coordinates and the term


is a source which drops out in the Galerkin projection procedure. This is due to the divergence-free nature of the POD basis functions and the fact that modes , , respectively. Now one takes the inner product of (23) with modes and inserts the expansion from equation (18). Due to the orthogonality of the POD modes, the following infinite-dimensional ODE system follows


The terms on the right hand side of (25) correspond to the production (), dissipation (), nonlinear transfer (), interaction with the mean flow (), and the dissipation due to the mean temperature field (), respectively. Closed forms for and have been obtained in Refs. sirovich_4 () and moehlis_2 (), but only for the velocity field. The general equations for the terms on the right hand side of (25) are as follows. The production term is given by


and the dissipation term by


The nonlinear mode coupling term is given by


where the coefficients are given by


The last two terms in Eq. (25), and , are calculated for a quasi-steady flow tarman_1 (). The term is related to the ensemble average of temperature, . However, from the ensemble average of (3) we obtain


where . In terms of the eigenvectors,


since is only a function of . However and therefore


Consequently, the equation for the term is


where is given by (33). For the low-dimensional description of free-shear-flow, Rajaee et al. mojtaba () keeps the original time dependence in (33), just replacing and considering this as a running-time-average factor, so the resulting term becomes cubic in the ODE system (25).

Finally, we have a dissipation term which is related to the mean temperature profile, . The term is nonzero only for the purely thermal and real modes . In the present case this corresponds to a maximum energy of . The term is given by


only if , otherwise . In Eq. (35), we used again . This completes the discussion of the different terms that arise due to the Galerkin projection of the Boussinesq equations onto the POD modes.

Model and % Energy
M1 210 76.894
M2 310 82.202
Total 15708 100
Table 1: Parameters of the two LDMs denoted as M1 and M2. We list the range of the horizontal wavenumbers and and the so-called quantum number . In the last two lines of the table, the total number of POD modes that has been calculated is given. It follows from the maximum range of horizontal wave and vertical quantum numbers. The horizontal wavenumber in -direction starts from zero due to symmetry of Eq. (16). The total number of POD modes from the snapshot analysis is thus .
Figure 2: (Color online) Reconstruction of the vertical profiles of mean convective heat flux (top) and temperature deviation from the linear profile (bottom). Data from DNS are compared with the two LDMs, M1 & M2, as well as with the complete set of POD modes obtained from the snapshot analysis.
Figure 3: Estimation of the maximum eddy viscosity and diffusivity following the procedure of Ref. 8. (Top) Ratio for the LDM model M2 with 310 POD modes. The dashed line marks the maximum of the ratio. (Bottom) and are separately shown for the same data.

iii.2 Truncation to a low-dimensional model

The integrals along the –axis contained in the coefficients of Eq. (25) are evaluated on the computational grid of the DNS. Since the integrands are discrete functions of the wave and quantum numbers, we include all degrees of freedom with , , and . This results in a maximum number of POD modes of . Out of this set of POD modes, we select small subsets of the most energetic POD modes which corresponds to the large-scale structures of the convection dynamics.

The choice of for our first LDM denoted as M1 follows from the restriction to modes with and as indicated in Tab.1. This level of truncation builds on experiences from similar studies by Aubry et al. aubry (), Moehlis et al. moehlis_3 (), Holmes et al. Holmes_1 (); Holmes_2 () and Podvin podvin_2 () in wall-bounded shear flows. They reproduced successfully the dynamics close to the walls for a range of values of the so-called Heisenberg parameter (which will be discussed below). Modes with in streamwise and in spanwise direction were taken. As explained by Holmes et al. Holmes_2 (), this choice is due to the fact that for the given spanwise domain length, the cross–stream interactions that contribute to the observed bursts of the velocity are well reproduced with at least five nonzero modes. Moreover, when we choose small quantum numbers, , the present LDM for convection yields solutions which decay monotonically with time. Our model and the resulting degeneracy restrictions for the average field–modes did not allow us to take . A significant improvement of the LDM dynamics is obtained for . This is due to the fact that the two additional average field–modes with have a degeneracy of 2 and will be incorporated. A second LDM called M2 was introduced with range of horizontal wave numbers (see Tab. 1). The latter LDM will be used for most of the following studies.

Figure 2 represents the Reynolds shear stress and the average temperature profile as reconstructed from the two LDMs. The calculation of both profiles is done by integration of Eqns. (33) and (30), respectively and explained in section III.1. As can be seen in the figure, the convergence to the vertical DNS profile for is very slow. Only the significant enhancement of the degrees of freedom up to results in an excellent agreement with the DNS profiles, for the present Rayleigh and Prandtl numbers.

The modes of the LDMs are however not the modes which contribute dominantly to thermal and kinetic energy dissipation. The missing couplings to the small-scale dissipating modes causes numerical stability problems of the LDM. Various methods have been proposed therefore to stabilize the truncated low-dimensional dynamical system as we have discussed in the introduction. Our studies showed that the model introduced by Cazemier et al. cazemier () worked best. This model introduces a closure based on the mean energy balance as derived from Eq. (25) which can be rewritten as


where is the highest quantum number in the LDM (). The three linear terms in Eq. (25) are summarized to


From Eq. (36), one can derive an equation for the total energy by multiplication with and summation over all quantum and wavenumbers. An additional linear damping term, , is then quantitatively determined from the requirement that the mean total energy of the extended dynamical system is in a statistically stationary state, i.e.,


Note, that in this real equation, the last term on the right hand side is nonzero only for the purely thermal modes with . Also, due to orthogonality, the last two indices of have to be equal.

In the so-called Heisenberg dissipation model by Aubry et al. aubry () the action of the neglected modes on the ones contained in the LDM is represented in an average sense, namely as a function of the dynamics of these coherent structures. An even simpler approach was chosen by Omurtag and Sirovich. omurtag () They simply introduced a constant empirical viscosity coefficient for turbulent channel flows, which has a similar effect as the Heisenberg eddy viscosity of Aubry and co-workers. First, we will apply here the constant eddy viscosity–diffusivity, but estimate the magnitude of by the method of Cazemier et al. cazemier () as described above. Since and thus , we also use the same amplitudes of for all fields. Equations (25) follow to


after truncation and addition of an eddy viscosity and diffusivity term


The constant eddy viscosity-diffusivity is determined by


Figure 3 shows the damping term for the POD modes of M2 and its ratio with the corresponding dissipation . In these plots represents the POD mode of the LDM. The modes are ordered by decreasing energy content (see also the second column of table 2). Note also that in Eq. (38) the dissipation term appears only when . In the top panel of the figure, the modes are grouped in accordance with their degeneracy as given in Tab. 2, the first two modes are and (the other two are obtained from complex conjugation), the second pair and , etc.

The order of magnitude of the eddy viscosity and diffusivity term is determined by the maximum value of the ratio , for . It is shown in the next section that all the regimes of interest, and their typical solutions, can be obtained by a variation of the real prefactor in the range for M1 and for M2. However, as pointed out by Kalb and Deane kalb (), this damping term can change sign in contrast to and can thus add as an additional production term.

It is also evident from Fig. 3 that the Heisenberg eddy viscosity-diffusivity introduces an overwhelming damping for the less energetic modes. This will cause stationary convection solutions for both LDMs. These limitations of the model with constant eddy viscosity–diffusivity make it necessary to extend the LDM to a mode–dependent or modal eddy viscosity–diffusivity coefficient, . Such a model matches the decreasing ratio .

order Modes in LDM Degeneracy Energy in%
in LDM
1 1, 2 6.1705 1.0000 4 30.879
2 3, 4 1.6710 0.2708 4 8.362
3 5, 6 0.5999 0.0972 4 3.002
4 7, 8 0.3788 0.0614 4 1.896
5 9, 10 0.3722 0.0603 2 0.932
6 11, 12 0.3583 0.0581 4 1.790
7 13, 14 0.2729 0.0442 8 2.690
15, 16
8 17, 18 0.2315 0.0375 4 1.158
9 19, 20 0.2274 0.0369 4 1.138
10 21, 22 0.2146 0.0348 4 1.078
11 23, 24 0.2078 0.0337 4 1.039
12 25, 26 0.1975 0.0320 4 0.971
13 27, 28