Reduced-order modeling of fully turbulent buoyancy-driven flows using the Green's function method

A one-dimensional (1D) reduced-order model (ROM) is developed for a 3D Rayleigh-Bénard convection system in the turbulent regime with Rayleigh number $\mathrm{Ra}={10}^6$. The state vector of the 1D ROM is horizontally averaged temperature. Using the Green's function (GRF) method, which involves applying many localized weak forcings to the system one at a time and calculating the responses using long-time averaged direct numerical simulations (DNSs), the system's linear response function (LRF) is computed. Another matrix, called the eddy flux matrix (EFM), which relates changes in the divergence of vertical eddy heat fluxes to changes in the state vector, is also calculated. Using various tests, it is shown that the LRF and EFM can accurately predict the time-mean responses of temperature and eddy heat flux to external forcings and that the LRF can predict well the forcing needed to change the mean flow in a specified way (inverse problem). The non-normality of the LRF is discussed and its eigenvectors or singular vectors are compared with the leading proper orthogonal decomposition modes of the DNS data. Furthermore, it is shown that if the LRF and EFM are simply scaled by the square root of the Rayleigh number, they perform equally well for flows at other $\mathrm{Ra}$, at least in the investigated range of $5\times {10}^5\le \mathrm{Ra}\le 1.25\times {10}^6$. The GRF method can be applied to develop 1D or 3D ROMs for any turbulent flow, and the calculated LRF and EFM can help with better analyzing and controlling the nonlinear system.

Figure 1

Schematic of the 3D Rayleigh-Bénard convection system. The temperature at the top and bottom walls is, respectively, maintained at constant values of $T_t$ and $T_b$, where $\mathrm{\Delta }T=T_b-T_t>0$. The velocity boundary condition at the walls is no slip, $\mathit{u}=0$. The horizontal directions ($x$ and $y$) are periodic. In this study, $\mathrm{Pr}=0.707$ and $5\times {10}^5\le \mathrm{Ra}\le 1.25\times {10}^6$.

Figure 2

Time series from DNS at $\mathrm{Ra}={10}^6$ for (a) anomalous temperature $\overline{T}-\langle \overline{T}\rangle$ at $z_0^{*}=-0.18$ and (b) the principal component of the leading POD (PC1). Both time series are normalized with their standard deviation $\sigma$ so that they have the same variability. Time is scaled by the advective timescale ${\tau }_{\text{adv}}$.

Figure 3

(a) Power spectral densities (PSDs) of the time series shown in Fig. 2. The PSDs are calculated by dividing all the data consisting of 110 000 samples ($\sim 12820{\tau }_{\text{adv}}$) into 500 windows with the same length ($\sim 26{\tau }_{\text{adv}}$) and carrying out a fast Fourier transform for each window. The results of all windows are then averaged to obtain the plotted PSDs. Frequency $\omega$ is normalized by the frequency of advective timescale $\tilde{\omega }=2\pi /{\tau }_{\text{adv}}$. (b) Autocorrelation $r_N$ of anomalous temperatures (solid line) and EHF (dashed line), shown as a function of time scaled by ${\tau }_{\text{adv}}$. The black lines show exponential fits.

Figure 4

(a) Singular values $s_i$ of different POD modes obtained by conducting singular-value decomposition on the longest available data set for horizontally averaged anomalous temperature. (b) Cumulative variance explained by the POD modes. The horizontal dashed and solid red lines mark 0.95 and 0.99, respectively.

Figure 5

Time-mean responses of temperature (left column) and eddy heat flux (right column) to forcings (a) and (b) $f_1$, (c) and (d) $f_3$, and (e) and (f) $f_4$ (see Table 1). Solid lines show the true response obtained from a long forced DNS and dashed lines show the predictions from $L$ and $E$ obtained using the GRF method. The red shading shows the uncertainty in the time-mean responses calculated from the DNS data (see the text for more details). For all cases, $\mathrm{Ra}={10}^6$.

Figure 6

Cases (a) and (b) C5 and (c) and (d) C6. (a) and (c) Solid lines show the target change in the mean flow ${\theta }_{\text{target}}$, while the dashed lines demonstrate the time-mean responses to forcings $f_5$ [shown in Fig. 7] and $f_6$ [shown in Fig. 7] obtained from a long forced DNS. (c) and (d) Changes in vertical eddy heat flux obtained from a long forced DNS (solid lines) or calculated by $E$ (dashed lines). As before, the red shading shows the uncertainty in the time-mean responses calculated from the DNS data. For all cases, $\mathrm{Ra}={10}^6$. More details are in Table 1.

Figure 7

Forcing that is predicted as $f=-L{\theta }_{\text{target}}$ to lead to the target time-mean responses ${\theta }_{\text{target}}$: (a) $f_5$ $[{\theta }_{\text{target}}$ is shown in Fig. 6] and (b) $f_6$ $[{\theta }_{\text{target}}$ is shown in Fig. 6].

Figure 8

Time-mean responses of temperature (left column) and eddy heat flux (right column) to forcings (a) and (b) $f_8$ ($\text{Ra}=5\times {10}^5$), (c) and (d) $f_9$ ($\text{Ra}=7.5\times {10}^5$), and (e) and (f) $f_{10}$ ($\text{Ra}=1.25\times {10}^6$) (see Table 1). Solid lines show the true response obtained from a long forced DNS, dotted lines show the predictions from $L\left({10}^6\right)$ and $E\left({10}^6\right)$, and dashed lines show the predictions obtained by the LRF and EFM scaled according to Eq. (27).

Figure 9

Scaling factors $c_f$ and $c_e$ demonstrated by blue circles and red crosses, respectively. The solid lines are the power fits to the corresponding discrete data whose functions are shown in the figure.

Figure 10

The first four leading (i.e., slowest-decaying) modes and neutral vector of $L$ and the leading POD modes. Solid blue (red) lines show the real (imaginary) part of $L$'s eigenvectors. The green dashed lines show the leading POD modes (time mean removed) obtained from the unforced DNS. The first four leading POD modes of unforced DNS explain around $48\%$, $20\%$, $15\%$, and $7.5\%$. In (a), the magenta dash-dotted and yellow dotted lined display, respectively, the neutral vector of $L$ and POD1 of the data obtained from the stochastic ODE (28); both lines coincide with the eigenvector. The eigenvalues of the shown eigenmodes are (a) $-0.058(1/{\tau }_{\text{adv}})$, (b) $-1.303(1/{\tau }_{\text{adv}})$, (c) $-1.506(1/{\tau }_{\text{adv}})$, and (d) $(-2.726+0.504i)(1/{\tau }_{\text{adv}})$ (see Fig. 11). All these modes are projected onto the grid space and normalized to have the magnitude of one.

Figure 11

(a) Eigenvalues of $L$ nondimensionalized using the advective timescale ${\tau }_{\text{adv}}$. (b) The $\epsilon$ pseudospectrum of $L$ [Eq. (29)] calculated following Trefethen and Embree [80], which shows the non-normality of $L$. The numbers on the isolines are $-log\left(\epsilon \right)$. In (a), only the first 15 eigenvalues (out of 25) are shown for better illustration and to keep the focus on the eigenvalues corresponding to the slowest-decaying modes. This has caused the large difference in the range of $x$ axis of the two panels.